基于多傳感器魯棒容積卡爾曼融合的半潛式起重平臺立柱快速壓載艙液位估計

肖業(yè)方，郭蘊華，高海波，胡 義

（武漢理工大學(xué)a.高性能艦船技術(shù)教育部重點實驗室；b.能源與動力工程學(xué)院，武漢 430063）

0 引 言

隨著世界上油井和海工平臺迎來新的報廢高峰期，半潛式起重平臺作為一種重要的海洋拆卸平臺得到了快速發(fā)展[1]。半潛式起重平臺的作業(yè)要求特殊，起吊過程中會產(chǎn)生巨大的傾覆力矩。為保持平臺穩(wěn)性，平臺排水量會急劇增加數(shù)千噸，此時需要采用空氣排載系統(tǒng)?？諝馀泡d系統(tǒng)利用高壓空氣將艙室內(nèi)的液體吹除，其排載速度可達常規(guī)離心泵的數(shù)十倍[2-5]。同時，正因為空氣排載系統(tǒng)的排載速度較快，其立柱快速壓載艙液位估計的可靠性和穩(wěn)定性對于半潛式起重平臺的安全作業(yè)就極為重要。為此，可在立柱快速壓載艙中配置多個液位傳感器，并通過濾波算法和信息融合技術(shù)提高液位估計的可靠性和穩(wěn)定性。

已經(jīng)有學(xué)者對濾波算法進行了深入研究，先后提出了卡爾曼濾波（Kalman Filter，KF）、拓展卡爾曼濾波（Extended Kalman Filtering，EKF）[6]、無跡卡爾曼濾波（Unscented Kalman Filter，UKF）[7]、容積卡爾曼濾波（Cubature Kalman Filters，CKF）[8]等方法。近年來，信息融合技術(shù)因為具有可擴展性、低通信量、對節(jié)點故障魯棒性強等優(yōu)點受到了越來越多的關(guān)注。Olfati-Saber等[9-11]針對離散線性高斯系統(tǒng)，利用平均一致性方法提出了分布式卡爾曼濾波（Distributed Kalman Filter，DKF）；Li等[12]基于一致性理論采用統(tǒng)計線性誤差傳播的方法，通過引入偽觀測矩陣提出了分布式UKF算法；丁家琳等[13]利用信息濾波和平均一致性理論提出一種分布式CKF 算法。近年來不少學(xué)者將這些技術(shù)用于液位測量系統(tǒng)中，丁雨淋等[14]以實時水位觀測為對象，通過水文變化語義知識約束，實現(xiàn)了卡爾曼濾波（Kalman Filter，KF）參數(shù)的自適應(yīng)調(diào)整，提高了水位測量的精度；劉經(jīng)宇等[15]以鍋爐的汽包水位測量為對象，將KF 算法用于單通道的汽包水位濾波，再將各通道濾波后的值進行求平均值方法融合，得到了更準(zhǔn)確的汽包水位；劉慶華等[16]以船舶高精度液位測量為對象，將EKF 算法用于船舶液位測量中，并將多傳感器濾波后的值進行加權(quán)融合，得到了更可靠的液位值。

美國“智能維護系統(tǒng)中心”研究表明：自動化系統(tǒng)40%以上的故障報警是由于傳感器系統(tǒng)自身的故障而產(chǎn)生的誤報[17]。傳感器的自身故障會導(dǎo)致觀測野值，進一步導(dǎo)致濾波估計的性能嚴(yán)重退化，并由錯誤的估計觸發(fā)錯誤的報警，嚴(yán)重干擾系統(tǒng)的運行和人員的判斷。為此，Huber 等[18]提出了廣義極大似然估計（M 估計），降低了觀測野值的權(quán)重。Tseng 等[19]將基于Huber 估計的CKF 濾波算法應(yīng)用于GPS 導(dǎo)航系統(tǒng)。彭美康等[20]在CKF 算法中引入了基于平方根逼近函數(shù)的修正因子以抑制野值的影響。這些魯棒估計方法對于含有較少野值的情形較為有效，一旦某一傳感器因為機械故障或者接觸不良處于半失效狀態(tài)，會導(dǎo)致觀測中存在高比例的野值。此時，已有融合算法的估計性能會嚴(yán)重退化，甚至出現(xiàn)融合估計反而不如單傳感器估計的情況。并且，液位估計方程不存在線性統(tǒng)計誤差傳播問題，現(xiàn)有的非線性分布式融合算法不適用于液位估計系統(tǒng)。針對以上問題，本文提出了一種多傳感器魯棒容積卡爾曼融合算法（MSRCKF），并將其應(yīng)用于半潛式起重平臺立柱快速壓載艙的液位估計中。在單一傳感器半失效的情況下，即使含大量野值，系統(tǒng)通過MSRCKF算法仍可以得到穩(wěn)定可靠的液位估計。

1 空氣排載數(shù)學(xué)模型

1.1 狀態(tài)方程

半潛起重平臺立柱快速壓載艙壓縮空氣排水過程分為兩個步驟，第一步是艙內(nèi)建壓，此時海底閥的閥門外處壓力大于內(nèi)處壓力而處于關(guān)閉狀態(tài)，艙內(nèi)開始充入高壓氣體；當(dāng)海底閥的閥門內(nèi)外處壓力相等時進入第二步，此時海底閥打開，艙內(nèi)持續(xù)充入高壓氣體，艙內(nèi)液體在高壓的作用下從海底閥排出?；诶硐霘怏w狀態(tài)方程與伯努利方程對其進行分析，以艙室狀態(tài)的真實變化建立系統(tǒng)的狀態(tài)方程：

式中，xk為k時刻下系統(tǒng)的狀態(tài)向量，f(xk-1)為k-1時刻的真實值對k時刻的預(yù)測，wk-1為k-1時刻系統(tǒng)的過程噪聲，Γ為過程噪聲分布矩陣。

xk是由艙內(nèi)水質(zhì)量（單位：t）、艙內(nèi)空氣體積（單位：m3）、艙內(nèi)空氣進入質(zhì)量（單位：kg）、艙內(nèi)氣壓（單位：Pa）、艙內(nèi)外氣壓差所對應(yīng)的水柱高度（單位：m）、艙內(nèi)液位（單位：m）、海底閥閥門處水流速（單位：m?s-1）和海底閥閥門處水體積流量（單位：m3?s-1）構(gòu)成的8維狀態(tài)向量，即系統(tǒng)的狀態(tài)向量維數(shù)n為8，xk可表示為

f(xk-1)的映射關(guān)系可用以下公式表示：

式中：Δt為時間步長，單位為s；ρwater為海水密度，單位為kg?m-3；QAir為充入壓縮氣體的質(zhì)量流量，單位為kg?s-1；TAir為艙內(nèi)氣體的開爾文溫度，單位為K；RAir為氣體常數(shù)，單位為J?(mol·K)-1；M為空氣摩爾質(zhì)量，單位為kg·mol-1；Ppri為艙內(nèi)初始?xì)鈮?，單位為Pa；Patm為標(biāo)準(zhǔn)大氣壓，單位為Pa；g為重力加速度，單位為m?s-2；fcap(?)為艙容表函數(shù)，由艙室的艙容表決定；H0為外界液位離艙底的距離，單位為m；kwater為海水阻力系數(shù)；S為海底閥的閥門總面積，單位為m2；Nwater為海底閥數(shù)量；r為海底閥閥門半徑，單位為m。

1.2 量測方程

以液位傳感器的測量建立系統(tǒng)的量測方程，即

式中，zk為k時刻下艙內(nèi)液位傳感器的測量值，H為量測轉(zhuǎn)換矩陣，vk為k時刻系統(tǒng)的測量噪聲。

vk與wk均為高斯白噪聲且互不相關(guān)，它們的統(tǒng)計特性為

2 魯棒CKF算法

2.1 CKF算法

CKF 算法的核心思想是通過球面徑向容積準(zhǔn)則進行相應(yīng)的數(shù)值積分運算，是當(dāng)前最接近貝葉斯濾波的近似算法，主要有時間更新和量測更新兩個步驟。假設(shè)已知系統(tǒng)在k-1 時刻的后驗密度函數(shù)為p（xk-1|zk-1）=N（x?k-1|k-1,Pk-1|k-1），CKF算法計算流程如下：

Step1——時間更新

對協(xié)方差矩陣Pk-1|k-1進行Cholesky分解得到Sk-1|k-1，即

式中，ξj為第j個容積點，j=1，2，…，2n，n為系統(tǒng)的狀態(tài)向量維數(shù)，In為n維單位矩陣，[?]j表示矩陣[?]的第j列。

通過非線性狀態(tài)函數(shù)對cubature點進行外推預(yù)測，有

2.2 Huber-CKF算法

Huber-CKF 算法與標(biāo)準(zhǔn)CKF 算法的區(qū)別主要在于引入了Huber 等價權(quán)函數(shù)，將濾波增益K中的Rk重構(gòu)為

式中，Ψz=diag[ψ(ei)]為穩(wěn)健因子，ψ(ei)為等價權(quán)函數(shù)，即

式中，ei為(zk-z?k k-1)(加權(quán)觀測殘差)的第i個分量，β為閾值。

2.3 魯棒CKF算法

對線性量測方程構(gòu)建線性回歸模型，即

式中，χk=yk-gkxk，為標(biāo)準(zhǔn)化殘差。且函數(shù)ρ(χk)的表達式為

事實上，因為真值xk不可知，在計算中用外推預(yù)測值x?k|k-1代替真值，故式（41）中δxk|k-1的值為0。同時，式（41）中的vk也可以視為觀測殘差，即vk=zk-z?k|k-1。為了與觀測殘差的經(jīng)典表示方法一致，令ek=vk，有

觀察式（42）可知，Φk是一個基于觀測殘差的無量綱數(shù)。構(gòu)建矩陣Ω為

式中，m為觀測向量維數(shù)。

根據(jù)矩陣Ω對誤差協(xié)方差矩陣進行重構(gòu)，有

由式（42）和式（43）可知，矩陣Ω只對協(xié)方差矩陣進行修正，故G?k也可表示為

式中，λk可視為k時刻下的一個自適應(yīng)因子，且λk=θ-1k。

從推導(dǎo)可知，魯棒CKF（R-CKF）算法的本質(zhì)為對觀測向量的誤差協(xié)方差進行修正，從而使損失函數(shù)最小，然后得到系統(tǒng)狀態(tài)的最佳估計。

R-CKF 算法借鑒了Huber等價權(quán)函數(shù)的思想，但又與之有著明顯區(qū)別：Huber-CKF 算法是對觀測值的每一維度分別進行檢測，當(dāng)某一維度含有野值時，降低這一維度的權(quán)重從而降低野值的影響;而R-CKF算法是根據(jù)觀測殘差的所有維度進行綜合評判，然后通過自適應(yīng)因子λk修正Rk。

當(dāng)傳感器測量值中的野值含量較高時，進行故障隔離可以保障算法的穩(wěn)定性和精確性。Huber-CKF 算法觀測值某一維度含有野值時，無法判斷其對整體的影響程度，沒法進行整體的取舍，而RCKF 算法可以使用自適應(yīng)因子λk作為綜合判據(jù)來判斷野值對整體的影響程度，當(dāng)λk大于閾值β2時，λk→∞，即對這一次的觀測值進行隔離。故自適應(yīng)因子λk的表達式為

對應(yīng)到CKF 算法，R-CKF 算法只在新息協(xié)方差矩陣?yán)锛尤肓俗赃m應(yīng)因子λk，令Rλ,k=λkRk，將式（24）修正為

當(dāng)含有大野值時，Rλ,k→∞導(dǎo)致Pzz,k|k-1→∞，則有Kk=0，從而對這一次的觀測值進行隔離。

3 魯棒容積卡爾曼融合算法

假定有N個傳感器同步對目標(biāo)進行測量，則k時刻傳感器i的量測方程為

類似地，將式（54）～（55）代入式（49），可得

所以，式（26）濾波增益可重寫為

式（68）～（69）中，Rf,λ,k=λf,kRf,k，且λf,k= diag(λ1,k,λ2,k,???,λN,k)。觀察魯棒因子λ的定義可知，對于融合中心來說，λi,k中觀測殘差減去的項應(yīng)該為融合中心的z?f,k|k-1，但由于融合系統(tǒng)不帶反饋，子系統(tǒng)得不到融合中心的信息，故這里用z?i,k|k-1來近似。

下面將式（68）～（69）進行變形使它可以利用各傳感器的局部預(yù)測與局部更新和中心預(yù)測x?k|k-1來表達。即融合中心不直接使用量測來進行全局估計，而僅僅聯(lián)合局部估計。將式（66）乘以式（67），有

整理可得

4 仿真實驗

4.1 仿真條件

本文的研究對象為招商重工（江蘇）有限公司建造的CM-163 半潛式起重拆解平臺。選取半潛式起重拆解平臺的單機起吊工中0～400 t起吊工況下的CSBT-MP-2艙室為仿真對象。CSBT-MP-2艙室的初始參數(shù)：艙室總?cè)莘e為2249 m3，排出水體積占總艙室容積的6%，艙內(nèi)溫度為293 K，外界液位離艙底的距離為12 m，海底閥閥門數(shù)量為2，半徑為0.495 m，空氣壓縮機質(zhì)量流量為5.612 kg?s-1，艙內(nèi)初始?xì)鈮簽?01272 Pa，初始艙內(nèi)液位為8.7 m，艙內(nèi)原始水質(zhì)量為1394.38 t，海水密度為1.025×103kg?m-3，海水阻力系數(shù)為4.2。

采用蒙特卡洛方法對Huber-CKF 算法、R-CKF 算法、Huber-CKF 融合算法和MSRCKF 算法進行仿真實驗，其中Huber-CKF 算法的閾值β為1.345，R-CKF 算法的閾值β1為0.98，β2為5。為了保證實驗的合理性，各算例的初始條件相同，蒙特卡洛次數(shù)為200次，時間步長為1 s，液位傳感器的量測噪聲標(biāo)準(zhǔn)差為0.05 m。在上述條件下，使用2 個液位傳感器對艙內(nèi)液位進行測量，其中1 號傳感器因為間歇性的機械故障處于半失效狀態(tài)，2號傳感器處于正常工作狀態(tài)。使用各算法的濾波估計誤差值來評價其性能。

算例1：在量測過程中，1號傳感器的觀測噪聲污染率為0.2，2號傳感器的觀測噪聲污染率為0，并將被污染觀測的觀測噪聲的標(biāo)準(zhǔn)差變?yōu)檎Ｖ档?0倍。

算例2：在量測過程中，1號傳感器的觀測噪聲污染率為0.4，2號傳感器的觀測噪聲污染率為0，并將被污染觀測的觀測噪聲的標(biāo)準(zhǔn)差變?yōu)檎Ｖ档?0倍。

4.2 結(jié)果分析

算例1～2 的液位誤差如圖1～2 所示，表1 給出了算例1 和算例2 下各算法的液位誤差平均值。其中，Huber-CKF-1和Huber-CKF-2分別對應(yīng)于1號傳感器和2號傳感器的Huber-CKF 算法；R-CKF-1和R-CKF-2分別對應(yīng)于1號傳感器和2號傳感器的R-CKF 算法；Huber-CKF 融合對應(yīng)于基于Huber-CKF的融合算法；MSRCKF對應(yīng)于本文提出的多傳感器魯棒容積卡爾曼融合算法。

表1 各算法的平均液位誤差Tab.1 Average level error of algorithms

圖1 算例1液位誤差圖Fig.1 Level error diagram of Example 1

從以上仿真實驗可以看出:

（1）如圖1和圖2所示，在觀測值中含有不同比例的野值的情況下，各算法都能隨著采樣次數(shù)的增加趨于收斂。在圖1 和圖2 中，對于不含野值的2 號傳感器觀測值，Huber-CKF 算法和R-CKF 算法的估計精度大致一樣；Huber-CKF 算法和R-CKF 算法都設(shè)有閾值，當(dāng)超過閾值時才會進行加權(quán)，在不含野值的情況下，大多不會超過其閾值，故兩者估計精度大致一樣。在圖1和圖2中，隨著1號傳感器觀測值中野值比例的提高，R-CKF算法較Huber-CKF算法在估計精度上的優(yōu)勢越來越明顯。

圖2 算例2液位誤差圖Fig.2 Level error diagram of Example 2

（2）如圖1 所示，在1 號傳感器觀測值中野值比例為0.2 時，Huber-CKF 融合與MSRCKF 這兩種融合算法均能取得比不含野值的2 號傳感器更小的估計誤差，這是因為融合中心獲得了比單個傳感器更多的信息量。如圖2 所示，隨著1 號傳感器觀測值中野值比例的提高，Huber-CKF 融合算法的估計誤差會增大，其估計誤差會比2號傳感器的估計誤差略大，這是因為其受到的大量不良數(shù)據(jù)的影響導(dǎo)致其性能出現(xiàn)退化，而MSRCKF算法由于其自適應(yīng)因子面對大的野值時會采取數(shù)據(jù)隔離，避免了不良數(shù)據(jù)的影響，從而得到比2號傳感器更小的估計誤差。

（3）表1 給出了2 個算例中各算法的所有采樣點的平均估計誤差。從平均估計誤差來看，MS?RCKF算法在2個算例中的平均估計誤差均是各算法中最小的，其性能更優(yōu)。

（4）如圖1 和圖2 所示，所有算法的誤差曲線在采樣50 次附近均出現(xiàn)了一個大的波動，這是由于運動方程中的狀態(tài)矢量發(fā)生較大改變引起的。此時，排載過程從第一步的艙內(nèi)建壓變?yōu)榈诙降拇蜷_海底閥吹除艙內(nèi)液體?？傮w而言，本文所提方法較Huber方法對狀態(tài)跳變更加不敏感。

由于在實際工程應(yīng)用中，量測噪聲基本不會服從純高斯分布，也難以保證某個傳感器不因發(fā)生故障而產(chǎn)生觀測野值，故使用2個或2個以上的傳感器通過MSRCKF算法可以使系統(tǒng)得到更加穩(wěn)定可靠的估計值。

5 結(jié) 語

為提高半潛式起重平臺快速壓載艙的液位估計的可靠性和穩(wěn)定性，本文提出了一種基于多傳感器魯棒容積卡爾曼融合的液位估計算法MSRCKF。針對故障傳感器測量含有觀測野值的問題，借鑒Huber估計，構(gòu)造了一個自適應(yīng)因子來抑制野值的影響，提高了系統(tǒng)的精度和魯棒性，并通過基于擴展信息濾波的矩陣變換完成了MSRCKF 算法的推導(dǎo)。MSRCKF 算法是一種分布式融合算法，它對融合中心的通信容量要求低，且工程上易于實現(xiàn)。仿真實驗證明，MSRCKF 算法在各種算例條件下，均能得到更加穩(wěn)定可靠的液位估計，具有一定的工程應(yīng)用價值。