浅海水声环境往往呈现水平非均匀性，其中海底底质声学特性的水平变化是其主要表现之一。由于直接采样测量海底底质声学特性的方法过程繁琐、成本高昂，许多研究聚焦于利用声学信号反演沉积层声学参数的间接方法。一种典型间接方法是基于单波束测深仪的沉积层地声参数反演方法^[1]。它基于海底高频反向散射的复合粗糙度模型^[2]，给出小入射角高频声波海底反向散射声信号包络时域建模方法^[3]，并通过海底反向散射信号包络进行地声参数反演^[4]。

地声参数反演问题通常看成一个非线性优化问题，将各种全局优化算法(模拟退火、遗传算法等)应用于求代价函数最优解问题^[5]。此类方法在寻优过程中需反复调用正演模型计算，普遍存在计算开销大、收敛速度慢，以及易陷入局部最优、超参数调优困难等问题。考虑到在实际的海洋测试中，海底底质特性在局部区域随距离变化的过程较为缓慢，这为地声参数的序贯估计提供了可行性。

卡尔曼滤波是一种典型序贯估计算法^[7]，可对线性或非线性系统中状态变量和参数进行动态估计^[8]，利用卡尔曼滤波框架进行序贯估计的方法已被有效应用于水声目标测距^[9]、声速剖面反演^[10]等领域。卡尔曼滤波过程分为预测和校正2个步骤，可根据正演模型预测与实测数据之间的误差反馈，自适应地更新对状态变量的激励，从而校正模型参数的预测值。

本文提出了一种利用卡尔曼滤波对海底反向散射回波进行沉积层参数序贯估计的方法，将模型应用于2018年8月采集的海试数据，获得测线上海底底质类型平均粒径分布，并与全局优化算法中性能较优的模拟退火算法进行了计算效率和反演精度的对比。

1 高频海底反向散射包络模型

海底反向散射包络计算需要考虑发射声脉冲、声源级、接收灵敏度、指向性、换能器姿态、球面扩散、海底界面的反向散射和沉积物体积的不均匀性等大量物理特征。假设海底界面的粗糙度各向同性并且满足高斯分布，其功率谱可以写成Von Karman谱形式W(k)=w₂k_b^γ，其中k_b为波数，w₂为谱强度，γ=3.5为本文中模型的谱指数。根据Jackson等^[2]提出的散射模型，海底反向散射回波包络强度分为海底表层散射和海底体积散射2部分，单波束声呐接收到的随时间t变化的强度可表示为：

$ I(t)=I_{i}(t)+I_{v}(t) $

(1)

式中：I_i为海底粗糙界面散射包络；I_v为海底介质体积散射包络，假设两者相互独立。I_i是一个亥姆霍兹衍射积分方程的Kirchhoff近似解，考虑了单波束声呐系统的安放位置、指向性、信号特征、水体扩展和吸收等，可以很好地描述沉积层界面反向散射的物理过程，特别是在小入射角条件下。体积散射信号包络I_v考虑了界面粗糙度对体积散射的影响。根据Sternlicht等^[3]的模型，散射回波包络声压最终表示为：

$ p_{a}\left(t, M_{z}\right)=\sqrt{\rho_{w} v_{w} I\left(t, M_{z}\right)} $

(2)

式中：ρ_w和v_w分别为水层的密度和声速；M_z为海底底质平均粒径，其他参数根据经验公式^[14]计算。

2 沉积层参数估计 2.1 全局优化方法

地声参数反演问题通常看作是非线性数值优化问题，以待反演参数为决策变量，确立代价函数和优化的可行域，即决策变量的搜索区间，并给定其须满足的约束条件，通过全局优化方法进行迭代计算，最终得到模型参数的全局最优解。

利用包络反演海底声学和散射特性参数通常采用最小二乘均方误差函数作为优化问题中的代价函数^[5]：

$ \boldsymbol{E}_{k}(\hat{\boldsymbol{X}})=\frac{\sum\limits_{n}\left[p_{a}(n, \hat{\boldsymbol{X}})-p_{k}(n)\right]^{2}}{\sum\limits_{n}\left[p_{a}(n, \hat{\boldsymbol{X}})\right]^{2}} $

(3)

式中：k为测量的散射回波包络的索引值；$ {\boldsymbol{\hat X}} $表示模型参数；n表示离散的采样点；p_a为正演模型输出的散射回波包络；p_k为实际测量获得散射回波包络。代价函数值越低，表示匹配度越高。

常用的非线性全局优化算法包括模拟退火算法^[11]、遗传算法^[12]、粒子群算法^[13]等。利用此类算法对海底底质平均粒径M_z进行搜索，通过多次迭代使式(3)中的代价函数最小化。在与序贯估计进行对比分析时，本文采用模拟退火作为反演模型的优化算法。

2.2 无迹卡尔曼滤波

无迹卡尔曼滤波(unscented Kalman filter，UKF)是通过无迹变换(unscented transform，UT)产生一系列确定样本来逼近状态的后验概率密度，既没有线性化忽略高阶项的过程，也不用求解雅可比矩阵^[7]。综合考虑滤波性能和运算效率，UKF是3种方法中的首选。

基于隐马尔可夫模型^[15]，卡尔曼滤波首先建立状态空间方程。假设非平稳随机过程中X_k为系统状态，状态空间方程表示为：

$ \boldsymbol{X}_{k}=\boldsymbol{f}_{k}\left(\boldsymbol{X}_{k-1}, \boldsymbol{V}_{k}\right) $

(4)

$ \boldsymbol{Y}_{k}=\boldsymbol{h}_{k}\left(\boldsymbol{X}_{k-1}, \boldsymbol{W}_{k}\right) $

(5)

式中：X_k和Y_k的方程表达式又分别称为系统的状态方程和观测方程；V_k和W_k分别为模型中状态噪声与观测噪声，两者之间互不相关。利用无迹变换选取2n+1个采样点作为Sigma点集：

$ \boldsymbol{X}_{k-11 k-1}^{(i)}=\left[\hat{\boldsymbol{X}}_{k-1 \mid k-1} \quad \hat{\boldsymbol{X}}_{k-1 \mid k-1}+\sqrt{(n+\lambda) \boldsymbol{P}_{k-1 \mid k-1}} \quad \hat{\boldsymbol{X}}_{k-1 \mid k-1}-\sqrt{(n+\lambda) \boldsymbol{P}_{k-1 \mid k-1}}\right] $

(6)

采样点的权值ω⁽ⁱ⁾为：

$ \left\{\begin{array}{l} \omega_{m}^{(0)}=\frac{\lambda}{n+\lambda} \\ \omega_{c}^{(0)}=\frac{\lambda}{n+\lambda}+\left(1-a^{2}+\beta\right) \\ \omega_{m}^{(i)}=\omega_{c}^{(i)}=\frac{\lambda}{2(n+\lambda)}, i=1,2, \cdots, 2 n \end{array}\right. $

(7)

UKF在此求出一组Sigma点的预测，并对他们加权求均值，得到系统状态和协方差的预测值：

$ \boldsymbol{X}_{k | k-1}^{(i)}=\boldsymbol{f}_{k}\left(k, \boldsymbol{X}_{k-1| k-1}^{(i)}\right) $

(8)

$ \hat{\boldsymbol{X}}_{k \mid k-1}=\sum\limits_{i=0}^{2 n} \omega^{(i)} \boldsymbol{X}_{k \mid k-1}^{(i)} $

(9)

$ \boldsymbol{P}_{k \mid k-1}=\sum\limits_{i=0}^{2 n} \omega^{(i)}\left[\hat{\boldsymbol{X}}_{k \mid k-1}-\boldsymbol{X}_{k \mid k-1}^{(i)}\right]\left[\hat{\boldsymbol{X}}_{k \mid k-1}-\boldsymbol{X}_{k \mid k-1}^{(i)}\right]^{\mathrm{T}}+\boldsymbol{Q} $

(10)

根据预测值，再次利用UT变换产生新的Sigma点集，并代入观测方程得到预测的观测量：

$ \boldsymbol{X}_{k \mid k-1}^{(i)}=\left[\hat{\boldsymbol{X}}_{k \mid k-1} \quad \hat{\boldsymbol{X}}_{k \mid k-1}+\sqrt{(n+\lambda) \boldsymbol{P}_{k \mid k-1}} \quad \hat{\boldsymbol{X}}_{k \mid k-1}-\sqrt{(n+\lambda) \boldsymbol{P}_{k \mid k-1}}\right] $

(11)

$ \boldsymbol{Z}_{k \mid k-1}^{(i)}=h\left(\boldsymbol{X}_{k \mid k-1}^{(i)}\right) $

(12)

接着，通过加权求和得到系统预测的均值及协方差，并计算卡尔曼增益和系统状态的协方差更新：

$ {\boldsymbol{\widehat Z}_{k\mid k - 1}} = \sum\limits_{i = 0}^{2n} {{\omega ^{(i)}}} \boldsymbol{Z}_{k\mid k - 1}^{(i)} $

(13)

$ \boldsymbol{P}_{z_{k}{z}_k}=\sum\limits_{i=0}^{2 n} \omega^{(i)}\left[\boldsymbol{Z}_{k \mid k-1}^{(i)}-\hat{\boldsymbol{Z}}_{k \mid k-1}\right]\left[\boldsymbol{Z}_{k \mid k-1}^{(i)}-\hat{\boldsymbol{Z}}_{k \mid k-1}\right]^{\mathrm{T}}+\boldsymbol{R} $

(14)

$ \boldsymbol{P}_{x_{k} z_{k}}=\sum\limits_{i=0}^{2 n} \omega^{(i)}\left[\boldsymbol{X}_{k \mid k-1}^{(i)}-\hat{\boldsymbol{X}}_{k \mid k-1}\right]\left[\boldsymbol{Z}_{k \mid k-1}^{(i)}-\hat{\boldsymbol{Z}}_{k \mid k-1}\right]^{\mathrm{T}} $

(15)

$ \boldsymbol{K}_{k}=\boldsymbol{P}_{x_{k}{z}_k} \boldsymbol{P}_{x_{k}{z}_k}^{-1} $

(16)

$ \boldsymbol{P}_{k \mid k}=\boldsymbol{P}_{k \mid k-1}-\boldsymbol{K}_{k} \boldsymbol{P}_{x_{k}{z}_k} \boldsymbol{K}_{k}^{\mathrm{T}} $

(17)

最后对系统状态量的预测值进行校正：

$ \left\{\begin{array}{l} \hat{\boldsymbol{X}}_{k \mid k}=\hat{\boldsymbol{X}}_{k \mid k-1}+\Delta \hat{\boldsymbol{X}}_{k} \\ \Delta \hat{\boldsymbol{X}}_{k}=\boldsymbol{K}_{k} \cdot \boldsymbol{E}_{k}\left(\hat{\boldsymbol{Z}}_{k \mid k-1}\right) \end{array}\right. $

(18)

式中E_k为误差函数，这里采用式(3)定义的最小二乘均方误差。式(18)在预测值$ {{\boldsymbol{\hat X}}_{k\left| {k - 1} \right.}} $的基础上加上预测误差的缩放量得到当前状态的估计值，卡尔曼增益K_k决定了对误差调整的大小，它和预测误差共同决定了系统中状态参数的更新量。

3 海上实验 3.1 实验概况

本试验于2018年8月在黄海进行，通过走航式测量海底散射回波，主要目的是对沉积物声学特征遥测的可行性进行研究。走航试验开始于试验海区E₁站位点(36°00.180 4′N 121°36.804 0′E)，利用收发合置换能器发射25 kHz单频信号和线性调频信号，于次日到达E₂站位点，全程约为130 km。

本次试验采用的换能器声源级和接收灵敏度较高，25 kHz声源级为212.72 dB，接收灵敏度为-186.1 dB，收发合置换能器具有较好的指向性，如图 1所示。图 2给出了收发合置换能器上罗经监测的换能器姿态，在实验的后半部分随海况恶化其姿态有所偏离，在反向散射包络的建模中引入了对指向性的纠正。

通过对发射的线性调频信号进行处理得到换能器距离海底深度，如图 3所示。水深在80 km处有相对较小的波动，其他较为平稳。图 4给出了E₁和E₂站位测量的水文声速剖面，在15~20 m海深处存在明显的跃变层，在此声速剖面基础上计算信号的传播时间。

3.2 数据预处理

实验中发射了频率为25 kHz，时间脉冲长度为2 ms的单频信号。信号采样频率为500 kHz，同时记录了功率放大器监视信号和收发合置换能器信号。散射信号包络的获取是通过对滤波后的原始信号进行Hilbert变换并取绝对值获得。图 5给出了第40组100个信号的海底散射回波包络图，可以看出由于海底的深度变化和换能器的上下浮动等因素，接收信号的包络到达时刻变化很大。由于收发合置换能器具有较高的源级，接收信号信噪比较高，故对信号包络可采用最小门限法^[5]进行对齐。图 6给出了信号的第40组包络对齐结果。在走航过程中海底的散射不同获得的信号包络差别较大，对获取的散射回波包络每100个样本进行平均，共得到80组信号包络。

基于最小门限法对齐的信号，对获取的信号回波包络进行80%的重叠平均处理，获得整个E₁-E₂测线海底散射回波包络，如图 7所示。利用获得的信号的80组包络分别进行海底声学参数全局优化反演和序贯估计。

3.3 序贯估计和全局优化反演

图 8给出了海底底质平均粒径值M_z分别取3.56、4.70、7.40时高频海底反向散射包络模型预报的散射包络。当M_z增大时，单波束声呐接收到的声强度增大，散射包络的面积增大。在该模型的基础上可利用预处理后的80组散射信号包络分别进行海底声学参数全局反演和序贯估计。

首先，利用全局优化法分别对80组散射包络进行地声参数反演，所用算法为SA，搜索区间为[0, 10]，初始温度T₀=100，温度变化率ΔT=0.98，设多个终止条件，当10次迭代内目标函数平均变化率δ < 10^-3或代价函数值小于0.02时终止迭代，且最大迭代次数N_max=50。其次，利用序贯估计法对80组散射包络依序进行地声参数估计，所用算法为UKF，为获取相对准确的初始状态值，将第1组包络的全局优化反演结果作为UKF初始状态，UKF初始状态协方差设为单位矩阵I。序列长度设为80，对应信号的80组散射包络。由于地声参数在短距离内的状态变化可视为缓变过程，令式(4)中状态方程为X_k=X_k-1+V_k。其中V_k表示k时刻过程中的随机扰动，服从均值为零、协方差矩阵为Q=qI的高斯分布；式(5)测量方程设为Y_k=h_k(X_k-1)+W_k，h_k为式(3)给出的函数，W_k表示测量噪声，服从均值为零、方差为r的高斯分布。根据状态转移过程和测量特性将q设为0.01，r试设为1和2。

图 9(a)和图 9(b)分别给出了SA和UKF预测包络与实测包络的最小二乘均方误差随距离变化的针状图。由于SA在误差收敛到一定阈值时终止搜索，其预测精度较为稳定。UKF根据预测与实测数据之间的误差反馈，自适应地调整待估参数的激励，由图可见在拐点附近UKF可及时调整状态参数的更新量以减小预测误差，但在海底底质剧变区测量误差增大，这是因为状态转移方程假定动态模型为缓变过程。SA和UKF的最小二乘均方误差均控制在0.05以下，表明预测结果可信度尚好。

利用序贯估计的海底声学参数预测海底散射包络随距离如图 10所示。通过对比图 7可知预测的海底散射包络与实测结果对应较好。

图 11给出了全局优化法和序贯估计法对海底底质平均粒径M_z的预测与海底底质分类结果。从图中可以看出UKF与SA的底质分类结果基本一致，从大约60 km处海底底质平均粒径变小，表明海底底质由软变硬。测量噪声方差r=2时相比r=1时测量结果曲线更加平滑，说明r的增大可使测量噪声对模型的影响增大，使卡尔曼增益K_k减小，从而减弱系统的校正权值，降低了系统的瞬态响应和稳态值。根据Wentworth分类标准，获得海底底质平均粒径为极细粉砂到中粉砂范围。

图 12给出了试验海区位置的海底底质类型分布，E₁-E₂测线为试验走航测线，在试验起始端底质类型为黏土质粉砂，与反演得到的平均粒径一致，如图 11中左侧方框所示；并很快过渡到砂质粉砂，后半段测线为粉砂质砂，与图 11中的右侧方框一致。序贯估计结果与试验预报结果符合很好。

利用UKF进行一步估计时，由于UKF无需同EKF一样计算非线性函数的雅可比矩阵，其计算复杂度仅来源于UT变换中产生Sigma点集时和预测校正时所调用的正演函数计算。因此，序贯估计法所需的时间代价同全局优化法一样由正演程序的计算次数决定。对于序贯估计法，每组包络固定执行2L+1次正演模型计算，L为待估参数个数；对于全局优化法，每组包络需执行的正演次数由寻优过程的迭代次数N决定，且根据预定的终止条件有所不同，每组反演所需的迭代次数根据具体收敛情况、伪峰和鞍点个数而不定，可达数十次到上百次不等。在Intel Core i5-9400F 2.90Ghz处理器的PC机环境下应用Matlab计算，二者所用计算时间对比如表 1所示。UKF在单次估计中平均计算时间相比SA法缩短83.4%，总体时间开销缩短79.6%。

4 结论

1) 序贯估计预测的海底散射回波包络与实测包络符合较好，全局优化反演和序贯估计法预测的最小二乘均方误差均可控制在0.05以下，预测精度尚好。

2) 海底底质平均粒径的预测值与海底底质分类分布对比具有较高的可信度。序贯估计的结果显示，海底底质类型在试验起始端为黏土质粉砂，并很快过渡到砂质粉砂，后半段测线为粉砂质砂，与试验预报结果符合很好。

3) 序贯估计法在保证和全局优化法精度相似的情况下有效提升算法效率，参数估计消耗的平均时间显著缩短，具有良好的研究和应用价值。