基于SBAS InSAR的新疆哈密砂墩子煤田開采沉陷監測與反演

沙永蓮，王曉文，劉國祥，張 瑞，張 波

(1.西南交通大學地球科學與環境工程學院,成都 611756；2.北京理工大學重慶創新中心,重慶 401120;3.高速鐵路運營安全空間信息技術國家地方聯合工程實驗室,成都 611756)

0 引言

煤炭是我國能源結構的主體，對煤炭資源開采活動進行合理的規劃和管理對我國經濟建設和生態環境保障具有重要意義[1]。我國煤炭資源多是通過地下開采，而這種開采模式往往會破壞地層內部原有的力學平衡狀態，導致采空區上方地表發生移動和變形[2]。煤田礦區地表不均勻變形可能導致地面形成裂縫或者塌陷，破壞原有地質環境條件并引發地質災害，從而嚴重影響地表相關建筑物等基礎設施的安全使用，給礦區當地經濟可持續發展造成隱患[3]。

針對礦區不均勻地表形變的時空演化特征進行監測，可為當地地表沉陷隱患治理，安全生產保障和當地生態環境保護等提供關鍵支撐信息。目前，水準測量、全球導航衛星系統(GNSS)和合成孔徑雷達干涉(InSAR)等多種大地測量方法已在礦區地表沉陷監測中得到廣泛應用[4-5]。相比于傳統大地測量手段，InSAR具有全天候觀測、覆蓋范圍廣和觀測效率高等優勢[6-7]。常規基于兩景SAR影像進行干涉處理測量地表形變的兩軌法DInSAR (Differential InSAR)，由于無法對形變觀測模型中的誤差項(如地形誤差、軌道誤差和大氣延遲誤差等)進行校正，常被用于瞬時、大幅度地表形變(如地震和火山)測量[8-9]。近十年來，眾多學者基于多景SAR影像提出多種時序InSAR分析方法以克服常規DInSAR相位失相干、大氣延遲等影響，并實現長時間序列地表形變監測[10-11]。其中，小基線集(small baseline subset，SBAS)時序InSAR方法能夠有效利用冗余干涉對進行干涉相位建模，可有效抑制形變恢復模型中的誤差項，已在城市地表沉降、活動斷層滑移和滑坡運動監測等領域中得到廣泛應用。

SBAS時序InSAR技術已在煤田礦區地表形變特征監測、采空區位置和邊界探測等方面展現了獨特的應用潛力[12]。趙偉穎等[13]應用SBAS InSAR對9景ALOS PALSAR數據進行處理，獲取了徐州市某煤礦采動區中村莊在2007—2011年的累計沉降量，表明了利用SBAS InSAR監測礦區地表移動邊界及建筑物區域地表移動臨界區域的可靠性；李達[14]基于13景TerraSAR-X 數據分析了陜西省榆林市某煤礦區地表形變的時空演化特征，發現不同工作面開采時間與地表沉陷速率有一定相關性；閻躍觀等[15]基于Sentinel-1A 影像采用SBAS InSAR 測量了河北某煤礦的地表形變，并根據形變結果計算了該煤礦綜采面走向邊界角。有學者將煤田地下綜采面簡化為一個矩形位錯模型，建立了綜采面參數(開采深度，采煤厚度)與位錯模型張裂分量之間的關系，王亞男[16]根據該模型反演了SBAS InSAR測量的內蒙古上灣煤礦地表形變，探討了采煤深度和厚度與地表形變量特征間的聯系。朱煜峰[17]反演了SBAS InSAR觀測的江西省豐城市某煤礦地表累積沉陷量，得到了InSAR觀測期間位錯模型張裂分量；然后利用反演的矩形位錯張裂量模擬了礦區地表沉陷，發現與水準測量數據一致，證明了簡化位錯模型在煤礦地下開采參數反演中的可行性。綜上可見，由InSAR觀測獲取的高分辨率形變場能夠反演采煤作業面深度等關鍵參數，但是基于InSAR地表形變為約束反演煤礦開采量的研究目前還比較少見。在反演得到綜采面參數基礎上，能否結合煤層視密度構建綜采面參數與開采量之間的關系仍值得進一步研究。

本文以新疆哈密砂墩子煤田為例研究了基于時序SBAS InSAR觀測形變反演地下煤礦綜采面參數及煤礦開采量的可行性。砂墩子煤田是國家重點煤礦，設計產能5.00 Mt/a，自2012年開始試生產以來，至今仍在開采中；該煤田處于半干旱地區，礦區表面及周邊植被覆蓋較少，非常有利于采用時序InSAR準確恢復其地表形變特征。本研究選取2018年9月—2019年10月共30景Sentinel-1A升軌影像，首先利用SBAS方法恢復了該礦區一年時間內的地表形變；其次，以SBAS InSAR測量的地表形變速率為約束，基于Okada矩形位錯模型反演了該煤礦綜采面參數，進而根據煤層視密度、綜采面參數和開采量的關系估算了砂墩子礦開采量，并與已有礦區資料進行了對比驗證。

1 研究區概況及數據源

砂墩子煤礦位于新疆哈密市三道嶺礦區西部，地理坐標為E92°26′09″～92°32′10″，N43°08′33″～43°12′26″，自西向東構成“L”形，北距312國道5 km，南距蘭新線6 km (圖1)。圖上，白色虛線為礦區邊界線，紅框表示InSAR分析區域，光學影像為研究區Sentinel-2光學影像，白色三角為主井位置。井田東西走向長約7 km，南北最寬約8 km，礦井分布三個井筒，主井位置為E92°28′05″，N43°10′04″。從煤炭資源賦存來看，該礦區資源儲量877.44 Mt，設計產能5.00 Mt/a，一期設計產能3.00 Mt/a[18]。該礦區水文地質條件簡單，屬簡單型礦床，主采煤層頂底板以粉砂巖-泥巖為主[19]。整個礦區氣候干旱，晝夜溫差大，屬于典型溫帶大陸性干旱氣候；降雨量少但集中，雨期多在夏季，易形成短時洪流。

圖1 砂墩子煤礦地理位置及其周邊地形圖Fig.1 Map of Shadunzi coal site

本文選取30景Sentinel-1A升軌SAR影像對研究區域進行SBAS InSAR形變監測分析，影像獲取時間為2018年9月7日—2019年10月8日。經過裁剪后監測區域面積約為66.30 km2，SAR影像空間分辨率約為13.95 m(方位向)×2.33 m(距離向)，SAR影像中心入射角為39.34°。

2 算法原理

2.1 SBAS InSAR 時序地表形變提取

SBAS InSAR是一種經典的時序InSAR方法，它通過選擇具有短時空基線的SAR影像對進行差分干涉處理，并針對相干目標點的相位進行建模和解算[20]，在一定程度上能夠克服DInSAR常常面臨的干涉相位失相干影響。此外，由于SBAS方法應用了盡可能多SAR干涉對參與解算，其能夠改正DInSAR測量中的地形殘差、軌道誤差和大氣延遲誤差等多項誤差，進而得到高分辨率和空間上連續的地表形變場。根據地表形變模型、干涉圖組合方法和選點方法等不同，眾多學者提出了不同的SBAS 時序InSAR處理策略。本研究采用了Hopper等[21-22]提出的StaMPS SBAS方法進行砂墩子煤礦地表形變測量。StaMPS結合振幅離差指數法和相位穩定性分析法選取濾波后相位微失相干 (slowly decorrelating filtered phase,SDFP)像素，能夠在容易失相干地區保障高密度的高相干目標點，以獲取連續可靠的形變測量結果[23]。同時，由于StaMPS SBAS方法沒有預先設定地表形變的時序變化特征(如假設地表形變隨時間發生一階線性變化)，因而更適宜于礦區等可能在時間序列上具有非線性高階變化特征的地表形變場提取。同時，StaMPS采用三維相位解纏方法恢復纏繞干涉圖的絕對相位，相比于經典二維解纏方法能夠有效提升相位解纏的精度和效率[24]。

在SAR數據處理過程中，首先根據SAR影像垂直基線和時間基線進行干涉對選取，其中最大垂直基線和時間基線分別是150 m和109 m，共組成100個干涉對。圖2展示了所有干涉圖基線組合情況。干涉圖生成過程中，設置相干系數閾值為0.4，對低相干區域像素進行掩模。設置幅度差分離差閾值0.4，利用StaMPS相位穩定性分析方法在整個研究區域獲取了SDFP點6 663個，平均監測密度為101點/km2，表明基于高分辨率Sentinel-1 SAR影像利用StaMPS SBAS方法能夠獲取高密度的監測點。

圖2 SBAS InSAR時空基線分布圖 Fig.2 The spatial-temporal baseline for SBAS InSAR analysis

2.2 基于矩形位錯模型的開采參數反演

(a)俯視圖 (b)正視圖

本研究采用貝葉斯后驗反演算法進行砂墩子煤礦綜采面參數反演。貝葉斯反演算法是基于統計理論的隨機反演算法，關鍵步驟是在反演過程中進行大量采樣，以獲得模型后驗概率分布，這里常用到馬爾科夫鏈條蒙特卡洛采樣方法(Markov chain Monte Carlo method，MCMC)，該方法可進行有效采樣并刻畫后驗概率分布[29]。在反演中主要以貝葉斯框架為理論基礎，結合MCMC方法采樣，通過獲取模型參數后驗概率密度函數，進而得到模型參數最優解。在貝葉斯算法中，后驗概率密度函數p(m|d)是模型參數m對形變觀測量數據d的條件概率函數，公式為:

(1)

式中，p(m)為模型參數的先驗概率密度函數。其中，先驗信息可來自資料、經驗等，常用假設先驗概率密度分布有均勻分布、高斯分布等。p(d)是與m無關的歸一化常量，p(d|m) 是給定形變數據d情況下模型參數m的似然函數，它能夠反映觀測數據與反演模型之間的匹配程度。似然函數可定量描述模型參數的不確定性，一般假設多為高斯分布，可表示為：

(2)

式中：N為點數；∑d為方差-協方差陣；G(m)為模型參數m的非線性模型函數。獲取先驗信息和似然函數后，利用貝葉斯公式獲取后驗概率密度函數p(m|d)，假設模型參數為n個，第i個參數的邊緣概率密度為：

(3)

本研究反演綜采面關鍵參數的具體流程如圖4所示,其中輸入的InSAR觀測形變數據d可為地表形變速率或累積形變量。在反演計算時為了減少不必要的數據冗余，一般首先需要對輸入數據進行降采樣處理，主要方法有均勻重采樣和四叉樹重采樣，同時構建觀測區InSAR形變方差-協方差矩陣∑d，用于計算似然函數；其次，根據先驗信息獲取先驗概率密度函數p(m)，并利用確定的初始模型參數mi(其中mi∈p(m))正演獲取模擬形變 (這里正演是指利用模型參數m和非線性模型函數G(m)計算得到地表形變結果)，然后基于以上結果利用式(2)計算似然函數；最后，根據 MCMC采樣方法和Metropolis-Hastings法則，自動選擇步長來獲取新的模型參數，從而計算新似然函數，并根據 MCMC和Metropolis-Hastings法則對兩個似然函數進行比較且保留對應的模型參數，以此來對模型參數進行大量有效的抽樣，最終獲得模型參數后驗概率密度函數，并取最大后驗概率解為參數最優值，即輸出模型參數的最優估值。

圖4 基于矩形位錯模型的貝葉斯參數反演流程圖Fig.4 Flow chart of the Bayesian parameter inversion based on rectangular dislocation model

3 實驗結果

3.1 砂墩子煤田地表形變場

由于SAR衛星的側視成像方式，利用SBAS InSAR測量得到的是沿衛星視線方向(line of sight，LOS)地表形變?？紤]煤田開采沉陷主要表現在垂直向發生位移，為更直觀地體現礦區沉陷特征，根據雷達波入射角將LOS向形變轉換至垂直向。圖5(a)展示了利用SBAS InSAR測量的砂墩子礦在2018年9月7日—2019年10月8日期間的年平均沉陷速率圖，圖中灰色虛線為過沉陷漏斗的剖面線，P1，P2，P3，P4為沉陷區選取的特征點，黑色虛線為礦區邊界。圖5(a)顯示煤田主井西北部存在一個顯著的沉陷漏斗，面積約為1.5 km2，其中紅色最深的區域表示沉陷中心，此處最大年均沉陷速率約為150 mm/a。通過提取穿越沉陷漏斗的兩條剖面線A1A2和B1B2上年均沉陷速率(圖5(b))，發現接近沉陷中心的兩條剖面線上的沉陷速率變化趨勢一致，均表現為從邊緣至中心急速持續性下沉趨勢，兩剖面線處最大年均沉陷速率約為140 mm/a。

(a)基于SBAS InSAR的砂墩子礦年均沉陷速率 (b)沿A1A2和B1B2剖面線地表沉陷速度變化

為直觀地表現砂墩子礦地表形變的時序變化特征，圖6顯示了在沉陷中心區域選取的特征點P1，P2，P3，P4在觀測時間段內的時序累計沉陷量變化，圖中紅色線條為形變序列的擬合曲線?？梢钥闯?，4個特征點在2018年9月—2019年6月期間均呈現持續性下沉趨勢；而在之后時間段里，4個特征點累計沉陷量逐漸趨于穩定。4個特征點時序形變特征較為一致，交叉驗證了SBAS InSAR恢復砂墩子煤田地表形變結果的可靠性。根據《砂墩子煤礦安全生產標準化煤礦申報表》顯示，砂墩子礦分多個采區并分階段開采，2019年6月后砂墩子礦累計沉陷量趨于穩定，可能是該沉陷區對應綜采面逐漸減小開采量或者停止開采的結果。

圖6 SBAS InSAR獲取的4個特征點時序累積沉降量Fig.6 Time series accumulative subsidence of the 4 points from SBAS InSAR analysis

此外，4個特征點在2019年6月—2019年8月發生相對抬升。根據水文地質調查，砂墩子礦地質條件富含較強的含水帶或導水帶，利于地表降雨和融雪水滲入地下。當地氣象資料及《潞新煤化工有限公司礦井改擴建工程報告》顯示，礦區降雨多集中于夏季(夏季月均降雨量約30 mm)，融雪水通常也在6至7月匯流而下，易形成短時洪流。另外，砂墩子主采煤層頂底板以粉砂巖-泥巖為主，遇水易膨脹。因此，SBAS InSAR觀測到的抬升現象可能是當地的地質條件與氣候環境因素綜合導致。值得注意的是，在礦區當地發生暴雨和山洪時期，井礦可能面臨突發性充水而更易發生事故，應采取積極有效預防措施，以免造成不必要的人員和財產損失。

3.2 Okada位錯模型反演結果

煤礦區地下開釆活動引發的地表形變幅度和形態取決于多方面因素，如開采深度、工作面尺寸等。本文基于Okada矩形位錯模型，以SBAS InSAR測量的2018—2019年地表年均形變速率為約束，反演了沉陷漏斗區對應綜采面的關鍵參數。注意由于輸入反演模型的觀測量是地表年均形變速率，故最終反演估計得到的開采量也代表年均值。反演所得參數結果會影響后續開采量的計算，因此驗證反演所得綜采面參數的準確性十分重要。本文通過以下3個方面來說明反演所得綜采面參數的準確性：①對比反演得到的綜采面參數與搜集的礦區真實資料；② 分析InSAR觀測形變與模擬形變之間的殘差大小和分布特征；③對比由模型計算得到的開采量和已有資料公布的礦井產能。

采用后驗貝葉斯反演方法獲取的砂墩子礦綜采面最優參數如表1所示。該綜采面中心投影至地表的平面坐標相對于坐標原點(E 92.44°，N 43.20°，即圖7左上角頂點)為 968.64 m 東，1 739.45 m 南，開采深度約349.89 m；綜采面走向與北方向夾角即走向角約177.41°，這與該煤田自北向南開采的實際情況相符；綜采面走向長(沿著煤層走向布置的工作面順槽長度)約1 001.27 m，傾向寬約211.80 m；綜采面與水平面的夾角即綜采面傾角約6.01°，根據礦區已有資料顯示，砂墩子礦主采4號煤層，煤層傾角5°～25°，平均傾角8°，該礦區采用綜采采煤工藝，綜采面大致沿煤層層面布置，也就是綜采面傾角基本反映該區域煤層傾角，即本研究反演得到的傾角基本符合實際情況。此外，本文將反演的綜采面參數結果與二維形變圖疊加展示(圖7)，可以看出反演的綜采面中心投影至地表后基本位于地表形變場中心位置，綜采面長度與地表形變場的空間尺度吻合，進一步表明了反演所得最優綜采面參數的可靠性。

表1 砂墩子礦綜采面最優擬合參數Tab.1 Optimal fitting parameters of fully mechanized mining face in Shadunzi Mine

圖7 反演參數與InSAR形變結果的疊加示意圖Fig.7 InSAR surface deformation map with the superposition of the inverted optimal parameters

圖8中(a)—(c)分別展示了SBAS InSAR觀測形變、模擬形變，以及由二者做差得到的殘差分布。從圖8(a)，(b)形變空間分布來看，InSAR觀測形變與模擬形變的分布特征基本一致；煤礦開采主要形變區域(圖8白色矩形框)對應的殘差值基本小于20 mm，較大殘差值主要分布于沉陷漏斗的北緣(圖8(c))；整個觀測區的殘差均值為-2.4 mm，說明殘差值相對InSAR觀測形變較小。殘差分布特征和大小表明基于Okada位錯模型能夠完好地模擬砂墩子煤礦地表形變場，反演得到的煤礦綜采面參數較為可靠。

(a)InSAR觀測形變 (b)模擬形變 (c)InSAR觀測形變與模擬形變差異

開拓煤量即開采量指在礦井可采儲量范圍內已完成設計規定的主井等開拓掘進工程構成的煤儲量，并減去開拓區內水文地質損失量、設計損失量和在可采期內不能回采的臨時煤柱及其他開采量。由于損失煤量通常占比較小可以忽略，本文根據煤田開采設計規范[30]將開拓煤量計算公式簡化為：

Q=L×W×M×D，

(4)

式中：L為煤層兩翼已開拓的走向長度；W為采區平均傾斜長；M為采區煤層平均厚度；D為煤的視密度，指單個煤粒的質量與外觀體積(包括煤的空隙)之比。

經反演獲得砂墩子礦綜采面走向長約1 001.27 m，傾向寬約211.80 m，根據《潞新煤化工有限公司礦井改擴建工程報告》顯示，砂墩子礦主采的4號煤層平均厚度約10.00 m，最大厚度32.48 m，煤層視密度為1.50 t/m3，將上述參數代入式(4)，可大致估算砂墩子礦年均開采量約3.18 Mt。根據新疆維吾爾自治區生態環境保護產業協會發布的《砂墩子礦環境影響評價》公告，砂墩子煤礦自2017年起實際生產規模達已到3.00 Mt/a，表明本文所估算的開采量與該礦實際產能基本相符。

4 結論

本文利用2018年9月7日—2019年10月8日期間覆蓋砂墩子礦區的30景Sentinel-1A升軌影像，基于時序SBAS InSAR方法對該礦區進行了地表沉陷監測。SBAS InSAR監測結果顯示砂墩子礦區主井西北側存在一個明顯的沉陷漏斗，對應該礦目前投產的一個綜采面，沉陷面積約1.5 km2，最大年均沉陷速率約150 mm/a。綜采面上方地表累積形變時序曲線顯示在2018年9月—2019年6月，沉陷漏斗區地表表現持續性下沉趨勢；2019年6月后累積沉陷量趨于穩定，可能是開采作業放緩或者停止的結果。

基于Okada矩形位錯模型，本文采用貝葉斯后驗反演方法利用SBAS InSAR觀測形變，反演得到砂墩子礦區綜采面最優參數解，其中采深約349.89 m，綜采面走向與北方向夾角約177.41°，走向長約1 001.27 m，傾向寬約211.80 m，傾角約6.01°，與已有礦區資料基本符合。將由最優參數解算的模擬形變與SBAS InSAR觀測形變結果對比，發現它們空間分布特征一致，兩者殘差較小。根據反演所得綜采面走向長與傾向寬參數，估算得到砂墩子礦在InSAR觀測期間的年均開采量約3.18 Mt，與已有礦區資料報道的礦區年產能3.00 Mt/a一致。本文研究結果表明了基于時序InSAR觀測形變反演和估算煤田綜采面關鍵參數的可行性。

志謝：感謝歐空局提供Sentinel-1A SAR影像及美國地質調查局(UGSG)提供的SRTM-1地形數據。