发明内容
本发明的目的在于针对现有技术的不足,提供了一种基于多波束测深的声速剖面及海底地形的联合反演方法和基于多波束测深的水温剖面的反演方法。
为实现上述目的,本发明所采取的其中一种技术方案是:本发明基于多波束测深的声速剖面及海底地形的联合反演方法包括如下步骤:
(1)通过多波束测深***的发射换能器阵向测量海域的海底发射多波束,并通过多波束测深***的接收换能器阵接收回波信号;
(2)多波束测深***根据接收到的回波信号,获得回波的到达角和到达时间;
(3)建立由式(1)所示的状态方程和式(2)所示的测量方程构成的状态空间模型:
式(1)和式(2)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,k≥1;zk表示第k个波束所到达的海底深度;当k>1时,zk-1表示第k-1个波束所到达的海底深度;当k=1时,zk-1表示海底深度的初始值;g表示声速梯度;μk表示第k个波束所到达的海底深度的过程噪声;yk表示由第k个波束回波的到达角的余弦值构成的测量变量;tk表示第k个波束回波的到达时间,c0表示多波束测深***的接收换能器阵所在位置的声速;vk表示第k个波束回波的到达角的余弦值的测量噪声。
(4)根据步骤(3)建立的状态空间模型及步骤(2)中接收到的回波的到达角和到达时间,利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值;根据声速梯度的反演值,利用式(3)计算获得测量海域的声速剖面的估计:
c(z)=c0+g(z-d0) (3)
式(3)中,c(z)表示测量海域的声速剖面;z表示海水深度;d0表示多波束测深***的接收换能器阵所在位置的深度;c0表示多波束测深***的接收换能器阵所在位置的声速;g表示声速梯度。
进一步地,本发明在所述步骤(4)中,所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。
本发明所采取的另一种技术方案是:本发明基于多波束测深的水温剖面的反演方法包括如下步骤:
(1)通过多波束测深***的发射换能器阵向测量海域的海底发射多波束,通过多波束测深***的接收换能器阵接收回波信号;
(2)多波束测深***根据接收到的回波信号,获得回波的到达角和到达时间;
(3)建立由式(4)所示的状态方程和式(5)所示的测量方程构成的状态空间模型:
式(4)和式(5)中,k表示发射换能器阵向测量海域的海底发射的多波束中的 第k个波束,k≥1;zk表示第k个波束所到达的海底深度;当k>1时,zk-1表示第k-1个波束所到达的海底深度;当k=1时,zk-1表示海底深度的初始值;g表示声速梯度;μk表示第k个波束所到达的海底深度的过程噪声;yk表示由第k个波束回波的到达角的余弦值构成的测量变量;tk表示第k个波束回波的到达时间,c0表示多波束测深***的接收换能器阵所在位置的声速;vk表示第k个波束回波的到达角的余弦值的测量噪声。
(4)根据步骤(3)建立的状态空间模型及步骤(2)中接收到的回波的到达角和到达时间,利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值;根据声速梯度的反演值,利用式(6)计算获得测量海域的声速剖面的估计:
c(z)=c0+g(z-d0) (6)
式(6)中,c(z)表示测量海域的声速剖面;z表示海水深度;d0表示多波束测深***的接收换能器阵所在位置的深度;c0表示多波束测深***的接收换能器阵所在位置的声速;g表示声速梯度;
(5)根据步骤(4)的声速剖面的估计值计算测量海域的水温剖面。
进一步地,本发明在所述步骤(4)中,所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。
进一步地,本发明所述步骤(5)中,采用Leroy声速公式、Dell Grosso声速公式、Wilson精确声速公式、Wilson简化声速公式、Mackenzie声速公式、Chen-Millero-Li声速公式或EM分层简化声速公式计算测量海域的水温剖面。
与现有技术相比,本发明的有益效果是:
(1)本发明采用多波束测深***能够同时获得数十甚至上百个窄测深波束,这样就能估计出数十甚至上百个多波束所到达的海底深度;
(2)本发明方法直接利用多波束测深***的测量数据实现测量海域声速剖面估计,进一步可以获得测量海域的水温剖面,相比于直接投放CTD的方法更为快速;
(3)本发明方法可以实现测量海域的海底地形的估计,在改善地形测量精度的同时,提供了一种显著提高环境参量(声速剖面及水温剖面)测量效率的新手段;
(4)本发明采用序贯滤波方法,鲁棒性高,能够消除由于测量噪声和环境因素变化带来的影响,能够准确地反演出测量海域的声速剖面及海底深度。
具体实施方式
下面结合附图和实施例对本发明作进一步具体描述:
如图1所示,目前多波束测深***通常包括:
——发射电路:用于发射电信号。
——发射换能器阵:用于将发射的电信号转化为声信号,并将声信号发射至测量海域。
——接收换能器阵:用于接收发射换能器阵发射的声信号的回波,并将回波声信号转换为电信号。
——接收放大电路:用于将接收换能器阵转换得到的电信号进行放大。
——***控制与信号处理子***:用于实现信号的高速采集和实时处理,并将采集的数据和处理结果通过千兆以太网传给显控与后处理平台,处理结果中包括回波的到达角和到达时间。
——显控与后处理平台:用于接收***控制与信号处理子***的处理结果,处理结果中包括回波的到达角和到达时间。
本发明可通过显控与后处理平台,利用所接收到的回波的到达角和到达时间以及所建立的状态空间模型,实现声速梯度及海底地形的联合反演,并利用声速梯度的反演值计算得到声速剖面的估计;进一步,可利用声速剖面的估计值计算测量海域的水温剖面。
——GPS、姿态仪等***传感器:用于实现母船位置、姿态及航向等的测定。
本发明基于多波束测深的声速剖面及海底地形的联合反演方法的流程如图2所示,包括如下步骤:
(一)通过多波束测深***的发射换能器阵向测量海域的海底发射多波束,通过多波束测深***的接收换能器阵接收回波信号。
(二)多波束测深***根据接收到的回波信号,获得回波的到达角和到达时间。
(三)建立由状态方程和测量方程构成的状态空间模型。其具体方法如下:
(1)状态方程建立
声速剖面用于表示声速随海水深度的变化关系。在声速剖面反演过程中,有效而简洁地表示声速剖面是必不可少的。由于海洋环境复杂多变,难以采用简单的连续函数表示声速剖面。一般可以近似地将声速剖面沿海水深度方向划分为等垂直间隔的一系列声速层,简单的方法是假设层内常声速,更精确的方法则是假设层内声速呈线性变化。当沿海水深度方向划分的层数很多时,以上两种假设都能很好地近似真实声速剖面,但是需要较多的参数。
而本发明采用线性声速剖面,只需要表面声速c0(多波束测深***的接收换能器阵所在位置的声速)和声速梯度g两个参数就可以描述声速剖面,更加简单有效。
声速剖面在测量海域呈线性变化,可表示为如式(7)所示的形式:
c(z)=c0+g(z-d0) (7)
式(7)中,c(z)表示声速剖面;z表示海水深度;d0表示多波束测深***的接收换能器阵所在位置的深度,在接收换能器阵处安装一个深度计即可获得;c0表示多波束测深***的接收换能器阵所在位置的声速,称之为表面声速,表面声速c0在实际中较为容易获得,在接收换能器阵处安装一个声速计即可获得表面声速c0;g表示声速梯度。因此,在采用线性声速剖面条件下,声速剖面可以用表面声速c0和声速梯度g描述。
由于多波束测深***发射的一帧多波束的回波信号经过同一个声速剖面的作用,则可假设在一帧多波束所经历的空间内声速梯度g为常数且为未知的参量。
将声速梯度g和海底深度zk这两个参量组成一个整体构成如式(8)所示的状态变量;此外,考虑海底地形的随机特性,将海底地形建模为如式(9)所示的形式。
zk=zk-1+μk (9)
式(8)和式(9)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,k≥1;zk表示第k个波束所到达的海底深度;当k>1时,zk-1表示第k-1个波束所到达的海底深度;当k=1时,zk-1(即z0)表示海底深度的初始值,它 可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器(例如CTD)测量确定;g表示声速梯度;xk表示声速梯度g和第k个波束所到达的海底深度共同构成的状态变量;μk表示第k个波束所到达的海底深度的过程噪声。
作为本发明的优选方案,过程噪声μk建模为服从均值为零、协方差为φk的高斯分布。在测量海域的海底地形起伏不大的情况下,φk一般取较小的值,按经验来说取为1-2即可;如果欲测量海域的海底地形起伏较大,φk一般取较大的值,按经验来说取为3-5即可。除此之外,过程噪声μk也可以建模为非高斯分布的形式。
建立如式(10)所示的状态方程:
xk=xk-1+wk (10)
式(10)中,xk表示声速梯度g和第k个波束所到达的海底深度共同构成的状态变量,xk-1表示声速梯度g和第k-1个波束所到达的海底深度共同构成的状态变量,wk表示状态变量xk的过程噪声(如式(11)所示)。
与过程噪声μk同理,过程噪声wk既可以建模为非高斯分布的形式,又可以建模为高斯分布。作为本发明的优选方案,过程噪声wk建模为服从均值为零、协方差矩阵为Qk的高斯分布。
式(12)中,在海底地形起伏不大的情况下,协方差φk一般取较小的值,按经验来说取为1-2即可;如果欲测量海域的海底地形起伏较大,φk一般取较大的值,按经验来说取为3-5即可。
由式(8)可知,xk-1可表达为:
综上,由式(8)、(9)、(11)和(13),可将式(10)所示的状态方程表达为如式(14)所示的形式:
式(14)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,k≥1;zk表示第k个波束所到达的海底深度;当k>1时,zk-1表示第k-1个波束所到达的海底深度;当k=1时,zk-1(即z0)表示海底深度的初始值,它可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器(例如CTD)测量确定;g表示声速梯度;μk表示第k个波束所到达的海底深度的过程噪声。
(2)测量方程建立
将回波的到达角的余弦值作为测量变量:
yk=cosθk (15)
式(15)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,θk表示第k个波束的回波的到达角,yk表示由第k个波束回波的到达角的余弦值构成的测量变量。
根据回波的到达角的余弦值与到达时间、声速梯度及海底深度的关系,测量方程可以表示为如式(16)所示的形式:
式(16)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,yk表示由第k个波束回波的到达角的余弦值构成的测量变量;tk表示第k个波束回波的到达时间,zk表示第k个波束所到达的海底深度,g表示声速梯度,c0表示多波束测深***的接收换能器阵所在位置的声速,称之为表面声速;vk表示第k个波束回波的到达角的余弦值的测量噪声。
与过程噪声μk同理,测量噪声vk既可以建模为非高斯分布的形式,又可以建模为高斯分布。作为本发明的优选方案,测量噪声vk建模为服从均值为零、协方差为Rk的高斯分布。
Rk的取值与具体采用的多波束测深***的***参数有关,按经验一般取为1×10-4-1×10-8,Rk的取值越小代表所采用的多波束测深***的角度分辨力越高。
至此,本发明建立了由式(14)所示的状态方程和式(16)所示的测量方程构成的状态空间模型,从而为实现序贯滤波提供了基础。
(四)根据上述建立的由状态方程和测量方程构成的状态空间模型,通过千兆以太网将回波的到达角与到达时间数据传给显控与后处理平台,利用序贯滤波方法实现测量海域的声速梯度及海底深度的反演。
通常情况下,可采用卡尔曼滤波(KF)、扩展卡尔曼滤波(EKF)、集合卡尔曼滤波(EnKF)、无味卡尔曼滤波(UKF)、质点滤波(PF)等序贯滤波方法实现 状态变量的估计。当状态方程和测量方程均是线性形式且过程噪声和测量噪声均服从高斯分布时,则以使用KF方法为较佳的选择;而当状态方程和测量方程中的一个为非线性形式,或者过程噪声和测量噪声中有一个不服从高斯分布时,KF方法不再适用;而EKF、EnKF、UKF或者PF方法在测量方程为线性形式和非线性形式时均可适用。
由于式(16)所示的测量方程是非线性形式的,所以本发明不能采用KF方法实现状态变量的估计,可采用EKF、EnKF、UKF或PF方法实现状态变量的估计。
以下以EKF方法为例来进一步说明本发明进行状态变量估计的过程。
EKF方法主要分为预测和更新两个过程:
(1)预测过程
预测过程中根据式(17)由第k-1个波束的后验状态变量估计值
预测得到第k个波束的先验状态变量估计值
并根据式(18)由第k-1个波束的后验误差协方差估计值
预测得到第k个波束的先验误差协方差估计值
式(17)和式(18)中,
表示第k-1个波束的后验状态变量估计值,
表示由第k-1个波束的后验状态变量估计值
预测出的第k个波束的先验状态变量估计值;
表示第k-1个波束的后验误差协方差估计值,
表示由第k-1个波束的后验误差协方差估计值
预测出的第k个波束的先验误差协方差估计值;Q
k表示过程噪声w
k的协方差矩阵;F
k表示状态过渡矩阵,
T表示矩阵转置操作,由式(14)可知F
k表示为:
(2)更新过程
更新过程中根据式(20)计算第k个波束的卡尔曼增益K
k,然后根据式(21)利用第k个波束的测量变量y
k优化在预测阶段中得到的第k个波束的先验状态变量估计值
获得第k个波束的一个新的更精确的后验状态变量估计值
最后根据式(22)利用第k个波束的先验误差协方差估计值
得到第k个波束的后验误差协方差估计值
Pk|k=(I-KkHk)Pk|k-1 (22)
式(20)-式(22)中,K
k表示卡尔曼增益,P
k|k-1表示第k个波束的先验误差协方差估计值,P
k|k表示第k个波束的后验误差协方差估计值,R
k表示测量噪声v
k的协方差,
-1表示矩阵求逆操作,
T表示矩阵转置操作,I表示单位矩阵;
表示第k个波束的先验状态变估计值,
表示当获得第k个波束的测量变量y
k后优化第k个波束的先验状态变量估计值
获得的一个新的更精确的后验状态变量估计值;H
k表示输出矩阵,由于式(16)所示的测量方程为非线性形式,不能直接得到H
k,利用式(16)对状态变量x
k求偏导,可得到H
k,如式(23)所示:
Hk=[H11,H12]
式(23)中,k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束,tk表示第k个波束回波的到达时间,zk表示第k个波束所到达的海底深度,g表示声速梯度,c0表示多波束测深***的接收换能器阵所在位置的声速,称之为表面声速。
(3)初值选取
使用EKF方法进行递归计算前,首先需要确定状态变量的初值
(即当k=1时,
的取值),表示如下:
式(24)中,
表示状态变量的初值;z
0表示海底深度的初始值,g
0表示声速梯度的初始值;z
0和g
0可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器(例如CTD)测量确定。
状态变量的初值
的选取对EKF算法的收敛速度影响有一定影响。如果初值
与真实值偏离较大,则EKF算法的收敛速度可能会较慢,而当初值的选取与真实值较为接近时,那么算法的收敛速度相对于任意选取时会变快。
通常,误差协方差矩阵的初值P0|0(即当k=1时,Pk-1|k-1的取值)设定为一个对角元素为A的对角矩阵,A是任意设定的一个大于0的值。
根据设定的初值,结合测量数据即多波束回波的到达角和到达时间,利用式(17)-式(18)的预测过程和式(20)-式(22)的更新过程就可以同时反演声速梯度g及海底深度zk,其效果参见图3和图4。由图3和图4可以看出,采用本发明的方法估计出的声速梯度与真实声速梯度十分接近,所估计的海底地形与真实海底地形也十分接近,可见本发明方法估计效果较好。
进一步,根据声速梯度g的反演值,利用式(7)就可以获得声速剖面的估计值c(z),其中,海水深度z的取值从多波束测深***的接收换能器阵所在位置的深度d0到海底深度zk等间隔取值,通常,间隔越小越好。
(五)根据步骤(四)的声速剖面的估计结果计算水温剖面
水温剖面用于表示海水温度随海水深度的变化关系。采用适用于测量海域的声速经验公式,根据步骤(四)中的声速剖面的估计结果,可进一步计算测量海域的水温剖面。声速可表示为海水温度(T)、海水盐度(S)和海水深度(z)的函数,或者表示为海水温度(T)、海水盐度(S)和海水压力(P)的函数。目前被普遍认可的声速经验公式主要有Leroy声速公式、Dell Grosso声速公式、Wilson精确声速公式、Wilson简化声速公式、Mackenzie声速公式、Chen-Millero-Li声速公式和EM分层简化声速公式这7种声速经验公式。其中,EM分层简化声速公式结构简单,计算精度高,且能够适用于海水盐度、海水深度或海水温度变化范围大的测量区域,是本发明的优选实施方式。
以下以EM分层简化声速公式为例来计算水温剖面。
EM分层简化声速公式如式(25)-式(28)所示:
表层声速计算公式:
(25)
淡水中深度达到200m、海水中深度达到1000m时的声速计算公式为:
C(T,z,S)=C(T,0,S)+16.5z (26)
淡水中深度达到2000m、海水中深度达到11000m时的声速计算公式为:
深度大于5000m时,声速计算公式应考虑纬度改正:
(28)
式(25)-式(28)中,T表示海水温度,z表示海水深度,S表示海水盐度,
表示所计算声速位置处的纬度,C(T,z,S)表示当海水温度为T、海水深度为z、海水盐度为S时的声速。
根据步骤(四)中反演得到的声速梯度g,利用声速梯度g的反演值根据式(4)就可以获得声速剖面的估计c(z)。
由步骤(四)中反演得到的海底深度zk即可知测量海域的海底深度,根据式(25)-式(28)的海水深度适用范围,选择其中合适的一个进行水温剖面的估计。从图4中可以看出,本发明中涉及的测量海域的海底深度在380m左右,所以可选取式(26)进行测量海域水温剖面的估计。
设海水盐度S为固定值,海水盐度S可以通过盐度计测量得到。将声速剖面c(z)、不同的海水深度z(海水深度z的取值从多波束测深***的接收换能器阵所在位置的深度d0到海底深度zk等间隔取值,通常,间隔越小越好)和海水盐度S代入式(26),即可计算出从多波束测深***的接收换能器阵所在位置的深度d0到海底深度zk的水温剖面T(参见图5)。