基于改進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。流量單位為m3/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%的流量相對精度的提高顯得必要且有意義。