D-InSAR和水準數據融合方法研究

楊 帆，趙增鵬，張 磊，張子文(遼寧工程技術大學測繪與地理科學學院，遼寧 阜新 123000)

水準測量是傳統的高程測量技術[1]，它具有精度高、成果可靠、儀器設備簡單、技術容易掌握等優點，但只能得到有限個監測點的變形值。InSAR監測技術是20世紀后期發展起來的新興交叉學科[2]。20世紀90年代，InSAR和D-InSAR技術得到了空前的發展并被廣泛地加以應用[3]。國內對InSAR技術應用于地表形變監測的研究始于近年，且主要集中于城市地面沉降監測[4-5]。但其監測結果受多種誤差源影響，沉降最大值附近精度較低。兩種數據有很好的互補性，因此為全面利用監測數據，更加準確地描述地面沉降區域的沉降現狀，可嘗試應用數據同化的方法將D-InSAR監測值和水準監測數據進行融合。

數據同化是指在考慮數據時空分布、觀測場和背景場誤差的基礎上，在數值模型的動態運行過程中融合新的觀測方法，對監測的動態模型軌跡進行一定的約束，從而提高變量的估計精度。作為數據同化的一種新方法，集合卡爾曼濾波(ensemble Kalman filter，EnKF)近年來取得了顯著的進展[6-7]，但是集合卡爾曼濾波在地面沉降監測多源數據融合中的應用在國內并不多見。本文提出一種基于集合卡爾曼濾波的D-InSAR和水準監測數據融合方法，使融合后的數據達到地面沉降監測數據的高空間分辨率與高高程變形精度的有效統一，使其更好地反映地面沉降區域的沉降規律，為地面沉降監測和預報提供數據上的支持。

1 融合方法

1.1 融合模型的建立

考慮水準監測數據比D-InSAR的處理結果更能真實地反映出地面沉降規律，在利用集合卡爾曼濾波進行D-InSAR值和水準監測數據融合的過程中，采用以水準監測數據為主的同化方式：根據觀測的水準值，利用最小二乘擬合方法擬合出工作面的區域沉降值(反演值)；將反演值與D-InSAR值分別作為觀測場和背景場代入到同化系統中，通過集合卡爾曼濾波同化算法融合D-InSAR值與反演值，并將同化結果與真實值相比較，從而判斷數據同化方法融合數據的準確性和可靠性，將D-InSAR值與反演值經集合卡爾曼濾波計算得出融合結果。

同化系統的流程如圖1所示。

圖1 同化系統流程

1.2 集合卡爾曼濾波原理

集合卡爾曼濾波由Evensen[8]于1994年最先提出并將其應用到海洋同化領域中。1998年，Houtekamer和Mitchell[9]首次在大氣資料同化領域中使用了EnKF方法，后被廣泛應用于大氣預報領域[10]，主要由預測和分析兩部分組成[11]。在集合卡爾曼濾波中，背景場由一組集合成員組成，通過對集合樣本的統計實現卡爾曼濾波方程組中背景誤差協方差矩陣的估計。它的基本流程為[12]：根據背景和觀測的誤差分布特征對背景場和觀測場加擾，再通過統計集合擾動場得到背景誤差協方差矩陣，進而得到分析方程中的增益矩陣，利用增益矩陣和分析方程，對背景集合進行分析，得到分析集合，繼續向前預報[13]。

集合卡爾曼濾波最顯著的特點為：誤差協方差是通過對集合成員的分析統計得到的，具體的分析過程如下：

首先定義t時刻的背景場為

(1)

集合平均為

(2)

集合成員的擾動量定義為

(3)

(4)

背景誤差協方差為

(5)

假設觀測是無偏的，定義觀測向量為

(6)

式中，N為觀測數；m為集合成員數。

基質由斜長石、橄欖石假象、單斜輝石組成，斜長石呈半自形板條狀，長徑一般0.03～0.2mm，雜亂似格架狀分布，表面相對較干凈，局部硅化；橄欖石、單斜輝石呈半自形—他形柱粒狀，粒徑一般<0.1mm，少數0.1～0.2mm，雜亂似填隙狀分布于斜長石格架間，顯間粒結構，橄欖石被皂石及少量硅質、碳酸鹽交代呈假象，單斜輝石局部皂石化、碳酸鹽化。巖內見少量碳酸鹽充填的顯微裂隙分布。

對觀測擾動量定義如下

(7)

(8)

觀測誤差協方差矩陣根據觀測擾動計算得到，定義如下

(9)

K為待定系數，也叫作Kalman增益矩陣，計算如下

K=BtHT(HBtHT+Ot)-1

(10)

在具體的計算中

(11)

(12)

式中，i為集合成員數。集合卡爾曼濾波計算流程如圖2所示。

2 實例驗證

2.1 試驗數據

本試驗采用EnKF方法，結合河北省某礦區開采沉降監測數據，在等間隔時間序列觀測值的基礎上建模，進行了開采區域沉降變形的預測。

選取1326工作面為研究對象，選取2014-03-06至2016-05-18的D-InSAR與水準累計沉降監測數據進行試驗。開采沉陷區概況和水準點布設情況如圖3所示，監測點累計沉降監測數據見表1。

背景場誤差均方差取0.01,觀測場誤差均方差取0.01。集合數m對預測值的精度有一定影響[14]。取值過大會影響計算速度，過小又會給系統帶來統計誤差[15]。為了兼顧計算精度和效率的平衡，本試驗經過模擬不同集合數m，最終取m等于100，可達到計算精度和效率的最大化。

圖2 集合卡爾曼濾波計算流程

圖3 工程概況

表1 監測點累計沉降量 mm

2.2 試驗結果分析

基于表1的數據建立EnKF預測模型，并采用Matlab語言編程實現。數據同化結果對比展示如圖4所示。

圖4 數據同化結果對比

各監測點處D-InSAR值和同化值相對于水準值的誤差如圖5所示。

圖5 誤差對比

為了精確評價同化效果，選取ME(平均誤差)、RMSE(均方根誤差)和MRE(平均相對誤差)3個指標對所得結果進行評價，計算公式為

(13)

(14)

(15)

精度評定結果見表2。

表2 監測點數據同化精度分析 mm

由圖4可知，4種數據總體變化趨勢一致，但呈現出一定的差別。反演值與水準值無論是在沉降盆地中心區域還是邊緣區域差別都比較小；D-InSAR值與水準值中心區域差別較大，邊緣區域差別較小；同化值結果處于兩者之間，與D-InSAR值相比有了很大的改善，同化值能夠較好地反映地面沉降區域的沉降特征。

由圖5可知，D-InSAR值與水準值的誤差基本在5 mm以上，而同化值與水準值的誤差大體都在5 mm以下，D-InSAR值與水準值的誤差遠大于同化值與水準值的誤差，從與水準值的誤差大小來看，同化值的效果好于D-InSAR值。

由表2可知，在ME、RMSE和MRE 3個指標下，同化值的效果均好于D-InSAR值。同化值的平均誤差為3.10 mm，均方根誤差為3.73 mm，平均相對誤差為0.07 mm。

3 結 語

在地面沉降監測中，D-InSAR監測值與水準監測數據總體變化趨勢相似，但在數值上呈現出一定的差別。在靠近沉降盆地中心區域，D-InSAR監測值與水準監測數據差別較大；在靠近沉降盆地邊緣區域，D-InSAR監測值與水準監測數據差別較小。本文利用集合卡爾曼濾波同化D-InSAR值和水準監測數據，使同化后的數據比D-InSAR監測值有了很大改善，更加符合地面沉降地區的沉降變形規律，為預測和防范地面沉降提供更為豐富的數據支持。

參考文獻：

[1] HU R L，YUE Z Q，WANG L C，et al.Review on Current Status and Challenging Issuse of Land Subsidence in China[J].Engineering Geology,2004,76(1):65-77.

[2] 劉國祥,張瑞,李陶,等.基于多衛星平臺永久散射體雷達干涉提取三維地表形變速度場[J].地球物理學報,2012,55(8):2598-2610.

[3] 舒寧.雷達影像干涉測量原理[M].武漢:武漢大學出版社,2003.

[4] 路旭,匡紹君,賈有良.用InSAR做地面沉降監測的試驗研究[J].大地測量與地球動力學,2002,22(4):66-70.

[5] ZHAO Chaoying,DING Xiaoli,ZHANG Qin,et al.Monitoring of Land Subsidence and Ground Fissures in Xi’an China 2005—2006:Mapped by SAR Interferometry[J].Environ Geol,2009,58(3):1533-1540.

[6] 劉國祥，丁曉利.沿海地區D-InSAR形變探測：精度與可應用性分析[J].測繪通報，2006(9)：9-13.

[7] 陳春，吳振森，孫樹計，等.集合卡爾曼濾波在電離層短期預報中的應用[J].空間科學學報，2010(2)：148-153.

[8] EVENSEN G.Sequential Data Assimilation with a Non-linear Quasi-geostrophic Model Using Monte Carlo Methods to Forecast Error Statistics[J].Journal of Geophysical Research Oceans,1994,99(C5)：10143-10162.

[9] HOUTEKAMER P L,MITCHELL H L.A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation[J].Monthly Weather Review,2001,129(1):123-137.

[10] GILLIJINS S,MENDOZA O B,CHANDRASEKAR J,et al.What Is the Ensemble Kalman Filter and How Well Does it Work?[J].American Control Conference,2006(6):4448-4453.

[11] 崔波，張家樹，楊宇.基于集合卡爾曼濾波的非線性目標跟蹤算法[J].計算機仿真，2013,30(4)：317-321.

[12] 董亞寧.基于集合卡爾曼濾波全球臭氧衛星觀測資料同化試驗研究[D].南京：南京信息工程大學，2016.

[13] 惠雪峰.基于集合卡爾曼濾波的森林資源動態預測[D].南京：南京林業大學，2012.

[14] 米麗倩，查劍鋒，王新.老采空區殘余沉降的集合卡爾曼濾波預測[J].金屬礦山，2012(8)：138-141.

[15] 李喜佳，肖志強，王錦地，等.雙集合卡爾曼濾波估算時間序列LAI[J].遙感學報，2014，18(1)：27-44.