基于卡爾曼濾波的D－InSAR和水準(zhǔn)監(jiān)測數(shù)據(jù)融合方法研究

李懷展 查劍鋒 米麗倩

1 中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，徐州市大學(xué)路1號，221008

2 山東省國土測繪院，濟(jì)南市經(jīng)十路2301號，250102

D－InSAR 技 術(shù)是 沉 降 監(jiān) 測 的 重 要 手 段［1－4］，其監(jiān)測結(jié)果受多種誤差源影響，沉降最大值附近精度較低；而水準(zhǔn)測量雖然精度較高，但只能得到有限個監(jiān)測點(diǎn)的變形值。兩種數(shù)據(jù)具有很好的互補(bǔ)性，可嘗試應(yīng)用數(shù)據(jù)同化的方法將D－InSAR 和水準(zhǔn)監(jiān)測數(shù)據(jù)融合。數(shù)據(jù)同化是指在考慮數(shù)據(jù)時空分布以及觀測場和背景場誤差的基礎(chǔ)上，在數(shù)值模型的動態(tài)運(yùn)行過程中融合新的觀測方法，對陸地、大氣、海洋的動態(tài)模型軌跡進(jìn)行一定的約束，從而提高變量的估計(jì)精度。作為數(shù)據(jù)同化的一種新方法，集合卡爾曼濾波近10a來取得很大進(jìn)展［4－8］。本文以高等級公路高精度沉降變形監(jiān)測為例，提出一種基于集合卡爾曼濾波同化的DInSAR－水準(zhǔn)數(shù)據(jù)融合的思路。

1 高等級公路多源數(shù)據(jù)對比分析

鄒濟(jì)公路在南屯煤礦井田內(nèi)長約8 400 m，煤礦開采對其有顯著影響。為確保公路的安全運(yùn)行，必須開展采動區(qū)高等級公路變形監(jiān)測及預(yù)報(bào)。對鄒濟(jì)公路進(jìn)行了水準(zhǔn)測量和InSAR 監(jiān)測，干涉雷達(dá)影像對及參數(shù)如表1。采用“二軌法”處理雷達(dá)影像［9－10］，并與實(shí)測水準(zhǔn)數(shù)據(jù)進(jìn)行對比。各期D－InSAR與水準(zhǔn)累加結(jié)果對比如圖1。

從圖1可以看出，疊加后的D－InSAR 監(jiān)測結(jié)果與水準(zhǔn)觀測值走勢上十分接近，能夠從整體上反映公路沉降變化的趨勢。但兩種數(shù)據(jù)在數(shù)值上呈現(xiàn)出一定的差異，D－InSAR 監(jiān)測值在整體上要小于水準(zhǔn)實(shí)測值，且兩種數(shù)據(jù)呈現(xiàn)出“中間差別大、兩邊差別小”的現(xiàn)象，即距離盆地中心越近，DInSAR 與水準(zhǔn)實(shí)測值的差別越大。這主要同雷達(dá)影像的空間分辨率和波長有關(guān)。D－InSAR 最大形變監(jiān)測梯度公式為：

表1 干涉雷達(dá)影像對及參數(shù)Tab.1 Interferometric SAR images and their parameters

其中，λ為波長，ps為像元大小。對于波長為3 cm 的TerraSAR－X，經(jīng)多視后像元大小為9m，干涉像元可探測到的最大形變梯度為1.6 mm／m。兩次觀測間隔內(nèi)地表的下沉差越大，形變的梯度則越大，當(dāng)形變梯度超過最大可探測梯度時，這些地表沉降值差別過大的點(diǎn)便無法準(zhǔn)確地在地表形變圖上體現(xiàn)出來，使得沉降越大的區(qū)域D－InSAR的誤差越大。

圖1 D－InSAR與水準(zhǔn)測量數(shù)據(jù)對比Fig.1 Comparison of D－InSAR and leveling data

綜合以上分析可知，D－InSAR 監(jiān)測能周期性地獲取監(jiān)測區(qū)域的面元變形信息，可達(dá)到cm級精度，但與常規(guī)水準(zhǔn)監(jiān)測相比精度不高；而水準(zhǔn)測量數(shù)據(jù)雖然精度較高，能達(dá)到mm級精度，但只能得到點(diǎn)元變形信息。所以，在實(shí)際應(yīng)用中，應(yīng)探索把這兩種數(shù)據(jù)進(jìn)行有效融合，達(dá)到數(shù)據(jù)的高空間分辨率與高精度的統(tǒng)一。

2 數(shù)據(jù)融合模型的建立

同化的本意是把不同的事物變得相近或相同［11］。設(shè)A為模型模擬結(jié)果，B為真實(shí)觀測值，則同化可有以下5層含義：第一，A不變，B變得和A相近或相同；第二，A基本不變，B向A同化；第三，A和B以相互靠攏的方式變得相近或相同；第四，B基本不變，A向B同化；第五，B不變，A變得和B相近或相同。如果以多源數(shù)據(jù)的融合集成利用為同化目標(biāo)，則應(yīng)該采用第三種同化方案。

考慮到基于點(diǎn)觀測的水準(zhǔn)實(shí)測數(shù)據(jù)比雷達(dá)影像的處理結(jié)果更為可靠、數(shù)值也更為準(zhǔn)確，利用集合卡爾曼濾波進(jìn)行衛(wèi)星影像監(jiān)測與水準(zhǔn)實(shí)測數(shù)據(jù)的融合時，采用以水準(zhǔn)實(shí)測數(shù)據(jù)（觀測）為主的同化方式。

由雷達(dá)影像處理得到的D－InSAR 監(jiān)測值反映的是區(qū)域面元信息，而路面水準(zhǔn)實(shí)測值獲取的是有限個監(jiān)測點(diǎn)元的變形，要將兩種數(shù)據(jù)進(jìn)行融合，必須對數(shù)據(jù)進(jìn)行有效的插值，實(shí)現(xiàn)二者的空間匹配。先對水準(zhǔn)數(shù)據(jù)的監(jiān)測點(diǎn)結(jié)果進(jìn)行插值，得到沉降實(shí)測估計(jì)值在整個區(qū)域的變形，再與區(qū)域的D－InSAR數(shù)據(jù)進(jìn)行融合計(jì)算。

同化系統(tǒng)如圖2所示，數(shù)據(jù)包括雷達(dá)影像數(shù)據(jù)、水準(zhǔn)實(shí)測數(shù)據(jù)，同化方法為Burgers等提出的擾動觀測集合卡爾曼濾波算法［12］，數(shù)據(jù)和同化方法組成了整個同化系統(tǒng)。

3 融合實(shí)驗(yàn)與結(jié)果分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文選取9308 工作面為研究對象，選取2012－01－05～04－02的D－InSAR 累積沉降結(jié)果進(jìn)行實(shí)驗(yàn)。由于9308工作面水準(zhǔn)監(jiān)測從2012－03－02開始，無法給出公路監(jiān)測點(diǎn)自01－05的水準(zhǔn)累積沉降結(jié)果，因此將2012－01－05～02－29的D－In－SAR 監(jiān)測數(shù)據(jù)作為初始累積沉降，與后期的水準(zhǔn)測量結(jié)果進(jìn)行累加，得到03－11～04－02的水準(zhǔn)累加值（水準(zhǔn)值）。

首先根據(jù)水準(zhǔn)值，利用遺傳算法反演出概率積分法預(yù)計(jì)參數(shù)，并在此參數(shù)基礎(chǔ)上，采用概率積分法預(yù)計(jì)出工作面的區(qū)域沉降值（反演值），作為觀測場引入到同化系統(tǒng)中，通過集合卡爾曼濾波同化算法融合D－InSAR值與反演值，并將同化結(jié)果與真實(shí)值相比較，從而判斷數(shù)據(jù)同化方法融合數(shù)據(jù)的準(zhǔn)確性和可靠性，將D－InSAR值與反演值經(jīng)集合卡爾曼濾波計(jì)算得出融合結(jié)果。

受9308工作面采動影響的鄒濟(jì)公路監(jiān)測點(diǎn)共24個，隨機(jī)選取其中11個作為驗(yàn)證點(diǎn)，其余為實(shí)驗(yàn)點(diǎn)，如圖3所示，其中“Δ”表示實(shí)驗(yàn)點(diǎn)，“○”表示驗(yàn)證點(diǎn)。

圖3 監(jiān)測點(diǎn)分布Fig.3 Location of bench mark

3.2 融合結(jié)果分析

將數(shù)據(jù)導(dǎo)入同化系統(tǒng)，利用集合卡爾曼濾波進(jìn)行整個面元的數(shù)據(jù)融合計(jì)算。將融合結(jié)果內(nèi)插后得出公路上各觀測站點(diǎn)的同化值，并與水準(zhǔn)值、反演值以及D－InSAR值進(jìn)行對比，如圖4所示。

圖4 驗(yàn)證點(diǎn)數(shù)據(jù)結(jié)果對比Fig.4 The comparison diagram of verification point’s data

從圖4可以看出，4種數(shù)據(jù)的變化趨勢總體上一致，但是數(shù)值大小呈現(xiàn)出一定的差異?？傮w而言，反演值與水準(zhǔn)值的中心區(qū)域差別較小，兩邊差別較大；D－InSAR 值與水準(zhǔn)值的中心區(qū)域差別較大，兩邊差別較?。煌祫t處于兩者之間，這與理論分析相吻合。實(shí)際上，圖4的結(jié)果驗(yàn)證了同化的含義，即作為分析場的同化值是將背景場與觀測場向中間靠近，其數(shù)值應(yīng)該介于作為背景場的D－InSAR值與作為觀測場的反演值之間。地表沉降的同化結(jié)果雖然與水準(zhǔn)值仍存在一定差距，但相比于反演值和D－InSAR結(jié)果有了很大改善。

為了更準(zhǔn)確地衡量系統(tǒng)的同化效果，在RMSE、r以及MAPE（平均絕對百分比誤差）3個指標(biāo)下比較了反演值、D－InSAR 值和同化值。其中，

式中，Yi表示第i個測點(diǎn)的水準(zhǔn)值，Xi為第i個測點(diǎn)的其他測量數(shù)據(jù)，N為測點(diǎn)數(shù)。

由表2 可得，在RMSE、r、MAPE 三個指標(biāo)下，大部分同化值結(jié)果要比D－InSAR 與反演值好，同化值的3個指標(biāo)各自有0、1、1個觀測時刻的數(shù)據(jù)效果位于D－InSAR 與反演值之間。各個時刻同化值的均方根誤差均小于22mm，滿足高等級公路變形監(jiān)測的精度要求。針對表2的驗(yàn)證數(shù)據(jù)，圖5進(jìn)一步從散點(diǎn)圖和誤差分布直方圖對同化結(jié)果與反演值以及D－InSAR值進(jìn)行比較。

綜合表2與圖5可得，與水準(zhǔn)值相比，同化值的效果最好，反演值次之，D－InSAR 值最差。目前高等級公路沉降觀測沒有相應(yīng)規(guī)范?，F(xiàn)場調(diào)研情況顯示，公路下沉200mm 以內(nèi)，采動變形對公路影響輕微，公路上產(chǎn)生不大于2mm 的裂縫，路面基本無隆起現(xiàn)象，行車無不適感，不用維修。根據(jù)國際測量工作者聯(lián)合會第13屆會議提出的精度要求，變形觀測中誤差均方根誤差應(yīng)小于允許變形值的1／10～1／20，高等級公路采動變形破壞測量精度應(yīng)達(dá)20 mm。同化值的總均方根誤差為17.7 mm，滿足高等級公路變形監(jiān)測的精度要求。

表2 公路驗(yàn)證點(diǎn)精度分析Tab.2 Precision analysis table of verification poing’s data

圖5 同化值與反演值及D－InSAR值的比較Fig.5 Comparison of the assimilation results and inversion，D－InSAR

4 結(jié) 語

D－InSAR監(jiān)測結(jié)果與水準(zhǔn)觀測值走勢十分接近，但兩種數(shù)據(jù)在數(shù)值上呈現(xiàn)出一定的差異，距離盆地中心越近，D－InSAR 與水準(zhǔn)實(shí)測值的差別越大?；诩峡柭鼮V波同化D－InSAR 監(jiān)測結(jié)果與水準(zhǔn)觀測值，實(shí)現(xiàn)高精度面狀監(jiān)測的方法是可行的，同化結(jié)果相比于反演值和D－InSAR 監(jiān)測結(jié)果有了很大改善。將該方法應(yīng)用于濟(jì)寧鄒濟(jì)公路，總均方根誤差為17.7 mm，滿足高等級公路變形監(jiān)測的精度要求。

［1］吳立新，高均海，葛大慶，等.工礦區(qū)地表沉陷D－InSAR監(jiān)測試驗(yàn)研究［J］.東北大學(xué)學(xué)報(bào)，2005（8）：778－782（Wu Lixin，Gao Junhai，Ge Daqing，et al.Experimental Study on Surfaee Subsidenee Monitoring with D－InSAR in Mining Area［J］.Journal of Northestern University，2005（8）：778－782）

［2］吳立新，高均海，葛大慶，等.基于D－InSAR 的煤礦區(qū)開采沉陷遙感監(jiān)測技術(shù)分析［J］.地理與地理信息科學(xué)，2004（2）：22－25（Wu Lixin，Gao Junhai，Ge Daqing，et al.Technical Analysis of the Remote Sensing Monitoring for Coal－Mining Subsidence Based on D－InSAR［J］.Geography and Geo－Information Science，2004（2）：22－25）

［3］王桂杰，謝謨文，邱騁，等.D－InSAR技術(shù)在大范圍滑坡監(jiān)測中的應(yīng)用［J］.巖土力學(xué)，2010（4）：1 337－1 344（Wang Guijie，Xie Mowen，Qiu Cheng，et al.Application of D－InSAR Technique to Landslide Monitoring［J］.Rock and Soil Mechanics，2010（4）：1 337－1 344）

［4］劉國祥，丁曉利.沿海地區(qū)D－InSAR 形變探測：精度與可應(yīng)用性分 析［J］.測 繪 通 報(bào)，2006（9）：9－13（Liu Guoxiang，Ding Xiaoli.Detecting of Surface Deformation in Coastlands by D－InSAR：Accuracy and Applicability Analysis［J］.Bulletin of Surveying and Mapping，2006（9）：9－13）

［5］費(fèi)劍鋒，韓月琪.集合卡爾曼濾波數(shù)據(jù)同化在一維波動方程中的應(yīng)用［J］.氣象科技，2005（2）：109－114（Fei Jianfeng，Han Yueqi.Application of Ensemble Kalman Filter Assimilation in One－Dimensional Wave Equation［J］.Meteorological Science and Technology，2005（2）：109－114）

［6］張學(xué)峰，黃大吉，章本照，等.集合數(shù)據(jù)同化方法的發(fā)展與應(yīng)用概述［J］.海洋學(xué)研究，2007（1）：88－94（Zhang Xuefeng，Huang Daji，Zhang Benzhao，et al.The Developments and Applications of Ensemble－Based Data Assimilation Methods［J］.Journal of Marine Science，2007（1）：88－94）

［7］陳春，吳振森，孫樹計(jì)，等.集合卡爾曼濾波在電離層短期預(yù)報(bào)中的應(yīng)用［J］.空間科學(xué)學(xué)報(bào)，2010（2）：148－153（Chen Chun，Wu Zhensen，Sun Shuji.Application of the Ensemble Kalman Filter in Short－Term Ionospheric Forecast［J］.Chinese Journal of Space Science，2010（02）：148－153）

［8］高山紅，吳增茂，謝紅琴.Kalman濾波在氣象數(shù)據(jù)同化中的發(fā)展 與 應(yīng) 用［J］.地 球 科 學(xué) 進(jìn) 展，2000（5）：571－575（Gao Shanhong，Wu Zengmao，Xie Hongqin.The Development and Applications of Kalman Filters in Meteorological Data Assimilation［J］.Advance in Earth Sciences，2000（5）：571－575）

［9］朱江，汪萍.集合卡爾曼平滑和集合卡爾曼濾波在污染源反演中的應(yīng)用［J］.大氣科學(xué)，2006（5）：871－882（Zhu Jiang，Wang Ping.Ensemble Kalman Smoother and Ensemble Kalman Filter Approaches to the Joint Air Quality State and Emission Estimation Problem［J］.Chines Journal of Atmospheric Sciences，2006（5）：871－882）

［10］何敏，陸曉燕，何秀鳳.利用D－InSAR 二軌法監(jiān)測徐州大屯中心區(qū) 地表形 變［J］.地 理 空 間 信 息，2011（5）：3－5（He Min，Lu Xiaoyan，He Xiufeng.Monitoring Ground Deformation in Datun Central District of Xuzhou Using Two Pass D－InSAR［J］.Geospatial Information，2011（5）：3－5）

［11］魏長婧，汪云甲，閆建偉.D－InSAR 技術(shù)二軌法監(jiān)測礦區(qū)地表沉陷的方法研究［J］.煤炭技術(shù)，2012（7）：121－122（Wei Changjing，Wang Yunjia，Yan Jianwei.Study of Application of D－InSAR Technology in Monitoring Mining Area Surface Settle by Two－Track Method［J］.Coal Technology，2012（7）：121－122）

［12］王躍山.數(shù)據(jù)同化——它的緣起、含義和主要方法［J］.海洋預(yù)報(bào)，1999（1）：12－21（Wang Yueshan.Data Assimilation－Its Cause，Its Meaning and Main Procedures［J］.Marine Forecasts，1999（1）：12－21）