CN102749622B - 基于多波束测深的声速剖面及海底地形的联合反演方法 - Google Patents

Abstract

本发明公开了一种基于多波束测深的声速剖面及海底地形的联合反演方法。它包括以下步骤：（1）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，通过多波束测深***的接收换能器阵接收回波信号；（2）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间；（3）建立由状态方程和测量方程构成的状态空间模型；（4）根据建立的状态空间模型及接收到的回波的到达角和到达时间，利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值，利用声速梯度的反演值获得测量海域的声速剖面的估计。进一步，利用声速剖面的估计值计算测量海域的水温剖面。本发明方法能够快速、准确地获得声速剖面及海底深度的估计。

Description

基于多波束测深的声速剖面及海底地形的联合反演方法

技术领域

本发明涉及的是水声测量技术领域，具体涉及一种基于多波束测深的声速剖面及海底地形的联合反演方法。 

背景技术

海底地形测量是人类对于海洋探索、海洋资源开发、海洋工程设计、海洋战场虚拟环境建设的一个重要手段。长期以来，传统单波束回声测深技术虽然在海底地形测量中发挥了重要作用，但由于其测量所需时间长、成本高、精度低、信息量少的缺陷，已经不能满足海洋经济发展、海洋权益争夺和海洋规划管理日益增长的需求。而多波束测深技术因其测量范围大、精度高、效率高、记录数字化和实时自动绘图的优点，逐渐成为目前一种非常重要且较为流行的海洋测量手段。多波束测深***是一种由多传感器组成的复杂***，是现代电子技术、先进信号处理技术、高分辨率显示技术、高精度导航定位技术以及其他相关高新技术等多种技术的高度集成。 

多波束测深***在测量过程中，利用特殊的发射和接收波束结构，在一个较宽的范围内形成许多较窄的测量波束，从而一次性获得一个较宽范围内水深信息。多波束测深***的发射换能器阵向海底发射纵向（沿航行方向）较窄横向（垂直航行方向）较宽的波束，形成一个测量条带；而接收波束是横向较窄但纵向较宽的波束，与发射波束相交形成许多测点（波束脚印）；最后利用各个脚印回波信号的到达角和到达时间，计算测点对应的深度值，实现海底地形的采样。 

然而，为了保证多波束测深***测深结果的精度，必须消除船在航行时的纵横摇、航偏及声速剖面的影响。其中声速剖面的正确与否，直接关系到测深结果的准确性。海水是不均匀的介质，声波在海水中传播时，声速受海水温度、压力和盐度等因素的影响而发生变化。不同声速在海水中构成一系列声速层，声线在不同声速层中传播时遵循Snell定律，使声线发生折射现象，导致声线弯曲。声线弯曲对多波束测深的影响主要反映在两个方面：（1）对声线传播距离的影响；（2）对水深测量值的影响。这两者导致波束脚印位置计算错误，从而使绘制的海底地形与真实的海底存在差异，形成虚假海底。所以必须要保证 声速剖面的及时更新，否则会影响测深数据的准确性。许多其它声纳***的测量结果也依赖于声速剖面的准确性，所以声速剖面的获取对于海洋研究具有重要的意义。 

目前获得声速剖面的方式主要有两种，即测量方法与反演方法。测量方法包括直接测量法和间接测量法两种，直接测量法利用声速仪，通过在固定的已知距离内发射与接收信号测量出声速值；间接测量法利用温盐深仪（CTD）测量海水的盐度、温度和深度，再根据经验公式计算出声速值。利用逐点测量的方法获取声速剖面往往需要母船停泊才能进行测量，需投入较多人力物力。而实际海洋环境条件复杂多变，江河入海口处淡水的流入、水团的运动及昼夜温差等因素，都使得海水声速剖面随时间和空间变化，这样就需要及时更新声速剖面，而逐点测量对于获取较大范围海域的声速剖面需要较长的时间，直接影响了多波束测深***的工作效率。 

近年来利用水声学方法反演声速剖面并利用声速与水文参数的关系反演相关水文参数的研究备受关注。用水声学方法通常可以同时反演出大范围海域的环境参数，并且可以实现长时间的连续监测，与传统方法相比有很大优势，提供了海洋环境测量一种比较方便、高效的途径。 

寻找能够准确、快速地获得随时间和空间演化的声速剖面的反演方法因而成为目前研究的重要点。 

发明内容

本发明的目的在于针对现有技术的不足，提供了一种基于多波束测深的声速剖面及海底地形的联合反演方法和基于多波束测深的水温剖面的反演方法。 

为实现上述目的，本发明所采取的其中一种技术方案是：本发明基于多波束测深的声速剖面及海底地形的联合反演方法包括如下步骤： 

（1）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，并通过多波束测深***的接收换能器阵接收回波信号； 

（2）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间； 

（3）建立由式（1）所示的状态方程和式（2）所示的测量方程构成的状态空间模型： 

$\left[\begin{array}{c}{z}_k\\ g\end{array}\right]=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]+\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(1\right)$

$y_k=\frac{z_k}{t_k}[\frac{1}{c_0}+g{z}_k(\frac{{t_k}^2}{{{2z}_k}^2}-\frac{1}{c_0^2})+{\left({\mathrm{gz}}_k\right)}^2(\frac{1}{c_0^3}-\frac{{t_k}^2}{6c_0{z_k}^2})]+v_k---\left(2\right)$

式（1）和式（2）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k＞1时，z_k-1表示第k-1个波束所到达的海底深度；当k=1时，z_k-1表示海底深度的初始值；g表示声速梯度；μ_k表示第k个波束所到达的海底深度的过程噪声；y_k表示由第k个波束回波的到达角的余弦值构成的测量变量；t_k表示第k个波束回波的到达时间，c₀表示多波束测深***的接收换能器阵所在位置的声速；v_k表示第k个波束回波的到达角的余弦值的测量噪声。 

（4）根据步骤（3）建立的状态空间模型及步骤（2）中接收到的回波的到达角和到达时间，利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值；根据声速梯度的反演值，利用式（3）计算获得测量海域的声速剖面的估计： 

c(z)=c₀+g(z-d₀)                    （3） 

式（3）中，c(z)表示测量海域的声速剖面；z表示海水深度；d₀表示多波束测深***的接收换能器阵所在位置的深度；c₀表示多波束测深***的接收换能器阵所在位置的声速；g表示声速梯度。 

进一步地，本发明在所述步骤（4）中，所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。 

本发明所采取的另一种技术方案是：本发明基于多波束测深的水温剖面的反演方法包括如下步骤： 

（1）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，通过多波束测深***的接收换能器阵接收回波信号； 

（2）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间； 

（3）建立由式（4）所示的状态方程和式（5）所示的测量方程构成的状态空间模型： 

$\left[\begin{array}{c}{z}_k\\ g\end{array}\right]=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]+\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(4\right)$

$y_k=\frac{z_k}{t_k}[\frac{1}{c_0}+g{z}_k(\frac{{t_k}^2}{{{2z}_k}^2}-\frac{1}{c_0^2})+{\left({\mathrm{gz}}_k\right)}^2(\frac{1}{c_0^3}-\frac{{t_k}^2}{6c_0{z_k}^2})]+v_k---\left(5\right)$

式（4）和式（5）中，k表示发射换能器阵向测量海域的海底发射的多波束中的 第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k＞1时，z_k-1表示第k-1个波束所到达的海底深度；当k＝1时，z_k-1表示海底深度的初始值；g表示声速梯度；μ_k表示第k个波束所到达的海底深度的过程噪声；y_k表示由第k个波束回波的到达角的余弦值构成的测量变量；t_k表示第k个波束回波的到达时间，c₀表示多波束测深***的接收换能器阵所在位置的声速；v_k表示第k个波束回波的到达角的余弦值的测量噪声。 

（4）根据步骤（3）建立的状态空间模型及步骤（2）中接收到的回波的到达角和到达时间，利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值；根据声速梯度的反演值，利用式（6）计算获得测量海域的声速剖面的估计： 

c(z)=c₀+g(z-d₀)                    （6） 

式（6）中，c(z)表示测量海域的声速剖面；z表示海水深度；d₀表示多波束测深***的接收换能器阵所在位置的深度；c₀表示多波束测深***的接收换能器阵所在位置的声速；g表示声速梯度； 

（5）根据步骤（4）的声速剖面的估计值计算测量海域的水温剖面。 

进一步地，本发明在所述步骤（4）中，所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。 

进一步地，本发明所述步骤（5）中，采用Leroy声速公式、Dell Grosso声速公式、Wilson精确声速公式、Wilson简化声速公式、Mackenzie声速公式、Chen-Millero-Li声速公式或EM分层简化声速公式计算测量海域的水温剖面。 

与现有技术相比，本发明的有益效果是： 

（1）本发明采用多波束测深***能够同时获得数十甚至上百个窄测深波束，这样就能估计出数十甚至上百个多波束所到达的海底深度； 

（2）本发明方法直接利用多波束测深***的测量数据实现测量海域声速剖面估计，进一步可以获得测量海域的水温剖面，相比于直接投放CTD的方法更为快速； 

（3）本发明方法可以实现测量海域的海底地形的估计，在改善地形测量精度的同时，提供了一种显著提高环境参量（声速剖面及水温剖面）测量效率的新手段； 

（4）本发明采用序贯滤波方法，鲁棒性高，能够消除由于测量噪声和环境因素变化带来的影响，能够准确地反演出测量海域的声速剖面及海底深度。 

附图说明

图1是多波束测深***的结构框图； 

图2是本发明联合反演方法的流程示意图； 

图3是利用本发明联合反演方法所获得的声速梯度估计结果与真实值的比较图； 

图4是利用本发明联合反演方法所获得的海底深度估计结果与真实值的比较图； 

图5是利用本发明联合反演方法所获得的水温剖面估计结果。 

具体实施方式

下面结合附图和实施例对本发明作进一步具体描述： 

如图1所示，目前多波束测深***通常包括： 

——发射电路：用于发射电信号。 

——发射换能器阵：用于将发射的电信号转化为声信号，并将声信号发射至测量海域。 

——接收换能器阵：用于接收发射换能器阵发射的声信号的回波，并将回波声信号转换为电信号。 

——接收放大电路：用于将接收换能器阵转换得到的电信号进行放大。 

——***控制与信号处理子***：用于实现信号的高速采集和实时处理，并将采集的数据和处理结果通过千兆以太网传给显控与后处理平台，处理结果中包括回波的到达角和到达时间。 

——显控与后处理平台：用于接收***控制与信号处理子***的处理结果，处理结果中包括回波的到达角和到达时间。 

本发明可通过显控与后处理平台，利用所接收到的回波的到达角和到达时间以及所建立的状态空间模型，实现声速梯度及海底地形的联合反演，并利用声速梯度的反演值计算得到声速剖面的估计；进一步，可利用声速剖面的估计值计算测量海域的水温剖面。 

——GPS、姿态仪等***传感器：用于实现母船位置、姿态及航向等的测定。 

本发明基于多波束测深的声速剖面及海底地形的联合反演方法的流程如图2所示，包括如下步骤： 

（一）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，通过多波束测深***的接收换能器阵接收回波信号。 

（二）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间。 

（三）建立由状态方程和测量方程构成的状态空间模型。其具体方法如下： 

（1）状态方程建立 

声速剖面用于表示声速随海水深度的变化关系。在声速剖面反演过程中，有效而简洁地表示声速剖面是必不可少的。由于海洋环境复杂多变，难以采用简单的连续函数表示声速剖面。一般可以近似地将声速剖面沿海水深度方向划分为等垂直间隔的一系列声速层，简单的方法是假设层内常声速，更精确的方法则是假设层内声速呈线性变化。当沿海水深度方向划分的层数很多时，以上两种假设都能很好地近似真实声速剖面，但是需要较多的参数。 

而本发明采用线性声速剖面，只需要表面声速c₀（多波束测深***的接收换能器阵所在位置的声速）和声速梯度g两个参数就可以描述声速剖面，更加简单有效。 

声速剖面在测量海域呈线性变化，可表示为如式（7）所示的形式： 

c(z)=c₀+g(z-d₀)                （7） 

式（7）中，c(z)表示声速剖面；z表示海水深度；d₀表示多波束测深***的接收换能器阵所在位置的深度，在接收换能器阵处安装一个深度计即可获得；c₀表示多波束测深***的接收换能器阵所在位置的声速，称之为表面声速，表面声速c₀在实际中较为容易获得，在接收换能器阵处安装一个声速计即可获得表面声速c₀；g表示声速梯度。因此，在采用线性声速剖面条件下，声速剖面可以用表面声速c₀和声速梯度g描述。 

由于多波束测深***发射的一帧多波束的回波信号经过同一个声速剖面的作用，则可假设在一帧多波束所经历的空间内声速梯度g为常数且为未知的参量。 

将声速梯度g和海底深度z_k这两个参量组成一个整体构成如式（8）所示的状态变量；此外，考虑海底地形的随机特性，将海底地形建模为如式（9）所示的形式。 

$x_k=\left[\begin{array}{c}{z}_k\\ g\end{array}\right]---\left(8\right)$

z_k=z_k-1+μ_k                            （9） 

式（8）和式（9）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k＞1时，z_k-1表示第k-1个波束所到达的海底深度；当k=1时，z_k-1（即z₀）表示海底深度的初始值，它 可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器（例如CTD）测量确定；g表示声速梯度；x_k表示声速梯度g和第k个波束所到达的海底深度共同构成的状态变量；μ_k表示第k个波束所到达的海底深度的过程噪声。 

作为本发明的优选方案，过程噪声μ_k建模为服从均值为零、协方差为φ_k的高斯分布。在测量海域的海底地形起伏不大的情况下，φ_k一般取较小的值，按经验来说取为1-2即可；如果欲测量海域的海底地形起伏较大，φ_k一般取较大的值，按经验来说取为3-5即可。除此之外，过程噪声μ_k也可以建模为非高斯分布的形式。 

建立如式（10）所示的状态方程： 

x_k=x_k-1+w_k                            （10） 

式（10）中，x_k表示声速梯度g和第k个波束所到达的海底深度共同构成的状态变量，x_k-1表示声速梯度g和第k-1个波束所到达的海底深度共同构成的状态变量，w_k表示状态变量x_k的过程噪声（如式（11）所示）。 

$w_k=\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(11\right)$

与过程噪声μ_k同理，过程噪声w_k既可以建模为非高斯分布的形式，又可以建模为高斯分布。作为本发明的优选方案，过程噪声w_k建模为服从均值为零、协方差矩阵为Q_k的高斯分布。 

$Q_k=\left[\begin{array}{cc}{\mathrm{\φ}}_k& 0\\ 0& 0\end{array}\right]---\left(12\right)$

式（12）中，在海底地形起伏不大的情况下，协方差φ_k一般取较小的值，按经验来说取为1-2即可；如果欲测量海域的海底地形起伏较大，φ_k一般取较大的值，按经验来说取为3-5即可。 

由式（8）可知，x_k-1可表达为： 

$x_{k-1}=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]---\left(13\right)$

综上，由式（8）、（9）、（11）和（13），可将式（10）所示的状态方程表达为如式（14）所示的形式： 

$\left[\begin{array}{c}{z}_k\\ g\end{array}\right]=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]+\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(14\right)$

式（14）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k＞1时，z_k-1表示第k-1个波束所到达的海底深度；当k＝1时，z_k-1（即z₀）表示海底深度的初始值，它可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器（例如CTD）测量确定；g表示声速梯度；μ_k表示第k个波束所到达的海底深度的过程噪声。 

（2）测量方程建立 

将回波的到达角的余弦值作为测量变量： 

y_k=cosθ_k                    （15） 

式（15）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，θ_k表示第k个波束的回波的到达角，y_k表示由第k个波束回波的到达角的余弦值构成的测量变量。 

根据回波的到达角的余弦值与到达时间、声速梯度及海底深度的关系，测量方程可以表示为如式（16）所示的形式： 

$y_k=\frac{z_k}{t_k}[\frac{1}{c_0}+g{z}_k(\frac{{t_k}^2}{{{2z}_k}^2}-\frac{1}{c_0^2})+{\left({\mathrm{gz}}_k\right)}^2(\frac{1}{c_0^3}-\frac{{t_k}^2}{6c_0{z_k}^2})]+v_k---\left(16\right)$

式（16）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，y_k表示由第k个波束回波的到达角的余弦值构成的测量变量；t_k表示第k个波束回波的到达时间，z_k表示第k个波束所到达的海底深度，g表示声速梯度，c₀表示多波束测深***的接收换能器阵所在位置的声速，称之为表面声速；v_k表示第k个波束回波的到达角的余弦值的测量噪声。 

与过程噪声μ_k同理，测量噪声v_k既可以建模为非高斯分布的形式，又可以建模为高斯分布。作为本发明的优选方案，测量噪声v_k建模为服从均值为零、协方差为R_k的高斯分布。 

R_k的取值与具体采用的多波束测深***的***参数有关，按经验一般取为1×10^-4-1×10^-8，R_k的取值越小代表所采用的多波束测深***的角度分辨力越高。 

至此，本发明建立了由式（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个波束的先验误差协方差估计值 

${\hat{x}}_{k|k-1}=F_k{\hat{x}}_{k-1|k-1}---\left(17\right)$

$P_{k|k-1}=F_k{P}_{k-1|k-1}{F}_k^T+Q_k---\left(18\right)$

式（17）和式（18）中， 

表示第k-1个波束的后验状态变量估计值， 表示由第k-1个波束的后验状态变量估计值 预测出的第k个波束的先验状态变量估计值； 表示第k-1个波束的后验误差协方差估计值， 

表示由第k-1个波束的后验误差协方差估计值 

预测出的第k个波束的先验误差协方差估计值；Q_k表示过程噪声w_k的协方差矩阵；F_k表示状态过渡矩阵，^T表示矩阵转置操作，由式（14）可知F_k表示为： 

$F_k=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]---\left(19\right)$

（2）更新过程 

更新过程中根据式（20）计算第k个波束的卡尔曼增益K_k，然后根据式（21）利用第k个波束的测量变量y_k优化在预测阶段中得到的第k个波束的先验状态变量估计值 

获得第k个波束的一个新的更精确的后验状态变量估计值 

最后根据式（22）利用第k个波束的先验误差协方差估计值 

得到第k个波束的后验误差协方差估计值 

$K_k=P_{k|k-1}{H}_k^T{(H_k{P}_{k|k-1}{H}_k^T+R_k)}^{-1}---\left(20\right)$

${\hat{x}}_{k|k}={\hat{x}}_{k|k-1}+K_k(y_k-H_k{\hat{x}}_{k|k-1})---\left(21\right)$

P_k|k=(I-K_kH_k)P_k|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）所示： 

H_k=[H₁₁,H₁₂] 

$H_{11}\≈\frac{1}{t_k{c}_0}-\frac{2{\mathrm{gz}}_k}{t_k{c}_0^2}+\frac{3{\left({\mathrm{gz}}_k\right)}^2}{t_k{c}_0^3}-\frac{g^2{t}_k}{6c_0}---\left(23\right)$

$H_{12}\≈\frac{t_k}{2}-\frac{z_k^2}{t_k{c}_0^2}+\frac{2{\mathrm{gz}}_k^3}{t_k{c}_0^3}-\frac{{\mathrm{gz}}_k{t}_k}{{3c}_0}$

式（23）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，t_k表示第k个波束回波的到达时间，z_k表示第k个波束所到达的海底深度，g表示声速梯度，c₀表示多波束测深***的接收换能器阵所在位置的声速，称之为表面声速。 

（3）初值选取 

使用EKF方法进行递归计算前，首先需要确定状态变量的初值 （即当k＝1时， 

的取值），表示如下： 

${\hat{x}}_{0|0}=\left[\begin{array}{c}{z}_0\\ g_0\end{array}\right]---\left(24\right)$

式（24）中， 

表示状态变量的初值；z₀表示海底深度的初始值，g₀表示声速梯度的初始值；z₀和g₀可以是任意设定的一个大于0的值,也可以由多波束测深***的***传感器（例如CTD）测量确定。 

状态变量的初值 的选取对EKF算法的收敛速度影响有一定影响。如果初值 

与真实值偏离较大，则EKF算法的收敛速度可能会较慢，而当初值的选取与真实值较为接近时，那么算法的收敛速度相对于任意选取时会变快。 

通常，误差协方差矩阵的初值P_0|0（即当k=1时，P_k-1|k-1的取值）设定为一个对角元素为A的对角矩阵，A是任意设定的一个大于0的值。 

根据设定的初值，结合测量数据即多波束回波的到达角和到达时间，利用式（17）-式（18）的预测过程和式（20）-式（22）的更新过程就可以同时反演声速梯度g及海底深度z_k，其效果参见图3和图4。由图3和图4可以看出，采用本发明的方法估计出的声速梯度与真实声速梯度十分接近，所估计的海底地形与真实海底地形也十分接近，可见本发明方法估计效果较好。 

进一步，根据声速梯度g的反演值，利用式（7）就可以获得声速剖面的估计值c(z)，其中，海水深度z的取值从多波束测深***的接收换能器阵所在位置的深度d₀到海底深度z_k等间隔取值，通常，间隔越小越好。 

（五）根据步骤（四）的声速剖面的估计结果计算水温剖面 

水温剖面用于表示海水温度随海水深度的变化关系。采用适用于测量海域的声速经验公式，根据步骤（四）中的声速剖面的估计结果，可进一步计算测量海域的水温剖面。声速可表示为海水温度（T）、海水盐度（S）和海水深度（z）的函数，或者表示为海水温度（T）、海水盐度（S）和海水压力（P）的函数。目前被普遍认可的声速经验公式主要有Leroy声速公式、Dell Grosso声速公式、Wilson精确声速公式、Wilson简化声速公式、Mackenzie声速公式、Chen-Millero-Li声速公式和EM分层简化声速公式这7种声速经验公式。其中，EM分层简化声速公式结构简单，计算精度高，且能够适用于海水盐度、海水深度或海水温度变化范围大的测量区域，是本发明的优选实施方式。 

以下以EM分层简化声速公式为例来计算水温剖面。 

EM分层简化声速公式如式（25）-式（28）所示： 

表层声速计算公式： 

$C(T,0,S)=1449.05+T[4.57-T(0.0521-0.00023T)]$ （25） 

$+[1.333-T(0.0126-0.00009T)](S-35)$

淡水中深度达到200m、海水中深度达到1000m时的声速计算公式为： 

C(T,z,S)=C(T,0,S)+16.5z            （26） 

淡水中深度达到2000m、海水中深度达到11000m时的声速计算公式为： 

$C(T,z,S)=C(T,0,S)+z[16.3+z(0.22-0.003z\sqrt{T+2})]---\left(27\right)$

深度大于5000m时，声速计算公式应考虑纬度改正： 

$C(T,z,S)=C(T,0,S)+z^{\′}[16.3+z^{\′}(0.22-0.003z^{\′}\sqrt{T+2})]$ （28） 

式（25）-式（28）中，T表示海水温度，z表示海水深度，S表示海水盐度， 

表示所计算声速位置处的纬度，C(T,z,S)表示当海水温度为T、海水深度为z、海水盐度为S时的声速。 

根据步骤（四）中反演得到的声速梯度g，利用声速梯度g的反演值根据式（4）就可以获得声速剖面的估计c(z)。 

由步骤（四）中反演得到的海底深度z_k即可知测量海域的海底深度，根据式（25）-式（28）的海水深度适用范围，选择其中合适的一个进行水温剖面的估计。从图4中可以看出，本发明中涉及的测量海域的海底深度在380m左右，所以可选取式（26）进行测量海域水温剖面的估计。 

设海水盐度S为固定值，海水盐度S可以通过盐度计测量得到。将声速剖面c(z)、不同的海水深度z（海水深度z的取值从多波束测深***的接收换能器阵所在位置的深度d₀到海底深度z_k等间隔取值，通常，间隔越小越好）和海水盐度S代入式（26），即可计算出从多波束测深***的接收换能器阵所在位置的深度d₀到海底深度z_k的水温剖面T（参见图5）。 

Claims (5)

1.一种基于多波束测深的声速剖面及海底地形的联合反演方法，其特征是，包括如下步骤：

（1）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，并通过多波束测深***的接收换能器阵接收回波信号；

（2）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间；

（3）建立由式（a）所示的状态方程和式（b）所示的测量方程构成的状态空间模型：

$\left[\begin{array}{c}{z}_k\\ g\end{array}\right]=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]+\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(a\right)$

$y_k=\frac{z_k}{t_k}[\frac{1}{c_0}+{\mathrm{gz}}_k(\frac{{t_k}^2}{{{2z}_k}^2}-\frac{1}{c_0^2})+{\left({\mathrm{gz}}_k\right)}^2(\frac{1}{c_0^3}-\frac{{t_k}^2}{{6c}_0{z_k}^2})]+v_k---\left(b\right)$

式（a）和式（b）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k>1时，z_k-1表示第k-1个波束所到达的海底深度；当k=1时，z_k-1表示海底深度的初始值；g表示声速梯度；μ_k表示第k个波束所到达的海底深度的过程噪声；y_k表示由第k个波束回波的到达角的余弦值构成的测量变量；t_k表示第k个波束回波的到达时间，c₀表示多波束测深***的接收换能器阵所在位置的声速；v_k表示第k个波束回波的到达角的余弦值的测量噪声；

（4）根据步骤（3）建立的状态空间模型以及步骤（2）中接收到的回波的到达角和到达时间，利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值；根据声速梯度的反演值，利用式（c）获得测量海域的声速剖面的估计：

c(z)=c₀+g(z-d₀)    （c）

式（c）中，c(z)表示测量海域的声速剖面；z表示海水深度；d₀表示多波束测深***的接收换能器阵所在位置的深度；c₀表示多波束测深***的接收换能器阵所在位置的声速；g表示声速梯度。

2.根据权利要求1所述的基于多波束测深的声速剖面及海底地形的联合反演方法，其特征是：在所述步骤（4）中，所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。

3.一种基于多波束测深的水温剖面的反演方法，其特征是，包括如下步骤：

（1）通过多波束测深***的发射换能器阵向测量海域的海底发射多波束，并通过多波束测深***的接收换能器阵接收回波信号；

（2）多波束测深***根据接收到的回波信号，获得回波的到达角和到达时间；

（3）建立由式（a）所示的状态方程和式（b）所示的测量方程构成的状态空间模型：

$\left[\begin{array}{c}{z}_k\\ g\end{array}\right]=\left[\begin{array}{c}{z}_{k-1}\\ g\end{array}\right]+\left[\begin{array}{c}{\mathrm{\μ}}_k\\ 0\end{array}\right]---\left(a\right)$

$y_k=\frac{z_k}{t_k}[\frac{1}{c_0}+{\mathrm{gz}}_k(\frac{{t_k}^2}{{{2z}_k}^2}-\frac{1}{c_0^2})+{\left({\mathrm{gz}}_k\right)}^2(\frac{1}{c_0^3}-\frac{{t_k}^2}{{6c}_0{z_k}^2})]+v_k---\left(b\right)$

式（a）和式（b）中，k表示发射换能器阵向测量海域的海底发射的多波束中的第k个波束，k≥1；z_k表示第k个波束所到达的海底深度；当k>1时，z_k-1表示第k-1个波束所到达的海底深度；当k=1时，z_k-1表示海底深度的初始值；g表示声速梯度；μ_k表示第k个波束所到达的海底深度的过程噪声；y_k表示由第k个波束回波的到达角的余弦值构成的测量变量；t_k表示第k个波束回波的到达时间，c₀表示多波束测深***的接收换能器阵所在位置的声速；v_k表示第k个波束回波的到达角的余弦值的测量噪声；

（4）根据步骤（3）建立的状态空间模型以及步骤（2）中接收到的回波的到达角和到达时间，利用序贯滤波方法获得测量海域的声速梯度及海底深度的反演值；根据声速梯度的反演值，利用式（c）计算获得测量海域的声速剖面的估计：

c(z)=c₀+g(z-d₀)    （c）

式（c）中，c(z)表示测量海域的声速剖面；z表示海水深度；d₀表示多波束测深***的接收换能器阵所在位置的深度；c₀表示多波束测深***的接收换能器阵所在位置的声速；g表示声速梯度；

（5）根据步骤（4）的声速剖面的估计值计算测量海域的水温剖面。

4.根据权利要求3所述的基于多波束测深的水温剖面的反演方法，其特征是：在所述步骤（4）中，所述序贯滤波方法为扩展卡尔曼滤波方法、集合卡尔曼滤波方法、无味卡尔曼滤波方法或质点滤波方法。

5.根据权利要求3或4所述的基于多波束测深的水温剖面的反演方法，其特征是：所述步骤（5）中，采用Leroy声速公式、Dell Grosso声速公式、Wilson精确声速公式、Wilson简化声速公式、Mackenzie声速公式、Chen-Millero-Li声速公式或EM分层简化声速公式计算测量海域的水温剖面。

Images