距离依赖的声速场反演与运动声源的跟踪定位

戴 淼， 李亚安

(西北工业大学 航海学院，西安 710072)

水下声传播环境是一种地区性、参数值多变的复杂环境，现场仪器测量难以大范围的精确获取所需要的海洋环境参数。声学反演技术将声传播模型与观测数据相结合，能够快速简便的获取所需要的参数。声速剖面 (Sound Speed Profile, SSP) 是影响海洋声传播的重要因素之一，其结构特征决定了水下声传播路径和传播损失[1]。经验正交函数 (Empirical Orthogonal Function, EOF) 是表示十分有效的基函数，通常只需要前几阶便可以精确的反演出SSP，大大减少了待反演参数[2]。已有研究成果表明：在水平变化环境下的平均SSP反演[3]和SSP特征量统计分布规律的获取[4]与同步观测数据相比，具有很好的一致性。在这里反演得到的是平均意义上的估计，SSP所具有的时空变化特性，尤其是声速快变区域的存在将影响反演精度。

序贯滤波方法基于系统的观测值对系统状态进行估计和更新，能够根据相关的物理统计模型预测未知海洋声学参数[5]。卡尔曼滤波器用于处理跟踪问题[6]，要求状态方程和观测方程均为线性。扩展卡尔曼滤波(Extended Kalman Filter, EKF)把非线性问题线性化，在海洋声学参数反演方面取得了一定效果[7]。Yardim等[8]通过EKF、无迹卡尔曼滤波(Unscented Kalman Filter, UKF)和粒子滤波(Particle Filter, PF)将地声反演转换为跟踪问题，对具有时空变化特性的环境参数进行估计。Carrière等[9]将声学测量数据同化问题建模为状态-空间模型，把沿海环境距离依赖声速场参数作为状态变量，利用EOF进行参数低维化，并首次应用集合卡尔曼滤波(Ensemble Kalman Filter, EnKF)估计了SSP随时间的演化[10]。Li等[11]对比多种序贯滤波算法分析了距离独立环境下地声参数反演跟踪情况。需要指出的是上述研究过程均采用固定收发装置，即声源与接收阵之间是相对静止的。

传统的固定收发装置得到的反演结果为声源和接收阵之间的平均SSP，空间分辨率存在局限性；在实际情况中，受季节、温度以及动态海洋(存在涡旋、内波、冷热水团等)的影响，收发装置之间存在的SSP会随着时间/距离的变化而变化，尤其是声速快变区域的存在，这就要求更多的待反演参数，故而反演算法的速度和精度都将受到制约。针对上述问题，由于SSP的时空演化特性可在一定条件下近似的建模为EOFs系数随时间/距离的演化，在已有序贯滤波算法基础上，本文提出借助运动声源的时变性，实现距离依赖快变声速场反演与运动声源位置的联合估计。通过引入运动声源在不同空间位置处发射信号，利用垂直阵进行接收，分别采用EKF和EnKF同时处理多个待反演参数。该方法相比固定收发布设装置增加了声场的空间信息，同时声场中包含的信息又可以逆推声源的位置，因而能够在重构真实声速场的同时对运动声源进行跟踪定位，并且根据估计结果的均方根误差对算法的估计性能进行分析、比较。研究表明：EnKF表现出较优越的估计性能，即使在阵元数目减少的局限条件下，该算法仍然能够很好的跟踪估计出运动声源参数和声速场在时间/距离维度上的演化轨迹。

1 理论与算法

1.1 距离依赖条件下的声传播模型

海洋环境复杂多变导致无法对其精确建模，一般采用量化水下各种物理过程对声传播影响的方法，运用合理的声传播模型最大程度的重构水下物理场。简正波模型通过求解与深度有关的特征方程，得到一组振动模式，振动的频率给出与模式传播相关联的水平波数，对加权的每个模式的贡献进行叠加来构成总声场。简正波声传播模型本质是利用简正波模式求解亥姆霍兹方程[12],频率为ω的空间任意点源s(r,z)在柱坐标下声压可表示为

(1)

式中：c(r,z)表示声速，在距离独立条下c(r,z)=c(z)，模态方程为

(2)

根据边界条件，声压正比于各阶模态的和与汉克尔函数的乘积

(3)

式中：φm(z)和km分别表示为深度方向上的模态和水平波数。

将距离独立环境推广至距离依赖简正波声传播理论：水平距离分成段，各段均近似视为距离独立环境，并对每一段以标准简正波求解声场，最后根据边界条件将各段的解连接在一起[10]：

(4)

(5)

(6)

声压在第l个界面上的连续性条件为

(7)

耦合简正波的数值求解过程非常复杂[13]，在实际计算中，通常采用绝热简正波理论[14]忽略模态之间的能量耦合项，即各阶模态从当前距离至下一距离是绝热耦合的。

1.2 非线性滤波算法

滤波算法将不同性质的测量数据有效结合起来，并且能够给出相关尺度范围内与观测数据吻合的动力过程的确定性描述，目的在于更好的估计所需要的参数，改进动力模型的性能。考虑到距离依赖SSP的动态演化特性，采用序贯滤波的思想对其进行预测估计，下面给出两种非线性滤波算法。

1.2.1 扩展卡尔曼滤波

扩展卡尔曼滤波 (Extended Kalman Filter, EKF)方法实质是将非线性模型在其状态变量均值的范围内进行泰勒级数展开，并在当前估计值处线性化，即对过程方程或测量方程求偏导，从而将非线性问题线性化。EKF算法可表示为以下递推形式：

首先给出两个非线性预设方程：

xk=f(xk-1,uk,wk)

(8)

yk=h(xk,vk)

(9)

函数f用k-1时刻的状态和k时刻的输入uk预测k时刻的状态，相应的，函数h以k时刻的状态预测k时刻的测量值，wk和vk为噪声向量。实际中并不知道每一时刻状态噪声和观测噪声的值，可以将其假设为零。

k时刻的预测值：

xk|k-1=f(xk-1|k-1,uk, 9)

(10)

yk|k-1=h(xk|k-1, 0)

(11)

状态协方差矩阵的预测值

(12)

k时刻的更新：

新息

rk=yk-yk|k-1

(13)

新息的协方差矩阵

(14)

卡尔曼增益

(15)

状态协方差矩阵

Pk|k=(I-GkFk)Pk|k-1

(16)

状态向量

xk|k=xk|k-1+Gkrk

(17)

上述函数用到了状态转换函数和观测函数的Jacobian:

(18)

(19)

式中：下标k|k-1表示由时刻的值预测得到k时刻的值，下标k|k表示使用k时刻新息后的更新值，下标k则表示k时刻的值。

1.2.2 集合卡尔曼滤波

集合卡尔曼滤波(Ensemble Kalman Filter, EnKF)基本思想:初始化产生一组系统状态的采样作为背景集合，根据集合的预测结果估计状态向量和观测向量的协方差，然后利用观测信息和协方差信息计算分析集合，通过系统模型传递分析集合并得到下一时刻的背景数据集合，递推向前进行预测。EnKF算法可表示为：

(20)

(21)

(22)

卡尔曼增益

(23)

式中：Rk为k时刻的观测误差协方差矩阵。

根据k时刻的观测信息，更新背景集合，得到k时刻的分析集合：

(24)

(25)

对比两种滤波算法可以看出：EnKF使用非线性集合估计真实统计值，省略了对模型算子和观测算子的线性化过程，在提高估计精度的同时降低了计算量。

1.3 声速剖面的经验正交函数表示

SSP反演作为一种多维优化问题，反演算法的复杂度和精度通常取决于待反演参数的个数。EOF用于表征海水的SSP，大大减少了描述声速垂直结构所需要的参数，其原理如下：

EOF将SSP随深度变化的关系参数化。通常可以由观测海域的实测SSP计算得到，不同时刻测得的SSP可表示为：

(26)

(27)

式中：fn为矩阵R的特征向量，即经验正交函数；λn为特征向量对应的特征值。协方差矩阵R表征了该海域声速起伏的不确定量。将上述的N个特征值λn按从大到小的顺序排列，选取前K个特征值对应的特征向量来表示SSP，则海域内任意一点的声速值可表示为：

(28)

式中：αk(x,y)为第k阶EOF系数。

由式(28)可知，EOF是SSP在深度维上的离散化。根据距离依赖的声传播理论，现将SSP在距离维进行离散化处理，即把声源和接收阵之间的水平距离分成段，相应的声速可表示为

(29)

ul(r)为门函数，定义为：

(30)

在第l个区域，有

(31)

式中：αk, l为相应距离处的EOFs系数。

在最小均方意义下，EOF用于表征SSP在保证反演精度的前提下，可降低算法复杂度及计算时间，从而提高算法性能。

1.4 运动声源的状态方程

(32)

式中：s为状态变量，Δt为状态更新的时间增量 (即两次测量的时间间隔) ，vz和va分别表示速度在深度方向的扰动和径向的加速度扰动，H为状态转移矩阵。

状态方程的第二部分用于描述SSP的状态，将EOFs系数作为状态参量，SSP的状态方程可由前n阶EOFs系数的状态方程代替：

(33)

式中：wk为EOFs系数在空间上演化的过程噪声。

在这里，将观测方程用于表征垂直线列阵 (VLA) 接收到的声压信号，根据简正波理论，合成声压可以表示为水平距离R，深度Z，频率ω，声速c以及海底边界条件参数(BCs)的非线性方程：

yk=h(Z,R,ω,c1, …，ck,BCs)+vk

(34)

式中：h为简正波声传播模型，ck为第k个距离上的声速剖面，vk为测量噪声。

运动声源的引入可增加水下声场的空间信息，同时声场中所包含的信息又能够逆推声源的位置，通过观测方程将状态向量 (待反演参数) 和观测向量 (VLA接收到的声压) 联系起来。

2 结果与分析

2.1 实测声速剖面的EOF反演

为证明EOFs系数表示SSP的有效性，将给出EOF反演SSP与实测声速剖面的对比。中美联合水声实验组于2001年6月在东中国海进行宽带信号传播实验(ASIAEX,ECS,2001)[15]，图1为现场CTD实测的声速样本和平均声速剖面，其中细线为实测的54组声速剖面，粗线为平均声速剖面。可见，实验期间温跃层 (30～60 m) 范围内声速随时间的起伏较大，变化最大峰值达到了11 m/s以上。

图1 东海不同时刻声速剖面样本Fig.1 Sound speed profiles in different time of the East China Sea

图2给出对该组声速剖面所提取到的EOFs：其中图2(a)为图(1)所示的声速起伏，并给出如图2(b)所示的前六阶EOFs值，由图2(c)可见，前三阶EOFs可表征近90%的声速起伏。利用图2提取的6阶EOFs对54条声速剖面中任意时刻CTD所测得的其中一条声速剖面进行拟合，并将反演得到的SSP与实测SSP对比于图3中。

验证结果表明：只需利用前几阶EOFs即较少EOFs系数便能够较为精确的表征该海域海水SSP的演变过程，此方法相比于直接将声速剖面表征为深度的函数大大减少了待反演参数的个数。在实际应用过程中，通常希望在保证一定精度的前提下用尽量少的参数来表示复杂的声场环境，EOF阶数的选取还应考虑计算量的大小。一般的，只需选取3～6阶EOFs便可较为准确的表征声速场的起伏情况。

(a) 声速剖面样本差异

(b) 前六阶EOFs

(c) 前10阶EOFs所占声速剖面起伏的能量比图2 提取的海水声速剖面经验正交函数Fig.2 EOFs of SSPs in the water volume

图3 某一时刻CTD实测声速剖面(实线)与EOF反演声速剖面(虚线)Fig.3 CTD measured SSP (solid line) and EOF invert SSP (dotted line)

2.2 距离依赖声速场重构和运动声源的跟踪定位

由于没有获取到同步的声压测量数据，因此本文中声压根据环境模型拟合得到，波导环境参考Benchmarking距离依赖浅海波导模型，海底参数根据东中国海(ASIAEX,ECS,2001)实测数据设置。假设仿真波导环境如图4所示：在平整海面条件下，海面边界为简单的压力释放边界。海深106 m,16元VLA等间距间隔4 m布设于距离海面15 m处，阵长60 m。运动声源初始深度为40 m，从水平方向初始位置1 km处以2.5 m/s的径向速度远离接收阵，并每隔20s发射一次信号，声源辐射频率为400 Hz，20 min内共移动60次。

图4 距离依赖的浅海环境参数模型Fig.4 Experimental configuration of the range dependent shallow water environment

为了减少待反演样本数目，拟采用前三阶EOFs系数作为声速参数来实现声速场的重构。假定过程噪声和测量噪声服从零均值的高斯分布，声场模型为绝热简正波模型[17]，声压作为观测值由声学软件KRAKENC根据环境参数合成。分别用EKF和EnKF滤波算法对运动声源参数(深度、距离和速度)以及EOFs前三阶系数同时进行跟踪，并将跟踪估计结果与真实值进行比较。其中，EnKF样本大小为50。

图5为EKF和EnKF的参数联合估计结果，其中，(a)～(c) 分别为运动声源深度、距离和速度随时间的演化，(d)～(f) 分别为EOFs前三阶系数随距离的演化。实线代表运动声源参数和前三阶EOFs系数随时间/距离的真实演化过程，点线和虚线分别是EKF和EnKF估计结果。可以看出，EKF除了运动声源的速度以外，其余参数均能够跟踪上，产生这一现象的主要原因是由于拟合的声压中并不包含声源的速度信息。EnKF对于所有反演参数均能够跟踪上，且估计结果相对稳定，波动较小，相比于EKF更能逼近待反演参数真实状态随时间/距离的演化轨迹。

图6给出上述两种滤波算法估计的运动声源参数在时间维的均方根误差，可以看出，EnKF在声源的深度、距离和速度随时间演化过程中均表现出较EKF优越的估计性能，几乎与声源的真实位置 (深度、距离) 相吻合，均方根误差在实时观测过程中接近于零。

图5 EKF和EnKF参数联合估计结果,图中实线为参数真实演化轨迹Fig.5 Simultaneous estimation results of EKF (dash-dotted line) and EnKF (dashed line). The true values trajectories are in solid line

图6 运动声源参数Fig.6 Root-mean-square-errors of the moving source parameters estimation in time dimension

由EKF和EnKF估计得到的声速参数重构的声速场和其相对误差如图7所示。由(a)，(c)可知，两种滤波算法均能够重构出真实的声速场，同时可以看到声速场个别区域存在较大的声速扰动现象。在滤波估计相对误差(b)，(d)中显示：EnKF估计效果整体优于EKF，相对误差均较小，尤其在声速变化较为剧烈区域，EnKF表现出更为理想的估计效果。

图8给出重构的声速场在距离维度上的均方根误差，两种滤波算法均能够大致跟踪上EOFs系数的真实演化趋势。其中，EnKF表现出更优越的估计性能，均方根误差波动非常小。经计算，EKF和EnKF在距离维度上的平均均方根误差分别为0.573 m/s和0.011 m/s。

图7 EKF和EnKF重构的距离依赖声速场Fig.7 Inversion of range dependent sound speed field

图8 滤波估计重构的声速场在距离维的均方根误差Fig.8. Root-mean-square-errors in range dimension of the invert sound speed fields from different filters estimation

需要指出的是：EnKF的估计性能整体优于EKF，主要是由于测量方程的非线性导致。EKF只适用于弱非线性系统，因其保留了非线性函数的一阶近似项，忽略了高阶项从而引起更大的误差；而EnKF不需要对模型算子和观测算子线性化，具有较高的处理效率和处理高维非线性系统的能力，方差和均值的估计精度至少可达2阶以上。

上述仿真均为16阵元的VLA得到的估计结果。但在实际应用中通常不能满足大量阵元数目的要求，因此需要验证EnKF滤波在减少阵元数目的局限条件下算法的宽容性。保留原有VLA的有效阵长60 m，将阵元个数减少为8个，阵间距增大至8.57 m, 由图9、10可知，在相同的时间/距离维度，EnKF仍然能够保持较为稳健的估计结果。

图9 不同阵元个数VLA的EnKF估计运动声源参数Fig.9 Root-mean-square-errors in time dimension of the moving source parameters EnKF estimation

图10 不同阵元个数VLA的EnKF估计重构声速场在距离维的均方根误差Fig.10 Root-mean-square-errors in range dimension of the EnKF estimate sound speed fields from different elements VLA

与16阵元VLA估计结果相比，8阵元VLA在相同时间/距离维度内的EnKF对运动声源参数/声速场估计结果的均方根误差均略有增大，经计算，在相同距离维度上重构的声速场平均均方根误差为0.024 m/s。阵元数目的减少降低了空间采样率，使估计结果的误差略有增大。但即使在阵元数目减少的局限条件下，EnKF仍然能够很好的跟踪估计出运动声源参数和EOFs系数即距离依赖的快变声速场的演化轨迹，体现了该算法的鲁棒性和宽容性。

3 结 论

受海洋内部复杂动力学的影响，在大范围海域或者声速扰动剧烈区域，收发装置之间存在的声速场会随着距离的变化而变化，因而要求更多的待反演参数，反演算法的速度和精度都将受到制约。本文针对距离依赖声速快变区域，对ASIAEX, ECS, 2001实验现场CTD测量的声速样本数据进行EOFs提取，明确只需前几阶EOFs便可较为准确的表示SSP，验证了具有时空变化特性的SSP在一定条件下可近似建模为经验正交函数系数随时间/距离演化的可行性，同时根据实验中观测的水文资料合成声压作为观测值。提出通过将运动声源参数(深度、距离和速度)和声速参数 (EOFs前三阶系数) 作为测量值放入状态方程，利用EKF和EnKF滤波算法实现了运动声源位置与声速场的联合估计。通过分析对比，表明本文提出的反演方法在重构真实声速场的同时可对运动声源进行跟踪定位。其中，EnKF表现出更优越的估计性能，并验证了算法的鲁棒性和宽容性。运动声源的引入，可改善利用传统的固定收发装置反演只能够得到声源和接收阵之间的平均SSP的局限性，使反演结果更接近于具有时空变化特性SSP的真实演化状态。同时，运动声源的灵活性可增加局部区域声压的采样密度，丰富声学和环境的观察信息，能够改善传统声层析法固定网络采样分辨率和观测范围的局限性，提高预报精度。

[ 1 ] CHIU L Y S, LIN Y T, CHEN C F, et al. Focused sound from three-dimensional sound propagation effects over a submarine canyon [J]. J. Acoust. Soc. Am., 2011, 129(6): 260-266.

[ 2 ] LEBLANC L R, MIDDLETON F H. An underwater acoustic sound velocity data model [J]. J. Acoust. Soc. Am., 1980, 67(6): 2055-2062.

[ 3 ] 何利，李整林，彭朝晖，等. 南海北部海水声速剖面声学反演[J]. 中国科学：物理学 力学 天文学, 2011, 41(1): 49-57.

HE Li, LI Zhenglin, PENG Chaohui, et al. Inversion for sound speed profiles in the northern of South China Sea [J]. Scientia Sinica Pysica,Mechanica & Astronomica, 2011, 41(1): 49-57.

[ 4 ] 李佳，杨坤德，雷波，等. 印度洋中北部声速剖面结构的时空变化及其物理机理研究[J]. 物理学报, 2012, 61(8): 282-264.

LI Jia, YANG Kunde, LEI Bo, et al. Research on the temporal-spatial distributions and the physical mechanisms for the sound speed profiles in north-central Indian Ocean [J]. Acta Physica Sinica, 2012, 61(8): 282-264.

[ 5 ] YARDIM C, MICHALOPOULOU Z H, GERSTOFT P. An overview of sequential Bayesian filtering in ocean acoustics [J]. IEEE J. Ocean Eng., 2011, 36(1): 71-89.

[ 6 ] CULVER R L, HODGKISS W S. Comparison of Kalman and least squares filters for locating autonomous very low frequency acoustic sensors [J]. IEEE J. Ocean. Eng.,1988, 13(4) :282-290.

[ 7 ] CANDY J V, SULLIVAN E J. Model-based environment inversion: A shallow water ocean application [J]. J. Acoust. Soc. Am.,1995, 98(3): 1446-1454.

[ 8 ] YARDIM C, GERSTOFT P, HODGKISS W S. Tracking of geoacoustic parameters using Kalman and particle filters [J]. J. Acoust. Soc. Am.,2009, 125(2): 746-760.

[11] LI J L, ZHOU H. Tracking of time-evolving sound speed profiles in shallow water using an ensemble Kalman-particle filter [J]. J. Acoust. Soc. Am., 2013, 133(3): 1377-1386.

[12] PEKERIS C L. Theory of propagation of explosive sound in shallow water [J]. Geologi11cal Society of America Memoris, 1948, 27: 1-116.

[13] EVANS R B. A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom [J]. J. Acoust. Soc. Amer.,1983, 74(1): 188-195.

[14] PORTER A D, REISS E L. Extension of the method of normal modes to sound propagation in an almost-stratified medium [J]. J. Acoust. Soc. Amer., 1965, 37(1): 19-27.

[15] DAHL P, ZHANG R H, MILLER J, et al. Overview of results from the Asian Seas International Acoustics Experiment in the East China Sea [J]. IEEE J. Oceanic Eng.,2004, 29(4): 920-928.

[16] Becker K M, Frisk G V. The impact of water column variability on horizontal wave number estimation and mode based geoacoustic inversion results [J]. J. Acoust. Soc. Amer., 2008, 123(2) : 658-666.

[17] PORTER M B. The KRAKEN normal mode program[C]. Technical Report, SACLANT Undersea Research Center, La Spezia, Italy, 1991.

[18] LUO X, MOROZ I M. Ensemble Kalman filter with the unscented transform [J]. Physica D: Nonlinear Phenomena, 2009, 238(5): 549-562.