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

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

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

0 引言

煤炭是我國能源結(jié)構(gòu)的主體，對煤炭資源開采活動進(jìn)行合理的規(guī)劃和管理對我國經(jīng)濟(jì)建設(shè)和生態(tài)環(huán)境保障具有重要意義[1]。我國煤炭資源多是通過地下開采，而這種開采模式往往會破壞地層內(nèi)部原有的力學(xué)平衡狀態(tài)，導(dǎo)致采空區(qū)上方地表發(fā)生移動和變形[2]。煤田礦區(qū)地表不均勻變形可能導(dǎo)致地面形成裂縫或者塌陷，破壞原有地質(zhì)環(huán)境條件并引發(fā)地質(zhì)災(zāi)害，從而嚴(yán)重影響地表相關(guān)建筑物等基礎(chǔ)設(shè)施的安全使用，給礦區(qū)當(dāng)?shù)亟?jīng)濟(jì)可持續(xù)發(fā)展造成隱患[3]。

針對礦區(qū)不均勻地表形變的時空演化特征進(jìn)行監(jiān)測，可為當(dāng)?shù)氐乇沓料蓦[患治理，安全生產(chǎn)保障和當(dāng)?shù)厣鷳B(tài)環(huán)境保護(hù)等提供關(guān)鍵支撐信息。目前，水準(zhǔn)測量、全球?qū)Ш叫l(wèi)星系統(tǒng)(GNSS)和合成孔徑雷達(dá)干涉(InSAR)等多種大地測量方法已在礦區(qū)地表沉陷監(jiān)測中得到廣泛應(yīng)用[4-5]。相比于傳統(tǒng)大地測量手段，InSAR具有全天候觀測、覆蓋范圍廣和觀測效率高等優(yōu)勢[6-7]。常規(guī)基于兩景SAR影像進(jìn)行干涉處理測量地表形變的兩軌法DInSAR (Differential InSAR)，由于無法對形變觀測模型中的誤差項(xiàng)(如地形誤差、軌道誤差和大氣延遲誤差等)進(jìn)行校正，常被用于瞬時、大幅度地表形變(如地震和火山)測量[8-9]。近十年來，眾多學(xué)者基于多景SAR影像提出多種時序InSAR分析方法以克服常規(guī)DInSAR相位失相干、大氣延遲等影響，并實(shí)現(xiàn)長時間序列地表形變監(jiān)測[10-11]。其中，小基線集(small baseline subset，SBAS)時序InSAR方法能夠有效利用冗余干涉對進(jìn)行干涉相位建模，可有效抑制形變恢復(fù)模型中的誤差項(xiàng)，已在城市地表沉降、活動斷層滑移和滑坡運(yùn)動監(jiān)測等領(lǐng)域中得到廣泛應(yīng)用。

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

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

1 研究區(qū)概況及數(shù)據(jù)源

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

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

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

2 算法原理

2.1 SBAS InSAR 時序地表形變提取

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

3 實(shí)驗(yàn)結(jié)果

3.1 砂墩子煤田地表形變場

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

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

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

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

此外，4個特征點(diǎn)在2019年6月—2019年8月發(fā)生相對抬升。根據(jù)水文地質(zhì)調(diào)查，砂墩子礦地質(zhì)條件富含較強(qiáng)的含水帶或?qū)畮В诘乇斫涤旰腿谘┧疂B入地下。當(dāng)?shù)貧庀筚Y料及《潞新煤化工有限公司礦井改擴(kuò)建工程報(bào)告》顯示，礦區(qū)降雨多集中于夏季(夏季月均降雨量約30 mm)，融雪水通常也在6至7月匯流而下，易形成短時洪流。另外，砂墩子主采煤層頂?shù)装逡苑凵皫r-泥巖為主，遇水易膨脹。因此，SBAS InSAR觀測到的抬升現(xiàn)象可能是當(dāng)?shù)氐牡刭|(zhì)條件與氣候環(huán)境因素綜合導(dǎo)致。值得注意的是，在礦區(qū)當(dāng)?shù)匕l(fā)生暴雨和山洪時期，井礦可能面臨突發(fā)性充水而更易發(fā)生事故，應(yīng)采取積極有效預(yù)防措施，以免造成不必要的人員和財(cái)產(chǎn)損失。

3.2 Okada位錯模型反演結(jié)果

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

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

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

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

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

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

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

Q=L×W×M×D，

(4)

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

經(jīng)反演獲得砂墩子礦綜采面走向長約1 001.27 m，傾向?qū)捈s211.80 m，根據(jù)《潞新煤化工有限公司礦井改擴(kuò)建工程報(bào)告》顯示，砂墩子礦主采的4號煤層平均厚度約10.00 m，最大厚度32.48 m，煤層視密度為1.50 t/m3，將上述參數(shù)代入式(4)，可大致估算砂墩子礦年均開采量約3.18 Mt。根據(jù)新疆維吾爾自治區(qū)生態(tài)環(huán)境保護(hù)產(chǎn)業(yè)協(xié)會發(fā)布的《砂墩子礦環(huán)境影響評價(jià)》公告，砂墩子煤礦自2017年起實(shí)際生產(chǎn)規(guī)模達(dá)已到3.00 Mt/a，表明本文所估算的開采量與該礦實(shí)際產(chǎn)能基本相符。

4 結(jié)論

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

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

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