肖業方,郭蘊華,高海波,胡 義
(武漢理工大學a.高性能艦船技術教育部重點實驗室;b.能源與動力工程學院,武漢 430063)
隨著世界上油井和海工平臺迎來新的報廢高峰期,半潛式起重平臺作為一種重要的海洋拆卸平臺得到了快速發展[1]。半潛式起重平臺的作業要求特殊,起吊過程中會產生巨大的傾覆力矩。為保持平臺穩性,平臺排水量會急劇增加數千噸,此時需要采用空氣排載系統??諝馀泡d系統利用高壓空氣將艙室內的液體吹除,其排載速度可達常規離心泵的數十倍[2-5]。同時,正因為空氣排載系統的排載速度較快,其立柱快速壓載艙液位估計的可靠性和穩定性對于半潛式起重平臺的安全作業就極為重要。為此,可在立柱快速壓載艙中配置多個液位傳感器,并通過濾波算法和信息融合技術提高液位估計的可靠性和穩定性。
已經有學者對濾波算法進行了深入研究,先后提出了卡爾曼濾波(Kalman Filter,KF)、拓展卡爾曼濾波(Extended Kalman Filtering,EKF)[6]、無跡卡爾曼濾波(Unscented Kalman Filter,UKF)[7]、容積卡爾曼濾波(Cubature Kalman Filters,CKF)[8]等方法。近年來,信息融合技術因為具有可擴展性、低通信量、對節點故障魯棒性強等優點受到了越來越多的關注。Olfati-Saber等[9-11]針對離散線性高斯系統,利用平均一致性方法提出了分布式卡爾曼濾波(Distributed Kalman Filter,DKF);Li等[12]基于一致性理論采用統計線性誤差傳播的方法,通過引入偽觀測矩陣提出了分布式UKF算法;丁家琳等[13]利用信息濾波和平均一致性理論提出一種分布式CKF 算法。近年來不少學者將這些技術用于液位測量系統中,丁雨淋等[14]以實時水位觀測為對象,通過水文變化語義知識約束,實現了卡爾曼濾波(Kalman Filter,KF)參數的自適應調整,提高了水位測量的精度;劉經宇等[15]以鍋爐的汽包水位測量為對象,將KF 算法用于單通道的汽包水位濾波,再將各通道濾波后的值進行求平均值方法融合,得到了更準確的汽包水位;劉慶華等[16]以船舶高精度液位測量為對象,將EKF 算法用于船舶液位測量中,并將多傳感器濾波后的值進行加權融合,得到了更可靠的液位值。
美國“智能維護系統中心”研究表明:自動化系統40%以上的故障報警是由于傳感器系統自身的故障而產生的誤報[17]。傳感器的自身故障會導致觀測野值,進一步導致濾波估計的性能嚴重退化,并由錯誤的估計觸發錯誤的報警,嚴重干擾系統的運行和人員的判斷。為此,Huber 等[18]提出了廣義極大似然估計(M 估計),降低了觀測野值的權重。Tseng 等[19]將基于Huber 估計的CKF 濾波算法應用于GPS 導航系統。彭美康等[20]在CKF 算法中引入了基于平方根逼近函數的修正因子以抑制野值的影響。這些魯棒估計方法對于含有較少野值的情形較為有效,一旦某一傳感器因為機械故障或者接觸不良處于半失效狀態,會導致觀測中存在高比例的野值。此時,已有融合算法的估計性能會嚴重退化,甚至出現融合估計反而不如單傳感器估計的情況。并且,液位估計方程不存在線性統計誤差傳播問題,現有的非線性分布式融合算法不適用于液位估計系統。針對以上問題,本文提出了一種多傳感器魯棒容積卡爾曼融合算法(MSRCKF),并將其應用于半潛式起重平臺立柱快速壓載艙的液位估計中。在單一傳感器半失效的情況下,即使含大量野值,系統通過MSRCKF算法仍可以得到穩定可靠的液位估計。
半潛起重平臺立柱快速壓載艙壓縮空氣排水過程分為兩個步驟,第一步是艙內建壓,此時海底閥的閥門外處壓力大于內處壓力而處于關閉狀態,艙內開始充入高壓氣體;當海底閥的閥門內外處壓力相等時進入第二步,此時海底閥打開,艙內持續充入高壓氣體,艙內液體在高壓的作用下從海底閥排出?;诶硐霘怏w狀態方程與伯努利方程對其進行分析,以艙室狀態的真實變化建立系統的狀態方程:

式中,xk為k時刻下系統的狀態向量,f(xk-1)為k-1時刻的真實值對k時刻的預測,wk-1為k-1時刻系統的過程噪聲,Γ為過程噪聲分布矩陣。
xk是由艙內水質量(單位:t)、艙內空氣體積(單位:m3)、艙內空氣進入質量(單位:kg)、艙內氣壓(單位:Pa)、艙內外氣壓差所對應的水柱高度(單位:m)、艙內液位(單位:m)、海底閥閥門處水流速(單位:m?s-1)和海底閥閥門處水體積流量(單位:m3?s-1)構成的8維狀態向量,即系統的狀態向量維數n為8,xk可表示為

f(xk-1)的映射關系可用以下公式表示:

式中:Δt為時間步長,單位為s;ρwater為海水密度,單位為kg?m-3;QAir為充入壓縮氣體的質量流量,單位為kg?s-1;TAir為艙內氣體的開爾文溫度,單位為K;RAir為氣體常數,單位為J?(mol·K)-1;M為空氣摩爾質量,單位為kg·mol-1;Ppri為艙內初始氣壓,單位為Pa;Patm為標準大氣壓,單位為Pa;g為重力加速度,單位為m?s-2;fcap(?)為艙容表函數,由艙室的艙容表決定;H0為外界液位離艙底的距離,單位為m;kwater為海水阻力系數;S為海底閥的閥門總面積,單位為m2;Nwater為海底閥數量;r為海底閥閥門半徑,單位為m。
以液位傳感器的測量建立系統的量測方程,即

式中,zk為k時刻下艙內液位傳感器的測量值,H為量測轉換矩陣,vk為k時刻系統的測量噪聲。
vk與wk均為高斯白噪聲且互不相關,它們的統計特性為

CKF 算法的核心思想是通過球面徑向容積準則進行相應的數值積分運算,是當前最接近貝葉斯濾波的近似算法,主要有時間更新和量測更新兩個步驟。假設已知系統在k-1 時刻的后驗密度函數為p(xk-1|zk-1)=N(x?k-1|k-1,Pk-1|k-1),CKF算法計算流程如下:
Step1——時間更新
對協方差矩陣Pk-1|k-1進行Cholesky分解得到Sk-1|k-1,即

式中,ξj為第j個容積點,j=1,2,…,2n,n為系統的狀態向量維數,In為n維單位矩陣,[?]j表示矩陣[?]的第j列。
通過非線性狀態函數對cubature點進行外推預測,有

Huber-CKF 算法與標準CKF 算法的區別主要在于引入了Huber 等價權函數,將濾波增益K中的Rk重構為

式中,Ψz=diag[ψ(ei)]為穩健因子,ψ(ei)為等價權函數,即

式中,ei為(zk-z?k k-1)(加權觀測殘差)的第i個分量,β為閾值。
對線性量測方程構建線性回歸模型,即

式中,χk=yk-gkxk,為標準化殘差。且函數ρ(χk)的表達式為

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

觀察式(42)可知,Φk是一個基于觀測殘差的無量綱數。構建矩陣Ω為

式中,m為觀測向量維數。
根據矩陣Ω對誤差協方差矩陣進行重構,有

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

式中,λk可視為k時刻下的一個自適應因子,且λk=θ-1k。
從推導可知,魯棒CKF(R-CKF)算法的本質為對觀測向量的誤差協方差進行修正,從而使損失函數最小,然后得到系統狀態的最佳估計。
R-CKF 算法借鑒了Huber等價權函數的思想,但又與之有著明顯區別:Huber-CKF 算法是對觀測值的每一維度分別進行檢測,當某一維度含有野值時,降低這一維度的權重從而降低野值的影響;而R-CKF算法是根據觀測殘差的所有維度進行綜合評判,然后通過自適應因子λk修正Rk。
當傳感器測量值中的野值含量較高時,進行故障隔離可以保障算法的穩定性和精確性。Huber-CKF 算法觀測值某一維度含有野值時,無法判斷其對整體的影響程度,沒法進行整體的取舍,而RCKF 算法可以使用自適應因子λk作為綜合判據來判斷野值對整體的影響程度,當λk大于閾值β2時,λk→∞,即對這一次的觀測值進行隔離。故自適應因子λk的表達式為

對應到CKF 算法,R-CKF 算法只在新息協方差矩陣里加入了自適應因子λk,令Rλ,k=λkRk,將式(24)修正為

當含有大野值時,Rλ,k→∞導致Pzz,k|k-1→∞,則有Kk=0,從而對這一次的觀測值進行隔離。
假定有N個傳感器同步對目標進行測量,則k時刻傳感器i的量測方程為


類似地,將式(54)~(55)代入式(49),可得

所以,式(26)濾波增益可重寫為

式(68)~(69)中,Rf,λ,k=λf,kRf,k,且λf,k= diag(λ1,k,λ2,k,???,λN,k)。觀察魯棒因子λ的定義可知,對于融合中心來說,λi,k中觀測殘差減去的項應該為融合中心的z?f,k|k-1,但由于融合系統不帶反饋,子系統得不到融合中心的信息,故這里用z?i,k|k-1來近似。
下面將式(68)~(69)進行變形使它可以利用各傳感器的局部預測與局部更新和中心預測x?k|k-1來表達。即融合中心不直接使用量測來進行全局估計,而僅僅聯合局部估計。將式(66)乘以式(67),有

整理可得

本文的研究對象為招商重工(江蘇)有限公司建造的CM-163 半潛式起重拆解平臺。選取半潛式起重拆解平臺的單機起吊工中0~400 t起吊工況下的CSBT-MP-2艙室為仿真對象。CSBT-MP-2艙室的初始參數:艙室總容積為2249 m3,排出水體積占總艙室容積的6%,艙內溫度為293 K,外界液位離艙底的距離為12 m,海底閥閥門數量為2,半徑為0.495 m,空氣壓縮機質量流量為5.612 kg?s-1,艙內初始氣壓為101272 Pa,初始艙內液位為8.7 m,艙內原始水質量為1394.38 t,海水密度為1.025×103kg?m-3,海水阻力系數為4.2。
采用蒙特卡洛方法對Huber-CKF 算法、R-CKF 算法、Huber-CKF 融合算法和MSRCKF 算法進行仿真實驗,其中Huber-CKF 算法的閾值β為1.345,R-CKF 算法的閾值β1為0.98,β2為5。為了保證實驗的合理性,各算例的初始條件相同,蒙特卡洛次數為200次,時間步長為1 s,液位傳感器的量測噪聲標準差為0.05 m。在上述條件下,使用2 個液位傳感器對艙內液位進行測量,其中1 號傳感器因為間歇性的機械故障處于半失效狀態,2號傳感器處于正常工作狀態。使用各算法的濾波估計誤差值來評價其性能。
算例1:在量測過程中,1號傳感器的觀測噪聲污染率為0.2,2號傳感器的觀測噪聲污染率為0,并將被污染觀測的觀測噪聲的標準差變為正常值的50倍。
算例2:在量測過程中,1號傳感器的觀測噪聲污染率為0.4,2號傳感器的觀測噪聲污染率為0,并將被污染觀測的觀測噪聲的標準差變為正常值的50倍。
算例1~2 的液位誤差如圖1~2 所示,表1 給出了算例1 和算例2 下各算法的液位誤差平均值。其中,Huber-CKF-1和Huber-CKF-2分別對應于1號傳感器和2號傳感器的Huber-CKF 算法;R-CKF-1和R-CKF-2分別對應于1號傳感器和2號傳感器的R-CKF 算法;Huber-CKF 融合對應于基于Huber-CKF的融合算法;MSRCKF對應于本文提出的多傳感器魯棒容積卡爾曼融合算法。

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

圖1 算例1液位誤差圖Fig.1 Level error diagram of Example 1
從以上仿真實驗可以看出:
(1)如圖1和圖2所示,在觀測值中含有不同比例的野值的情況下,各算法都能隨著采樣次數的增加趨于收斂。在圖1 和圖2 中,對于不含野值的2 號傳感器觀測值,Huber-CKF 算法和R-CKF 算法的估計精度大致一樣;Huber-CKF 算法和R-CKF 算法都設有閾值,當超過閾值時才會進行加權,在不含野值的情況下,大多不會超過其閾值,故兩者估計精度大致一樣。在圖1和圖2中,隨著1號傳感器觀測值中野值比例的提高,R-CKF算法較Huber-CKF算法在估計精度上的優勢越來越明顯。

圖2 算例2液位誤差圖Fig.2 Level error diagram of Example 2
(2)如圖1 所示,在1 號傳感器觀測值中野值比例為0.2 時,Huber-CKF 融合與MSRCKF 這兩種融合算法均能取得比不含野值的2 號傳感器更小的估計誤差,這是因為融合中心獲得了比單個傳感器更多的信息量。如圖2 所示,隨著1 號傳感器觀測值中野值比例的提高,Huber-CKF 融合算法的估計誤差會增大,其估計誤差會比2號傳感器的估計誤差略大,這是因為其受到的大量不良數據的影響導致其性能出現退化,而MSRCKF算法由于其自適應因子面對大的野值時會采取數據隔離,避免了不良數據的影響,從而得到比2號傳感器更小的估計誤差。
(3)表1 給出了2 個算例中各算法的所有采樣點的平均估計誤差。從平均估計誤差來看,MS?RCKF算法在2個算例中的平均估計誤差均是各算法中最小的,其性能更優。
(4)如圖1 和圖2 所示,所有算法的誤差曲線在采樣50 次附近均出現了一個大的波動,這是由于運動方程中的狀態矢量發生較大改變引起的。此時,排載過程從第一步的艙內建壓變為第二步的打開海底閥吹除艙內液體??傮w而言,本文所提方法較Huber方法對狀態跳變更加不敏感。
由于在實際工程應用中,量測噪聲基本不會服從純高斯分布,也難以保證某個傳感器不因發生故障而產生觀測野值,故使用2個或2個以上的傳感器通過MSRCKF算法可以使系統得到更加穩定可靠的估計值。
為提高半潛式起重平臺快速壓載艙的液位估計的可靠性和穩定性,本文提出了一種基于多傳感器魯棒容積卡爾曼融合的液位估計算法MSRCKF。針對故障傳感器測量含有觀測野值的問題,借鑒Huber估計,構造了一個自適應因子來抑制野值的影響,提高了系統的精度和魯棒性,并通過基于擴展信息濾波的矩陣變換完成了MSRCKF 算法的推導。MSRCKF 算法是一種分布式融合算法,它對融合中心的通信容量要求低,且工程上易于實現。仿真實驗證明,MSRCKF 算法在各種算例條件下,均能得到更加穩定可靠的液位估計,具有一定的工程應用價值。