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

沙永蓮，王曉文，劉國(guó)祥，張 瑞，張 波

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

0 引言

煤炭是我國(guó)能源結(jié)構(gòu)的主體，對(duì)煤炭資源開采活動(dòng)進(jìn)行合理的規(guī)劃和管理對(duì)我國(guó)經(jīng)濟(jì)建設(shè)和生態(tài)環(huán)境保障具有重要意義[1]。我國(guó)煤炭資源多是通過地下開采，而這種開采模式往往會(huì)破壞地層內(nèi)部原有的力學(xué)平衡狀態(tài)，導(dǎo)致采空區(qū)上方地表發(fā)生移動(dòng)和變形[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]。

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

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

本文以新疆哈密砂墩子煤田為例研究了基于時(shí)序SBAS InSAR觀測(cè)形變反演地下煤礦綜采面參數(shù)及煤礦開采量的可行性。砂墩子煤田是國(guó)家重點(diǎn)煤礦，設(shè)計(jì)產(chǎn)能5.00 Mt/a，自2012年開始試生產(chǎn)以來，至今仍在開采中；該煤田處于半干旱地區(qū)，礦區(qū)表面及周邊植被覆蓋較少，非常有利于采用時(shí)序InSAR準(zhǔn)確恢復(fù)其地表形變特征。本研究選取2018年9月—2019年10月共30景Sentinel-1A升軌影像，首先利用SBAS方法恢復(fù)了該礦區(qū)一年時(shí)間內(nèi)的地表形變；其次，以SBAS InSAR測(cè)量的地表形變速率為約束，基于Okada矩形位錯(cuò)模型反演了該煤礦綜采面參數(shù)，進(jìn)而根據(jù)煤層視密度、綜采面參數(shù)和開采量的關(guān)系估算了砂墩子礦開采量，并與已有礦區(qū)資料進(jìn)行了對(duì)比驗(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國(guó)道5 km，南距蘭新線6 km (圖1)。圖上，白色虛線為礦區(qū)邊界線，紅框表示InSAR分析區(qū)域，光學(xué)影像為研究區(qū)Sentinel-2光學(xué)影像，白色三角為主井位置。井田東西走向長(zhǎng)約7 km，南北最寬約8 km，礦井分布三個(gè)井筒，主井位置為E92°28′05″，N43°10′04″。從煤炭資源賦存來看，該礦區(qū)資源儲(chǔ)量877.44 Mt，設(shè)計(jì)產(chǎn)能5.00 Mt/a，一期設(shè)計(jì)產(chǎn)能3.00 Mt/a[18]。該礦區(qū)水文地質(zhì)條件簡(jiǎn)單，屬簡(jiǎn)單型礦床，主采煤層頂?shù)装逡苑凵皫r-泥巖為主[19]。整個(gè)礦區(qū)氣候干旱，晝夜溫差大，屬于典型溫帶大陸性干旱氣候；降雨量少但集中，雨期多在夏季，易形成短時(shí)洪流。

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

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

2 算法原理

2.1 SBAS InSAR 時(shí)序地表形變提取

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

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

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

2.2 基于矩形位錯(cuò)模型的開采參數(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對(duì)形變觀測(cè)量數(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ù)，它能夠反映觀測(cè)數(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個(gè)，第i個(gè)參數(shù)的邊緣概率密度為：

(3)

本研究反演綜采面關(guān)鍵參數(shù)的具體流程如圖4所示,其中輸入的InSAR觀測(cè)形變數(shù)據(jù)d可為地表形變速率或累積形變量。在反演計(jì)算時(shí)為了減少不必要的數(shù)據(jù)冗余，一般首先需要對(duì)輸入數(shù)據(jù)進(jìn)行降采樣處理，主要方法有均勻重采樣和四叉樹重采樣，同時(shí)構(gòu)建觀測(cè)區(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法則，自動(dòng)選擇步長(zhǎng)來獲取新的模型參數(shù)，從而計(jì)算新似然函數(shù)，并根據(jù) MCMC和Metropolis-Hastings法則對(duì)兩個(gè)似然函數(shù)進(jìn)行比較且保留對(duì)應(yīng)的模型參數(shù)，以此來對(duì)模型參數(shù)進(jìn)行大量有效的抽樣，最終獲得模型參數(shù)后驗(yàn)概率密度函數(shù)，并取最大后驗(yàn)概率解為參數(shù)最優(yōu)值，即輸出模型參數(shù)的最優(yōu)估值。

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

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

3.1 砂墩子煤田地表形變場(chǎng)

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

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

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

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

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

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

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

采用后驗(yàn)貝葉斯反演方法獲取的砂墩子礦綜采面最優(yōu)參數(shù)如表1所示。該綜采面中心投影至地表的平面坐標(biāo)相對(duì)于坐標(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í)際情況相符；綜采面走向長(zhǎng)(沿著煤層走向布置的工作面順槽長(zhǎng)度)約1 001.27 m，傾向?qū)捈s211.80 m；綜采面與水平面的夾角即綜采面傾角約6.01°，根據(jù)礦區(qū)已有資料顯示，砂墩子礦主采4號(hào)煤層，煤層傾角5°～25°，平均傾角8°，該礦區(qū)采用綜采采煤工藝，綜采面大致沿煤層層面布置，也就是綜采面傾角基本反映該區(qū)域煤層傾角，即本研究反演得到的傾角基本符合實(shí)際情況。此外，本文將反演的綜采面參數(shù)結(jié)果與二維形變圖疊加展示(圖7)，可以看出反演的綜采面中心投影至地表后基本位于地表形變場(chǎng)中心位置，綜采面長(zhǎng)度與地表形變場(chǎng)的空間尺度吻合，進(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觀測(cè)形變、模擬形變，以及由二者做差得到的殘差分布。從圖8(a)，(b)形變空間分布來看，InSAR觀測(cè)形變與模擬形變的分布特征基本一致；煤礦開采主要形變區(qū)域(圖8白色矩形框)對(duì)應(yīng)的殘差值基本小于20 mm，較大殘差值主要分布于沉陷漏斗的北緣(圖8(c))；整個(gè)觀測(cè)區(qū)的殘差均值為-2.4 mm，說明殘差值相對(duì)InSAR觀測(cè)形變較小。殘差分布特征和大小表明基于Okada位錯(cuò)模型能夠完好地模擬砂墩子煤礦地表形變場(chǎng)，反演得到的煤礦綜采面參數(shù)較為可靠。

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

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

Q=L×W×M×D，

(4)

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

經(jīng)反演獲得砂墩子礦綜采面走向長(zhǎng)約1 001.27 m，傾向?qū)捈s211.80 m，根據(jù)《潞新煤化工有限公司礦井改擴(kuò)建工程報(bào)告》顯示，砂墩子礦主采的4號(hào)煤層平均厚度約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é)會(huì)發(fā)布的《砂墩子礦環(huán)境影響評(píng)價(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升軌影像，基于時(shí)序SBAS InSAR方法對(duì)該礦區(qū)進(jìn)行了地表沉陷監(jiān)測(cè)。SBAS InSAR監(jiān)測(cè)結(jié)果顯示砂墩子礦區(qū)主井西北側(cè)存在一個(gè)明顯的沉陷漏斗，對(duì)應(yīng)該礦目前投產(chǎn)的一個(gè)綜采面，沉陷面積約1.5 km2，最大年均沉陷速率約150 mm/a。綜采面上方地表累積形變時(shí)序曲線顯示在2018年9月—2019年6月，沉陷漏斗區(qū)地表表現(xiàn)持續(xù)性下沉趨勢(shì)；2019年6月后累積沉陷量趨于穩(wěn)定，可能是開采作業(yè)放緩或者停止的結(jié)果。

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

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