靳　鵬,楊理踐,高松巍

(沈陽工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院, 沈陽 110870)

?

管道慣性測(cè)繪原始數(shù)據(jù)迭代去噪算法

靳鵬,楊理踐,高松巍

(沈陽工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院, 沈陽 110870)

管道慣性測(cè)繪內(nèi)檢測(cè)通常以搭載在清管器(PIG)內(nèi)部的三維正交陀螺儀、加速度計(jì)，再結(jié)合里程計(jì)數(shù)據(jù)，解算出管道地理坐標(biāo)軌跡。通常使用迭代卡爾曼算法進(jìn)行原始數(shù)據(jù)去噪計(jì)算，但迭代的終止條件難以判斷。利用管道發(fā)球筒及初始段管道平直，PIG初始姿態(tài)可控、可測(cè)的特點(diǎn)，結(jié)合靜基座對(duì)準(zhǔn)效果的評(píng)估，判斷濾波效果，優(yōu)化去噪過程，使初始對(duì)準(zhǔn)及后續(xù)管道地理坐標(biāo)的解算精度明顯提高。

管道慣性測(cè)繪內(nèi)檢測(cè)；慣性測(cè)量單元；卡爾曼濾波；初始對(duì)準(zhǔn)

我國油氣主要通過管道運(yùn)輸，管道的完整性管理是確保管道安全和經(jīng)濟(jì)運(yùn)行的重要手段[1]?；趹T性技術(shù)的管道測(cè)繪內(nèi)檢測(cè)是當(dāng)前國際上主要的管道工程測(cè)繪方法。研究人員可在管道正常運(yùn)行狀態(tài)下，使用慣性測(cè)量器件(Inertial Measurement Unit,IMU)測(cè)繪管道的三維坐標(biāo)，如果以地面為參考點(diǎn)進(jìn)行修正，能精確描繪管道中心線三維走向圖[2]。系統(tǒng)主要采用捷聯(lián)慣性導(dǎo)航系統(tǒng)與里程計(jì)組成組合導(dǎo)航系統(tǒng)進(jìn)行航位推算和定位[3]。

管道組合導(dǎo)航系統(tǒng)采用Kalman算法對(duì)IMU原始數(shù)據(jù)去噪。并采用迭代濾波的方法，但過分濾波有可能導(dǎo)致原始數(shù)據(jù)失真，因此何時(shí)終止迭代成為新的研究方向。

在管道慣性測(cè)繪內(nèi)檢測(cè)的工程應(yīng)用中，有很多因素有助于解決問題：① 發(fā)球時(shí)要求PIG(清管器)保持較長時(shí)間靜止，從統(tǒng)計(jì)角度上可看成穩(wěn)態(tài)大樣本，能夠以實(shí)測(cè)數(shù)據(jù)通過統(tǒng)計(jì)方法計(jì)算出當(dāng)前試驗(yàn)狀態(tài)下的兩個(gè)協(xié)方差；其次，靜置狀態(tài)下的加速度計(jì)和陀螺儀通常是已知的。② 發(fā)球筒及初始管道平直，靜置其中的IMU的載體PIG的姿態(tài)可控、可測(cè)。針對(duì)油氣管道發(fā)球工藝，利用其靜置時(shí)間長，靜置狀態(tài)下傳感器輸出數(shù)據(jù)部分可知，靜置狀態(tài)下IMU載體姿態(tài)部分具有可控、可測(cè)，初始運(yùn)動(dòng)段管道平直等特點(diǎn)，為了對(duì)陀螺儀和加速度計(jì)原始數(shù)據(jù)進(jìn)行降噪處理，設(shè)計(jì)了迭代Kalman算法及終止迭代的判斷方法。試驗(yàn)表明，隨著迭代濾波的進(jìn)行，初始對(duì)準(zhǔn)會(huì)經(jīng)歷一個(gè)不準(zhǔn)到準(zhǔn)，再到不準(zhǔn)的過程。據(jù)此可對(duì)迭代過程進(jìn)行控制，對(duì)提高濾波器濾波效果有直接的幫助。

1　迭代濾波過程及濾波效果的監(jiān)控

1.1Kalman濾波

定義一個(gè)離散時(shí)間過程的狀態(tài)變量x∈Rn，xk是x在k時(shí)刻的n維狀態(tài)，用差分方程的方法進(jìn)行表示如下

(1)

式中：u為控制向量；令B為0;A為差分方程的增益矩陣，表示前一時(shí)刻到后一時(shí)刻的系統(tǒng)狀態(tài)轉(zhuǎn)移方陣。

(2)式中：H為測(cè)量系統(tǒng)的參數(shù)；wk和vk分別為系統(tǒng)噪聲和觀測(cè)噪聲。

假設(shè)為白噪聲，得遞歸Kalman方程為

(3)

(4)

(5)

(6)

(7)

式中：xk,k-1為利用k-1狀態(tài)進(jìn)行估計(jì)得到的當(dāng)前k狀態(tài)(非最優(yōu)估計(jì))；xk-1為k-1狀態(tài)計(jì)算得到的最優(yōu)估計(jì)，xk定義與其類似；Pk,k-1為xk,k-1的估計(jì)誤差協(xié)方差；Pk-1為xk-1的估計(jì)誤差協(xié)方差(已經(jīng)最小化，達(dá)到最優(yōu)估計(jì))，Pk定義與其類似；Kgk為Kalman增益；Zk為對(duì)x實(shí)測(cè)得到的k時(shí)刻數(shù)據(jù)；Q為狀態(tài)噪聲協(xié)方差陣;R為觀測(cè)噪聲協(xié)方差陣，可通過wk和vk統(tǒng)計(jì)得到。

輸入初始狀態(tài)，并在每次迭代計(jì)算時(shí)輸入當(dāng)前待濾波的測(cè)量數(shù)據(jù)Zk，根據(jù)式(3)(7)可以求得每個(gè)時(shí)刻的最優(yōu)估計(jì)xk。

1.2觀測(cè)噪聲的提取

先處理vk，靜置時(shí)式(2)中xk=xsk,常數(shù)矩陣xsk，令xk為加速度矢量，不考慮擾動(dòng)重力,則

(8)

式中：g為試驗(yàn)現(xiàn)場(chǎng)的重力加速度。

對(duì)于SINS(捷聯(lián)慣性導(dǎo)航系統(tǒng))來說，此時(shí)H也是常數(shù)陣，是靜置時(shí)的姿態(tài)矩陣(通過靜基座初始對(duì)準(zhǔn)算法獲得)，記為Hs，則Hsxsk為重力加速度在三個(gè)坐標(biāo)軸上的投影，而此時(shí)測(cè)得的靜置觀測(cè)數(shù)據(jù)Zsk為三個(gè)坐標(biāo)軸上加速度計(jì)的讀數(shù)，代入式(2)有：

(9)

當(dāng)xk為角速度矢量，則有

(10)

式中：ωie為地球自轉(zhuǎn)角速度值；φ為當(dāng)?shù)鼐暥龋部梢酝ㄟ^式(9)計(jì)算得到vk。

1.3狀態(tài)噪聲的提取和迭代去除

從式(1)(忽略u(píng))可知，當(dāng)靜置狀態(tài)時(shí)，有xsk=xsk-1，代入式(8)或式(10)，有

(11)

但xsk的真值應(yīng)該為式(8)或式(10)的常量，則無法用式(11)求出xsk的wk，自然也無法統(tǒng)計(jì)出Q。如果使用觀測(cè)值代替xsk的真值參加計(jì)算，則直接引入了觀測(cè)噪聲，將觀測(cè)噪聲和狀態(tài)噪聲混雜在一起計(jì)算，使式(11)的精度受到影響，此時(shí)可以使用濾波后的更高精度的觀測(cè)值代替xsk的真值，并形成迭代計(jì)算，從而進(jìn)一步減少觀測(cè)噪聲的影響，以求得滿足工程需要的wk和Q。算法步驟如下：

(1) 使用儀器標(biāo)定的參數(shù)計(jì)算得到當(dāng)前的Q0和R0及其他初始條件，利用加速度計(jì)和陀螺儀數(shù)據(jù)進(jìn)行靜基座初始對(duì)準(zhǔn)，獲得Hs。

(2) 使用式(8)(9)統(tǒng)計(jì)得到R。

(5) 如Vx和Vs同時(shí)滿足預(yù)設(shè)精度，則此時(shí)求得的Q可作為算法的輸出，否則把當(dāng)前Q引入下次迭代，并利用濾波后的加速度計(jì)和陀螺儀數(shù)據(jù)重新進(jìn)行靜基座初始對(duì)準(zhǔn)，獲得Hs，回到步驟(2)繼續(xù)。

1.4中止迭代過程的判定方法

步驟(5)中，多次迭代的濾波過程有可能導(dǎo)致信號(hào)失真。這種現(xiàn)象通過試驗(yàn)可以觀察到：迭代計(jì)算得到更小的Vx和Vs，使用相對(duì)應(yīng)的Q和R進(jìn)行濾波，有時(shí)卻導(dǎo)致對(duì)準(zhǔn)和位置解算失敗。

而評(píng)估初始對(duì)準(zhǔn)的仰俯角和方向角精度時(shí)，可以利用管道發(fā)球系統(tǒng)的工藝要求,即發(fā)球管道初始段相對(duì)平直的特點(diǎn)，使靜置時(shí)仰俯角β和滾轉(zhuǎn)角γ為趨于0的小量，建立下述模型。

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

由設(shè)定，上述公式各項(xiàng)中sinεz、sinεx、sinγ、sinβ均為小量，忽略包含兩個(gè)小量相乘的項(xiàng)，再考慮cosγ、cosεx和cosβ近似為1，上述公式化簡(jiǎn)為：

(20)

(21)

(22)

由式(13)，(14)，可計(jì)算b′系航向角α′和仰俯角β′，即當(dāng)前靜置數(shù)據(jù)對(duì)準(zhǔn)后航向角和仰俯角的真值為：

(23)

(24)

(25)

(26)

將式(23)和(24)分別代入式(25)和(26)，αi、αj、βi、βj是兩次對(duì)準(zhǔn)計(jì)算的結(jié)果，α′和β′在試驗(yàn)和工程現(xiàn)場(chǎng)可測(cè)，求得σij>0，則第j次濾波得到的陀螺儀數(shù)據(jù)精度更高；τij>0，則第j次濾波得到的加速度計(jì)數(shù)據(jù)精度更高，反之亦然。

這樣，可在迭代中始終監(jiān)控陀螺儀和加速度計(jì)數(shù)據(jù)質(zhì)量，在幾次迭代計(jì)算中選擇數(shù)據(jù)質(zhì)量相對(duì)最好的一組，及時(shí)跳出迭代計(jì)算，試驗(yàn)數(shù)據(jù)中優(yōu)化后的Kalman濾波參數(shù)Q和R是由這些數(shù)據(jù)統(tǒng)計(jì)得到的。

2　試驗(yàn)結(jié)果

2.1靜態(tài)數(shù)據(jù)的濾波效果

以最簡(jiǎn)單的加速度信號(hào)預(yù)處理為例，設(shè)計(jì)靜態(tài)數(shù)據(jù)處理試驗(yàn)。采用車載PIG構(gòu)成試驗(yàn)平臺(tái)，不考慮安裝誤差，靜置前調(diào)整IMU滾轉(zhuǎn)角γ近似為0，然后通過調(diào)平系統(tǒng)保證β′為0(對(duì)γ有微小影響，但不影響其值為一個(gè)小量)，并確定試驗(yàn)用IMU系統(tǒng)軸線的方向角α′(令其為0，不影響算法的理論證明過程)，然后靜置。根據(jù)IMU參數(shù)計(jì)算得到Q0為9.610 9×10-9和R0為3.459 9×10-7。因?yàn)槭菍?duì)數(shù)據(jù)進(jìn)行簡(jiǎn)單的初始降噪，令A(yù)為單位陣，B為0，Hs由初始對(duì)準(zhǔn)模塊求出。

此外，以Z軸加速度計(jì)的局部數(shù)據(jù)(第4組數(shù)據(jù)中的前8 000個(gè))為例，圖1給出了使用Q0和R0以及Q=0.002 5和R=0.024 8的濾波效果對(duì)比。

表1　使用Q0、R0及優(yōu)化后Q、R的Kalman濾波效果對(duì)比

圖1　Z軸加速度計(jì)第4組部分試驗(yàn)數(shù)據(jù)優(yōu)化濾波的效果對(duì)比

圖2　載車運(yùn)動(dòng)軌跡及其解算效果示意

2.2同等條件下軌跡解算的質(zhì)量對(duì)比

圖2是載車運(yùn)動(dòng)軌跡及其解算效果的圖示。試驗(yàn)過程如下：先延靜置方向角α′的直行軌跡，長度為79.9 m，而后轉(zhuǎn)向繼續(xù)運(yùn)動(dòng)，全程1 854 m(里程計(jì)數(shù)據(jù))回到出發(fā)點(diǎn)。因采用相對(duì)坐標(biāo)解算試驗(yàn)軌跡，設(shè)靜止?fàn)顟B(tài)的α′=0，則圖中初始段軌跡應(yīng)該是垂直指向正上方的直線，以此驗(yàn)證對(duì)準(zhǔn)的效果。又因?yàn)檎麄€(gè)運(yùn)動(dòng)軌跡理論上應(yīng)該封閉，通過觀查軌跡未封閉區(qū)間的大小，可直觀反映SINS位置解算的精度。

圖2中包括過渡濾波數(shù)據(jù)的解算軌跡，以及優(yōu)化濾波前后數(shù)據(jù)的解算軌跡。三次計(jì)算除了數(shù)據(jù)預(yù)處理部分，對(duì)準(zhǔn)和軌跡解算的算法完全一樣?？梢灾庇^地看到，過度濾波得到的軌跡從對(duì)準(zhǔn)開始就失真，使整個(gè)軌跡有明顯順時(shí)針旋轉(zhuǎn)的趨勢(shì)，未封閉區(qū)間達(dá)到72.08 m；而使用Q0、R0濾波數(shù)據(jù)計(jì)算的軌跡，因剩余噪聲的積累越來越大，初期有一定精度，后期明顯失真，未封閉區(qū)間67.85 m；優(yōu)化算法后得到Q=0.002 2，R=0.024 9，濾波后數(shù)據(jù)解算的軌跡最接近真值，未封閉區(qū)間35.64 m。證明其對(duì)原始數(shù)據(jù)的濾噪效果更好，算法對(duì)提高系統(tǒng)位置解算精度有一定幫助。

3　結(jié)語

(1) 靜置狀態(tài)下該算法可顯著去除陀螺儀和加速度計(jì)噪聲，使加速度計(jì)測(cè)得的重力加速度更接近真值。

(2) 通過評(píng)估靜基座對(duì)準(zhǔn)精度，該算法可防止過度濾波導(dǎo)致信號(hào)損失。

(3) 算法可提高管道慣性測(cè)繪內(nèi)檢測(cè)系統(tǒng)的初始對(duì)準(zhǔn)精度和位置解算精度。

(4) 算法獲得了更加優(yōu)化的Kalman濾波參數(shù)Q和R，實(shí)質(zhì)上提高了IMU器件的標(biāo)稱精度。

[1]王富祥，馮慶善，楊建新，等.油氣管道慣性測(cè)繪內(nèi)檢測(cè)及其應(yīng)用[J].油氣儲(chǔ)運(yùn)，2012，31(5)：372-375.

[2]岳步江，唐雅琴，張永江，等.組合導(dǎo)航技術(shù)在油氣管道測(cè)繪系統(tǒng)中的應(yīng)用[J]. 中國慣性技術(shù)學(xué)報(bào)，2010，16(6)：712-716.

[3]EUN H S，NASER E S.Navigation Kalman filter design for pipeline pigging[J].The Journal of Navigation，2005，58: 283-289.

[4]譚紅力，黃新生，岳冬雪. 等效轉(zhuǎn)動(dòng)矢量在捷聯(lián)慣導(dǎo)系統(tǒng)二位置對(duì)準(zhǔn)中的應(yīng)用[J].中國慣性技術(shù)學(xué)報(bào), 2007,15(3)：269-272.

The Arithmetic of Pipeline Geometry Mapping Original Data Iterative Denoising

JIN Peng, YANG Li-jian, GAO Song-wei

(School of Information Science and Engineering, Shenyang University of Technology, Shenyang 110870, China)

Orthogonal gyroscopes, accelerometers, with odometer, are usually mounted inside the PIG for mapping by pipeline geometry PIG. The solution calculates the geographical coordinates of the pipeline track. An iterative Kalman algorithm was usually used to restrain the original data noise, but iteration termination conditions were difficult to judge. In this paper, because the initial segment of the pipeline and launch tube are straight, and PIG initial attitude is controllable and measurable, combined with the stationary base alignment effect evaluation, the filter performance can be determined and the de-noising process can be optimized, so that the initial alignment and subsequent pipeline geographical coordinates calculation accuracy has been improved significantly.

Pipeline geometry PIG; Inertial measurement unit; Kalman filter; Initial alignment

2015-09-08

“十二五”國家科技部支撐計(jì)劃資助項(xiàng)目(2011BAK06B01-03)；國家高技術(shù)研究發(fā)展計(jì)劃“863”資助項(xiàng)目(2012AA040104)

靳鵬(1974-)，男，講師，主要從事長輸油管道內(nèi)檢測(cè)技術(shù)及相關(guān)理論、無損檢測(cè)方面的工作。

10.11973/wsjc201603005

TG115.28

A

1000-6656(2016)03-0014-04