融合多種差分干涉測量的礦區(qū)地表形變監(jiān)測

樂穎， 夏元平*， 錢文龍

(1.東華理工大學測繪工程學院， 南昌 330013； 2.福州市勘測院， 福州 350108)

中國煤礦資源非常豐富，在開采過程中往往會出現高強度開采現象。而長期過度開采煤礦資源，將導致開采區(qū)域地表發(fā)生沉降，甚至可能對周邊的居民區(qū)、道路和農田等造成嚴重影響[1-2]。目前，隨著社會科技的進步，監(jiān)測技術手段越來越多。合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)技術由于其全天時、全天候的特點，已經廣泛應用于地表形變監(jiān)測[3-4]、滑坡監(jiān)測[5-7]、地形測繪、冰川測量[8]和地質災害監(jiān)測[9]等領域。PS-InSAR技術對數據數量要求較多，且監(jiān)測結果與研究區(qū)域地貌密切相關，在植被較少的區(qū)域可取得較好的效果，而當地區(qū)覆蓋較多植被時，取得的監(jiān)測效果并不理想。煤礦開采區(qū)域往往位于山區(qū)，植被覆蓋較多，且具有較大的開采量與較為明顯的形變，易引起失相干問題。因此，PS-InSAR技術難以準確估算煤礦沉降中心形變值，但當PS-InSAR技術中的第一步反演所使用的控制點大多選擇在形變小且相干性高的區(qū)域，符合SBAS-InSAR中參考點的選擇要求，可將其直接應用于SABS-InSAR中作為GCP點從而提高軌道精煉的精度。而D-InSAR、SBAS-InSAR這兩種技術方法在礦山形變監(jiān)測領域應用較多。

為了避免單一方法的局限性，結合多種技術方法進行監(jiān)測已經成為一種新的趨勢。徐小波等[10]聯(lián)合D-InSAR與偏移量追蹤技術(offset-tracking)技術對高強度開采沉陷區(qū)進行了監(jiān)測，并取得了良好效果。陳艷青等[11]利用D-InSAR和SBAS-InSAR技術監(jiān)測大橋歷史形變趨勢，發(fā)掘雷達衛(wèi)星的監(jiān)測形變能力，為以后的監(jiān)測預警提供了寶貴的經驗。王勇等[12]融合D-InSAR與精密水準測量數據，利用線性回歸構建校正模型，并應用于滑坡形變監(jiān)測，從而提高了形變結果的精度。付波霖等[13]利用時序InSAR技術提取了西藏一滑坡點的形變信息，并對形變結果進行了交叉驗證。柴華彬等[14]為解決因植被旺盛引起的InSAR失相干問題，提出了一種融合實測數據的SBAS-InSAR方法，通過對失相干區(qū)域采用反距離加權插值法，從而獲取地表沉降信息。但是，大部分的學者更多的是根據不同技術得到的結果進行對比，或在缺少水準數據時通過交叉驗證說明結果的準確性。并且，傳統(tǒng)D-InSAR方法會受到時間失相關、空間失相關和大氣效應等因素影響，從而限制了其在礦山形變中的應用；SBAS-InSAR可以彌補這些因素的制約，將兩種技術融合或可取得更優(yōu)的監(jiān)測精度。鑒于此，現提出一種融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變監(jiān)測，將D-InSAR結果進行Kriging插值以補充SBAS-InSAR在礦區(qū)沉陷中心空缺部分，使得最終的礦區(qū)形變結果更加完整，并利用水準數據進行驗證。

1 時序InSAR技術

1.1 基于相干性系數的 D-InSAR 地表形變信息提取方法

D-InSAR作為一門極具潛力的技術，基本思想是通過干涉相干性分析提取地表變化信息，并利用二次差分方法從干涉相位圖中去除地形與其他因素的影響，從而達到提取形變信息的目的?！岸壏ā辈罘指缮婕夹g表達式[15-16]為

φ=φflat+φtopo+φdefo+φatm+φnoise

(1)

式(1)中：φflat為平地相位；φtopo為地形起伏造成的地形相位；φdefo為地表發(fā)生形變；φatm為由于大氣發(fā)生延遲產生的大氣相位；φnoise為觀測時伴隨的噪聲相位。

總之，該技術最終目的就是將式(1)中φflat、φtopo、φatm與φnoise等相位信息剔除，最后只保留φdefo形變相位的技術。

相干系數于1993年最早由普拉蒂提出，該系數取值范圍為[0,1]，取值越大說明相干性越高，理論模型表達式[17]為

(2)

式(2)中：E[·]為數學期望；u*為共軛復數；u1、u2為主、副圖像的信號。然而，在實際的InSAR處理過程中，同一像素內很難獲取計算所需數量的采樣值。因此，根據式(2)的理論模型，基于SAR影像復數數據計算獲得相干系數的標準表達式[18]為

(3)

式(3)中：u1(n,m)、u2(n,m)分別為主、副影像數據塊內某個坐標(n,m)處的復數值；|·|2為對應的二階范數；M、N為計算相干性的數據塊尺寸大??；m、n為數據塊內對應的行列號。

1.2 結合PS-InSAR的SBAS-InSAR技術

SBAS-InSAR干涉技術是由Ferretti等[19]提出的一種用于獲取地表形變信息的時間序列分析方法，其基本原理是對相干目標進行相位分析，最終獲取時序形變。通過選擇合適的時空基線閾值構成差分干涉對，選定相干目標點，并利用線性相位變化構建模型與解算，通過時空濾波方法去除大氣延遲，同時減弱了高程、大氣誤差以及D-InSAR處理過程中存在的失相關問題，從而獲取地表的時間序列形變信息。

若SAR影像數據集中包含N+1景覆蓋同一區(qū)域的影像，則影像獲取時間序列[20]為

t=[t0,t1,…,tN]T

(4)

根據三基線和最小法選定其中一幅影像為主影像進行配準后，在全部自由組合差分干涉對中選取合適的時間基線和空間基線閾值干涉對，若獲得M幅差分干涉圖，則有

(5)

假設t0為參考時刻，則任意時刻ti(i=1,2,…,N)相對于t0時刻的差分相位φ(ti)為未知數，數據處理過程中所獲取的差分干涉相位δφ(tk)(i=1,2,…,M)為觀測量。若忽略失相關、高程誤差以及大氣延遲等因素影響，則第i幅(i=1,2,…,M)差分干涉圖中像元(r,x)相位計算公式為

δφi(r,x)=φ(tA,r,x)-φ(tB,r,x)

(6)

式(6)中：λ為雷達波長；d(tA,r,x)和d(tB,r,x)分別為像元在主影像成像時刻tA和副影像成像時刻tB沿雷達視線方向的形變。若假設d(tA,r,x)=0，那么時間序列上的相位值為

(7)

式(7)中：d(ti,r,x)為任意時刻ti(i=1,2，…,N)相對于t0時刻沿雷達視線方向的形變。

圖1 研究區(qū)域地理位置分布圖Fig.1 Geographical distribution of the study area

2 融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變監(jiān)測

2.1 研究區(qū)域概況

峰峰礦區(qū)位于114°3′E～114°16′E、36°20′N～36°34′N[21]。峰峰礦區(qū)包括九龍礦、梧桐礦、新三礦、薛村礦、孫莊礦、牛兒莊礦、萬年礦、通二礦、四礦等礦區(qū)，圖1繪制了研究區(qū)域具體的地理位置分布情況。

現主要研究九龍礦區(qū)和萬年礦區(qū)，從圖1可以看出，這兩大礦區(qū)的具體地理位置分別位于峰峰礦區(qū)的新坡鎮(zhèn)和武安市的伯延鎮(zhèn)。九龍礦隸屬于冀中能源峰峰集團，是該集團有限公司名下的主力礦井之一。九龍礦位于太行山與華北平原的過渡地帶，屬于緩傾斜山前沖洪積地貌類型。礦區(qū)內地勢西低東高，南低北高，總體形態(tài)為四周高中間低的盆地。萬年礦區(qū)所在區(qū)域，四季分明，冬季寒冷干燥，夏季炎熱多雨，地勢低平，板塊主體由南向北傾斜。值得注意的是，萬年礦區(qū)和九龍礦區(qū)均位于農田下方，為了避免農田因測量而遭到破壞，同時考慮到觀測站的穩(wěn)定性，將水準測量觀測站布設在附近的道路上[22]。

2.2 研究方法

融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變監(jiān)測，其目的主要是對D-InSAR和SBAS-InSAR這兩種方法所得結果進行Kriging值，從而使得監(jiān)測結果更加完整，以彌補SBAS-InSAR技術在監(jiān)測礦區(qū)沉降中心存在空缺值的不足。在雷達影像數據處理過程中需要借助外部DEM，用來提供參考地形或者參考地理坐標系，常見的SRTM系統(tǒng)主要是30 m分辨率的SRTM1和90 m分辨率的SRTM3(https://srtm.csi.cgiar.org/)，本文選取的是SRTM 30 m數據。圖2為了融合D-InSAR和SBAS-InSAR技術的礦區(qū)形變監(jiān)測流程圖。本文研究方法主要改進部分及具體操作步驟如下。

圖2 融合D-InSAR和SBAS-InSAR技術的礦區(qū)形變監(jiān)測Fig.2 Deformation monitoring in mining area by integrating DInSAR and SBAS-InSAR technology

(1)D-InSAR方法中的軌道精煉和重去平：由于軌道精煉時缺少形變先驗認識及地面控制點信息，盲目選擇控制點可能會引入誤差。本文研究在選取D-InSAR方法中軌道精煉控制點時，引入相干系數法，該方法通過設置相干性系數閾值以剔除相干性低的像元點，從而保證選取的控制點相干性很高；再結合相位解纏結果，避免選取殘余地形條紋上并遠離形變區(qū)。

(2)SBAS-InSAR方法中主影像的選?。褐饔跋襁x取對于差分干涉、相位解纏、估計平均速率等都有影響。因此，為了選取效果更好的主影像，本文研究采用三基線和最小法來選取主影像[23]。通過計算每景影像的時間基線、空間基線和多普勒基線，將影像集中這三種基線和最小的影像作為主影像，表達式為

(8)

(3)SBAS-InSAR方法中的軌道精煉：軌道精煉過程中需要選取若干參考點，該步驟主要由人為主觀判斷選擇，并需要遠離形變區(qū)和位于沒有殘余地形條紋上。若盲目選擇控制點可能會引入誤差，故此步驟對于操作人員的專業(yè)知識和經驗能力要求極高。由于PS-InSAR技術中的第一步反演所使用的控制點大多選擇在形變小且相干性高的區(qū)域，符合SBAS-InSAR中參考點的選擇要求，因此可將其直接應用于SABS-InSAR中作為軌道精煉的控制點。利用地理編碼將PS-InSAR技術中的參考點文件從地理坐標系轉換為SAR坐標系，并將其作為SBAS-InSAR技術的軌道精煉控制點引入數據處理中，用于去除軌道數據所引起的誤差。

表1 雷達影像數據的各個基線參數Table 1 Baseline parameters of radar image data

(4)D-InSAR結果和SBAS-InSAR結果的融合：本文采用的是Kriging插值法，將DInSAR方法得到的形變結果用于填補SBAS-InSAR結果在礦區(qū)中的空缺值，充分利用兩種方法的互補性，以此獲取最終的礦區(qū)地表形變結果。此外，Kriging法可對周圍的測量值進行加權，從而得出未測量位置的預測，其表達式由數據的加權總和組成：

(9)

式(9)中：Z(si)為第i個位置的測量值；λi為第i個位置測量值的未知權重；S0為預測位置；N為測量值數。需要注意的是，在使用Kriging方法時，權重不僅取決于測量點之間的距離和預測位置，還取決于基于測量點的整體空間排列。若需要在權重中使用空間排列，必須量化空間自相關。因此，在普通Kriging法中，權重λi取決于測量點和預測位置的距離以及預測位置周圍測量值之間空間關系的擬合模型。

2.3 監(jiān)測結果評價指標

本文研究以標準差和均方根差作為精度指標，對研究區(qū)域形變結果精度進行定量評價[24]。標準差(standard deviation，SD)，是離均差平方的算術平均數的算術平方根，是測繪領域經常用來統(tǒng)計分布程度測量依據，用σ表示，σ越小表明數據離散程度越小，穩(wěn)定性越高，其計算公式為

圖3 2015-06-17—2015-10-03間的DInSAR形變結果Fig.3 DInSAR deformation results between 2015-06-17 and 2015-10-03

(10)

均方根差(root mean square error，RMSE)是一種常用的測量數值之間差異的量度，值越小說明觀測值與真值之間誤差越小，所得結果精度越高，其計算公式為

(11)

許諾是在電視上看到海島發(fā)生地震的新聞后，匆忙放下手頭的工作趕過來的，他提著一顆心飛了一路，看到丁小慧毫發(fā)無損地站在他面前，他才松了口氣，緊緊地抱了抱她，“你是想嚇死我??？”

3 形變結果分析

3.1 融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變監(jiān)測結果

3.1.1 基于相干性系數的D-InSAR地表形變結果

圖4 2015-06-29—2015-10-27間的DInSAR形變結果Fig.4 D-InSAR deformation results between 2015-06-29 and 2015-10-27

以2015年6月17日—10月3日和2015年6月29日—10月27日數據為例，將25景影像分別進行基于相干性系數的D-InSAR處理，圖3和圖4分別為兩個時間段的數據所得形變結果圖。為了更清晰地體現研究區(qū)域的形變結果，將萬年礦區(qū)與九龍礦區(qū)形變結果進行局部放大，兩者局部放大形變結果分別見圖3和圖4右側的上方和下方。從圖3和圖4中可以明顯看出，2015年6月17日—10月3日期間的形變量變化范圍為-81.51～72.88 mm；2015年6月29日—10月27日期間的形變量變化范圍為-130.46～77.06 mm。此外，通過對比分析圖3和圖4，能夠發(fā)現兩者地表形變量變化范圍具有細微差異，主要是由于兩者使用不同時間段的影像所致，但最終形變結果的沉降分布情況仍基本保持一致。從圖3和圖4的局部放大圖中可以明顯看出，萬年礦區(qū)有兩處沉降區(qū)域，九龍礦區(qū)有一處沉降區(qū)域。

圖5 礦區(qū)的形變速率圖Fig.5 Deformation rate diagram of the mining area

3.1.2 結合PS-InSAR的SBAS-InSAR技術的礦區(qū)形變結果

SBAS-InSAR技術相比于D-InSAR技術具有時序分析的優(yōu)勢，不僅可獲得地表形變量結果，還可得到年均形變速率。將PS-InSAR技術與SBAS-InSAR技術相結合，并對25景SAR數據進行處理，可得到研究區(qū)域2015年6月—2017年3月間的年均形變速率，圖5繪制了研究區(qū)域的形變速率結果。

從圖5的左側整體圖可以看出，峰峰礦區(qū)較多煤礦開采區(qū)域存在地表沉降現象，年均形變速率最高達到了35 mm/a；而圖5右側局部放大區(qū)域分別繪制的是萬年礦區(qū)和九龍礦區(qū)的形變速率變化情況，能夠明顯看出由煤礦開采引起的沉降特征呈漏斗形狀，位于沉陷中心位置的沉降量最大，沉降速率也最快，且年均沉降速率最大可達-103.82 mm/a，符合礦區(qū)的沉降特征。

3.1.3 融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變結果

圖6 萬年礦區(qū)形變時序結果Fig.6 Time series results of mine deformation in Wannian

將25景Sentinel-1影像根據3.1.1節(jié)與3.1.2節(jié)中的兩種方法分別得到礦區(qū)形變結果，兩者形變結果中研究區(qū)域的沉降趨勢變化和沉降位置分布基本一致。利用Kriging插值法將這兩種方法的形變結果進行融合得到礦區(qū)最終的地表形變結果，綜合分析可反映兩個礦區(qū)的沉降分布情況和沉降變化特征，圖6和圖7分別繪制了萬年礦區(qū)和九龍礦區(qū)在2015年6月—2017年3月間的形變時序結果，其中正值代表地表向上發(fā)生抬升，負值代表地表向下發(fā)生沉降。從圖6和圖7中可以看出，地表形變在空間上的分布情況是隨時間變化的，且累積沉降量隨時間逐漸增加，并逐漸呈現漏斗形狀，最大累積沉降量可達176.39 mm?？傮w而言，萬年礦區(qū)和九龍礦區(qū)的地表沉降變化特征符合煤礦開采的變化特征，且融合后的累積沉降量和SBAS-InSAR得到的年均形變速率中的地表沉降分布情況吻合。

圖7 九龍礦區(qū)的形變時序結果Fig.7 Time series results of mine deformation in Jiulon

3.2 融合形變結果分析

相干系數圖反映了兩幅影像的相關程度，可以用于評價SAR影像精確配準后所得到干涉條紋圖的質量，主要是通過計算干涉條紋圖的相干系數實現。并有研究表明，相干系數圖質量的好壞還與季節(jié)的變化有關，主要是由于植被覆蓋的影響，夏季的相干系數圖質量普遍偏低，冬天相干系數圖質量通常較高。根據1.1節(jié)中相干系數理論知識可知，相干系數的值越接近1，則兩幅影像的相干性越高。礦區(qū)開采往往伴隨著大范圍的地表沉降變化，可以將其與光學影像相結合進行露天礦區(qū)的開采識別[25]。圖8為研究區(qū)域的相干系數圖，從圖8中可以看出，除了部分水體和礦區(qū)大量級開采區(qū)域的相干性較低，其余整體的相干性都較高。從而說明，相干系數圖的質量較優(yōu)，所得地表形變結果準確性較高。

圖8 研究區(qū)域的相干性系數圖Fig.8 Coherence coefficient diagram of the study area

3.2.2 典型區(qū)域形變結果剖面分析

(1)萬年礦區(qū)的剖面分析。為了研究礦區(qū)沉降漏斗的演化過程，將萬年礦區(qū)沿著開采面的橫縱剖面分別進行形變結果分析，圖9繪制了地表形變橫、縱剖面曲線結果。從圖9中可以明顯看出萬年礦區(qū)橫縱剖面的沉降情況，2016年10月9日、2016年12月20日和2017年3月2日這三景SAR數據的沉降趨勢高度一致，在132610工作面開采區(qū)出現了明顯的沉降漏斗。其中，圖9(a)中，橫向剖面結果顯示有一個沉降中心，最大沉降量達到了88 mm左右；圖9(b)中，縱向剖面結果顯示有兩個沉降中心，這與時序形變結果顯示的沉降漏斗基本一致，最大沉降量達到了98 mm左右。

(2)九龍礦區(qū)的剖面分析。與萬年礦區(qū)研究沉降漏斗的演化過程一樣，將九龍礦區(qū)沿著開采面的橫縱剖面分別進行形變結果分析，圖10繪制了地表形變橫、縱剖面曲線結果。從圖10中可以觀察到九龍礦區(qū)橫縱剖面的沉降情況，2016年10月9日、2016年12月20日和2017年3月2日這三景SAR數據的沉降趨勢高度一致，在15235工作面開采區(qū)出現了明顯的沉降漏斗。其中，圖10(a)中，橫向剖面結果顯示有一個沉降中心，最大沉降量達到了126 mm左右；圖10(b)中，縱向剖面結果顯示有一個沉降中心，最大沉降量達到了110 mm左右。

圖9 萬年礦區(qū)地表形變剖面曲線圖Fig.9 Surface deformation profile curve of Wannian mining area

圖10 九龍礦區(qū)地表形變剖面曲線圖Fig.10 Surface deformation profile curve of Jiulong mining area

3.3 融合形變結果評估

表2 水準數據和形變結果誤差對比Table 2 Comparison of errors between leveling data and deformation results

4 結論

為了確保結果的準確性，采用了多種InSAR技術對研究區(qū)域進行形變監(jiān)測。首先根據PS-InSAR方法的GCP點改進了SBAS-InSAR方法中軌道精煉控制點選取的問題，并利用Kriging插值法將基于相干系數的D-InSAR方法處理結果填補SBAS-InSAR技術中的空缺值，使得最終得到的礦區(qū)形變結果更加全面、有效。

(1)通過基于相干性系數的D-InSAR地表形變結果對礦區(qū)進行形變監(jiān)測，結果表明，2015年6月17日—10月3日期間的形變量變化范圍為-81.51～72.88 mm；2015年6月29日—10月27日期間的形變量變化范圍為-130.46～77.06 mm。并測得萬年礦區(qū)有兩處沉降區(qū)域，九龍礦區(qū)有一處沉降區(qū)域。

(2)利用Kriging插值法融合D-InSAR和SBAS-InSAR技術的礦區(qū)地表形變結果，綜合分析可反映萬年礦區(qū)和九龍礦區(qū)這兩個典型礦區(qū)的沉降分布情況和沉降變化特征。結果表明，地表形變在空間上的分布情況是隨時間變化的，且累積沉降量隨時間逐漸增加，并逐漸呈現漏斗形狀，最大累積沉降量可達176.39 mm。

(3)將本文結果與水準測量數據進行對比分析。結果表明，通過本文方法得到的形變結果和水準處理數據基本一致，且誤差均小于13 cm，監(jiān)測精度滿足煤礦開采相關測量技術要求。

融合D-InSAR方法和SBAS-InSAR方法可以避免單一方法所產生的偏差，綜合利用兩種方法的優(yōu)勢可彌補另一種方法中的不足，從而提高形變監(jiān)測結果的準確性。但是對于大形變量沉降的礦區(qū)監(jiān)測，由于C波段可監(jiān)測形變梯度的限制，將在今后的研究利用其他波段的衛(wèi)星影像對礦區(qū)進行形變監(jiān)測。