基于DS-InSAR的烏達煤田火區長時序地表形變監測與分析

李 柱， 范洪冬， 高彥濤， 許耀宗

(1.中國礦業大學礦山生態修復教育部工程研究中心，徐州 221116； 2.河南省地質礦產勘查開發局測繪地理信息院，鄭州 450006)

0 引言

煤火一般是指地下煤層火災、近地表煤層火災、煤矸石火災和煤堆火災[1]。大多數煤火是由于不適當的采礦造成的，也有一些是因處理不當的礦山廢棄物、閃電和森林火災引起的[2]。作為一種全球性災害，煤火對人類健康、煤礦安全和基礎設施造成了嚴重的負面影響。煤火的燃燒不僅消耗了大量寶貴的煤炭資源，同時引發地面沉降、裂縫、裂隙等地質災害，釋放的大量有毒氣體(CO，NO2，SO2，H2S)和微量元素(Hg，F，As，Se)還會造成嚴重的環境問題[3]。因此，全面、詳細以及準確地監測煤火，對煤火的演變趨勢、自燃規律和燃燒狀態進行研究以及開展煤火治理活動具有重要的意義。

遙感方法具有獲取信息速度快、監測區域面積大、周期短、全天時、全天候和經濟效益高的特點，這些獨特的優勢可以彌補傳統煤火探測與監測技術的不足[4]。利用遙感方法進行地面沉降監測和溫度反演，已經成為煤火探測與監測領域的一個重要方向。國內外學者利用熱紅外遙感技術對烏達煤田煤火區進行了大量研究[5-8]，可以較高精度地探測出煤火區及其演變趨勢，但熱紅外遙感技術很難實現對深部煤火、煤火燃燒階段以及煤火引起的地表形變的監測，這使得熱紅外遙感技術監測煤火時存在一定的不足。

差分合成孔徑雷達干涉測量(differential interferometric synthetic aperture Radar，D-InSAR)技術作為一種先進的對地觀測手段已在各種變形監測領域得到了廣泛的應用，但是該技術易受時空失相關及大氣延遲等因素影響，難以在長時序上得到高精度的地表形變結果[9]，而近年來發展的分布式目標合成孔徑雷達干涉測量(distributed scatterer interfero-metric synthetic aperture Radar，DS-InSAR)方法則有效地解決了長時序上低相干地區監測點位密度不足的問題[10]。目前，國內外一些學者利用合成孔徑雷達干涉測量(interferometric synthetic aperture Radar，InSAR)技術對煤火區的探測與監測進行了相關研究，大多數研究側重于將InSAR監測的地表形變結果作為輔助信息與地表熱異常信息相結合進行煤火區域的探測[11-14]，而對煤火引起的地表形變規律研究較少[15-16]； 另外，對于烏達煤田煤火的研究多數是利用熱紅外遙感技術進行探測與監測，只有少數研究采用InSAR技術[17-19]，并且缺乏近些年的地表形變監測結果。

針對上述問題以及煤田火區的地表形變復雜、失相干嚴重等現象，本文研究了一種基于DS-InSAR技術的煤田火區監測方法，利用該方法對烏達煤田火區2017年3月—2019年4月的63景Sentinel-1A影像進行了處理，通過監測結果對烏達煤田及其煤火區的地表形變進行了分析，驗證了該方法的適用性。

1 研究區域概況及數據源

1.1 研究區域概況

烏達煤田位于中國內蒙古自治區西南部的烏海市烏達區，如圖1所示。該區東鄰黃河與鄂爾多斯黃土高原，南部是烏達工業園區，西鄰賀蘭山脈北部邊緣，北鄰烏蘭布和沙漠。煤田呈南北向斜狀，面積約為55 km2，包括五虎山、黃白茨和蘇海圖3個煤礦，海拔從1 100～1 300 m不等，平均高程為1 200 m。氣候為典型的溫帶大陸性氣候，年降水量為168 mm，蒸發量為3 500 mm左右，常年干旱少雨。

圖1 研究區域Fig.1 The study area

1.2 數據源

Sentinel-1A衛星由歐州航天局于2014年4月3日發射，采用近極地太陽同步軌道，軌道高度為693 km，重訪周期為12 d，搭載了一個C波段的合成孔徑雷達(synthetic aperture Radar，SAR)，能提供超寬幅、干涉寬幅、條帶以及波浪4種成像模式。本文選取2017年3月—2019年4月期間的63景干涉寬幅模式成像下的SAR影像為數據源，其距離向和方位向的像元尺寸分別為2.3 m和14.0 m。為去除干涉圖中的地形相位，借助空間分辨率為90 m的SRTM DEM數據作為外部數字高程模型(digital elevation model，DEM)進行差分干涉處理。

2 DS-InSAR方法原理

Ferretti等[20]于2011年提出了第二代永久散射體技術SqueeSAR。拉開了DS-InSAR方法研究的序幕。與傳統的時序InSAR方法不同的是，DS-InSAR方法在其基礎上增加了同質點選取和相位優化2個步驟，有效地解決了傳統時序方法相干點選取數目少和空間分布不均勻的問題。

2.1 同質點選取

(1)

在勻質區域的單視復數SAR影像服從瑞利分布，其變異系數CV為常數，即

(2)

假設像元P后向散射性穩定，則式(1)可以轉換為：

(3)

2.2 特征值分解相位優化

(4)

式中: 向量y=[y1y2,…,yN]為分布式目標的同質點在N景SAR影像上的復數觀測向量x=[x1,x2,…,xN]經過振幅歸一化處理得到；NSHPs為分布式目標同質點的數目； H為Hermitian共軛轉置。

(5)

DS-InSAR數據處理具體流程如下： ①根據設置的時空基線閾值選擇干涉對，并生成干涉圖； ②利用FaSHPS方法對強度影像進行同質點識別，并通過特征值分解法對干涉圖進行相位優化； ③利用外部DEM去除干涉圖中的地形相位生成差分干涉圖，并利用優化后的相位計算同質點的時空相干性，設定相干性閾值選取分布式目標(distributed scatterer，DS)點； ④利用差分干涉圖對DS點上的相位進行解纏，并建立解纏相位與地表形變速率、DEM誤差和殘余相位之間的觀測方程； ⑤利用奇異值分解法估算地表形變相位和DEM殘余相位，并對殘余相位進行時空濾波分離出非線性形變相位、大氣相位和噪聲相位； ⑥將線性形變相位與非線性形變相位相加解算研究區域視線向(line of sight，LOS)的地表時序形變。數據處理流程如圖2所示。

圖2 數據處理流程Fig.2 Data processing flow chart

3 監測結果與分析

3.1 監測結果對比

分別采用TCP-InSAR與DS-InSAR方法對覆蓋烏達煤田的63景Sentinel-1A數據進行處理，獲得了2017年3月—2019年4月時間段內的地表時序形變信息，年平均形變速率如圖3所示。從圖3可以看出，二者雖在空間上形變趨勢具有較好的一致性，但DS-InSAR監測結果能夠更好地反映出煤田的地表形變情況，并對研究區域形變明顯的8個區域進行A—H編號。

(a) TCP-InSAR (b) DS-InSAR

TCP-InSAR方法共選出了156 853個監測點； DS-InSAR方法共選出了194 585個DS點。從選取的監測點數量而言，DS-InSAR能夠選出更多的DS點，監測點位密度比TCP-InSAR提高1.24倍，并且在低相干區域也可以獲得足夠數量的監測點； 從監測點位的空間分布來看，DS-InSAR方法選出的DS點在保證點密度的同時分布更均勻，在一定程度上提高了形變監測結果的解算精度，為研究區提供了更為詳細準確的形變信息，有效地彌補了在低相干區域TCP-InSAR等傳統時序形變監測方法中存在的缺陷。

為驗證DS-InSAR監測結果的可靠性，通過對比分析TCP-InSAR和DS-InSAR獲得的形變速率值來進行交叉驗證。本文以TCP點為基準，選取最鄰近的DS點作為同名點，共篩選出95 760對同名點，通過同名點對的形變速率值及其差異的統計分析，繪制了形變速率相關性和差異分布直方圖，如圖4所示。通過形變速率相關性(圖4(a))可知，DS-InSAR與TCP-InSAR方法獲得的地表形變速率之間具有較強的相關性，利用Pearson計算的相關系數R為0.84。由形變速率差異分布直方圖(圖4(b))可知，2種方法獲得的形變速率差異較小，兩者之間的標準差為5.83 mm/a。因此，在缺少實測水準數據的情況下，與TCP-InSAR方法監測結果的交叉驗證，表明了DS-InSAR地表形變監測結果是較為可靠的。

3.2 烏達煤田地表形變分析

為了能夠全面分析烏達煤田地表形變特點，本文利用DS-InSAR對烏達煤田進行了長時序的地表形變監測，時序形變累積如圖5所示。圖5中的14幅圖是每隔一段時間相對于起始時間2017年3月17日的累積形變量。從圖5中可以看出，在監測時間段內烏達煤田存在嚴重的地表形變現象。在空間上，沉降區域分布不均勻。在時間序列上，隨著時間的推移研究區域內逐漸出現了形變量級與沉降范圍大小不同的多個沉降區域。其中，A，B，C這3個區域的形變現象最為嚴重，沉降范圍最廣，這3個區域的地表形變在整個監測時間段內持續存在，地表形變量與沉降范圍越來越大，導致B和C區域從2018年8月15日之后逐漸連接形成一個大的沉降區域； G區域的地表形變是從2017年5月16日開始產生，前期形變變化較小，到后期形變變化加快； D和F區域從2017年8月20日之后地表開始產生形變，并逐漸形成多個沉降中心； H區域的地表形變從2017年12月18日開始發生地表形變； E區域產生地表形變的時間最晚，但地表形變變化較快，從2018年10月14日—2019年4月12日半年的時間里產生了與D區域相同量級的形變。為更詳細具體地分析地表形變情況，利用圖3(b)地表形變速率圖對其進行定量分析。從圖5中可知，地表形變最為嚴重的A，B，C這3個區域的最大形變速率分別為-178.4 mm/a，-166.3 mm/a和-215 mm/a，在這3個沉降區域的沉降中心，由于形變梯度大，地表形變量在監測時間段內超過時序InSAR的監測能力，失相干現象嚴重，無法得到足夠密度的DS點； 另外，F和G這2個區域同樣存在較為嚴重的形變，最大形變速率分別為-96.5 mm/a和-100.9 mm/a； 相比之下，D，E和H區域的形變量相對較小，最大形變速率分別為-46.5 mm/a，-53.9 mm/a和-48.5 mm/a。

圖5 烏達煤田地表形變時序累積圖

3.3 煤田火區地表形變分析

烏達煤田嚴重的地表形變現象主要與煤火燃燒及人工采挖有關[17,19]，根據文獻[7]中2018年1月底實地調查獲取的部分區域的明顯煙點和煤火點，這些實測的煤火點附近地區伴隨著地表形變現象。為了研究煤火燃燒引起的地表形變，對實測煤火區的地表形變進行重點分析。

3.3.1 煤田火區地表形變時序分析

為分析煤田火區的時序地表形變規律，將實測煤火點編號為H1—H6(圖3(b))，對應D，A和G這3個沉降區域。選取3個區域中實測煤火點附近且累積形變量最大的5個DS點作為特征點分析地表形變，編號記為P1—P5(圖3(b))。為了驗證特征點時序形變的可靠性，從TCP-InSAR監測結果中選取距離特征點最近的TCP點作為驗證點對其進行驗證，由于TCP-InSAR選點的原因，只在P3和P5點旁選了P3′和P5′點。圖6(a)—(c)分別為D，A和G這3個沉降區域的特征點及驗證點時序形變圖，特征點與驗證點之間的時序形變趨勢基本吻合，具有很好的一致性，2點間的時序形變值整體存在偏大或偏小的問題，但考慮到2點之間具有一定的距離，屬于合理范圍，所以可以利用選取的特征點進行地表時序形變分析。

由圖6可知，3個沉降區域的形變量級與時序形變趨勢各不相同，不同時間段地表形變變化快慢也不相同。P1從2018年5月底才開始發生形變并持續到監測時間結束，這段時間內地表形變變化快，煤火燃燒劇烈。處于同一沉降區域的P2，P3和P4開始發生形變的時間各不相同，P2點最晚，從2017年9月開始，并且中間一段時間因煤火處于休眠階段未發生明顯的地表形變現象； P3點地表形變變化緩慢，但在整個監測時間段內形變持續發生； 相比之下，P4點從2017年6月開始出現形變現象，且在之后時間段內出現了3個未發生形變的時期。P5點的時序形變趨勢更加復雜，多個時間段內出現形變停滯現象，但在發生形變的時期內形變變化較快。煤火區的這些復雜地表時序形變現象歸因于煤火燃燒階段的復雜性、燃燒強度的非均質性，并且可以從這些特征點時序上的地表形變變化快慢推斷煤火的燃燒強度及燃燒階段。從整體上分析這5個特征點時序地表形變規律，可以發現秋冬季節的地表形變變化相對較快，這是由于烏達煤田處于干旱半干旱地區，秋冬季節干旱少雨，煤層易發生火災[23]，導致地表形變變化較快。

3.3.2 典型煤火區地表形變特征分析

為分析煤田火區的地表形變特征，選取3個沉降區域中形變最嚴重的A區作為典型煤火區。為清晰直觀地表征區域A地表形變的空間形態，將地表累積形變圖(圖7(a))進行克里金插值，并將插值之后的圖轉換成三維形變圖(圖7(b))。從圖7(b)可知，區域A的地表形變在空間上變化較快，整體形變規律復雜，邊緣形狀不規則，且具有多個發育程度不同的沉降中心。從三維形變圖的形變等值線中可知，這些沉降中心附近的等值線密集，向四周逐漸變得稀疏，對應的沉降中心凹陷較深，向四周逐漸平緩，邊緣近似長橢圓形，并且這些沉降中心的形變延伸方向隨著煤火蔓延和燃燒強度而各不相同。

(a) 剖面線示意圖 (b) 三維形變圖

為定量分析典型煤火區的地表形變特征，在圖7(a)區域A的沉降中心附近給出了2條剖面線a1和a2，并沿剖面線提取地表累積形變量，結果如圖8所示。由圖8可知，剖面線a1穿過了3個形變量級與沉降范圍大小不同的沉降中心，最嚴重的沉降中心累積形變量達到-311.6 mm； 剖面線a2經過2個發育程度相似的沉降中心，最大累積形變量從左到右分別為-369.5 mm和-355.2 mm； 這些沉降中心到邊緣的形變量變化快，形變梯度較大。

圖8 典型煤火區剖面累積形變圖

4 結論

本研究利用DS-InSAR對烏達煤田火區2017年3月—2019年4月采集的63幅Sentinel-1A影像進行了處理，分析了長時序地表形變，得到如下結論：

1)與TCP-InSAR技術相比，DS-InSAR監測點密度提高1.24倍，且分布均勻，能夠更好地反映出研究區域的地表形變情況。研究區域存在多個形變量級與沉降范圍大小不同的沉降區域，且在空間上分布不均勻，最大地表形變速率約為-215 mm/a。

2)監測時間段內烏達煤田存在因煤火燃燒導致的地表形變。煤火區的地表時序形變規律復雜，秋冬季節地表形變變化相對較快，并且利用地表時序形變可以推斷出煤火的燃燒狀態，形變嚴重的煤火區出現多個形變延伸方向且發育程度不同的沉降中心。

3)本研究僅利用DS-InSAR獲取的結果對烏達煤田及其煤火區的地表形變進行了分析，將來還需要與現場實測資料以及開采沉陷理論知識相結合，進一步研究煤火燃燒導致的地表形變機理，最終實現對煤田火區地表形變全面詳細準確的監測。