基于改进Sage-Husa自适应滤波的ADCP船速基准精确获取

陈志高，柳锐，吴子豪，陈志平，吴自银

（1.东华理工大学 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室，江西 南昌 330013；2.东华理工大学 测绘与空间信息工程学院，江西 南昌 330013；3.自然资源部海底科学重点实验室，浙江 杭州 310012）

流速测量是开发和利用水（海洋）资源、开展水利（海洋）工程建设所面临的基础问题之一［1］。声学多普勒流速剖面仪（acoustic Doppler current profiler，ADCP）借助水跟踪（water tracking，WT）获得的相对流速剖面加上底跟踪（bottom track‐ing，BT）获得的船速即可实时获取三维流速剖面，具有速度快、精度高、不干扰流场等优势，被国际海委会定为4 种先进的海洋观测仪器之一［2］，在流速观测和流量测验中被广泛采用［3-6］。然而，受复杂测量环境影响，ADCP 在实际应用中仍存在底质流动或悬沙浓度较高时底跟踪船速失效，进而直接导致流速测量结果出现偏差或者数据丢失等问题，严重影响了测量结果的可靠性，限制了ADCP 的作业环境。

针对底质流动或高浓度悬沙情况下ADCP 底跟踪船速失效的情况，国内外学者一般采用全球定位系统（global positioning system，GPS）前后历元的GGA 定位信息计算出速度，并对底跟踪船速进行替换［7-9］。研究结果表明，采用GPS 船速替换方法可以较好地解决底跟踪船速失效的问题。然而，一方面，受岸边树木遮挡导致可用卫星数目较少、风浪导致船体姿态剧烈变化等影响，GPS 定位解易出现异常［10］，这将导致GPS 船速精度降低，进而引起流速计算出现偏差；另一方面，底跟踪船速精度通常优于1 cm/s，高于利用GPS 定位解计算的船速精度［11］，直接的船速替换将有损底跟踪数据有效时绝对流速的测量精度。因此，采用船速直接替换未能根本解决ADCP 船速基准的准确获取问题。

卡尔曼（Kalman）滤波在航迹估计、姿态融合解算等方面广泛使用［12-13］。Kalman 滤波模型包括状态方程和观测方程，其核心思想为预测+测量反馈。在Kalman 滤波计算过程中，误差协方差阵对滤波系统的精度起着重要的作用。因此，系统噪声协方差阵Q和观测误差协方差R的确定尤为关键。在实际的测量过程中，噪声的随机特性主要是受测量仪器以及周围环境的影响，误差协方差阵Q和观测误差协方差阵R不可能先验已知，因而单一固定的先验方差可能会严重地影响滤波精度。在针对不确定的误差协方差阵Q和R的算法中，Sage-Husa 算法具有一定的代表性［14］。该算法依据最小均方误差准则，利用观测数据在递推求解的同时，通过时变噪声统计估值器，实时估计和修正系统噪声和量测噪声的统计特性参数。该算法基本忽略了初值Q和R，这有助于实时估计噪声的统计特性。动态模型中的不确定性会影响滤波过程，并且GPS 异常解及底跟踪失效等异常数据将直接影响船速滤波结果，进而导致流量计算结果产生偏差。同时应用Sage-Husa 算法时，容易出现R失去正定性或Q失去非负定性的问题，使滤波发散。

基于上述问题，本文结合自适应抗差滤波原理提出一种改进的Sage-Husa 自适应滤波算法。通过状态不符值统计量构造自适应因子ak并基于丹麦法构建抗差等价权因子w，二者可分别抑制动力学模型扰动异常（波浪、潮汐对船速的扰动）和观测信息异常（GPS 和底跟踪观测异常值），提高滤波的自适应性和鲁棒性；为防止在更新时出现R失去正定性或Q失去非负定性导致滤波发散的问题，对Q和R增加Kalman 滤波的约束条件。通过以上2 步，实现底跟踪船速和GPS 船速的有机融合：既可在底质流动或高浓度悬沙情况下获取稳健的GPS 船速，又可保留正常情况下高精度的底跟踪船速，进而增强ADCP 船速基准的准确性与稳定性，提高ADCP 走航观测的普适性。

1 滤波模型

1.1 Kalman滤波

由ADCP 底跟踪船速VBT与外部传感器GPS 获得的船速VGPS组成的观测值矩阵可表示为：

式中：k为时间步长；T 为矩阵转置；下标E和N分别表示东分量和北分量。

状态向量为：

式中VE、VN分别是船速的东分量和北分量。

Kalman 滤波模型的观测方程和状态方程可以表示为：

式中：Hk为观测矩阵；Φk为状态转移矩阵。相应的，系统噪声协方差阵Q和量测噪声协方差阵R分别为：

式中：SP为真实的船舶动态位置的方差；σGPS、σBT为测量误差。

然而，在实际情况中，噪声的随机特性主要是受测量仪器以及周围环境的影响，Q阵和R阵不可能先验已知，因而单一固定的先验方差可能会严重地影响滤波精度。

1.2 Sage-Husa算法

在衍生于Kalman 滤波的自适应算法中，Sage-Husa 算法是针对不确定的Q阵和R阵的经典算法。该算法依据最小均方误差准则，利用观测数据在递推求解的同时，通过时变噪声统计估值器，实时估计和修正系统噪声和量测噪声的统计特性参数，达到降低模型误差、抑制滤波发散、提高滤波精度的目的［15］。具体的调节方式为：

式中：vk为观测残差向量；dk为加权因子；b为遗忘因子，一般取0.95～0.99。

然而，虽然通过时变噪声统计估值器可以自适应更新Q和R，提高Kalman 滤波的性能，但是动态模型中的不确定性会影响滤波过程；同时，应用Sage-Husa 算法时,容易出现R失去正定性或Q失去非负定性的问题。

1.3 改进Sage-Husa算法

针对上述Sage-Husa算法的不足，本文根据自适应抗差滤波理论，一方面，构建自适应因子ak和抗差等价权因子w，抑制动力学模型扰动异常和观测信息异常，以提高滤波的自适应性和鲁棒性；另一方面，对Q和R矩阵增加约束条件，以保证Sage-Husa算法的稳定性。本文改进算法的基本流程如图1所示。

图1 改进的Sage-Husa算法流程Fig.1 Flow chart of improved Sage-Husa algorithm

1.3.1 基本流程

依然假设量测噪声Vk与系统噪声Wk是统计特性已知的零均值高斯白噪声且互不相关，Vk和Wk均与初始状态X0无关，改进算法的递推方程如下。

预测方程为：

预测误差协方差阵为：

观测残差向量为：

滤波增益矩阵为：

滤波解为：

协方差阵更新为：

系统噪声协方差阵Q阵和量测噪声协方差阵R阵更新为：

式中：ak为自适应因子；Pt为观测信息的抗差等价权矩阵；w为抗差等价权因子，其对角元素wi可采用观测残差判别统计量并基于丹麦法进行构造［16］。

与常规Sage-Husa 算法不同，本文改进算法根据自适应抗差滤波原理，首先基于丹麦法构造抗差因子w，计算观测信息的抗差等价权矩阵Pt，当观测信息中含有粗差时，抗差估计通过抗差因子w淘汰粗差；然后，根据状态不符值统计量，基于三段函数法构造自适应因子ak，自适应因子有助于减少预测状态误差模型中的不确定性［17］；最后，为了防止在更新时，Rk+1失去正定性或Qk+1失去非负定性的问题，对Qk+1和Rk+1增加卡尔曼滤波的约束条件，将式（14）和（15）的等式右边的第2 项处理结果的对角线元素取绝对值，非对角线元素取零［14］。

本文重点介绍自适应因子ak和抗差等价权因子w的具体构造过程。

1.3.2 自适应因子ak的构造

对于自适应因子的构造，杨元喜等［18］给出了4种构造自适应因子ak的误差判别统计量，包括状态不符值统计量、预测残差统计量、方差分量比统计量和速度不符值统计量。为避免观测信息不可靠时对ak构造的影响，本文选择基于状态不符值统计量并采用三段函数法构建ak，表示为［19-20］：

1.3.3 抗差等价权因子w的构造

ADCP底跟踪船速失效及GPS定位解易出现异常均会直接导致流速测量结果出现偏差，进而影响流量计算结果的可靠性。因此需要运用抗差估计原理抵制粗差干扰。丹麦法构建抗差等价权因子可避免发生矩阵奇异，当残差的可信度处在模糊界限时，该方法能够较有效地抑制粗差对滤波精度的影响。丹麦法抗差等价权函数为［21］：

式中：wi为第i个观测向量的等价权因子；vi为第i个观测向量的标准化残差；v（ki）和分别为残差的第i个分量和其标准差矩阵的第i个分量；c2一般取值为1.5～2.0。

2 船速滤波实验结果及分析

2.1 实验数据

为验证各种滤波方法的有效性，在长江口水域（图2）开展了走航式ADCP 流速测量实验。实验水域为潮汐河段，底质稳定；流速测量采用RDI公司生产的300 kHz 瑞江系列4 波束ADCP，Leica SR530用于导航定位，THALES 3011 GPS 罗经为其提供外部方位（定向精度优于0.5°），并采用实时动态量测技术（real time kinematic，RTK）获取船速，其精度优于±0.05 m/s。实验设计的走航断面CD如图2所示，其长度约为1 200 m；对CD 断面进行往返观测，共2个测回，4个测次。

图2 位于长江口的实验区域Fig.2 Located in the experimental area of the Yangtze River estuary

为测试不同底质流动速度对观测流速及船速滤波结果的影响，实验将模拟的底质流动速度叠加到实测的底跟踪船速中。模拟的底质流动速度包括移动偏差（moving bed，MB）与随机噪声2 部分。移动偏差的强度fMB由以下函数确定［22-23］：

式中：t为观测时刻；T为观测总时间；P为每个观测点相对于断面中心点的位置；e 为自然常数；F为最大强度；fMB为每个历元的强度。由式（20）可以生成一个在断面边缘附近为0，然后沿断面航迹逐渐增大，直至断面中心达到最大值的移动偏差序列。最后，在每个历元中加入N（0，5 cm/s）的高斯白噪声（随机噪声），即可得到最终模拟的底质流动速度。如图3 所示，实验共模拟了移动偏差fMB为5 cm/s（MB5），10 cm/s（MB10）和15 cm/s（MB15）的3 种不同强度底质流动速度。

图3 实测流速与模拟的3种底质流动速度Fig.3 Measured water velocity and simulated three kinds of moving bed velocity

2.2 抗差性能比较

由于测量船在走航过程中存在障碍物遮挡或GPS 卫星星座几何形状较差的情况，因此GPS（RTK）定位信息可能存在粗差。如图4 所示，为了展现不同的滤波方法的抗差性能，实验在fMB=5 cm/s的条件下，对第1测次的GPS数据在131～133 s、151～155 s 之间在东、北分量加入大小为0.5 m/s 的粗差。由图4 可知，Kalman 滤波、Sage-Husa 算法受粗差影响较大，导致滤波后船速曲线产生明显“隆起”，所以需要在滤波前事先将其剔除；而本文改进算法得到的滤波结果受粗差影响很小，所以在滤波前可以不作滤波剔除工作。

图4 加入粗差后不同滤波方法的抗差结果Fig.4 Results of different methods after adding gross er‐ror

2.3 滤波精度分析

图5 为在模拟移动偏差fMB为5 cm/s 的条件下，由标准Kalman 滤波、Sage-Husa 算法和本文改进Sage-Husa 算法这3 种滤波方法获得的不同测次的船速结果。由图5 可知，与Kalman 滤波和Sage-Hu‐sa 算法所得结果相比，本文方法的滤波结果更接近用于参考的GPS 船速，表明在存在底质流动的情况下，采用本文改进滤波方法可以提高船速精度。

图5 不同滤波方法得到的船速Fig.5 Ship velocity obtained by different filtering methods

图6 所示为在fMB为5 cm/s 的条件下，由3 种滤波方法（下标AKF表示Kalman滤波，下标SHA表示Sage-Husa 算法，下标ISH 表示本文改进方法）获得的船速与GPS 船速之间的东、北分量的均方根误差（root mean squared error，RMSE）计算公式为：

图6 不同滤波方法的船速东、北分量误差Fig.6 Error of ship velocity east and north components with different filtering methods

式中：Vfit为滤波后的船速；Vdata为实测的GPS 船速；n为该测次的总历元。

从图6 可以看出，本文方法与Kalman 滤波和Sage-Husa算法相比，其东、北分量误差曲线更接近0；其RMSE 在4 个测次取均值，可以求得本文方法在东分量的RMSE 为7.6 cm/s，而Kalman 滤波和Sage-Husa 算法分别为10.0 cm/s 和9.0 cm/s，分别相对提高了2.4 cm/s 和1.4 cm/s；本文方法在北分量的RMSE 为8.6 cm/s，而Kalman 滤波和Sage-Husa算法分别为12.1 cm/s和9.3 cm/s，分别提高了3.5 cm/s和0.7 cm/s。

为进一步验证不同底质流动情况下3 种滤波方法的效果，实验采用矢量标准差σ和均方根误差RMSE 作为精度指标进行评价。其中，，σe和σn分别为东和北分量的标准差；RMSE计算公式与式（21）一致。表1和表2分别为3种滤波方法在不同底质流动下的矢量标准差σ和均方根误差RMSE计算结果。由表1和表2可知：

表1 3种滤波方法下的船速矢量标准差σTable1 Standard deviation σ of boat velocity under 3 filtering methods cm/s

表2 3种滤波方法下的船速RMSETable 2 Root mean square error of velocity under 3 filtering methods cm/s

1）Sage-Husa 算法在底质流动MB5、MB10、MB15下的σ对比Kalman滤波分别提高了1.6 cm/s、1.4 cm/s 和1.2 cm/s，RMSE 则分别提高了1.1 cm/s、1.0 cm/s 和0.7 cm/s。究其原因，受波浪、风、泥沙浓度等影响ADCP 走航测量环境较为复杂，单一固定的Q阵和R阵会影响滤波效果，而Sage-Husa算法通过时变噪声统计估值器可以自适应更新Q和R，提高Kalman滤波的性能。

2）本文改进方法在底质流动MB5、MB10、MB15下的σ相对Sage-Husa算法又分别提高了1.0 cm/s、1.2 cm/s 和1.5 cm/s，RMSE 则分别提高了1.0 cm/s、1.0 cm/s 和1.3 cm/s。本文改进方法相对Sage-Husa算法精度得到了进一步提升，主要原因是本文改进方法不仅通过时变噪声统计估值器更新Q和R，还通过构造自适应因子控制动力学模型误差的影响。

2.4 流量计算

由ADCP 流量测量原理可知，走航式ADCP 测量所得流量是由各个ping（历元）获得的微断面流量累加得到，断面总流量计算公式为［24］：

式中：Qm为实测流量；Vx/y，WT为ADCP 水跟踪获得的相对流速的东/北分量；Vx y,BV为GPS 或底跟踪获得的船速的东/北分量；N为微断面总个数；i为第微断面序号；hi为第i个微断面的深度；Δti为相邻微断面i和i-1之间的采样时间间隔。

图7所示为模拟底质流动速度MB为5 cm/s的条件下，利用实测GPS 船速与3 种滤波方法得到的船速计算得到的绝对流速及其对应的流量计算结果。图7 的各个子图中，从左至右分别为以实测GPS 速度、传统Kalman滤波、Sage-Husa滤波以及本文改进方法滤波后的速度作为船速参考计算的绝对流速，相应的流量计算结果分别记为QGPS、QAKF、QSHA及QISH。流量单位为m³/s，负号表示涨潮，正号表示落潮。为便于比较，表3 还列出了fMB为10 cm/s 和15 cm/s 的条件下，由GPS 船速与3 种滤波方法船速作为参考计算出来的流量。由图7及表3可得出以下结论：

表3 3种滤波方法下计算的断面流量Table 3 Section discharge calculation under three filtering methods m3·s-1

图7 不同滤波方法的流量计算结果Fig.7 Discharge calculation of different filtering methods

1）随着底质流动速度fMB的增大，计算得到的断面流量结果变小。究其原因，由于底质流动速度与实际流速方向基本一致，受底质流动影响，实测的底跟踪船速产生和底质流动速度一样的矢量偏差，而计算流量的式（22）实际上是由船速矢量和流速矢量相乘转化而来，因此该公式计算的流量结果必然偏小；

2）随着底质流动速度的增大，本文改进算法相对于传统Kalman 滤波方法其流量计算精度明显提高。在fMB为5 cm/s的条件下，本文方法的精度相对于Kalman 滤波和Sage-Husa 算法分别提升了0.9%和0.6%，在fMB为10 cm/s 的条件下分别提升了3.0%和0.4%，在fMB为15 cm/s 的条件下则分别提升了6.0%和1.2%。尽管本文方法对船速精度的提高并不明显（约为2 cm/s），但对于年均流量为2.8×104m3/s的长江口这种大型感潮河段来说，本文方法1%～6%的流量相对精度的提高显得必要且有意义。

3 结论

1）传统Kalman滤波结果受粗差影响较大，需事先将粗差剔除；本文改进算法引入的抗差因子具有较强的抗差性，滤波前可不作粗差剔除工作。

2）本文改进算法相较Kalman 滤波方法和Sage-Husa算法，其船速滤波精度提高了1.0～2.7 cm/s。主要原因是本文改进方法不仅通过时变噪声统计估值器更新误差协方差阵Q和R，还通过构造自适应因子控制动力学模型误差的影响。

3）本文方法的流量计算精度相对于传统Kal‐man滤波和Sage-Husa算法在底质流动速度为5 cm/s、10 cm/s 与15 cm/s 的情况下分别提升了0.9%～3.0%和0.4%～1.2%。尽管本文方法对船速精度的提高并不明显，但对于年均流量为2.8×104m3/s的长江口这种大型感潮河段来说，本文方法1%～6%的流量相对精度的提高显得必要且有意义。