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

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

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

1 融合方法

1.1 融合模型的建立

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

同化系統(tǒng)的流程如圖1所示。

圖1 同化系統(tǒng)流程

1.2 集合卡爾曼濾波原理

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

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

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

(1)

集合平均為

(2)

集合成員的擾動量定義為

(3)

(4)

背景誤差協(xié)方差為

(5)

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

(6)

式中，N為觀測數(shù)；m為集合成員數(shù)。

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

對觀測擾動量定義如下

(7)

(8)

觀測誤差協(xié)方差矩陣根據(jù)觀測擾動計算得到，定義如下

(9)

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

K=BtHT(HBtHT+Ot)-1

(10)

在具體的計算中

(11)

(12)

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

2 實例驗證

2.1 試驗數(shù)據(jù)

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

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

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

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

圖3 工程概況

表1 監(jiān)測點累計沉降量 mm

2.2 試驗結果分析

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

圖4 數(shù)據(jù)同化結果對比

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

圖5 誤差對比

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

(13)

(14)

(15)

精度評定結果見表2。

表2 監(jiān)測點數(shù)據(jù)同化精度分析 mm

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

由圖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 結 語

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

參考文獻：

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

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

[4] 路旭,匡紹君,賈有良.用InSAR做地面沉降監(jiān)測的試驗研究[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] 劉國祥，丁曉利.沿海地區(qū)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] 董亞寧.基于集合卡爾曼濾波全球臭氧衛(wèi)星觀測資料同化試驗研究[D].南京：南京信息工程大學，2016.

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

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

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