杜建濤，閆 麗，趙超英

(1.長安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054；2.中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710077)

我國是一個(gè)煤炭能源大國，已探明的煤儲(chǔ)量位居世界前列，而礦產(chǎn)資源的過度開采容易破壞礦區(qū)地質(zhì)結(jié)構(gòu)，進(jìn)而引發(fā)地表開裂、塌陷甚至滑坡等地質(zhì)災(zāi)害[1]。因此，礦區(qū)沉陷監(jiān)測(cè)是煤礦安全開采和可持續(xù)發(fā)展的重要工作，對(duì)預(yù)防潛在地質(zhì)災(zāi)害具有重要意義[2]。合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Rader，InSAR)因其高精度、大覆蓋范圍及可重復(fù)觀測(cè)等技術(shù)優(yōu)勢(shì)在采礦沉陷監(jiān)測(cè)中得到廣泛應(yīng)用，極大地克服了傳統(tǒng)測(cè)量手段空間分辨率低、周期長、成本高等局限。Zhao Chaoying等[3]、Yang Chengsheng 等[4]對(duì)山西大同礦區(qū)大量地表塌陷進(jìn)行InSAR 監(jiān)測(cè)；Zheng Meinan 等[5]、Du Sen等[6]基于多種先進(jìn)InSAR 技術(shù)分析了河北峰峰礦區(qū)的沉陷機(jī)理；Yang Zefa 等[7]、Ma Chao 等[8]、賀衛(wèi)中等[9]針對(duì)陜西榆林礦區(qū)地面塌陷采用InSAR 技術(shù)進(jìn)行大范圍、多維形變監(jiān)測(cè)。然而，目前針對(duì)采礦沉陷的InSAR 時(shí)空演化分析和精細(xì)化監(jiān)測(cè)研究不足，且該技術(shù)在河北蔚縣礦區(qū)沉陷監(jiān)測(cè)的應(yīng)用也較少。為了揭示該礦區(qū)的地表形變演化特征，本文利用Sentinel-1A/B 數(shù)據(jù)對(duì)采礦沉陷開展InSAR 監(jiān)測(cè)，獲取了礦區(qū)升降軌數(shù)據(jù)沿視線向(LOS)形變監(jiān)測(cè)結(jié)果，并在此基礎(chǔ)上重點(diǎn)分析了西細(xì)莊礦近1 a 內(nèi)地表二維形變的時(shí)空演化特征。

1 研究區(qū)概況

蔚縣位于河北省張家口市最南邊，境內(nèi)南部為高山地區(qū)；中部地勢(shì)平坦，壺流河蜿蜒而過；北部為丘陵區(qū)，溝壑密布，地下蘊(yùn)藏著豐富的煤炭資源。蔚縣礦區(qū)總面積264 km2，距離蔚縣縣城約10 km，是河北省最大的地下采煤區(qū)之一。礦區(qū)井田分布如圖1 所示，崔家寨井田、單侯井田、南留莊井田和西細(xì)莊井田為目前主要開采區(qū)域；鄭溝灣礦、興源礦及周邊部分礦井已開采殆盡。受小煤窯破壞，崔家寨井田西翼四采區(qū)目前處于停采狀態(tài)。礦區(qū)含煤地層為侏羅系下花園組，分為上、下兩段，上段為雜色含煤地層，可靠性較差，均不可采；下段含煤8 層，位于底部的1 號(hào)煤層(西細(xì)莊礦主采層，平均厚度3.16 m)和中部的5 號(hào)、6 號(hào)煤層(崔家寨礦主采層，平均厚度分別為0.9 m 和2.04 m)，層位穩(wěn)定，厚度和分布面積均較大，為礦區(qū)主采煤層，埋藏深度為143.6～336.5 m。煤層頂板巖性為泥巖、粉砂巖、細(xì)砂巖，底板巖性為粉砂巖、細(xì)砂巖及泥巖。地表第四系松散層厚140～245 m，基巖厚度30～55 m，松散層下段砂礫石含水層厚7.5～34.3 m，巖性為粗砂、礫石及礫卵石[10]。

圖1 河北蔚縣地理位置及SAR 影像覆蓋范圍Fig.1 Location of Yuxian County and SAR image coverage

2 數(shù)據(jù)采集與處理

本文采用61景Sentinel-1A/B干涉寬幅(Interferometric Wide swath，IW)模式數(shù)據(jù)進(jìn)行礦區(qū)形變監(jiān)測(cè)，包括33 景降軌數(shù)據(jù)(軌道號(hào)T120，時(shí)間為2017 年7 月至2018 年9 月)和28 景升軌數(shù)據(jù)(軌道號(hào)T142，時(shí)間為2017 年8 月至2018 年8 月)，影像具體參數(shù)如表1 所示。此外，本文還利用精密軌道數(shù)據(jù)AUX_POEORB文件來精化影像初始軌道狀態(tài)矢量，利用30 m 分辨率的SRTM DEM 數(shù)據(jù)來模擬和去除地形相位。

表1 Sentinel-1A/B SAR 數(shù)據(jù)參數(shù)Table 1 SAR data and their imaging parameters

對(duì)獲取的61 景升降軌Sentinel-1A/B 數(shù)據(jù)分別進(jìn)行小基線集(Small Baseline Subsets，SBAS)干涉處理[11]，然后采用融合升降軌數(shù)據(jù)的多維形變時(shí)序估計(jì)方法，生成垂向和東西向的二維形變時(shí)間序列[12]。關(guān)鍵處理步驟如下：

a.Sentinel-1 預(yù)處理 由于Sentinel-1 數(shù)據(jù)的特殊成像模式，首先需對(duì)原始影像進(jìn)行預(yù)處理，將其轉(zhuǎn)化成可以干涉的單視復(fù)數(shù)(SLC)格式，具體包括采用AUX_POEORB 軌道文件進(jìn)行基線精化，對(duì)主、從影像進(jìn)行去斜和配準(zhǔn)，使方位向配準(zhǔn)精度達(dá)到千分之一像元。

b.SLC 影像干涉 設(shè)置時(shí)間和空間基線閾值生成小基線差分干涉圖。為提高干涉圖的相干性，并抑制噪聲，對(duì)差分干涉圖進(jìn)行自適應(yīng)頻率域?yàn)V波處理，采用基于Delaunay 三角網(wǎng)的最小費(fèi)用流算法進(jìn)行相位解纏。為抑制時(shí)空去相干的影像，本文時(shí)間基線閾值取36 d，空間基線閾值取200 m。干涉組合的時(shí)空基線分布如圖2 所示。

圖2 干涉組合時(shí)空基線分布Fig.2 Baseline of InSAR combination

c.生成形變圖 挑選高質(zhì)量干涉組合，估計(jì)并去除趨勢(shì)項(xiàng)和大氣延遲相位，構(gòu)建相干像元的形變相位解算模型，采用最小二乘法進(jìn)行參數(shù)求解，得到SAR 影像獲取時(shí)刻間的平均形變速率。

d.計(jì)算二維形變時(shí)間序列 融合多軌道SAR 數(shù)據(jù)的多維形變時(shí)序估計(jì)技術(shù)，利用奇異值分解法來計(jì)算多源SAR 數(shù)據(jù)在重疊時(shí)間段內(nèi)的二維形變速率和形變時(shí)間序列。為了避免觀測(cè)方程的病態(tài)，需要對(duì)設(shè)計(jì)矩陣進(jìn)行正則化處理[13]。觀測(cè)方程表示為式(1)。

式中A=(-cosαsinθΔt,cosθΔt)為設(shè)計(jì)矩陣；α為衛(wèi)星飛行方向的方位角；θ為雷達(dá)波入射角；Δt為相鄰SAR影像時(shí)間間隔；ξ為正則化參數(shù)；I為單位矩陣；V=(VE,VU)T，VE和VU為待求參數(shù)，分別表示相鄰SAR 影像獲取時(shí)刻的東西向和垂向形變速率；D為形變相位矩陣。待求參數(shù)估計(jì)值可表示為：

最后，對(duì)相鄰SAR 獲取時(shí)刻間的形變速率在時(shí)間域積分，即可得到地表東西向和垂向的形變時(shí)間序列：

式中d為累積形變時(shí)間序列；n為不同軌道影像時(shí)間跨度重合的總數(shù)量。

3 結(jié)果分析

3.1 礦區(qū)沉陷監(jiān)測(cè)數(shù)據(jù)

將獲取的升降軌SAR 數(shù)據(jù)按上述流程處理，生成蔚縣礦區(qū)LOS 向年平均形變速率圖，重采樣到地理坐標(biāo)系下，如圖3 所示，其中，圖3a 為降軌形變監(jiān)測(cè)結(jié)果，圖3b 為升軌形變監(jiān)測(cè)結(jié)果。

圖3 河北蔚縣礦區(qū)2017—2018 年地表形變場(chǎng)Fig.3 Surface deformation site of Yuxian mining area of Hebei Province from 2017 to 2018

從圖3 可以看出，升、降軌數(shù)據(jù)獲取的礦區(qū)形變場(chǎng)基本一致，整個(gè)礦區(qū)在監(jiān)測(cè)期間均出現(xiàn)不同程度的地表形變現(xiàn)象，且沉陷區(qū)普遍呈近橢圓形的漏斗狀分布。礦區(qū)崔家寨井田、單侯井田和西細(xì)莊井田內(nèi)采煤造成的地面沉陷較為嚴(yán)重，南留莊井田在監(jiān)測(cè)期間形變相對(duì)不明顯。其中，崔家寨井田沉陷面積最大，且沉降速率達(dá)-30 cm/a。資料顯示，礦區(qū)在2015 年調(diào)整井田邊界和優(yōu)化資源配置，提高了各礦井生產(chǎn)能力，其中，崔家寨煤礦生產(chǎn)能力180 萬t/a，多個(gè)立井開拓；單侯煤礦生產(chǎn)能力150 萬t/a，3 個(gè)立井開拓；西細(xì)莊礦生產(chǎn)能力60 萬t/a，一對(duì)立井開拓；南留莊礦生產(chǎn)能力30 萬t/a，一對(duì)立井開拓[14]。地下煤炭資源的長期高強(qiáng)度開采破壞開采區(qū)周圍原始應(yīng)力平衡，從而導(dǎo)致覆巖移動(dòng)并向上波及地表，使地表發(fā)生形變和沉陷。

為了定量分析礦區(qū)地面沉陷災(zāi)害情況，本文對(duì)不同量級(jí)年沉陷速率區(qū)域的面積進(jìn)行統(tǒng)計(jì)分析。圖4 為礦區(qū)沉陷等值線圖，等值線由內(nèi)向外依次為-20、-15、-10、-5 cm/a，從圖4 可以看出，崔家寨、麥子破、回回墓等多個(gè)村莊受到采礦沉陷影響。表2 為不同量級(jí)年沉陷速率區(qū)域面積統(tǒng)計(jì)結(jié)果，其中年沉陷速率超過-10 cm/a 的區(qū)域達(dá)到2.16 km2，年沉陷速率超過-20 cm/a 的區(qū)域達(dá)到0.10 km2。大面積的地表沉陷對(duì)礦區(qū)的基礎(chǔ)建設(shè)及地表生態(tài)環(huán)境造成了嚴(yán)重的影響[15]。

圖4 河北蔚縣礦區(qū)沉陷等值線(單位：cm/a)Fig.4 Subsidence contours in Yuxian mining area in Hebei Province

表2 河北蔚縣礦區(qū)沉陷面積統(tǒng)計(jì)Table 2 Subsidence area statistics in Yuxian mining area of Hebei Province

對(duì)比圖3a、圖3b 可以發(fā)現(xiàn)，升、降軌數(shù)據(jù)獲取的監(jiān)測(cè)結(jié)果中部分區(qū)域在形變范圍和量級(jí)上存在一些差異。為此，本文提取跨越西細(xì)莊礦塌陷區(qū)的剖線AA′，對(duì)其形變空間特征進(jìn)行分析(圖5)。從圖5 可以看出，升、降軌數(shù)據(jù)獲取的西細(xì)莊礦最大年形變速率相差約5 cm/a，且沉陷中心位置也具有一定偏差。究其原因，可能是由于升、降軌數(shù)據(jù)具有不同的軌道方向和入射角度，其對(duì)地表真實(shí)形變的敏感度不同。另外，采礦沉陷的LOS 向形變中可能包含水平形變的影響。

3.2 二維形變時(shí)序

圖5 蔚縣西細(xì)莊煤礦年形變速率剖線Fig.5 Annual deformation rate along profile A—A′ of Xixizhuang mine in Yuxian County

由于西細(xì)莊井田的干涉圖整體相干性較高，便于時(shí)序分析，因此，在升、降軌一維形變分析基礎(chǔ)上，仍以該礦為試驗(yàn)區(qū)，采用融合多軌道SAR 數(shù)據(jù)的多維形變估算法，研究礦區(qū)二維形變時(shí)空演化特征。將統(tǒng)一到相同地理坐標(biāo)系下的升降軌形變圖，按式(1)構(gòu)建方程組，求解不同時(shí)間段內(nèi)的二維形變速率，然后依據(jù)式(3)，獲取西細(xì)莊礦2017 年8月至2018 年8 月期間垂向與東西向的形變時(shí)間序列，如圖6 所示。分析圖6 可以發(fā)現(xiàn)，西細(xì)莊礦地表在監(jiān)測(cè)期間以垂向形變?yōu)橹?，也具有明顯的東西向水平形變，且東、西水平形變方向相反，量級(jí)也存在一定差異。垂向形變?cè)诳臻g上呈規(guī)則的橢圓形，東西向水平形變則表現(xiàn)為兩個(gè)半圓形，二者形變范圍均隨時(shí)間逐漸擴(kuò)大，形變量級(jí)也相應(yīng)增加。

為定量分析礦區(qū)地表形變的時(shí)空演化特征，共選取4 個(gè)特征點(diǎn)P1、P2、P3、P4(圖6)，進(jìn)行二維形變時(shí)間序列分析(圖7)，其中P1、P2 點(diǎn)分別位于沉陷中心區(qū)域及采礦工作面傾向邊緣，P3、P4 點(diǎn)位于采礦工作面走向邊緣。

由圖7 可以看出，在2017 年8 月至2018 年8月期間，P1 點(diǎn)的垂向累積形變最大，達(dá)–30 cm，且具有持續(xù)下沉的趨勢(shì)；P3、P4 點(diǎn)位于東、西向水平形變最大的區(qū)域，形變量分別達(dá)到–10 cm 和7 cm，其形變方向相反，其中正值表示向東移動(dòng)，負(fù)值表示向西移動(dòng)，且形變都指向沉陷漏斗中心，同時(shí)二者均具有較大的垂向形變，形變量分別達(dá)–22 cm 和–13 cm?？梢钥闯?，P3 點(diǎn)的垂向和東西向形變量都較點(diǎn)P4 大，這與工作面的開采方向有關(guān)，P2 點(diǎn)的垂向形變達(dá)–20 cm，但和P1 點(diǎn)都具有水平方向保持不變的特點(diǎn)。此外，不同于城市地面沉降的線性特征，采礦沉降的形變過程與時(shí)間呈現(xiàn)出非線性關(guān)系，從圖中P1 點(diǎn)的垂向形變時(shí)間序列曲線(圖7a)可以看出，該點(diǎn)處的地表形變具有先緩慢下沉后加速沉降最后趨于穩(wěn)定的特征，這一形變演化過程對(duì)應(yīng)于礦區(qū)地表單點(diǎn)沉降的初始沉降、主要沉降和殘余沉降3 個(gè)時(shí)期，與“S”型增長的Logistic 時(shí)間函數(shù)模型能夠很好地吻合[16-17]，更能反映采礦沉陷的真實(shí)規(guī)律，因此，利用InSAR 數(shù)據(jù)進(jìn)行多維形變估算獲得的地表沉陷值，與實(shí)際情況較為吻合，可為沉陷的分析、治理提供參考數(shù)據(jù)。

圖6 蔚縣西細(xì)莊礦二維形變時(shí)間序列Fig.6 Time sequence of 2D deformation of Xixizhuang mine in Yuxian County

圖7 蔚縣西細(xì)莊礦特征點(diǎn)二維形變時(shí)間序列Fig.7 Time sequence of 2D deformation of feature points of Xixizhuang mine in Yuxian County

4 結(jié)論

a.采用InSAR 技術(shù)獲取了河北蔚縣采礦區(qū)大范圍地面沉陷的空間分布特征，且部分區(qū)域年形變速率達(dá)到-20 cm/a。該結(jié)果顯示出InSAR 技術(shù)用于礦區(qū)地表形變監(jiān)測(cè)的優(yōu)越性，為礦區(qū)安全開采和生態(tài)恢復(fù)提供科學(xué)依據(jù)。

b.通過對(duì)研究區(qū)西細(xì)莊井田地表形變的精細(xì)化監(jiān)測(cè)與分析，獲得采礦沉陷的時(shí)空演化特征，即礦區(qū)地表動(dòng)態(tài)沉降過程符合Logistic 時(shí)間函數(shù)模型，該特征揭示了真實(shí)的采礦沉陷規(guī)律。

c.基于時(shí)間函數(shù)模型參數(shù)的采礦沉陷預(yù)測(cè)和礦區(qū)形變?nèi)S分解將是下一步研究計(jì)劃。