杜建濤,閆 麗,趙超英
(1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院,陜西 西安 710054;2.中煤科工集團(tuán)西安研究院有限公司,陜西 西安 710077)
我國(guó)是一個(gè)煤炭能源大國(guó),已探明的煤儲(chǔ)量位居世界前列,而礦產(chǎn)資源的過(guò)度開采容易破壞礦區(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è)量手段空間分辨率低、周期長(zhǎng)、成本高等局限。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í)空演化特征。
蔚縣位于河北省張家口市最南邊,境內(nèi)南部為高山地區(qū);中部地勢(shì)平坦,壺流河蜿蜒而過(guò);北部為丘陵區(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
本文采用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文件來(lái)精化影像初始軌道狀態(tài)矢量,利用30 m 分辨率的SRTM DEM 數(shù)據(jù)來(lái)模擬和去除地形相位。
表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ù),利用奇異值分解法來(lái)計(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ù)量。
將獲取的升降軌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 萬(wàn)t/a,多個(gè)立井開拓;單侯煤礦生產(chǎn)能力150 萬(wàn)t/a,3 個(gè)立井開拓;西細(xì)莊礦生產(chǎn)能力60 萬(wàn)t/a,一對(duì)立井開拓;南留莊礦生產(chǎn)能力30 萬(wàn)t/a,一對(duì)立井開拓[14]。地下煤炭資源的長(zhǎng)期高強(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é)果,其中年沉陷速率超過(guò)-10 cm/a 的區(qū)域達(dá)到2.16 km2,年沉陷速率超過(guò)-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 向形變中可能包含水平形變的影響。
圖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)。此外,不同于城市地面沉降的線性特征,采礦沉降的形變過(guò)程與時(shí)間呈現(xiàn)出非線性關(guān)系,從圖中P1 點(diǎn)的垂向形變時(shí)間序列曲線(圖7a)可以看出,該點(diǎn)處的地表形變具有先緩慢下沉后加速沉降最后趨于穩(wěn)定的特征,這一形變演化過(guò)程對(duì)應(yīng)于礦區(qū)地表單點(diǎn)沉降的初始沉降、主要沉降和殘余沉降3 個(gè)時(shí)期,與“S”型增長(zhǎng)的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
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.通過(guò)對(duì)研究區(qū)西細(xì)莊井田地表形變的精細(xì)化監(jiān)測(cè)與分析,獲得采礦沉陷的時(shí)空演化特征,即礦區(qū)地表動(dòng)態(tài)沉降過(guò)程符合Logistic 時(shí)間函數(shù)模型,該特征揭示了真實(shí)的采礦沉陷規(guī)律。
c.基于時(shí)間函數(shù)模型參數(shù)的采礦沉陷預(yù)測(cè)和礦區(qū)形變?nèi)S分解將是下一步研究計(jì)劃。