融合三次樣條插值的D-InSAR沉陷變形監測技術

余禮仁，徐良驥，龐 會，王振兵

(安徽理工大學，安徽 淮南 232000)

融合三次樣條插值的D-InSAR沉陷變形監測技術

余禮仁，徐良驥，龐 會，王振兵

(安徽理工大學，安徽 淮南 232000)

傳統煤礦開采沉陷形變監測只能獲得離散監測點的地表形變量，無法獲取沉陷區整體的形變趨勢，而D-InSAR技術能夠獲取整體地表形變信息，但D-InSAR技術需要大量遙感影像數據，否則導致時間失相干性而無法獲得連續的形變量。本文針對D-InSAR影像數據解譯過程中影像數據較少的問題，首先利用D-InSAR技術對淮北礦業集團袁店二礦7225、7226工作面進行了監測分析，獲取了幾個時間段內整體的下沉形變場；然后提出了采用三次樣條插值與D-InSAR技術相結合的開采沉陷監測方法，基于D-InSAR影像上監測點監測值，利用三次樣條插值建立了內插反演函數，在已建立函數的基礎上得出了其他雷達衛星重訪周期內的下沉值；最后將內插反演下沉值與實測水準數據進行了對比分析。結果表明：D-InSAR監測結果能夠有效地反映開采沉陷的影響范圍，能夠較為準確地提取區域地表的形變信息；結合三次樣條插值的D-InSAR技術監測結果最大誤差和最大相對誤差分別為31.5 mm和17%，該方法能有效地解決D-InSAR影像數據缺少的問題。

礦區沉陷；變形監測；D-InSAR；三次樣條插值

Abstract： Traditional coal mining subsidence deformation monitoring can only obtain the surface shape of the discrete monitoring variables,unable to get deformation trend of subsidence area as a whole,while the D-InSAR technique can obtain the overall surface deformation information.But D-InSAR technique requires a large number of remote sensing image data,otherwise it will lead to loss of coherence time and won’t be able to obtain continuous variables.Based on D-InSAR image data interpretation in the process of image data with fewer problems,this paper first using the D-InSAR technology yuen shop two mine of huaibei mining industry group in 7225,7226 working face monitoring analysis, obtained several period subsidence deformation field of the whole,and then put forward combined cubic spline interpolation with D-InSAR technique mining subsidence monitoring method,based on D-InSAR imaging of monitoring value on monitoring,using the cubic spline interpolation to establish interpolation inverse function, the function has been established on the basis of other radar satellite revisit cycle sinking value.Finally the interpolation inversion subsidence value and the measured level data were analyzed.The results show that D-InSAR technique can effectively reflect the influence range of mining subsidence,and can extract the deformation information of regional surface accurately.Combining with the cubic spline interpolation of D-InSAR technology monitoring results maximum error and the maximum relative error is 31.5 mm and 17% respectively.The method can effectively solve the problem of the lack of D-InSAR imaging data.

Keywords： mining subsidence;deformation monitoring;D-InSAR;cubic spline interpolation

傳統的礦區開采沉陷變形監測常采用精密水準測量、全站儀三角測量和GPS測量等方法[1-4]。但常規的監測方法存在費用高、需要大量人力物力、監測周期長等缺點，而D-InSAR技術作為在SAR基礎上發展而來的新型空間對地觀測技術，具有連續空間覆蓋能力、高度自動化和高精度的優勢，近年來在礦區開采沉陷中得到廣泛的應用[5-6]。D-InSAR技術主要是獲取監測區域內地表形變信息，地表形變是一個長期持續的過程，大多數研究方法是采用多幅影像進行差分干涉處理，得到地表時序形變信息。在此基礎上提出了干涉疊加技術，主要包括永久散射體干涉測量(PS-InSAR)和小基線集技術(SBAS-InSAR)[7-10]。兩種監測方法都要求獲得連續時間序列的SAR影像數據，因此雖然D-InSAR技術在礦區開采沉陷變形監測方面具有很多成功且典型的工程實例，但在具體的實踐中和后續數據處理中受到很多限制性約束，對于礦山開采沉陷監測而言，SAR影像數據不足的情況下主要影響因素是時間失相干，從而不能獲得連續時間段的地表下沉值，因此SAR影像源數據的獲取和影像費用高的問題成為影響該技術的發展前景之一[11]。為了克服這一缺陷，本文提出將合成孔徑雷達差分干涉測量技術(D-InSAR)與三次樣條插值法相融合的開采沉陷監測方法。該方法利用D-InSAR影像上對應監測點在已知雷達衛星重訪周期內的監測值作為三次樣條插值的節點，通過三次樣條插值法來內插反演其他雷達重訪周期內監測點的下沉值。

1 D-InSAR沉陷監測原理與內插反演模型的建立

1.1 D-InSAR測量技術原理

合成孔徑雷達差分干涉測量(D-InSAR)技術以合成孔徑雷達(SAR)復數據所提供的相位及強度信息為信息源，利用同一目標區域的兩幅或多幅干涉紋圖獲取SAR的相位信息，從而提取地表三維信息和高程變化信息，然后通過差分處理(除去地球表面、地形起伏等因素)來獲取地表微小形變[12-13]。

目前，D-InSAR獲取地表形變信息的實現方法主要有3種：二軌法、三軌法和四軌法[14-17]。二軌法是獲取同一地區形變發生前后的兩幅SAR影像，將其進行干涉處理，干涉結果與已有DEM數據進行差分，消除地形因素而獲取地表形變信息的方法；三軌法需要3幅SAR影像，其中兩幅進行差分干涉處理生成干涉圖，作為DEM信息，再將另一幅影像與生成的干涉圖進行干涉處理，最終得出地表形變信息；四軌法基本思想是獲取形變前后的兩對SAR影像，分別進行干涉處理，形成形變前后的兩幅干涉相位圖，形變前的主要用來生成DEM，形變后的主要用來反映形變信息，將二者進行差分，得到形變相位，再將相位轉換為斜距從而計算出雷達視線向的形變量。

1.2 融合D-InSAR監測值的三次樣條插值沉陷內插反演函數的建立

三次樣條插值函數是最常用的插值曲線擬合函數，其良好的收斂性、可靠的穩定性，以及具有二階光滑度等優點使其在函數逼近、微積分和微分方程等科學計算中應用廣泛。在礦區開采沉陷中，由于沉降是一個隨時間平滑變化的過程，因此，三次樣條插值在開采沉陷數據處理中起到非常重要的作用。

本文利用D-InSAR技術獲得(礦區沉降面上)一系列與時間相關的離散數據{xj}={x0,x1,…,xn}，xj表示雷達衛星重訪周期。使用三次樣條插值進行沉陷監測，建立已有觀測數據和未知預測數據之間的函數關系。以區間[xj,xj+1]為例，S(x)的二階導數為S″(x)=Mj(j=0,1,…,n)，則該區間上S(x)的三次樣條插值函數為

綜上所述，支氣管哮喘患兒血清IgE、IL-4水平升高，RBP水平下降，MP感染能夠加重血清IgE、IL-4升高程度，降低RBP水平。

(1)

(2)

線性方程組的解為M，代入式(2)中可以得到各個小區間上的三次樣條插值函數式[18-19]。采用三次樣條插值建立開采沉陷監測函數步驟如下：

(1) 為構建三次樣條插值所需要的已知觀測數據，首先對已有的雷達影像進行差分干涉處理，得到一系列影像周期內監測點的下沉值。以與首幅影像日期(2012-11-01)間隔24 d為第一個周期開始起算，將得到的下沉值作為三次樣條插值函數的已知數據點。

(2) 運用MATLAB軟件進行編程，使用三次樣條插值函數spline，建立內插模型

x=[x1,…,xi,…,xn];

y=[y1,…,yi,…,yn];

cs=spline(x,y);

xx=linspace(0,38,38);

plot(x,y,′o′,xx,ppval(cs,xx),′-′)

2 D-InSAR數據處理和結果分析

2.1 研究區概況及相關數據

本次試驗采用淮北礦業集團袁店二礦作為研究區域。在2012年11月—2013年12月對7226工作面進行了開采工作，7226工作面長約768 m，工作面平均寬為158.9 m，工作面標高為-453～-416 m，煤層平均厚度3.89 m，煤層傾角為3°～14°。由于該工作面上方有村莊，進一步的地下采煤會引起地表建筑物的損壞，因此在該工作面采用注漿充填開采方式，從2013年7月開始進行注漿開采。2013年10月—2015年5月對7225工作面進行了開采工作，7225工作面長約960 m，工作面平均寬為176 m，工作面平均標高為-450 m，煤層平均厚度為4.07 m，煤層傾角為6°～19°。

本文采用9景C波段RADARSAT-2影像數據對試驗礦區進行研究分析，影像分辨率為3 m，雷達衛星重訪周期為24 d。本文選用的外部DEM是SRTM3 DEM數據，是美國奮進號航天飛機在2000年2月11—22日進行的為期11 d的航天飛機雷達地形測繪任務SRTM所獲得的。為減少時間失相關，盡可能選取時間間隔最小的兩幅影像進行干涉處理，9景影像共組成5組干涉對，其相關參數見表1。

表1 研究區SAR影像干涉對參數表

2.2 D-InSAR數據處理及結果分析

為了得到礦區地表下沉形變場，采用SARscape軟件中D-InSAR Displacement Workflow工作流分別對5組干涉對進行二軌差分處理，經過多視和濾波抑制斑點噪聲，相位解纏后進行地理編碼和輻射定標得到形變圖，將得到的5組形變圖和在7225、7226工作面設立的地表移動觀測站Cass圖導入ArcMap中進行疊加分析，得到5個時間段內工作面地表形變圖，如圖1所示。

圖1 7225、7226工作面不同時間段內沉降分布

(1) 由圖1(a)可知，形成的最大下沉區域位于距離7226工作面340 m的7221工作面，該工作面開采起止時間是2011年1月—12月，最大下沉區域是由于該工作面開采殘余變形的影響，并非7226工作面開采引起，而7226工作面在開采初期沉降值和影響范圍很小。

(2) 由圖1(b)可知，注漿減沉開采的影響導致注漿周邊區域下沉分布不均勻，在2013年11月—12月期間，地表最大下沉區域和影響范圍隨著7226工作面開采逐漸向前推進，形成了明顯的形變區沉降漏斗。

(3) 由圖1(c)可知，由于7225工作面開采與7226工作面殘余變形疊加影響，造成下沉影響范圍較大的區域分別出現在7226工作面的尾端和7225工作面的開切眼處，而在7226工作面注漿區域地表下沉值依然分布不均勻。

(4) 由圖1(d)可知，隨著7225工作面采煤的推進，下沉影響范圍增加，而下沉值最大區域不在7225工作面正上方，而是在兩個工作面之間，這是由于受到7226工作面開采的殘余變形的影響。

(5) 由圖1(e)可知，隨著兩個工作面采煤工作的結束，地表受殘余變形的影響，在216 d內下沉范圍增廣明顯，但下沉速率開始減緩，下沉值大于10 mm等值線的區域平均沉降速率為0.23 mm/d，屬于下沉衰退期[20]。

3 結合三次樣條插值函數的開采沉陷監測

為了克服SAR影像數據不足而存在的時間失相干現象，只能得到7225、7226工作面5個時間段內的下沉值和影響范圍，不能獲取2012年11月—2015年12月時間內的連續下沉影響范圍。本文采用三次樣條插值方法，將D-InSAR監測值作為三次樣條插值函數的初始值。根據開采沉陷相關知識可知，開采工作面的走向線和傾向線方向的下沉值最能代表地表的下沉規律，因此選取走向線上L系列監測點中的ML02、ML06、ML09、ML11和傾向線S系列監測點中的MS01、MS04、MS08、MS12為研究對象，如圖1所示。以雷達衛星一個重訪周期的時間間隔為三次樣條插值函數的自變量步長，將1號和4號干涉對中兩個周期內的下沉值作平均值處理，所選取的8個監測點下沉信息見表2。

表2 不同雷達衛星重訪周期內所選監測點下沉值 mm

利用三次樣條插值法求出2012年11月—2015年5月上述8個水準點每個雷達衛星重訪周期內的下沉值，運用MATLAB軟件繪制地表下沉曲線圖，如圖2所示。由圖2可知，雷達衛星重訪周期內該區域地表下沉趨勢略有不同，符合實際地表下沉規律。最大下沉值出現在靠近中心盆地監測點ML11處，單個周期內最大下沉量為42.4 mm。

為了進一步驗證精度，實地觀測了2013年7月23日(對應SAR重訪周期第11周)—2014年7月18日(對應SAR重訪周期第26周)7225、7226工作面16個周期(共360 d)內的水準數據，將內插反演值與實測值對比，結果見表3。由表3可知，插值解算累計下沉值與實測結果基本一致，其中最大誤差出現在監測點ML11處，為31.5 mm，其余監測點誤差均在20 mm以內。最大相對誤差出現在監測點MS01處，為17%，其余監測點處相對誤差都不大于8.2%。

表3 11～26雷達衛星重訪周期內累計反演值與實測值比較

圖2 第11～26個雷達衛星重訪周期內所選監測點下沉值

4 結 論

(1) 本文采用D-InSAR技術對淮北礦業集團袁店二礦7225、7226工作面開采沉陷進行了分析，獲取了5個時期內地表沉陷的下沉影響范圍，結果與實際礦區開采導致地表下沉影響范圍符合得很好。

(2) 融合D-InSAR技術與三次樣條插值的方法有效地解決了在SAR影像數據不足的情況下對袁店二礦沉陷區進行監測分析的問題。應用結果表明：該方法內插解算反演結果較實測水準監測結果最大誤差和最大相對誤差分別為31.5 mm和17%，內插反演結果與實測值比較吻合，從而有效地說明了該方法的可行性。

[1] 岳建平，田林亞．變形監測技術與應用[M].北京:國防工業出版社，2007．

[2] WEBLEY P W,BINGLEY R M,DODSON A H，et al.Atmospheric Water Vapour Correction to InSAR Surface Motion Measurements on Mountains:Results from a Dense GPS Network on Mount Etna[J].Physics and Chemistry of the Earth，2002，27(4):363-370.

[3] 陳炳乾,鄧喀中,范洪冬.基于D-InSAR技術和SVR算法的開采沉陷監測與預計[J].中國礦業大學學報,2014,43(5):880-886.

[4] CHEN Bingqian,DENG Kazhong,FAN Hongdong,et al.Large-scale Deformation Monitoring in Mining Area by D-InSAR and 3D Laser Scanning Technology Integration[J].International Journal of Mining Science and Technology,2013,23(4):555-561.

[5] BAYUAJI L,SUMANTYO J T S,KUZE H.ALOS PALSAR D-InSAR for Land Subsidence Mapping in Jakarta,Indonesia[J].Canadian Journal of Remote Sensing,2010,36(1):1-8.

[6] PEI Liang,LI Wenjie,TAN Yang.Study of Monitoring Mining Subsidence in Coal Mining Area by D-InSAR Technology[J].Journal of Coal Science and Engineering,2008,14(4):591-593.

[7] 楊成生,劉媛媛,敖萌.基于SBAS時序分析的大同地面沉降與地下水活動研究[J].國土資源遙感,2015,27(1):127-132.

[8] 侯安業,張景發,劉斌,等.PS-InSAR與SBAS-InSAR監測地表沉降的比較研究[J].大地測量與地球動力學,2012,34(4):125-128.

[9] 李永生,張景發,羅毅,等.利用高分辨率聚束模式TerraSAR-X影像的PSInSAR監測地表變形[J].武漢大學學報(信息科學版),2012,27(12):1452-1455,1514.

[10] 許才軍,何平,溫揚茂.利用PSInSAR研究意大利Etna火山的地表形變[J].武漢大學學報(信息科學版),2011,26(9):1012-1016.

[11] 吳立新,高均海,葛大慶,等.基于D-InSAR的煤礦區開采沉陷遙感監測技術分析[J].地理與地理信息科學,2004,20(2):22-25.

[12] 范洪冬,鄧喀中,祝傳廣,等.基于時序SAR技術的采空區上方高速公路變形監測及預測方法[J].煤炭學報,2012,37(11):1841-1846.

[13] ALLAN M L,MUTHAMA M N,KINYANJUI S Z.Time Series Analysis Model for the Rate of Influx of Refugees in Kenya[J].The International Journal of Engineering and Science,2013,2(9):7-18.

[14] 魏長婧,汪云甲,閆建偉.D-InSAR技術二軌法監測礦區地表沉陷的方法研究[J].煤炭技術,2012,31(7):129-130.

[15] 王小兵.基于DInSAR技術的礦山開采沉陷監測研究現狀[J].金屬礦山,2015(S1):65-71.

[16] SHARIFIKIA M.A Comprehensive Interferometric Process for Monitoring Land Deformation Using ASAR and PALSAR Satellite Interferometric Data[J].GIScience & Remote Sensing,2015,52(1):58-77.

[17] 劉國林,郝華東,陶秋香.卡爾曼濾波相位解纏及其與其他方法的對比分析[J].武漢大學學報(信息科學版),2010,25(10):1174-1178.

[18] 陳文略,王子羊.三次樣條插值在工程擬合中的應用[J].華中師范大學學報(自然版),2004,38(4):418-422.

[19] 張威，楊月婷.數值分析[M].5版.北京：清華大學出版社，2010:41-44.

[20] 高延法,賈君瑩,李冰,等.地表下沉衰減函數與塌陷區穩定性分析[J].煤炭學報,2009,34(7):892-896.

D-InSARDeformationMonitoringTechnologyBasedonCubicSplineInterpolation

YU Liren，XU Liangji，PANG Hui，WANG Zhenbing

(Anhui University of Science and Technology,Huainan 232000,China)

P23

A

0494-0911(2017)09-0051-05

2017-01-03

安徽省對外科技合作計劃(1503062020)；淮北礦業集團有限責任公司2014年研發計劃

余禮仁(1993—)，男，碩士生，研究方向為礦山開采沉陷及D-InSAR技術。E-mail：yu_liren@qq.com

徐良驥

余禮仁，徐良驥，龐會，等.融合三次樣條插值的D-InSAR沉陷變形監測技術[J].測繪通報，2017(9)：51-55.

10.13474/j.cnki.11-2246.2017.0286.