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

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

融合三次樣條插值的D-InSAR沉陷變形監(jiān)測技術(shù)

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

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

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

礦區(qū)沉陷；變形監(jiān)測；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

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

1 D-InSAR沉陷監(jiān)測原理與內(nèi)插反演模型的建立

1.1 D-InSAR測量技術(shù)原理

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

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

1.2 融合D-InSAR監(jiān)測值的三次樣條插值沉陷內(nèi)插反演函數(shù)的建立

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

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

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

(1)

(2)

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

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

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

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數(shù)據(jù)處理和結(jié)果分析

2.1 研究區(qū)概況及相關(guān)數(shù)據(jù)

本次試驗采用淮北礦業(yè)集團袁店二礦作為研究區(qū)域。在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影像數(shù)據(jù)對試驗礦區(qū)進行研究分析，影像分辨率為3 m，雷達衛(wèi)星重訪周期為24 d。本文選用的外部DEM是SRTM3 DEM數(shù)據(jù)，是美國奮進號航天飛機在2000年2月11—22日進行的為期11 d的航天飛機雷達地形測繪任務SRTM所獲得的。為減少時間失相關(guān)，盡可能選取時間間隔最小的兩幅影像進行干涉處理，9景影像共組成5組干涉對，其相關(guān)參數(shù)見表1。

表1 研究區(qū)SAR影像干涉對參數(shù)表

2.2 D-InSAR數(shù)據(jù)處理及結(jié)果分析

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

圖1 7225、7226工作面不同時間段內(nèi)沉降分布

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

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

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

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

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

3 結(jié)合三次樣條插值函數(shù)的開采沉陷監(jiān)測

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

表2 不同雷達衛(wèi)星重訪周期內(nèi)所選監(jiān)測點下沉值 mm

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

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

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

圖2 第11～26個雷達衛(wèi)星重訪周期內(nèi)所選監(jiān)測點下沉值

4 結(jié) 論

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

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

[1] 岳建平，田林亞．變形監(jiān)測技術(shù)與應用[M].北京:國防工業(yè)出版社，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技術(shù)和SVR算法的開采沉陷監(jiān)測與預計[J].中國礦業(yè)大學學報,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] 侯安業(yè),張景發(fā),劉斌,等.PS-InSAR與SBAS-InSAR監(jiān)測地表沉降的比較研究[J].大地測量與地球動力學,2012,34(4):125-128.

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

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

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

[12] 范洪冬,鄧喀中,祝傳廣,等.基于時序SAR技術(shù)的采空區(qū)上方高速公路變形監(jiān)測及預測方法[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技術(shù)二軌法監(jiān)測礦區(qū)地表沉陷的方法研究[J].煤炭技術(shù),2012,31(7):129-130.

[15] 王小兵.基于DInSAR技術(shù)的礦山開采沉陷監(jiān)測研究現(xiàn)狀[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] 張威，楊月婷.數(shù)值分析[M].5版.北京：清華大學出版社，2010:41-44.

[20] 高延法,賈君瑩,李冰,等.地表下沉衰減函數(shù)與塌陷區(qū)穩(wěn)定性分析[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)；淮北礦業(yè)集團有限責任公司2014年研發(fā)計劃

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

徐良驥

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

10.13474/j.cnki.11-2246.2017.0286.