基于DS-InSAR的西部礦區地表時序沉降監測與分析

王 波 譚志祥,2 鄧喀中,2

(1.中國礦業大學環境與測繪學院,江蘇 徐州 221116;2.自然資源部國土環境與災害監測重點實驗室,江蘇 徐州 221116)

多年來對東部礦區煤炭資源的高強度開發,使得 東部礦區煤炭產量不斷下降。為滿足國民經濟發展對于煤炭資源的需求,我國加大了對西部礦區煤炭資源的開發利用力度。與東部礦區相比,西部礦區煤炭資源儲量更大,而且具有煤層厚度大、傾角小等優點,因此有著更大的生產潛力[1]。然而由于歷史原因,西部礦區煤礦開采時間比東部礦區要短,因此針對西部礦區開采引起的地表形變的監測工作開展較少,相應的研究成果也比較匱乏。而且由于東西部礦區之間地質采礦條件的巨大差異,導致以東部礦區為研究對象所得到的相關研究成果在西部礦區適用性不強[2],這給西部礦區地表形變預計、村莊保護煤柱留設以及各類建筑物保護帶來了困難,不利于煤礦生產工作安全高效進行。因此有必要針對西部礦區特定的地質采礦條件,開展礦區地表形變監測工作,并以此為基礎進行西部礦區地表形變規律和機理的研究,為后續煤礦生產工作提供技術依據。

水準測量等傳統的地表形變監測方法在西部礦區的應用具有較大的局限性,主要體現在受西部礦區高原堆積型沙丘地貌的影響,地面水準點以及地表特征點的選取和觀測困難。而且傳統的測量方法往往費時費力,難以進行大范圍面狀地表形變監測[3]。合成孔徑雷達差分干涉測量技術(Differential InSAR,D-InSAR)能夠不依賴地面特征點來獲取地表形變信息,因此被廣泛應用于地表形變的監測與分析中[4]。然而該技術極易受到時空去相干因素的影響,難以獲取礦區開采引起的地表形變范圍。針對這一問題,科研工作者將注意力集中在了那些能夠在長時間序列合成孔徑雷達(Synthetic Aperture Radar,SAR)影像中保持高相干性的點目標上,提出了永久散射體干涉測量技術(Permanent Scatterer InSAR,PS-InSAR)[5],在此基礎上進行了相干點目標識別[6]、三維相位解纏算法[7]等研究,將D-InSAR技術推進到了時間序列InSAR的新階段。以PS-InSAR和小基線集方法(Small Baseline Subsets InSAR,SBAS-InSAR)[8]為代表的常規時序InSAR技術在干涉對優化的基礎上,提取分析了影像中受時空去相干影響較小而表現出相對穩定散射特性的高相干像元在整個監測時段內的變化,有效克服了D-InSAR技術在應用過程中的失相干問題。然而礦區及周邊地表的覆蓋多為植被和沙地,導致常規的時序InSAR技術在礦區及其周圍難以獲得足夠數量、分布均勻的高相干點,地表形變監測精度受到了很大限制。

近些年發展起來的融合分布式目標的時序InSAR技術[9]將研究的重點從散射特性相對穩定的高相干點目標擴大到了散射特性中等的分布式目標,可同時利用這兩種目標進行地表形變監測,從而極大地提高了監測區域點目標的密度,有效克服了常規時序InSAR技術在低相干區域由于無法獲取足夠數量的點目標而無法對地表形變進行準確監測的不足。該技術先是在地震形變監測、山區滑坡災害監測、冰川運動監測等方面進行應用并取得了良好的效果[10],證明了該方法與常規時序InSAR技術相比具有優越性。之后部分學者開始嘗試將該技術應用到礦區地表形變監測中并取得了一定的成果[11-18]。但是上述研究的試驗區域大多位于東部礦區,對于該技術在西部礦區地表形變監測方面的應用研究涉及較少,缺乏針對性的研究。為此,本研究在現有成果的基礎上,利用一種針對西部礦區高原堆積型沙丘地貌的基于DS-InSAR的時序地表形變監測方法,獲取了西部礦區內蒙古石拉烏素煤礦地表形變及其時空分布信息,并驗證了該方法的可靠性,豐富了DS-InSAR技術在西部礦區地表形變監測應用方面的研究。

1 研究區概況與數據源

1.1 研究區概況

石拉烏素煤礦位于內蒙古自治區鄂爾多斯市伊金霍洛旗東勝煤田內,礦區面積約70 km2,設計生產能力為1 000萬t/a。礦區位于鄂爾多斯高原的中南部,具有典型的高原堆積型沙丘地貌特征,地表大部分被第四系風積沙所覆蓋,局部有白堊下統地層露出,植被稀疏,為沙漠—半沙漠地區。本文研究區域所在的221采區位于石拉烏素煤礦的中部(圖1),采區地形總體趨勢是東高西低、北高南低,地面標高+1 332～+1 411 m,采區內煤層埋深近700 m,上覆大厚度的中生代砂巖,上覆巖層性質具體如表1所示,其地質采礦條件與東部礦區相比存在明顯差異。在監測時段內221采區有17、18、06A 3個工作面開采,采區北部的17、18工作面采用綜合機械化一次采全高工藝開采,平均采厚約5.40 m,南部的06A工作面采用綜合機械化放頂煤工藝開采,平均采厚約9.78 m。各工作面的參數取值見表2。

表1 煤層上覆巖層性質Table 1 Properties of overlying strata of coal seam

表2 221采區工作面參數Table 2 Parameters of working face in 221mining area

圖1 研究區地理位置Fig.1 Geographical locationof the study area

1.2 數據源

Sentinel-1衛星由歐空局發射,并免費向全球提供衛星影像下載服務,為SAR相關研究提供了極大的便利。本研究選用了2016年9月—2018年6月共52景覆蓋研究區域的IW模式下VV極化的Sentinel-1A升軌影像。Sentinel-1衛星影像的分辨率為5 m×20 m,衛星的重訪周期為12 d。此外本研究采用30m分辨率的航天飛機雷達測圖計劃(ShuttleRadar Topography Mission,SRTM)數字高程模型(Digital Elevation Model,DEM)從干涉圖中去除地形相位的影響和地理編碼。

2 監測方法選擇及其原理

2.1 高相干點目標選取

從時序SAR影像中準確地識別高相干點目標是時間序列InSAR分析的重要環節,也是實現聯合高相干點目標與分布式目標精確提取地表形變信息的前提條件。受西部礦區高原堆積型沙丘地貌的影響,研究區域內能夠在長時間序列中保持相位穩定的永久散射體目標的數量稀少,但存在數量較多在短時間內經濾波后相位失相關性變小而保持高相干性的點目標。本研究對所選SAR影像進行處理得到濾波后的差分干涉圖集,通過幅度差分離差法[19]選取出了此類點目標作為高相干點目標。幅度差分離差法原理可用下式表示

式中,DΔA為幅度差分離差值;σΔA為主、輔影像幅度差的標準差;μA為差分干涉圖幅度的均值。通過設定幅度差分離差閾值,便可實現高相干點目標的選取。

2.2 分布式目標識別與相位優化

2.2.1 FaSHPS同質像元識別

受礦區特定地質條件的影響,研究區域內高相干點目標的空間分布極不均勻,大部分高相干點目標集中分布在地表裸露巖石較多的區域,而在地表大范圍的積沙區域內此類點目標數量十分稀少。為實現對研究區完整地表形變信息的提取,必須在高相干點目標的基礎上充分利用廣大積沙區域內的分布式目標。參數假設檢驗方法中的快速同質點選取方法(Fast Statistically Homogeneous Pixel Selection,FaSHPS)[20]具有同質像元選取效率高且提取出的同質像元之間的異質性小的優點,為此本研究選用該方法對研究區域內的同質像元進行識別,并通過設置時間相干性閾值實現對分布式目標的提取。FaSHPS方法的思路為,假設時序SAR影像數據集內任一像元時間維上振幅的均值服從高斯分布,從而可以將傳統的假設檢驗方法轉化為置信區間的估計問題,只要利用簡單的邏輯計算來判斷兩個像元服從的分布函數是否相同,便可實現同質像元的識別。

設由N景時序SAR影像組成數據集,其中任一像元L在時間維上的振幅均值為(L),由中心極限定律可知,隨著樣本數目N加大,振幅均值(L)趨于高斯分布,此時其區間估計可表示為

式中,P{·}表示概率;μ(L)為像元L的數學期望;z1-α/2為標準正態分布中置信度為α時對應的分位點;Var(A(L))為像元L時間維上振幅的真實方差。

當樣本的數量足夠大時,單視復數SAR影像的時序振幅信息近似服從瑞利分布,此時變異系數CV為常量,其值約為0.52。若像元在時間維上的后向散射特性穩定,則式(2)可表示為

當給定置信度α時,式(3)便給定了一個區間,若待估像元在時間維上的平均強度值在目標像元所對應的置信區間內,則待估像元與目標像元為同質像元。

2.2.2 特征值分解相位優化

分布式目標像元內地物勻質分布[21],其回波信號由像元內各子散射體相干疊加而成,易受時空失相干影響,相位穩定性較差。因此需要基于分布式目標像元散射特性的整體統計特征來去除噪聲相位,實現對相位的優化。特征值分解法在相位優化的過程中無需迭代計算,因此其優化效率較高,同時該方法的相位優化質量也較好。故本研究選用該方法對分布式目標像元相位進行優化。特征值分解法的核心思想為從相干矩陣中分離,估計出多元散射機制并得到其對應的相位分量,通過分離出主導散射機制的相位分量來達到優化相位的目的[22-23]。

式中,Λ為N階對角矩陣,對角線上的元素為相干矩陣的N個特征值,用λi表示;每個不同的特征值分別對應著不同的特征向量,用μi表示;特征值不同,其對應的散射機制也不同。特征值越大,在像元整體散射信號中其對應的子散射體所占比重就越大。因此將最大特征值對應的特征向量作為主散射體,其余的特征向量作為噪聲,就可以實現同質像元相位的優化,此時相干矩陣可以表示為

在完成同質像元識別以及相位優化之后,通過計算相位優化前后干涉相位之間的擬合度便可對相位優化的質量進行評價。相位優化質量評估完成之后,設置時間相干性閾值,從分布式目標中篩選出相干性較好、信噪比較高的部分。將其與篩選出的高相干點目標一并融入時序InSAR處理流程中,便可對整個研究區域地表形變進行高精度監測。針對西部礦區高原堆積型沙丘地貌的基于DS-InSAR的時序地表形變監測方法的主要流程如2所示。

3 監測結果與可靠性分析

3.1 SBAS-InSAR與DS-InSAR監測結果對比

本研究分別利用傳統的SBAS-InSAR方法和圖2所示的基于DS-InSAR的時序地表形變監測方法提取了2016年9月—2018年6月研究區地表沉降的時空分布信息。兩種方法獲取的研究區地表衛星視線向的年平均沉降速率如圖3所示。傳統的SBASInSAR方法共提取到了57 803個監測點,監測點數目較少,且空間分布極不均勻。SBAS-InSAR方法監測到研究區地表存在一處較為明顯的下沉盆地,最大年平均沉降速率出現在該下沉盆地的中心位置,為139 mm/a,在該下沉盆地的南部同樣監測到一處較為明顯的沉降區域,但由于在該沉降區域內獲取的監測點數量稀少,只能大致反映該區域的外圍輪廓,無法得到其整體的沉降信息。DS-InSAR方法共提取到了239 750個監測點,與傳統的SBAS-InSAR方法相比該方法提取到的監測點數量顯著增加,約為前者的4倍,且監測點空間分布均勻。DS-InSAR方法監測到研究區地表存在南北方向相鄰的2個明顯的下沉盆地,且在2個下沉盆地中北部下沉盆地的沉降更為明顯,沉降范圍也更大。在北部下沉盆地的中心位置監測到了研究區域內的最大年平均沉降速率,其值為179 mm/a。

圖2 基于DS-InSAR的時序地表形變監測方法數據處理流程Fig.2 Data processing flow of time series surface deformation monitoring method based on DS-InSAR

圖3 研究區地表衛星視線向年平均沉降速率(單位:mm/a)Fig.3 Annual average subsidence rate of LOS direction in study area

結合衛星影像底圖對比分析兩種方法的監測結果可知,受西部礦區高原堆積型沙丘地貌的影響,傳統的SBAS-InSAR方法只能在礦區地表裸露巖石較多的區域獲取足夠數量的高相干點目標,在積沙區域由于地物的后向散射特性較弱獲取的觀測點數量稀少,監測結果難以準確反映礦區地表整體的沉降信息,不能對沉降區域范圍進行準確探測。本研究基于DS-InSAR的時序地表形變監測方法有效克服了傳統SBAS-InSAR方法在西部礦區地表沉降監測中的不足,在充分利用裸露巖石較多區域的高相干點目標的基礎上,通過提取研究區域內的分布式目標顯著增加了積沙區域觀測點的數量,從而能夠對礦區地表形變進行更高精度地監測,獲取礦區地表整體豐富的沉降細節信息。

3.2 DS-InSAR時序監測結果

以SAR影像數據集中首期影像的日期為開始時間,計算其余影像對應時間相對于開始時間的累計沉降值,將觀測結果從衛星的視線方向轉換為垂直方向,即可得到礦區地表的時序累積沉降。為更加直觀地顯示整個監測時段內研究區地表沉降區域的形成與地下工作面開采之間的關系,將地下工作面的位置疊加到時序累積沉降圖上,并提取出研究區域內受采動影響地表發生沉降的區域,即可得到沉降區域的累計沉降及沉降細節信息。研究區地表沉降區域范圍如圖4所示,沉降區域內的累計沉降及沉降細節信息如圖5所示。

圖4 研究區地表沉降區域范圍Fig.4 Scope of ground subsidence in the study area

圖5 研究區地表累計沉降及沉降細節信息Fig.5 Surface cumulative subsidence and sunsidence details in the study area

由圖5可知:地表發生沉降的區域與地下工作面的位置一致。17、18兩個工作面的開采導致研究區地表北部下沉盆地形成,本次試驗監測到的最大累計沉降出現在該下沉盆地的中心位置,為402 mm。受06A工作面開采的影響,研究區地表形成了南部下沉盆地,該盆地內監測到的累計沉降較小,其中心位置的最大累計沉降約220 mm。監測結果表明,針對西部礦區高原堆積型沙丘地貌的基于DS-InSAR的時序地表形變監測方法在研究區域內獲取了更多數量、空間分布均勻的監測點目標,能夠對礦區地表開采引起的沉降區域的形成和發展過程進行直觀地描述,并準確地顯示地表沉降區域的范圍。

3.3 DS-InSAR監測結果可靠性分析

本研究選取了17、18工作面上方2條觀測線在4個不同時段內的水準實測數據,通過將觀測線上DS-InSAR監測結果與水準實測數據進行對比分析來驗證其可靠性。觀測線與工作面的位置關系如圖6所示。因水準實測數據的日期與SAR影像日期并非完全一致,故將DS-InSAR監測結果按時間進行插值使其與水準實測數據的日期保持一致。本研究在2條觀測線上共選取了93個監測點目標,觀測線上DS-InSAR沉降監測結果與水準實測數據對比如圖7所示。

圖6 觀測線與工作面位置關系Fig.6 Position relationship between observation lines and working faces

由圖7可知:DS-InSAR監測結果與水準實測數據擬合較好,除極個別點該方法監測結果與水準實測數據偏差略大外,各觀測線上二者得到的累計沉降值及沿觀測線的沉降趨勢基本一致。為更加直觀地體現兩種方法監測結果間的相關性與誤差分布,繪制了如圖8所示的相關性及誤差分布圖,并通過計算每條觀測線上的誤差均值、絕對誤差均值及標準差,對監測結果的精度進行評價,結果見表3。

圖7 DS-InSAR監測結果與水準實測數據對比Fig.7 Comparison of DS-InSAR monitoring results and leveling measured data

圖8 DS-InSAR監測結果與水準實測數據相關性及誤差分布Fig.8 Correlation and error distribution between DS-InSAR monitoring results and leveling measured data

由圖8和表3可知:兩種方法監測結果間的相關性較高,皮爾遜相關系數為0.97,二者間的誤差較小且整體分布合理,證明針對西部礦區高原堆積型沙丘地貌的基于DS-InSAR的時序地表沉降監測方法在研究區地表非充分采動情況下的監測結果具有較高的可靠性。

表3 DS-InSAR監測結果精度評價Table 3 Accuracy evaluation of DS-InSAR monitoring results

4 沉降監測結果分析

4.1 地表沉降與地下開采關聯性分析

結合采工圖和DS-InSAR時序監測結果可知,2016年9月17工作面開始掘進,但此后較長一段時間內地表并未監測到沉降,這是由于工作面的采深大且上方覆有大厚度的中生代砂巖,工作面開采影響尚未傳播到地表。地表開始發生沉降時,工作面推進距離已達采動距,之后隨著工作面繼續推進,地表沉降區域的范圍和累計沉降值不斷增大。DS-InSAR監測結果顯示,2016年11月份左右研究區地表開始出現沉降,12月份地表已經形成一處較為明顯的沉降區域,2017年5月份17工作面開采結束后,沉降區域發展成為一個明顯的下沉盆地,其形狀近似為圓形,盆地內各點的下沉值隨著與盆地中心距離的增加而減小。2017年6月,17工作面西鄰的18工作面開始回采,由于兩個工作面之間相距很近,18工作面的開采并沒有在礦區地表形成新的下沉盆地,而是使原有下沉盆地的范圍和累計沉降值進一步增大,但受兩個工作面之間留設的120 m寬保護煤柱以及上覆巖層應力狀態的影響,地表仍未達到充分采動且下沉極不充分。2017年12月份18工作面回采結束,其平均采厚達5.47 m,但直至2018年6月監測結束該下沉盆地內監測到的最大沉降也僅為402 mm。

06A工作面從2017年9月開始掘進,因此根據2017年12月地表沉降監測結果可知,除了在開采已經結束的17、18工作面上方存在一處下沉盆地外,在其南部也出現了一處較為明顯的沉降區域。2018年4月06A工作面推進結束,到2018年6月該工作面上方已經形成了一處明顯的下沉盆地,且其沉降范圍已經與北部下沉盆地相當。通過比較南北兩個下沉盆地的形成和發展過程不難發現,從地表開始發生沉降到形成明顯的下沉盆地過程中南部下沉盆地所用時間要明顯短于北部下沉盆地。這是因為雖然3個工作面的采深和上覆巖層的性質十分接近,但06A工作面采用綜合機械化放頂煤工藝開采,其采厚達9.78 m,遠大于17、18兩個工作面,因此其開采影響傳播到地表所用的時間更短,但受該礦區特殊地質采礦條件的影響,直至監測時段結束該工作面上方地表同樣未達到充分采動。

4.2 南部沉降區地表動態變形分析

為進一步驗證基于DS-InSAR的時序地表沉降監測方法的性能,本研究利用該方法獲取了圖9所示06A工作面開采期間南部沉降區地表沿該工作面走向和傾向主斷面的動態沉降曲線,對該區域缺少的水準實測數據起到了很好的補充作用。為充分反映工作面推進過程中地表的動態變形情況,將沉降監測結果進一步處理得到了沿06A工作面走向主斷面的動態傾斜與動態曲率曲線,分別如圖10和圖11所示。由圖10、圖11可知:在06A工作面推進過程中,沿該工作面走向的傾斜和曲率變形不斷增大,至2018年4月底開始趨于穩定,傾斜在20號點附近(采空區中心)近似為0,在8號點和26號點附近(采空區邊界)達到最大,為0.41 mm/m;地表曲率變形較小,均小于0.02 mm/m2。

圖9 06A工作面走向與傾向主斷面動態沉降曲線Fig.9 Dynamic subsidence curves of main section of strike and incline of 06A working face

圖10 06A工作面走向主斷面動態傾斜曲線Fig.10 Dynamic inclination curves of main section of strike of 06A working face

圖11 06A工作面走向主斷面動態曲率曲線Fig.11 Dynamic curvature curves of main section of strike of 06A working face

5 結 論

本研究針對西部礦區高原堆積型沙丘地貌,以內蒙古石拉烏素煤礦為例,采用基于DS-InSAR的時序地表形變監測方法對2016—2018年間覆蓋研究區域的52景Sentinel-1A影像數據進行處理,獲取了研究區域的地表時序沉降信息,并將監測結果結合井下開采資料進行了分析。所得結論如下:

(1)與傳統的SBAS-InSAR方法相比,本研究方法在充分利用高相干點目標的基礎上,依靠FaSHPS方法對研究區域內的分布式目標進行提取,并通過特征值分解對其相位進行優化,之后將二者聯合進行時序處理來獲取地表形變信息,可有效提升研究區地表監測點密度和形變監測精度。

(2)在監測時段內,該方法監測到研究區地表形成了南北相鄰的2處明顯的沉降區域,最大沉降約400 mm。與井下開采資料對比分析表明,該方法監測到的地表沉降區域的位置和范圍較為準確,且沉降區域的形成及發展能與地下工作面開采進程相對應,可為開采沉陷影響邊界確定及其機理研究提供參考。

(3)驗證了基于DS-InSAR的時序地表形變監測方法的可靠性,并利用該方法獲取了礦區06A工作面開采期間南部沉降區域沿06A工作面走向和傾向主斷面的動態沉降曲線,對該區域缺少的水準實測數據起到了很好的補充作用,并將沉降監測結果進一步處理,得到了沿06A工作面走向主斷面的動態傾斜和曲率變形曲線,可為地表動態變形規律研究及形變控制提供依據。