基于改進Sage-Husa自適應濾波的ADCP船速基準精確獲取

陳志高，柳銳，吳子豪，陳志平，吳自銀

（1.東華理工大學 自然資源部環(huán)鄱陽湖區(qū)域礦山環(huán)境監(jiān)測與治理重點實驗室，江西 南昌 330013；2.東華理工大學 測繪與空間信息工程學院，江西 南昌 330013；3.自然資源部海底科學重點實驗室，浙江 杭州 310012）

流速測量是開發(fā)和利用水（海洋）資源、開展水利（海洋）工程建設(shè)所面臨的基礎(chǔ)問題之一［1］。聲學多普勒流速剖面儀（acoustic Doppler current profiler，ADCP）借助水跟蹤（water tracking，WT）獲得的相對流速剖面加上底跟蹤（bottom track‐ing，BT）獲得的船速即可實時獲取三維流速剖面，具有速度快、精度高、不干擾流場等優(yōu)勢，被國際海委會定為4 種先進的海洋觀測儀器之一［2］，在流速觀測和流量測驗中被廣泛采用［3-6］。然而，受復雜測量環(huán)境影響，ADCP 在實際應用中仍存在底質(zhì)流動或懸沙濃度較高時底跟蹤船速失效，進而直接導致流速測量結(jié)果出現(xiàn)偏差或者數(shù)據(jù)丟失等問題，嚴重影響了測量結(jié)果的可靠性，限制了ADCP 的作業(yè)環(huán)境。

針對底質(zhì)流動或高濃度懸沙情況下ADCP 底跟蹤船速失效的情況，國內(nèi)外學者一般采用全球定位系統(tǒng)（global positioning system，GPS）前后歷元的GGA 定位信息計算出速度，并對底跟蹤船速進行替換［7-9］。研究結(jié)果表明，采用GPS 船速替換方法可以較好地解決底跟蹤船速失效的問題。然而，一方面，受岸邊樹木遮擋導致可用衛(wèi)星數(shù)目較少、風浪導致船體姿態(tài)劇烈變化等影響，GPS 定位解易出現(xiàn)異常［10］，這將導致GPS 船速精度降低，進而引起流速計算出現(xiàn)偏差；另一方面，底跟蹤船速精度通常優(yōu)于1 cm/s，高于利用GPS 定位解計算的船速精度［11］，直接的船速替換將有損底跟蹤數(shù)據(jù)有效時絕對流速的測量精度。因此，采用船速直接替換未能根本解決ADCP 船速基準的準確獲取問題。

卡爾曼（Kalman）濾波在航跡估計、姿態(tài)融合解算等方面廣泛使用［12-13］。Kalman 濾波模型包括狀態(tài)方程和觀測方程，其核心思想為預測+測量反饋。在Kalman 濾波計算過程中，誤差協(xié)方差陣對濾波系統(tǒng)的精度起著重要的作用。因此，系統(tǒng)噪聲協(xié)方差陣Q和觀測誤差協(xié)方差R的確定尤為關(guān)鍵。在實際的測量過程中，噪聲的隨機特性主要是受測量儀器以及周圍環(huán)境的影響，誤差協(xié)方差陣Q和觀測誤差協(xié)方差陣R不可能先驗已知，因而單一固定的先驗方差可能會嚴重地影響濾波精度。在針對不確定的誤差協(xié)方差陣Q和R的算法中，Sage-Husa 算法具有一定的代表性［14］。該算法依據(jù)最小均方誤差準則，利用觀測數(shù)據(jù)在遞推求解的同時，通過時變噪聲統(tǒng)計估值器，實時估計和修正系統(tǒng)噪聲和量測噪聲的統(tǒng)計特性參數(shù)。該算法基本忽略了初值Q和R，這有助于實時估計噪聲的統(tǒng)計特性。動態(tài)模型中的不確定性會影響濾波過程，并且GPS 異常解及底跟蹤失效等異常數(shù)據(jù)將直接影響船速濾波結(jié)果，進而導致流量計算結(jié)果產(chǎn)生偏差。同時應用Sage-Husa 算法時，容易出現(xiàn)R失去正定性或Q失去非負定性的問題，使濾波發(fā)散。

基于上述問題，本文結(jié)合自適應抗差濾波原理提出一種改進的Sage-Husa 自適應濾波算法。通過狀態(tài)不符值統(tǒng)計量構(gòu)造自適應因子ak并基于丹麥法構(gòu)建抗差等價權(quán)因子w，二者可分別抑制動力學模型擾動異常（波浪、潮汐對船速的擾動）和觀測信息異常（GPS 和底跟蹤觀測異常值），提高濾波的自適應性和魯棒性；為防止在更新時出現(xiàn)R失去正定性或Q失去非負定性導致濾波發(fā)散的問題，對Q和R增加Kalman 濾波的約束條件。通過以上2 步，實現(xiàn)底跟蹤船速和GPS 船速的有機融合：既可在底質(zhì)流動或高濃度懸沙情況下獲取穩(wěn)健的GPS 船速，又可保留正常情況下高精度的底跟蹤船速，進而增強ADCP 船速基準的準確性與穩(wěn)定性，提高ADCP 走航觀測的普適性。

1 濾波模型

1.1 Kalman濾波

由ADCP 底跟蹤船速VBT與外部傳感器GPS 獲得的船速VGPS組成的觀測值矩陣可表示為：

式中：k為時間步長；T 為矩陣轉(zhuǎn)置；下標E和N分別表示東分量和北分量。

狀態(tài)向量為：

式中VE、VN分別是船速的東分量和北分量。

Kalman 濾波模型的觀測方程和狀態(tài)方程可以表示為：

式中：Hk為觀測矩陣；Φk為狀態(tài)轉(zhuǎn)移矩陣。相應的，系統(tǒng)噪聲協(xié)方差陣Q和量測噪聲協(xié)方差陣R分別為：

式中：SP為真實的船舶動態(tài)位置的方差；σGPS、σBT為測量誤差。

然而，在實際情況中，噪聲的隨機特性主要是受測量儀器以及周圍環(huán)境的影響，Q陣和R陣不可能先驗已知，因而單一固定的先驗方差可能會嚴重地影響濾波精度。

1.2 Sage-Husa算法

在衍生于Kalman 濾波的自適應算法中，Sage-Husa 算法是針對不確定的Q陣和R陣的經(jīng)典算法。該算法依據(jù)最小均方誤差準則，利用觀測數(shù)據(jù)在遞推求解的同時，通過時變噪聲統(tǒng)計估值器，實時估計和修正系統(tǒng)噪聲和量測噪聲的統(tǒng)計特性參數(shù)，達到降低模型誤差、抑制濾波發(fā)散、提高濾波精度的目的［15］。具體的調(diào)節(jié)方式為：

式中：vk為觀測殘差向量；dk為加權(quán)因子；b為遺忘因子，一般取0.95～0.99。

然而，雖然通過時變噪聲統(tǒng)計估值器可以自適應更新Q和R，提高Kalman 濾波的性能，但是動態(tài)模型中的不確定性會影響濾波過程；同時，應用Sage-Husa 算法時,容易出現(xiàn)R失去正定性或Q失去非負定性的問題。

1.3 改進Sage-Husa算法

針對上述Sage-Husa算法的不足，本文根據(jù)自適應抗差濾波理論，一方面，構(gòu)建自適應因子ak和抗差等價權(quán)因子w，抑制動力學模型擾動異常和觀測信息異常，以提高濾波的自適應性和魯棒性；另一方面，對Q和R矩陣增加約束條件，以保證Sage-Husa算法的穩(wěn)定性。本文改進算法的基本流程如圖1所示。

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

1.3.1 基本流程

依然假設(shè)量測噪聲Vk與系統(tǒng)噪聲Wk是統(tǒng)計特性已知的零均值高斯白噪聲且互不相關(guān)，Vk和Wk均與初始狀態(tài)X0無關(guān)，改進算法的遞推方程如下。

預測方程為：

預測誤差協(xié)方差陣為：

觀測殘差向量為：

濾波增益矩陣為：

濾波解為：

協(xié)方差陣更新為：

系統(tǒng)噪聲協(xié)方差陣Q陣和量測噪聲協(xié)方差陣R陣更新為：

式中：ak為自適應因子；Pt為觀測信息的抗差等價權(quán)矩陣；w為抗差等價權(quán)因子，其對角元素wi可采用觀測殘差判別統(tǒng)計量并基于丹麥法進行構(gòu)造［16］。

與常規(guī)Sage-Husa 算法不同，本文改進算法根據(jù)自適應抗差濾波原理，首先基于丹麥法構(gòu)造抗差因子w，計算觀測信息的抗差等價權(quán)矩陣Pt，當觀測信息中含有粗差時，抗差估計通過抗差因子w淘汰粗差；然后，根據(jù)狀態(tài)不符值統(tǒng)計量，基于三段函數(shù)法構(gòu)造自適應因子ak，自適應因子有助于減少預測狀態(tài)誤差模型中的不確定性［17］；最后，為了防止在更新時，Rk+1失去正定性或Qk+1失去非負定性的問題，對Qk+1和Rk+1增加卡爾曼濾波的約束條件，將式（14）和（15）的等式右邊的第2 項處理結(jié)果的對角線元素取絕對值，非對角線元素取零［14］。

本文重點介紹自適應因子ak和抗差等價權(quán)因子w的具體構(gòu)造過程。

1.3.2 自適應因子ak的構(gòu)造

對于自適應因子的構(gòu)造，楊元喜等［18］給出了4種構(gòu)造自適應因子ak的誤差判別統(tǒng)計量，包括狀態(tài)不符值統(tǒng)計量、預測殘差統(tǒng)計量、方差分量比統(tǒng)計量和速度不符值統(tǒng)計量。為避免觀測信息不可靠時對ak構(gòu)造的影響，本文選擇基于狀態(tài)不符值統(tǒng)計量并采用三段函數(shù)法構(gòu)建ak，表示為［19-20］：

1.3.3 抗差等價權(quán)因子w的構(gòu)造

ADCP底跟蹤船速失效及GPS定位解易出現(xiàn)異常均會直接導致流速測量結(jié)果出現(xiàn)偏差，進而影響流量計算結(jié)果的可靠性。因此需要運用抗差估計原理抵制粗差干擾。丹麥法構(gòu)建抗差等價權(quán)因子可避免發(fā)生矩陣奇異，當殘差的可信度處在模糊界限時，該方法能夠較有效地抑制粗差對濾波精度的影響。丹麥法抗差等價權(quán)函數(shù)為［21］：

式中：wi為第i個觀測向量的等價權(quán)因子；vi為第i個觀測向量的標準化殘差；v（ki）和分別為殘差的第i個分量和其標準差矩陣的第i個分量；c2一般取值為1.5～2.0。

2 船速濾波實驗結(jié)果及分析

2.1 實驗數(shù)據(jù)

為驗證各種濾波方法的有效性，在長江口水域（圖2）開展了走航式ADCP 流速測量實驗。實驗水域為潮汐河段，底質(zhì)穩(wěn)定；流速測量采用RDI公司生產(chǎn)的300 kHz 瑞江系列4 波束ADCP，Leica SR530用于導航定位，THALES 3011 GPS 羅經(jīng)為其提供外部方位（定向精度優(yōu)于0.5°），并采用實時動態(tài)量測技術(shù)（real time kinematic，RTK）獲取船速，其精度優(yōu)于±0.05 m/s。實驗設(shè)計的走航斷面CD如圖2所示，其長度約為1 200 m；對CD 斷面進行往返觀測，共2個測回，4個測次。

圖2 位于長江口的實驗區(qū)域Fig.2 Located in the experimental area of the Yangtze River estuary

為測試不同底質(zhì)流動速度對觀測流速及船速濾波結(jié)果的影響，實驗將模擬的底質(zhì)流動速度疊加到實測的底跟蹤船速中。模擬的底質(zhì)流動速度包括移動偏差（moving bed，MB）與隨機噪聲2 部分。移動偏差的強度fMB由以下函數(shù)確定［22-23］：

式中：t為觀測時刻；T為觀測總時間；P為每個觀測點相對于斷面中心點的位置；e 為自然常數(shù)；F為最大強度；fMB為每個歷元的強度。由式（20）可以生成一個在斷面邊緣附近為0，然后沿斷面航跡逐漸增大，直至斷面中心達到最大值的移動偏差序列。最后，在每個歷元中加入N（0，5 cm/s）的高斯白噪聲（隨機噪聲），即可得到最終模擬的底質(zhì)流動速度。如圖3 所示，實驗共模擬了移動偏差fMB為5 cm/s（MB5），10 cm/s（MB10）和15 cm/s（MB15）的3 種不同強度底質(zhì)流動速度。

圖3 實測流速與模擬的3種底質(zhì)流動速度Fig.3 Measured water velocity and simulated three kinds of moving bed velocity

2.2 抗差性能比較

由于測量船在走航過程中存在障礙物遮擋或GPS 衛(wèi)星星座幾何形狀較差的情況，因此GPS（RTK）定位信息可能存在粗差。如圖4 所示，為了展現(xiàn)不同的濾波方法的抗差性能，實驗在fMB=5 cm/s的條件下，對第1測次的GPS數(shù)據(jù)在131～133 s、151～155 s 之間在東、北分量加入大小為0.5 m/s 的粗差。由圖4 可知，Kalman 濾波、Sage-Husa 算法受粗差影響較大，導致濾波后船速曲線產(chǎn)生明顯“隆起”，所以需要在濾波前事先將其剔除；而本文改進算法得到的濾波結(jié)果受粗差影響很小，所以在濾波前可以不作濾波剔除工作。

圖4 加入粗差后不同濾波方法的抗差結(jié)果Fig.4 Results of different methods after adding gross er‐ror

2.3 濾波精度分析

圖5 為在模擬移動偏差fMB為5 cm/s 的條件下，由標準Kalman 濾波、Sage-Husa 算法和本文改進Sage-Husa 算法這3 種濾波方法獲得的不同測次的船速結(jié)果。由圖5 可知，與Kalman 濾波和Sage-Hu‐sa 算法所得結(jié)果相比，本文方法的濾波結(jié)果更接近用于參考的GPS 船速，表明在存在底質(zhì)流動的情況下，采用本文改進濾波方法可以提高船速精度。

圖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。

為進一步驗證不同底質(zhì)流動情況下3 種濾波方法的效果，實驗采用矢量標準差σ和均方根誤差RMSE 作為精度指標進行評價。其中，，σe和σn分別為東和北分量的標準差；RMSE計算公式與式（21）一致。表1和表2分別為3種濾波方法在不同底質(zhì)流動下的矢量標準差σ和均方根誤差RMSE計算結(jié)果。由表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 算法在底質(zhì)流動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 走航測量環(huán)境較為復雜，單一固定的Q陣和R陣會影響濾波效果，而Sage-Husa算法通過時變噪聲統(tǒng)計估值器可以自適應更新Q和R，提高Kalman濾波的性能。

2）本文改進方法在底質(zhì)流動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算法精度得到了進一步提升，主要原因是本文改進方法不僅通過時變噪聲統(tǒng)計估值器更新Q和R，還通過構(gòu)造自適應因子控制動力學模型誤差的影響。

2.4 流量計算

由ADCP 流量測量原理可知，走航式ADCP 測量所得流量是由各個ping（歷元）獲得的微斷面流量累加得到，斷面總流量計算公式為［24］：

式中：Qm為實測流量；Vx/y，WT為ADCP 水跟蹤獲得的相對流速的東/北分量；Vx y,BV為GPS 或底跟蹤獲得的船速的東/北分量；N為微斷面總個數(shù)；i為第微斷面序號；hi為第i個微斷面的深度；Δti為相鄰微斷面i和i-1之間的采樣時間間隔。

圖7所示為模擬底質(zhì)流動速度MB為5 cm/s的條件下，利用實測GPS 船速與3 種濾波方法得到的船速計算得到的絕對流速及其對應的流量計算結(jié)果。圖7 的各個子圖中，從左至右分別為以實測GPS 速度、傳統(tǒng)Kalman濾波、Sage-Husa濾波以及本文改進方法濾波后的速度作為船速參考計算的絕對流速，相應的流量計算結(jié)果分別記為QGPS、QAKF、QSHA及QISH。流量單位為m3/s，負號表示漲潮，正號表示落潮。為便于比較，表3 還列出了fMB為10 cm/s 和15 cm/s 的條件下，由GPS 船速與3 種濾波方法船速作為參考計算出來的流量。由圖7及表3可得出以下結(jié)論：

表3 3種濾波方法下計算的斷面流量Table 3 Section discharge calculation under three filtering methods m3·s-1

圖7 不同濾波方法的流量計算結(jié)果Fig.7 Discharge calculation of different filtering methods

1）隨著底質(zhì)流動速度fMB的增大，計算得到的斷面流量結(jié)果變小。究其原因，由于底質(zhì)流動速度與實際流速方向基本一致，受底質(zhì)流動影響，實測的底跟蹤船速產(chǎn)生和底質(zhì)流動速度一樣的矢量偏差，而計算流量的式（22）實際上是由船速矢量和流速矢量相乘轉(zhuǎn)化而來，因此該公式計算的流量結(jié)果必然偏??；

2）隨著底質(zhì)流動速度的增大，本文改進算法相對于傳統(tǒng)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 結(jié)論

1）傳統(tǒng)Kalman濾波結(jié)果受粗差影響較大，需事先將粗差剔除；本文改進算法引入的抗差因子具有較強的抗差性，濾波前可不作粗差剔除工作。

2）本文改進算法相較Kalman 濾波方法和Sage-Husa算法，其船速濾波精度提高了1.0～2.7 cm/s。主要原因是本文改進方法不僅通過時變噪聲統(tǒng)計估值器更新誤差協(xié)方差陣Q和R，還通過構(gòu)造自適應因子控制動力學模型誤差的影響。

3）本文方法的流量計算精度相對于傳統(tǒng)Kal‐man濾波和Sage-Husa算法在底質(zhì)流動速度為5 cm/s、10 cm/s 與15 cm/s 的情況下分別提升了0.9%～3.0%和0.4%～1.2%。盡管本文方法對船速精度的提高并不明顯，但對于年均流量為2.8×104m3/s的長江口這種大型感潮河段來說，本文方法1%～6%的流量相對精度的提高顯得必要且有意義。