楊 利，薛東劍，王海方，付 林，張 婷

(1.成都理工大學(xué)地球科學(xué)學(xué)院，四川 成都 610059；2.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756)

礦產(chǎn)資源的開(kāi)發(fā)利用是我國(guó)經(jīng)濟(jì)的重要組成部分。然而礦產(chǎn)資源開(kāi)發(fā)引發(fā)的環(huán)境地質(zhì)問(wèn)題日益突出，特別是礦山開(kāi)采易造成崩塌、滑坡、礦震、地面沉陷等地質(zhì)災(zāi)害，嚴(yán)重威脅到礦區(qū)周邊人民的生產(chǎn)生命安全以及國(guó)家的財(cái)產(chǎn)安全[1]。對(duì)礦區(qū)地表進(jìn)行動(dòng)態(tài)監(jiān)測(cè)可以有效獲取地表移動(dòng)變化規(guī)律，最大限度地降低地表沉陷帶來(lái)的損失[2]。因此對(duì)礦區(qū)地表沉降進(jìn)行有效監(jiān)測(cè)，分析沉降特征，對(duì)開(kāi)采沉陷預(yù)計(jì)與治理具有重要意義[3]。

合成孔徑雷達(dá)(synthetic aperture radar，SAR)具有全天時(shí)、全天候、多波段、多極化、高分辨率、不易受氣象條件制約、對(duì)地表具有一定的穿透力等優(yōu)勢(shì)[4]。 干涉合成孔徑雷達(dá)(interferometric synthetic aperture radar，InSAR)是利用目標(biāo)與兩天線的位置幾何關(guān)系，通過(guò)計(jì)算兩次過(guò)境時(shí)SAR影像的相位差來(lái)獲得數(shù)字高程模型的技術(shù)[5]。GABRIEL等[6]提出的差分干涉合成孔徑雷達(dá)(differential InSAR，D-InSAR)技術(shù)則是引入外部DEM或三軌差分、四軌差分實(shí)現(xiàn)對(duì)地表厘米級(jí)的微小變形的監(jiān)測(cè)，高精度、高時(shí)效、低成本、大范圍是該技術(shù)的監(jiān)測(cè)優(yōu)勢(shì)。然而時(shí)空去相干及大氣延遲等因素使得干涉條紋無(wú)法生成或不連續(xù)，極大制約著該技術(shù)監(jiān)測(cè)的精度和運(yùn)用，較難完成高精度、長(zhǎng)時(shí)間間隔地表監(jiān)測(cè)[7]。為克服這些因素的影響，BERARDINO等[8]提出了小基線集方法(small baseline subset，SBAS)。 SBAS-InSAR技術(shù)通過(guò)影像的自由組合得到較多干涉，選擇空間和時(shí)間基線較小、具有高相干性的干涉圖，以減小地形對(duì)差分的影響，提高變形監(jiān)測(cè)的準(zhǔn)確性，能夠更好地消除D-InSAR技術(shù)中的時(shí)空失相關(guān)、大氣延遲相位以及地形相位誤差等影響，獲取研究區(qū)域毫米級(jí)的形變信息[9-10]。

本文針對(duì)礦區(qū)的地表沉降問(wèn)題，結(jié)合地形考慮Sentinel-1A影像的可視情況，采用SBAS-InSAR技術(shù)獲取研究區(qū)的地表形變信息。利用地學(xué)統(tǒng)計(jì)方法，分析沉降區(qū)在空間上的分布，結(jié)合形變速率和累計(jì)形變量，對(duì)工作區(qū)范圍內(nèi)的典型沉降區(qū)進(jìn)行時(shí)序分析。

1 研究區(qū)概況及數(shù)據(jù)準(zhǔn)備

1.1 研究區(qū)概況

研究區(qū)處于螢石溝-紅旗溝脆韌性剪切帶及其所形成的斷裂構(gòu)造集中帶之中東段，是五龍溝礦區(qū)的一個(gè)重要構(gòu)造發(fā)育區(qū)，隸屬五龍溝礦區(qū)三大主要控成礦構(gòu)造區(qū)帶之一。研究區(qū)內(nèi)出露地層主要以下元古界金水口群、上元古界青白口系丘吉東溝組和下古生界奧陶系祁曼塔格群變火山巖組為主，中元古界長(zhǎng)城紀(jì)小廟組次之，溝谷和山前有大面積的第四系沖洪積物分布。 研究區(qū)中心坐標(biāo)：95°56′4.14″E，36°12′24.2″N，面積約5 144 377 m2，總體構(gòu)造線呈北西-南東向展布。礦石采礦均采用地下平硐開(kāi)拓方式，隨著地下礦石開(kāi)采量的增加，研究區(qū)水閘東溝-黃龍溝方向(北西-南東向)形成了大范圍的地下采空區(qū)，導(dǎo)致研究區(qū)及其附近存在崩塌、泥石流、地面塌陷等地質(zhì)災(zāi)害。

圖1 研究區(qū)地質(zhì)構(gòu)造簡(jiǎn)圖

1.2 數(shù)據(jù)源

本文采用SRTMDEM 90M分辨率原始高程數(shù)據(jù)作為參考DEM為SBAS-InSAR提供地形信息和地理位置參考，利用從歐洲航天局質(zhì)量控制中心下載的精密軌道數(shù)據(jù)對(duì)軌道信息進(jìn)行修正，去除因軌道誤差引起的系統(tǒng)性誤差，以保證結(jié)果的準(zhǔn)確性。由于觀測(cè)幾何的限制，不同方向形變的貢獻(xiàn)值不相同。結(jié)合地形考慮影像的可視性，采用的Sentinel-1A衛(wèi)星以IW模式獲取的30期降軌單視復(fù)數(shù)產(chǎn)品進(jìn)行時(shí)序反演，該數(shù)據(jù)輻射精度為1 dB，相位誤差為5°，地面分辨率為5 m×20 m(距離向×方位向)，時(shí)間跨度為2017年11月18日—2019年9月21日。

2 礦區(qū)形變監(jiān)測(cè)

SBAS技術(shù)是將輸入影像自由組合，以多主影像的干涉對(duì)為基礎(chǔ)，基于高相干點(diǎn)恢復(fù)研究區(qū)域的時(shí)間序列形變信息[12]。原理是根據(jù)每個(gè)圖像中基線的時(shí)間和位置，將研究區(qū)中的SAR圖像組合成小的基線子集，進(jìn)行差分干涉處理，得到較多的干涉像對(duì)，保證了較高的空間密度和時(shí)間采樣[13]。然后利用奇異值分解算法聯(lián)合多個(gè)小基線子集求解地表形變相位最小范數(shù)意義下的最小二乘解[14]。該算法有效解決了總體法方程秩虧和單個(gè)集合時(shí)間采樣率較低的問(wèn)題[15]。本文利用SBAS-InSAR技術(shù)進(jìn)行數(shù)據(jù)處理的流程如圖2所示。

圖2 SBAS-InSAR技術(shù)流程圖

基于青海省的地形、氣候和研究區(qū)的植被茂密程度分析，將時(shí)間基線閾值設(shè)為550 d，以避免出現(xiàn)無(wú)相干數(shù)據(jù)對(duì)；為避免完全空間失相關(guān)，將空間基線閾值設(shè)為45%，并進(jìn)行了3D解纏，以便估算相干性低的區(qū)域、建立新的連接。利用參考DEM模擬并去除干涉圖中的地形相位，采用自適應(yīng)濾波方法平滑相干斑噪聲，去除大氣延遲相位，提高干涉條紋的清晰度。選擇Delaunay MCF算法對(duì)干涉圖進(jìn)行相位解纏，并剔除解纏結(jié)果不理想的像對(duì)。選取相位穩(wěn)定、未發(fā)生形變的68個(gè)地面控制點(diǎn)對(duì)干涉結(jié)果進(jìn)行軌道精煉和重去平，消除了斜坡效應(yīng)，對(duì)相位偏移進(jìn)行修正，估算和去除殘余的恒定相位。去除地形殘余相位、大氣延遲相位、噪聲相位是SBAS技術(shù)的核心，常分成兩步進(jìn)行。第一次形變反演是建立線性模型計(jì)算出所有像對(duì)的殘余地形及地表形變速率，利用奇異值分解法解算相關(guān)參數(shù)，得出干涉像對(duì)平均形變速率和高程誤差；第二次形變反演是通過(guò)大氣濾波估算和去除大氣相位和地形殘余相位，獲取更加準(zhǔn)確的時(shí)間序列最終位移圖。經(jīng)過(guò)SBAS兩次反演得到的地表形變信息是在SAR坐標(biāo)系中，將反演結(jié)果通過(guò)地理編碼轉(zhuǎn)化為地理坐標(biāo)，最終獲得研究區(qū)的地表形變圖。

3 實(shí)驗(yàn)結(jié)果分析與討論

在雷達(dá)圖像干涉處理的基礎(chǔ)上，利用GIS工具，結(jié)合研究區(qū)SBAS反演的年平均沉降速率、累計(jì)形變等信息，進(jìn)而分析工作區(qū)的沉降分布、沉降量及典型區(qū)域沉降中心的形變時(shí)間序列，并對(duì)沉降區(qū)的剖面進(jìn)行時(shí)序形變分析。

3.1 礦區(qū)整體形變分析

圖3 研究區(qū)形變速率圖

在監(jiān)測(cè)時(shí)間范圍內(nèi)，獲得的五龍溝礦區(qū)年平均速率如圖3所示，該范圍內(nèi)地表移動(dòng)速率相對(duì)來(lái)說(shuō)比較穩(wěn)定，但沉降區(qū)表現(xiàn)明顯。 工作區(qū)的平均沉降速率為-1.40 mm/a，最大沉降速率為-11.74 mm/a。如圖4(a)所示，-2～0 mm/a為主要的形變速率區(qū)間，占工作區(qū)面積的72.35%，而形變速率為-11.8～-10 mm/a的強(qiáng)烈形變區(qū)域僅占工作區(qū)面積的0.54%。

五龍溝礦區(qū)以2017年11月18日為基準(zhǔn)，在2018-05-05期開(kāi)始出現(xiàn)較明顯的沉降區(qū)域，隨著時(shí)間的推移，采礦活動(dòng)的繼續(xù)開(kāi)展及礦區(qū)自身重力的影響，沉降區(qū)域沉降量增加，沉降范圍逐漸擴(kuò)大，到2019年9月21日，五龍溝礦區(qū)范圍內(nèi)最大沉降量可達(dá)-21.68 mm。如圖4(b)所示，工作區(qū)范圍內(nèi)的整體型變量不大，主要集中在-3～0 mm，占該工作區(qū)面積的61.57%，其中，沉降量超過(guò)-6 mm的面積約為364 344 m2，占工作區(qū)面積的7.09%。

圖4 形變面積占比圖

3.2 典型區(qū)形變分析

目前工作區(qū)有原礦堆存場(chǎng)、選礦廠、尾礦庫(kù)、地下采礦區(qū)(水閘東溝采區(qū)和黃龍溝采區(qū))、礦區(qū)道路、硐口工業(yè)場(chǎng)地等設(shè)施。如圖5所示，工作區(qū)范圍內(nèi)主要存在兩個(gè)明顯沉降區(qū)，且分布在工作區(qū)的西部，這與實(shí)際情況相符。在此基礎(chǔ)上圈定兩個(gè)典型沉降區(qū)域，分別為沉降區(qū)A和沉降區(qū)B。

圖5 典型沉降區(qū)分布圖

沉降區(qū)A的總面積約為124 575 m2，占工作區(qū)面積的2.42%，該沉降區(qū)形變較小，主要集中在-7～-10 mm之間，其所占面積是沉降區(qū)A面積的53.67%；沉降區(qū)B的總面積為155 311 m2，占工作區(qū)面積的3.02%，形變量明顯大于沉降區(qū)A，其形變量在-22～-19 mm范圍內(nèi)的面積大約是沉降區(qū)A的7倍。工作區(qū)內(nèi)兩沉降區(qū)的形變量面積占比情況見(jiàn)表1。

表1 典型沉降區(qū)形變量面積占比

以2017年11月18日為基準(zhǔn)，在監(jiān)測(cè)時(shí)間范圍內(nèi)，地表監(jiān)測(cè)點(diǎn)的沉降沒(méi)有出現(xiàn)突變現(xiàn)象，整體沉降位移量最大不超過(guò)30 mm。沉降區(qū)A的最大形變點(diǎn)位于沉降區(qū)中部偏左下方，如圖6(a)所示，其最大沉降值為-20.27 mm。該點(diǎn)以近似穩(wěn)定的速率呈現(xiàn)線性沉降的趨勢(shì)。沉降區(qū)B的最大形變點(diǎn)位于沉降區(qū)的左上方，如圖6(b)所示，其最大沉降值達(dá)到-21.68 mm。該點(diǎn)的沉降呈現(xiàn)非線性特征變化。沉降區(qū)A主要采用分段空?qǐng)龇?、無(wú)底柱分段崩落采礦法回采。在2017-11-18—2019-04-06時(shí)間段內(nèi)，該區(qū)域空區(qū)已形成，隨著礦區(qū)地下開(kāi)采的擴(kuò)展，采空區(qū)的暴露面積、圍巖所承受的應(yīng)力以及圍巖蠕變量相應(yīng)變大，誘發(fā)頂板出現(xiàn)冒落，造成地表微小形變。2019-04-06期后較平穩(wěn)的原因是該區(qū)通過(guò)自然冒落或強(qiáng)制崩落圍巖方式進(jìn)行處理。沉降區(qū)B采礦方式與沉降區(qū)A相同。 在監(jiān)測(cè)段內(nèi)正開(kāi)采3 505 m中斷，平衡拱承受的應(yīng)力沒(méi)有達(dá)到破壞極限，頂板未出現(xiàn)冒落，因此地表呈現(xiàn)的沉降變化不明顯。隨著工作的推進(jìn)，平衡拱受空區(qū)傾向和走向跨度擴(kuò)展的影響，拱的受壓逐漸接近極限，超過(guò)極限后平衡拱被破壞，且由于礦區(qū)自身重力的影響，冒落高度增大，頂板巖層移動(dòng)變形造成地表微小沉降。相較于沉降區(qū)A，沉降區(qū)B在局部地段未支護(hù)，其穩(wěn)固性較差，地表沉降趨勢(shì)較為明顯，其中2019-08-04處形變起伏最為顯著。

地下空區(qū)形成后，采礦區(qū)會(huì)存在對(duì)應(yīng)的沉陷中心。豎直方向上的變形將導(dǎo)致地面路基、土質(zhì)的松弛，影響其承受能力，繼而引起和加大水平方向的傾斜和拉伸變形。沉降區(qū)A以最大形變點(diǎn)位為中心逐漸向四周擴(kuò)散，其中向東北向擴(kuò)散的趨勢(shì)較明顯；沉降區(qū)B以最大形變點(diǎn)位為中心呈帶狀沿西北-東南方向逐漸向四周擴(kuò)散。在沉降區(qū)擴(kuò)散趨勢(shì)范圍內(nèi)分別繪制沉降區(qū)A和沉降區(qū)B的剖面線，通過(guò)時(shí)序分析發(fā)現(xiàn)，兩沉降區(qū)中心沉降(圖7)均呈現(xiàn)逐年增強(qiáng)趨勢(shì)。沉降區(qū)A呈漏斗狀逐漸向四周擴(kuò)散，且向四周擴(kuò)散趨勢(shì)明顯。沉降區(qū)B則是呈現(xiàn)盆狀逐漸向外圍擴(kuò)散，擴(kuò)散趨勢(shì)較為平穩(wěn)。

圖6 沉降區(qū)最大形變點(diǎn)形變序列圖

4 結(jié) 論

1) 利用SBAS-InSAR方法得到了青海省五龍溝礦區(qū)2017年11月至2019年9月期間各時(shí)間段的地表微小緩慢形變，獲取了礦區(qū)地面沉降的時(shí)間序列圖及其空間展布。監(jiān)測(cè)結(jié)果顯示工作區(qū)整體形變不大，但存在兩個(gè)明顯沉降區(qū)。

2) 形變干涉結(jié)果表明五龍溝礦區(qū)總體比較穩(wěn)定。在監(jiān)測(cè)時(shí)間范圍內(nèi)，工作區(qū)地表整體平均沉降速率為-1.40 mm/a，其中形變速率為-2～2 mm/a的區(qū)域占實(shí)驗(yàn)區(qū)面積的81.54%；地表形變量為-3～3 mm的區(qū)域占工作區(qū)面積的79.02%。工作區(qū)內(nèi)最大形變速率為-11.8 mm/a，最大沉降量為-21.68 mm，其中沉降量超過(guò)-6 mm的面積為364 344.35 m2，僅占工作區(qū)面積的7.09%。

3) 沉降區(qū)A呈現(xiàn)線性特征，以最大沉降點(diǎn)為中心向四周擴(kuò)散，且向東北向擴(kuò)散的趨勢(shì)較明顯。沉降區(qū)B呈現(xiàn)非線性特征，以最大沉降點(diǎn)位為中心，沉降區(qū)呈帶狀沿西北-東南方向逐漸向四周擴(kuò)散。