深海多波束声线精确跟踪方法与流程

本发明涉及声线精确跟踪方法，特别涉及深海多波束声线精确跟踪方法。

背景技术：

众所周知，水下目标探测是海底资源调查研究、海洋环境开发和海洋工程设计的工作基础，而多波束测深作为目前水下目标高精度探测的主要手段，可借助其实现海底地形地貌高精度探测、海底资源精准探测开发、水下基准点布设和观测、水下辅助导航、海底滑坡及垮塌等地质灾害的预警、跨海大桥建设辅助勘测等。目前，影响多波束测深精度的主要因素为潮位变化、换能器吃水变化、涌浪变化及声线跟踪。其中潮位、换能器吃水改正和涌浪变化可借助预报或实测数据直接进行深度改正，而声线跟踪与波束初始入射角、声速剖面、跟踪模型密切相关，因此，声线跟踪成为多波束测深数据处理工作中最关键的步骤。

多波束声线跟踪的三个影响因素共同制约着多波束测深精度，具体体现为：

声线精确跟踪首要的前提条件是波束初始入射角的准确计算，而波束初始入射角决定着声线实际传播剖面，故影响到后期在声速剖面中的声线跟踪，最终影响到波束海底投射点坐标的计算精度。

由于水体声速主要受温度、盐度和压力三个因素影响，因此，在不同水团中的声速具有时空变化的特点。而在数据处理过程中，常用实测点位的声速剖面或根据声速经验公式计算的存在较大误差的声速剖面替代未实测点位的声速剖面，常会造成海底地形变形(两边上翘或下塌)(图1)，边缘水深值存在约15％的误差而难以被接受^[1]，特殊情况下，其影响会更大。

声线跟踪模型构建作为声线跟踪的重点工作，目前普遍认为层内常梯度声线跟踪模型最能接近声线实际传播轨迹，但在计算每个水柱的旅行时间中，用垂直方向的平均声速替代了真实弧线上的平均声速，导致声线在每个水柱的旅行时间存在误差，从而影响声线整个旅行时间的精确计算，最终导致波束的归位计算存在较大误差，尤其在边缘波束，如图2所示，左、右图分别为声速剖面改正前后的海底地形，改正前图像存在明显的皱褶，而改正后除了边缘波束外中间区域地形较平滑。但由于声线跟踪不够精准，导致边缘波束的皱褶仍然消除不掉^[2]。其次，在数千米的深海声线跟踪过程中，由于原始声速剖面文件的水柱层数过多，声线跟踪耗时太长，严重影响了波束归位计算效率。

综上所述，传统的声线跟踪方法常造成多波束测深精度低下，尤其在数千米的深海地区，传统方法常给多波束最终测量成果带来非常显著的影响，甚至常出现超过水深值5％的垂直误差而难以被接受。本发明在全面分析传统多波束声线跟踪模型缺陷的基础上，提出了一种精确反映波束实际传播轨迹的思想和算法，最大限度地削弱了上述因素的影响，研究成果将较大地提高深海地区多波束测深精度，使之能够广泛应用于高精度的深海地形地貌测量工程中。

多波束测深精度在较大程度上取决于实际声线的跟踪，而声线的精确跟踪又分别与波束的初始入射角、声速剖面和声线跟踪模型三者具有密切关系，下面分别就以上三者的研究现状进行分析。

在波束初始入射角计算方面，国内外文献极少详细阐述，目前多波束初始入射角的获取方式有两种方式：1)直接取自换能器提供的各个波束分配角^[3]；2)仅顾及了横摇角，认为波束初始入射角是多波束分配角与横摇角之和^[4]。以上两种初始入射角获取方法均基于一个假设，即多波束Ping测深断面与测量船航迹方向正交，而实际多波束换能器在横摇、纵摇姿态角作用下，Ping测深断面与航迹方向已不再正交，而前者未顾及这一影响，后者虽顾及了横摇的影响，但忽视了纵摇的影响；此外，在对姿态影响的处理方面，两种方法均先通过声线跟踪获得波束点坐标，尔后再借助姿态角构建旋转矩阵进行强制旋转变换，获得波束海底测深点在理想船体坐标系下的坐标。不同的是，前者借助纵摇、横摇角构建的旋转矩阵来实施变换，而后者仅借助纵摇角构建的旋转矩阵进行变化。

在声速剖面构建方面，主要涉及声速和水层界面两个因素。水声速是一个受海洋环境影响的复杂量，它是由水体温度、盐度、静压力及空气泡含量等共同决定，并且具有时空变化特性。水层界面即为垂直声速变化较快的界面，当两个相邻水层的声速梯度大于一定的边界值，则认为这两个水层为不同的水层；反之，可将相邻的两个水层合并为同一个水层。最早的声速测量起步于1827年，Colladon与Strum在日内瓦湖上进行了人类历史上第一次声速测量，虽然测量的结果与Laplace理论预测的结果基本一致，但是由于实验过程较简单，只能测得湖水的平均声速，而无法获取声速剖面^[5]。Munk和Wunsch最早提出声层析方法是微扰法，其思想为使用一个假定的声速背景模型，声速相对背景模型的偏移和由背景模型所计算出的声波旅行时间与其测量值的偏差成比例，而所采用的传播模型是基于射线理论的^[6]。基于之前的研究基础，这两位学者又提出了基于Abel变换的非扰动方法^[7]，之后这个方法又被Jones推广到非对称领域，阐述了声道与距离变化的关系^[8-10]。加拿大学者Dinn和法国学者Helene等通过大量的实验数据研究验证了声速误差主要对多波束系统的波束指向角和波束路径产生较大影响，同时给出了具体的量化分析结果^[11]。Gonharov等运用匹配场处理方法将挪威海声层析实验获取的数据反演声速剖面，反演得到的结果与实测的声速剖面一致，验证了该反演方法的有效性^[12]。唐俊峰基于射线理论构建三维声场模型，并运用经验正交函数表示声速剖面，构造出声速剖面反演的代价函数，同时将遗传算法应用于声速剖面反演计算中，通过实验证明，反演的声速剖面与实验测量的结果基本一致。但是在复杂的海洋环境中，声速也会随深度变化存在明显的水平梯度，在唐俊峰的声速剖面反演研究中并未考虑到声速水平梯度变化的情况，以及抑制在此情况下可能引入的传播时间误差，这些误差也会降低多波束的测深精度^[13]。丁继胜等提出一种基于实测声速剖面的EOF表示方法^[14]，该方法利用经验正交函数建立了声速剖面场的数学模型，获得了实测区任意点位声速剖面，尤其提高了边缘波束脚印跟踪的精度。

在声线跟踪模型方面，Hamid等(2013)利用高斯牛顿和扩展卡尔曼滤波解提高了常梯度声速剖面中声线节点定位和跟踪精度，相比基于直线波传播的跟踪模型精度较高^[15]。吴德明(1992)研究了一种应用于长基线直角坐标阵形双曲定位系统的声线修正迭代法^[16]，其基本思想为声速在分层中呈现等梯度分布，并由迭代法求出合理的声线和定位点。之后，李迎春和吴德明分别通过内插和差分方程获取更合理的声线和修正量，从而使得该迭代方法得到更好的改进^[17]。该思想的局限性是在定位模型不变的前提下通过平均声速、近似函数或迭代进行逼近，而忽略了实际声速是随水深变化的现实，从而产生了较大的测深误差。姜薇(2005)提出一种三维声线追踪的正三棱锥前向伸展算法^[18]，通过实验证明该方法具有较高的运算速度和精确度，可应用于三维超声层析图像重建的散射修正方面。兰华林(2007)针对深海应答器定位导航系统中的声线弯曲问题进行了较全面的分析，并提出了基于距离误差最小的球面交汇迭代修正方法^[19]。孙万卿(2007)提出一种基于有穷状态自动机的二维的浅海声线跟踪方法^[20]，解决了声线的反射、折射、全反射等问题并快速地找出声线的路径，并与特定声速剖面情况的数值解和解析解进行比较，验证了此方法的可行性，数值解精度达到10-4。陆秀平(2012)针对常梯度声线跟踪方法中的层内平均声速取自声线圆弧对应弦段的平均声速问题，认为传统方法不够严密，故提出一种沿声线传播圆弧路径进行积分求取层内平均声速的改进方法^[21]，实验表明此方法在深海多波束测深中，可显著提高边缘波束脚印坐标的计算精度。

由此可见，现有技术存在下列技术缺陷：

(1)目前的温盐深场模型分辨率较低，并未在声线跟踪中得到较好的应用；

(2)波束初始入射角计算未考虑船体姿态的影响，从而导致波束脚印坐标计算精度较低；

(3)在数据处理过程中，常用实测点位的声速剖面或根据声速经验公式计算的存在较大误差的声速剖面替代未实测点位的声速剖面，常会造成海底地形变形(两边上翘或下塌)，边缘水深值存在约15％的误差而难以被接受。

(4)传统的层内常梯度声线跟踪模型在计算每个水柱的旅行时间中，用垂直方向的平均声速替代了真实弧线上的平均声速，导致声线在每个水柱的旅行时间存在误差，从而影响声线整个旅行时间的精确计算，最终导致波束的归位计算存在较大误差。其次，在数千米的深海声线跟踪过程中，未过滤合并梯度较小的水柱层，导致声线跟踪耗时太长，严重影响了波束归位计算效率。

技术实现要素：

本发明要解决的技术问题是在全面分析水声信道传播特性、综合多源海洋物理水文观测数据的基础上，建立一种深海多波束声线精确跟踪方法，以满足深海多波束高精度测深的要求。

为达到上述目的，本发明的技术方案如下：

深海多波束声线精确跟踪方法，主要包括以下步骤：

(1)在全面分析海洋表层温度、盐度易受各种因素影响的基础上，并综合海洋卫星和Argo浮标多源海洋物理水文观测数据，建立一种时空变化的海洋温盐场模型；

(2)通过分析船体瞬时姿态对波束入射角的影响，顾及船体姿态精确计算每个波束初始入射角；

(3)在声速剖面未实测区域，分别通过时空温盐场模型计算声速和基于经验正交函数反演声速，在相应点位取两种方法获取声速的均值作为该点位的声速，以此计算三维声速剖面模型；

(4)构建一种高效的常梯度声线精确跟踪模型；

(5)提供声线跟踪精度评估方法。

在本发明的一个实施例中，建立一种时空变化的海洋温盐场数学模型，主要包括以下步骤：

1)通过海洋卫星监测获取海洋次表层物理水文观测数据；

2)通过Argo浮标获取测量剖面的物理水文观测数据；

3)在综合海洋卫星和Argo浮标多源观测数据的基础上，运用海洋动力学数值模拟模型(FVCOM)建立时空温盐场模型；

4)在关键的网格节点上运用CTD进行实测，获取节点的实测温盐深数据；

5)最后对建立的时空温盐场模型进行内外符合精度计算。

在本发明的一个实施例中，计算三维声速剖面模型，主要包括以下步骤：

1)在测区未实施声速剖面实测的点位，利用经验正交函数反演该点位的声速剖面；

2)同时，在该未实测点位，运用温盐场模型和声速经验公式获取该点位的声速剖面；

3)在同一点位垂直方向不同深度上，取以上两种方法获得声速值的均值作为相应深度上的声速值；

4)以此类推，计算网格各个节点的三维声速剖面数据；

5)在关键的网格节点上运用CTD进行实测，获取节点的实测温盐深数据；

6)对构建的三维声速剖面模型进行内外符合精度计算。

在本发明的一个实施例中，各种因素为空气、太阳、洋流以及潮汐因素。

本发明要解决的技术问题主要包括如下几个方面：

(1)顾及海洋水文要素相关性及时空变化的温盐场模型构建。

(2)研究船体瞬时姿态与波束实际投射方向的几何关系，精确计算波束初始入射角。

(3)在未实施声速剖面测量区域，综合基于时空温盐场模型声速计算方法和基于经验正交函数声速计算方法，构建高精度的声速剖面。

(4)构建一种高效的常梯度声线精确跟踪模型，既要过滤合并声速梯度变化较小的水柱层，又要准确计算声线在每个水柱层的旅行时间，实现声线精确跟踪。

通过上述技术方案，本发明的有益效果是：

(1)充分顾及了船体姿态，提出了一种顾及姿态的波束初始入射角精确计算方法，能够较大提高波束脚印坐标的精度。

(2)综合海洋卫星和Argo浮标多源海洋物理水文观测数据(海洋次表层和Argo浮标实测剖面观测数据)，利用非结构网格和有限体积法构建时空温盐场模型，可较大提高温盐深剖面的分辨率和精度。

(3)在部署的未实测的网格节点上，利用实测声速剖面数据和经验正交函数反演其他未实测点位的声速，同时通过温盐深场数据代入声速经验公式获取相应点位的声速，取两者声速的均值作为相应点位的声速值，进而构建未实测点位的声速剖面，可较大提高声速剖面的分辨率和精度。

(4)构建了一种高效的常梯度声线精确跟踪模型，在不失精度的情况下，显著提高深海传统声线跟踪的效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是错误声速剖面引起的地形变形；

图2是声速剖面改正前后的海底地形效果图；

图3是本发明的总体框图；

图4是本发明的测区时空温盐场模型构建及评估；

图5是本发明的直角坐标系图解；

图6是本发明的换能器基阵旋转角度模型；

图7是本发明的波束点空间旋转示意图；

图8是本发明的测区三维声速剖面构建与评估；

图9是本发明的自适应分层原理图；

图10是本发明的常梯度声线跟踪示意图；

图11是本发明的波束海底投射点坐标内符合精度评估示意图；

图12是本发明的波束海底投射点坐标外符合精度评估示意图。

具体实施方式

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体图示，进一步阐述本发明。

首先，本发明涉及如下技术术语：

FVCOM模型

FVCOM模型全称是非结构网格有限体积法海洋模式(The Finite-Volume Coastal Ocean Model),FVCOM的原始方程主要包含动量方程、质量连续方程以及温度、盐度和密度方程；在物理和数学上用Mellor-Yamada2.5阶垂向湍流闭合模型及Smagorinsky水平湍流闭合模型对方程组进行闭合。在垂向上适用坐标系或通用垂向坐标系对不规则底部地形进行拟合，在水平上利用无结构三角形网格对水平计算区域进行空间离散；在数值计算上，利用对水平三角形控制体进行通量有限体积积分的方式对控制方程进行离散求解；该有限体积积分方法结合了有限元方法的自由几何拟合特性和有限差分方法的离散结构简单及计算高效的特性，利用有限体积积分方法，能够更好地保证复杂几何结构的河口海湾及海洋计算中的质量、动量、盐度、温度及热量的守恒性。

经验正交函数

经验正交函数分析方法(Empirical Orthogonal Function,缩写为EOF)，也称特征向量分析(Eigenvector Analysis)，或者主成分分析(Principal Component Analysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法；EOF分析方法能够把随时间变化的变量场分解为不随时间变化的空间函数部分以及只依赖时间变化的时间函数部分。空间函数部分概括场的地域分布特点，而时间函数部分则是由场的空间点的变量线性组合所构成，称为主要分量；这些分量的头几个占有原场内空间点所有变量的总方差的很大部分，这就相当于把原来场的主要信息浓缩在几个主要分量上，因而研究主要分量随时间变化的规律就可以代替场的时间变化研究，且可以通过这一分析得出的结果来解释场的物理变化特征。

层内常梯度声线跟踪模型

由于海水介质温度、盐度和压力的变化，导致声波传播的速度分布不均匀；声波传播产生折射，折射效应使得声线的传播路径不再是直线，而是一条连续的曲线，其弯曲程度(声线的曲率半径)与声速分布相关。声波从换能器发射点到接收点时延大于声波直线传播时延，为得到多波束测深系统波束脚印的准确位置，需要沿波束的传播路径逐层追踪声线，而认为每个水柱层内声速的梯度是恒定的，已知波束的初始入射角和单程旅行时间，则可根据Snell法则，进行声线追踪，获得波束在海底投射点的船体坐标。

本发明深海多波束声线精确跟踪方法主要包括以下步骤：

(1)在全面分析海洋表层温度、盐度易受各种因素(空气、太阳以及潮汐等)影响的基础上，并综合海洋卫星和Argo浮标多源海洋物理水文观测数据，建立一种时空变化的海洋温盐场模型；一种时空变化的海洋温盐场数学模型主要包括以下步骤：

1)通过海洋卫星监测获取海洋次表层物理水文观测数据；

2)通过Argo浮标获取测量剖面的物理水文观测数据；

3)在综合海洋卫星和Argo浮标多源观测数据的基础上，运用海洋动力学数值模拟模型(FVCOM)建立时空温盐场模型；

4)在关键的网格节点上运用CTD进行实测，获取节点的实测温盐深数据；

5)最后对建立的时空温盐场模型进行内外符合精度计算。

(2)通过分析船体瞬时姿态对波束入射角的影响，顾及船体姿态精确计算每个波束初始入射角；

(3)在声速剖面未实测区域，分别通过时空温盐场模型计算声速和基于经验正交函数反演声速，在相应点位取两种方法获取声速的均值作为该点位的声速，以此计算三维声速剖面模型；计算三维声速剖面模型主要包括以下步骤：

1)在测区未实施声速剖面实测的点位，利用经验正交函数反演该点位的声速剖面；

2)同时，在该未实测点位，运用温盐场模型和声速经验公式获取该点位的声速剖面；

3)在同一点位垂直方向不同深度上，取以上两种方法获得声速值的均值作为相应深度上的声速值；

4)以此类推，计算网格各个节点的三维声速剖面数据；

5)在关键的网格节点上运用CTD进行实测，获取节点的温盐深数据；

6)对构建的三维声速剖面模型进行内外符合精度计算。

(4)构建一种高效的常梯度声线精确跟踪模型；

(5)提供声线跟踪精度评估方法。

参见图3至图12所示，本发明的具体实施方式如下：

(1)总体技术路线

第一步，在测区按一定间距部署声速剖面网格；第二步，在网格关键节点上实施声速剖面测量，并利用经验正交函数反演未实测节点的声速剖面；同时，利用海洋卫星和Argo浮标获取测区海洋次表层和Argo浮标实测剖面的时空温盐场数据，运用声速经验公式和插值算法获取相应点位的声速值；第三步，在网格节点不同深度的点位上，取以上两种方法获取声速值的均值作为相应深度的声速值，从而构建高精度的三维声速剖面；最终，在高精度的三维声速剖面基础上，结合精确的波束初始入射角，利用改进的常梯度声线跟踪方法实现波束海底投射点坐标的精确计算。

(2)海洋时空温盐场模型建立及评估

由于海洋温盐场建模的数据主要来源于实际观测点采集的数据。这些数据大都是有限的、离散的且分布不规则，而海洋温盐场的实际数据又是连续的。充分利用海洋卫星获取海洋次表层的物理水文观测数据，同时通过Argo浮标获取测量剖面的物理水文观测数据，在综合海洋卫星和Argo浮标多源观测数据的基础上，充分利用这些有限的、散乱的海洋温盐信息，来客观地、全面地、合理地恢复海洋温盐场分布形态，是海洋温盐场建模的必要环节与研究重点。因此，基于海洋动力学方程并结合实测数据的数值模拟方法成为一种有效的获取海洋环境要素数据的方法。较为常用的海洋动力学数值模拟模型有POM(princeton ocean mode1)、FVCOM(an unstructured grid，finite-volume coastal ocean model)、HAMSOM(hamburg shelf ocean model)、HYCOM(hybrid coordinate ocean model)等；而FVCOM模型是基于非结构网格和有限体积法，在高分辨率及小尺度计算方面上优势明显，故采用此模型进行时空温盐场建模，模型构建之后需要在关键的网格节点上运用CTD进行实测，并对构建的时空温盐场数据进行内、外符合精度计算，以此评估构建模型的精度。时空温盐场模型构建及评估过程如图4所示。

下面式子是其在笛卡尔坐标下作为核心的动量方程。

FVCOM模型采用σ坐标系(见图5)，其与Z坐标系的变换关系如下：

式中σ是t坐标系的时间自变量；H是水深，ζ是潮高。

通过式(2)我们可以得出，在海底，σ＝-1，在海表面，σ＝0。FVCOM模型在σ坐标系下的控制方程分别为：

式中，T代表温度，S代表盐度；t代表时间；x，y，z分别表示东向、北向和垂直方向；u，v，w分别为x，y和z方向上的速度分量；K_h为热量垂直旋转扩散系数；F_T，F_S分别为水平方向上热量与盐度的扩散项。

在边界条件上：

σ＝0时，

σ＝-1时，

式中，Q_n(x,y,t)其中为表面净热量通量，包括四部分：向下的短波，长波辐射，显通量和潜通量；SW(x,y,ζ,t)为海表面的短波通量；c_p为海水比热；A_H为水平热量扩散系数；α为底面地形；ρ为密度；其中分别为降水率和蒸发率。

(3)顾及姿态的波束初始入射角计算

由于船体坐标系的中心通常是以换能器为中心，因此以换能器理想水平状态作为基准面来分析船体姿态对波束初始入射角的影响，如图6所示的换能器基阵坐标系中，水平状态的换能器基准面位于OABC平面内，O为换能器中心，OA为基准面纵轴正方向，OC为基准面横轴正方向。设OA长度为a，OC长度为c，A、B两点的坐标分别为(a,0,0)和(0,c,0)。在某一姿态(横摇、纵摇角分别为r和p)影响下基准面变化为OA₁B₁C₁，即基阵面由水平先绕OX轴旋转角度α(α≠r)，再绕OY轴旋转角度β形成。A点、C点经过两次旋转后分别转到A₁和C₁位置，A₁、C₁两点在水平面OXY上的投影分别为A₂和C₂。在此状态下，OA₁与水平面夹角∠A₁OA₂即为纵摇角p，OC₁与水平面夹角∠C₁OC₂即为横摇角r。根据纵摇角、横摇角和旋转角定义，r和α符号一致，p和β符号一致。

由以上过程可知，基准面OABC经过α和β两次旋转得到OA₁B₁C₁，则有：

则旋转后A₁点坐标为：

旋转后C₁点坐标为：

由式(7)得到旋转后A₁点坐标再根据三角形正弦定理可计算出基阵面的纵摇角p(即∠A₁OA₂)：

式(9)中Z_A1为A₁点在Z轴上的坐标，根据β与p符号一致可得：

β＝p (10)

类似地，由式(8)得到旋转后C₁点坐标再根据三角形正弦定理可计算出基阵面的横摇角r(即∠C₁OC₂)：

式(11)中Z_C1为C₁点在Z轴上的坐标，r和α符号一致，并将β＝p带入式(11)得：

sinr＝sinαcosp (12)

则有：

α＝arcsin(sinr/cosp) (13)

由式(10)和式(13)可知，在旋转变换中，绕OY轴旋转角β等于纵摇角p，而绕OX轴旋转角α并不等于横摇角r。因此，在目前较为精细的声线跟踪计算中，即使顾及姿态影响，定义初始入射角为θ₀+r显然是不正确的。

为了得到姿态(r、p)影响下真实的波束入射角，下面推导实际波束初始入射角θ′₀的计算模型。

由以上推导可知，实际声线可由理想状态下的声线经α、β旋转变换R后得到。设理想状态下，第i个波束分配初始入射角为θ_i，在不失精度的情况下，假设经历第一个水层以常声速传播，传播距离为R_i，则波束在第一水层下界的落点P_i坐标为(0，R_isinθ_i,R_icosθ_i)，而在姿态影响下的实际坐标(x_i,y_i,z_i)为：

式(14)可借助图7解释。假设换能器基阵水平时，第i号波束的波束角度为θ_i，斜距为R，则点A坐标为(0，R_isinθ_i,R_icosθ_i)，换能器基阵在横摇r和纵摇p的影响下，A点旋转到了B点，第i号波束的实际入射角度为θ′_i(即∠BOD)，定义经旋转后第i号波束的水平角度即为波束横距BD与OY轴的夹角，其表达式为：

由式(15)可获得姿态影响下的波束实际初始入射角，之后根据改进的精确声线跟踪方法进行三维声线精确跟踪(参见图7)。

(4)三维声速剖面构建及评估

由于多波束实测区域较大，不可能每个点位进行声速剖面测量，只能在一定距离实施声速剖面测量，其他点位的声速剖面通过以下步骤获得：①利用经验正交函数反演声速剖面；②利用温盐深场数据代入声速经验公式获得声速剖面；③取相同位置上两种方法获取的声速的均值作为该位置的声速，由此构建整个声速剖面。三维声速剖面构建与评估如图8所示。①利用经验正交函数反演声速剖面

以网格形式布施M个声速剖面进行非等间隔采样，在声速剖面变化剧烈的深度上采样较密集，变化缓慢的深度上采样较稀疏。则M个散射后的样本声速剖面c₁(z₁)，c₂(z₂)，…，c_M(z_j)的平均声速剖面为：

上式中，下角标i和j分别表示第i个声速剖面和第j个深度。根据M个样本声速剖面和平均声速剖面定义协方差矩阵R为：

上式中N是深度上的采样点数，Δc_i(z_j)是第i个样本声速剖面与平均声速剖面之差，即

将协方差矩阵R进行特征分解得到特征值λ_n和对应的特征向量f_n，则R可以表示为：

选取前K个较大的特征值所对应的特征向量作为经验正交函数。

已知平均声速和各阶经验正交函数f_k(z)后，待反演的声速剖面可以表示为：

上式中，α_k(x,y)是经验正交函数系数，为待反演参数，通常是水平坐标x,y的函数，当声速在水平方向上变化不是很大时可以将α_k(x,y)近似认为是常数。通过式(21)，声速剖面反演这样一个复杂的问题就转化为经验正交函数系数求解的形式。

②利用声速经验公式获取声速剖面

利用上述获得的温盐深数据和声速经验公式，得到相应位置的声速，从而构建出声速剖面。一些学者通过分析7种声速经验公式的适用范围，普通认为EM分层简化声速公式在1-12000m整个深度层中具有良好的计算精度，因此本文采用该公式计算声速。

水面声速：

C(0,T,S)＝1449.05+T[4.57-T(0.0521-0.00023T)]+[1.333-T(0.0126-0.00009T)](S-35) (22)

海水中水面到1000m深度的声速：

C(T,D,S)＝C(0,T,S)+16.5D (23)

海水中深度1000m到11000m的声速：

式(22)、(23)和(24)中的T，D，S分别表示温度、深度和盐度。

③最终计算的声速剖面

分别经过声速剖面反演和声速经验公式计算获得c^P(z)、c^T(z)后，取相同位置声速的均值为该位置最终的声速值，进而构建声速剖面。

(5)高效的常梯度声线精确跟踪模型

由于数千米深海的声线跟踪过程耗时太长，且传统常梯度声线跟踪算法在计算每个水层旅行时间的不精确导致较大误差，为此，提出一种自适应分层的常梯度声线精确跟踪算法。

首先利用自适应分层方法对原始声速剖面进行抽稀，该方法思想是通过对声速剖面数据进行合理抽稀并保留信息量多的声速层，在不失跟踪精度的基础上减少跟踪时间。声速剖面自适应分层原理如图9所示。

在抽稀过程中，首先需要设定合理的拟合窗口大小n,n取值过小会造成拟合多项式的参数不准确，n取值过大会造成拟合后层数过少。其次，从顶层开始选取连续的n个声速值Ci(i＝1,2,…,n)进行曲线拟合。如果曲率小于等于设定的阈值，则将此n层水柱合并为一层，声速值取此n层的均值；若曲率大于设定的阈值(表明此数据中声速剖面变化较大)，则此n层不能合并为一层，需要释放n个点中的第一个点，下移一个点继续选取n个点声速值进行拟合。通过此方法将声速剖面中声速梯度变化小的水柱层过滤，适当抽稀获取声速梯度变化较大的水柱层，进而减少声线跟踪时间。

在抽稀后的声速剖面中进行声线精确跟踪，海水中声速与海水温度、盐度和静压力相关，随着深度的变化而变化。由于很难获得声速随深度变化的函数，通常只能借助声速剖面仪获得一定深度间隔的声速剖面，并根据声速变化，借助Snell法则，沿声线跟踪，获得波束海底投射点在船体坐标系下的坐标。在每个水层中，由于只知道该层上界和下界声速，因此通常假设声速在该层以常梯度g传播，采用类似处理方法处理其他各水层，并沿着声线跟踪到海底，即实现常梯度声线跟踪和波束海底坐标计算。

假设声线从换能器发射，经过N个水层，每个水层上界Z_i声速为C_i，下界Z_i+1声速为C_i+1，层内声速以常梯度g_i变化，则层i内声速函数为：

Cⁱ(z)＝C_i+g_i(z-z_i) (26)

如图10所示，在常梯度g_i声速变化下，声线在第i层内的传播轨迹为一连续的、带有一定曲率半径为R_i的弧段，若Snell常数为p,则R_i为：

R_i＝-1/pg_i (27)

在第i层内，声线经过的弧线微元为ds,所需要的时间为dt＝ds/C，而垂向微元dz与ds的关系为ds＝dz/cosθ,θ则为波束在弧线微元的入射角，声线在此层内传播的时间t_i和水平位移Δy_i分别为：

在分层声线跟踪时，除了计算整层的垂直位移、水平位移和传播时间外，还需要依据传播剩余时间计算剩余层的垂直位移和水平位移。假设声线在第i层内传播时，声线在该层内r点处(如图10)结束，此时剩余时间t_r等于波束单程旅行时间t_all减去第i层以前累计的传播时间，则声线在剩余层内的垂直位移Δz_r和水平位移Δy_r为：

Δy_r＝R_i(sinα_i-sinα_r)

则声线传播总的垂直位移z和水平位移y为：

(6)声线跟踪精度评估

建立一种顾及姿态的高效声线精确跟踪模型后，需要对该模型的精度进行评估，而对声线的跟踪精度的评估归根结底是对波束海底投射点坐标精度的评估。为此，可通过以下方法对波束海底投射点坐标进行内外符合精度评估，为水下导航定位服务提供依据。

①内符合精度计算

对波束海底投射点坐标内符合精度评价可选择3条平行的相互间有50％重叠度的相邻条带以及中央波束与前3条正交的间距较小的2个条带，由于多波束中央波束相对于边缘波束受姿态影响小得多，在声速剖面正确情况下可视为实际深度，因此以后2个条带的中央波束作为检测线，作为前3个条带所对应平面位置测深点深度的精度评估。如图11所示，条带①、②、③为相邻相互重叠度为50％的待评估条带，以航行方向相反的条带④、⑤的中央波束航迹线作为内符合精度计算的检测线。

以条带④、⑤的检测线的平面坐标为基准，分别在条带①、②、③区域插值出相应位置的深度值，并与检测线上相应位置的深度值作内符合精度计算。设m为重复测线的数目，n为重新测线公共段数据点的数目，则对于所有重复测线观测值，总的内符合精度计算公式为：

δ_ij＝Δg_ij-Δg_i (33)

其中，δ_ij为第j条重新测线公共段上的第i个波束海底投射点观测值Δg_ij与该点各重复测线观测的平均值Δg_i之差。

②外符合精度计算

外符合精度计算可通过已经布设好的水下大地基准点(坐标已知)，当船体位于不同位置时分别测量水下不同大地基准点的斜距，然后根据斜距、姿态可推算出不同水下大地基准点的平面坐标和水深值，取重复测量的各个水下大地基准点的平面坐标均方差和深度值均方差，以此评估声线跟踪的外符合精度。如图12所示，当船在①位置时，分别测量出3个水下大地基准点(a)、(b)、(c)的斜距，根据船体姿态和斜距可计算出各个基准点的平面坐标和水深值，以此类推，可计算出船在②、③位置时相对水下基准点的斜距，进而计算出各个基准点的平面坐标和水深值。以3次观测的水面坐标和水深值的均方差作为声线跟踪平面坐标和深度值的外符合精度。

结合图7，假设船在①位置时，根据声线跟踪计算出船到水下大地基准点(a)的斜距R₁^a，再根据船体姿态计算波束入射角(式(15))和方位角(式(16))，计算水下基准点(a)的水深值、X坐标和Y坐标。

同理，可得当船分别在②、③位置时，以此类推可得其中i＝b,c；j＝2,3。

在水下基准点(a)处声线跟踪三个轴向的外符合度分别为：

同理，可得水下基准点(b)及(c)处声速跟踪三个轴向的外符合度。

由此可见，本发明中解决了下列技术问题：

(1)顾及海洋水文要素相关性及时空变化的温盐场模型构建。

(2)研究船体瞬时姿态与波束实际投射方向的几何关系，精确计算波束初始入射角。

(3)在未实施声速剖面测量区域，综合基于时空温盐场模型声速计算方法和基于经验正交函数声速计算方法，构建高精度的声速剖面。

(4)构建一种高效的常梯度声线精确跟踪模型，既要过滤合并声速梯度变化较小的水柱层，又要准确计算声线在每个水柱层的旅行时间，实现声线精确跟踪。

另外，本发明的技术特点如下：

(1)顾及海洋水文要素相关性及时空变化的温盐场模型构建方法。

(2)顾及姿态的波束初始入射角精确计算方法。

(3)基于时空温盐深模型的高分辨率高精度声速剖面模型构建方法。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

参考文献：

[1]董庆亮,韩红旗,方兆宝等.声速剖面改正对多波束测深的影响[J].海洋测绘,2007,27(2):56-58.

[2]黄承义,肖春桥,田丰等.基于EM系列多波束水深测量的声速剖面校正方法及技巧[J].测绘通报,2012,(s1):698-699.

[3]赵建虎.多波束测深及图像数据处理[M].武汉大学出版社,2008.

[4]Alexandrou D.,De Moustier C.Adaptive noise canceling applied to Sea Beam sidelobe interference rejection[J].IEEE Journal of Oceanic Engineering,1988,13(2):70-76.

[5]Lasky M.Review of undersea acoustics to 1950[J].Journal of the Acoustical Society of America,1977,61(2):283-297.

[6]Munk W.,Wunsch C.Ocean acoustic tomography:a scheme for large scale monitoring[J].Deep Sea Research Part A Oceanographic Research Papers,1979,26(2):123-161.

[7]Munk W.,Wunsch C.Ocean acoustic tomography:Rays and modes[J].Journal of Trauma,2003,55(21):1004-1005.

[8]Jones R.M.,Shang E.C.,Georges T.M.Nonperturbative modal tomography inversion.Part I.Theory[J].Journal of the Acoustical Society of America,1993,94(4):2296-2302.

[9]Jones R.M.,Georges T.M.Nonperturbative ocean acoustic tomography inversion[J].Journal of the Acoustical Society of America,1994,96(1):439-451.

[10]Jones R.M.,Howe B.M.,Mercer J.A.,et al.Nonperturbative ocean acoustic tomography inversion of 1000-km pulse propagation in the Pacific Ocean[J].Journal of the Acoustical Society of America,1994,96(5):439-451.

[11]Dinn D.F.,Loncarevic B.D.,Costello G.The effect of sound velocity errors on multi-beam sonar depth accuracy[A].In Oceans'95Mts/ieee Challenges of Our Changing Global Environment Conference Proceedings[C],1995:1001-1010vol.1002.

[12]Goncharov V.V.,Voronovich A.G.An experiment on matched-field acoustic tomography with continuous wave signals in the Norway Sea[J].Journal of the Acoustical Society of America,1993,93(4):1873-1881.

[13]唐俊峰.海水中的声速剖面反演研究[D].哈尔滨工程大学,2006.

[14]丁继胜,周兴华,唐秋华等.多波束勘测声速剖面场的EOF表示方法[J].武汉大学学报:信息科学版,2007,32(5):446-449.

[15]Ramezani H.,Jamali-Rad H.,Leus G.Target Localization and Tracking for an Isogradient Sound Speed Profile[J].IEEE Transactions on Signal Processing,2013,61(6):1434-1446.

[16]吴德明.一种用于声线修正的迭代法[J].声学学报:中文版,1992,(2):104-110.

[17]王恕铨,Farmer D.求三维本征声线的一种新方法[J].声学学报:中文版,1992,(2):155-157.

[18]姜薇,李太宝.三维声线追踪的正三棱锥前向伸展算法[J].声学学报:中文版,2005,(5):404-408.

[19]兰华林.深海水声应答器定位导航技术研究[D].哈尔滨工程大学,2007.

[20]孙万卿,葛亮.一种基于有穷状态自动机的二维浅海声线追踪方法[J].兵工学报,2007,28(11):1382-1387.

[21]陆秀平,边少锋,黄谟涛等.常梯度声线跟踪中平均声速的改进算法[J].武汉大学学报:信息科学版,2012,37(5):590-593。