基于SBAS-InSAR技術(shù)的礦區(qū)地表形變監(jiān)測(cè)研究

梁 芳，楊維芳*，李蓉蓉

（1.蘭州交通大學(xué) 測(cè)繪與地理信息學(xué)院，甘肅 蘭州 730070；2.地理國(guó)情監(jiān)測(cè)技術(shù)應(yīng)用國(guó)家地方聯(lián)合工程研究中心，甘肅 蘭州 730070；3.甘肅省地理國(guó)情監(jiān)測(cè)工程實(shí)驗(yàn)室，甘肅 蘭州 730070）

煤炭是人類(lèi)使用的主要能源，對(duì)國(guó)家經(jīng)濟(jì)發(fā)展具有重要推動(dòng)作用；然而進(jìn)行大面積地下煤炭開(kāi)采活動(dòng)易引起地面形變，破壞地表結(jié)構(gòu)，誘發(fā)滑坡、塌陷、地裂縫等地質(zhì)災(zāi)害，嚴(yán)重威脅區(qū)域的可持續(xù)發(fā)展，給人們的生命財(cái)產(chǎn)安全帶來(lái)隱患[1-2]，因此及時(shí)有效地對(duì)礦區(qū)地面進(jìn)行形變監(jiān)測(cè)具有重要的研究意義。傳統(tǒng)水準(zhǔn)測(cè)量、GPS等測(cè)量技術(shù)存在測(cè)站稀疏、工作量大、周期長(zhǎng)、監(jiān)測(cè)范圍小、點(diǎn)位保存困難等缺點(diǎn)[3]，不能適應(yīng)礦區(qū)形變監(jiān)測(cè)的要求。

近年來(lái)，衛(wèi)星遙感領(lǐng)域不斷發(fā)展，隨著大量星載SAR數(shù)據(jù)的出現(xiàn)，合成孔徑雷達(dá)應(yīng)用技術(shù)也在不斷完善，差分雷達(dá)干涉測(cè)量技術(shù)（D-InSAR）具有連續(xù)覆蓋、高精度和高度自動(dòng)化的能力，成為獲取地表形變信息的有效手段[4]。然而，D-InSAR技術(shù)存在時(shí)空基線和大氣效應(yīng)等誤差，影響了其測(cè)量精度，很難獲得連續(xù)的形變結(jié)果[5-7]。SBAS-InSAR技術(shù)通過(guò)影像的自由組合得到較多干涉圖，選擇空間和時(shí)間基線較小，具有高相干性的干涉圖可減小地形對(duì)差分的影響，提高變形監(jiān)測(cè)的準(zhǔn)確性，能更好地消除技術(shù)中的時(shí)空失相關(guān)、大氣延遲相位以及地形相位誤差等影響，獲取研究區(qū)域毫米級(jí)的形變信息[8-9]。因此，本文采用SBAS-InSAR技術(shù)對(duì)29景Sentinel-1A影像數(shù)據(jù)進(jìn)行處理，獲得了研究區(qū)的形變信息，并從形變特征和形變機(jī)理等方面對(duì)形變信息進(jìn)行了監(jiān)測(cè)研究。

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

本文選取的研究區(qū)位于陜西省榆林市榆陽(yáng)區(qū)東北部27 km的毛烏素沙漠南沿、屬北部草灘區(qū)東部，與神木縣大保當(dāng)鎮(zhèn)接壤，西南與榆陽(yáng)區(qū)牛家梁鎮(zhèn)接壤，北部與榆陽(yáng)區(qū)孟家灣鄉(xiāng)相連，神延鐵路、榆府公路、110和330 kV線路均穿鎮(zhèn)而過(guò)；以沙漠灘地為主，地形較平緩，屬低緩的半沙漠緩坡丘陵地貌，總體地形西部稍低、東部稍高；地下水、煤炭資源儲(chǔ)存量大，地質(zhì)資源量為1 872.32 Mt，可采儲(chǔ)量為995.78 Mt，井田地質(zhì)構(gòu)造簡(jiǎn)單，煤層賦存穩(wěn)定，煤質(zhì)優(yōu)良，可作為多種用途資源；礦井外部建設(shè)條件良好，交通四通八達(dá)。

1.2 數(shù)據(jù)來(lái)源

本文下載覆蓋陜西省榆林市榆陽(yáng)區(qū)2019-01-01—2019-12-15的29景歐洲航天局Sentinel-1A升軌影像數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，數(shù)據(jù)類(lèi)型為干涉寬幅模式（IW）的單視復(fù)數(shù)（SLC）數(shù)據(jù)，地面分辨率為5 m×20 m，極化方式為VV同向極化。實(shí)驗(yàn)還利用30 m分辨率的SRTM1高程數(shù)據(jù)消除地形相位，利用各影像對(duì)應(yīng)的精密定軌星歷數(shù)據(jù)減少軌道誤差。

2 數(shù)據(jù)處理流程與結(jié)果

2.1 SBAS-InSAR數(shù)據(jù)處理流程

SBAS又稱(chēng)短基線集，是由Berardino P等在2002年提出的一種基于分布式目標(biāo)算法的時(shí)序InSAR技術(shù)[10-11]。其基本思想是將所有覆蓋同一地區(qū)的SAR影像進(jìn)行組合，形成若干個(gè)集合，再通過(guò)有效合并得到所有可用的小基線干涉圖。利用奇異值分解（SVD）獲取影像序列間地表形變速率的最小范數(shù)最小二乘解，能有效克服時(shí)間和空間失相干的影響，提高時(shí)間采樣頻率和空間覆蓋率，且能在一定程度上消除大氣相位和限制地形誤差的影響[12-14]。SBAS-InSAR方法監(jiān)測(cè)地表形變的技術(shù)流程如圖1所示。

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

實(shí)驗(yàn)基于ENVI軟件的SARscape，利用SBAS-In-SAR技術(shù)對(duì)礦區(qū)進(jìn)行形變監(jiān)測(cè)：①設(shè)置最大時(shí)間基線閾值為180 d，空間基線閾值為臨界基線的45%，生成影像連接圖；②利用Delaunay三角網(wǎng)中的最小費(fèi)用流法和Goldstein濾波對(duì)連接好的像對(duì)進(jìn)行干涉處理，并需移除相干性低和解纏效果不好的像對(duì)；③在相位解纏效果好的圖中選擇50個(gè)高相干的控制點(diǎn)來(lái)估算和去除殘余的恒定相位和解纏后仍存在的相位坡道，并基于所選的控制點(diǎn)進(jìn)行軌道精煉和重去平；④第一次估算形變速率和殘余地形；⑤計(jì)算時(shí)間序列上的位移，即對(duì)第一次得到的形變速率進(jìn)行定制的大氣濾波，用以估算和去除大氣相位，得到時(shí)間序列上的最終形變結(jié)果，并將生成的所有LOS方向的形變量都投影到垂直方向，通過(guò)ArcGIS軟件輸出。

2.2 SBAS-InSAR數(shù)據(jù)處理結(jié)果

經(jīng)實(shí)驗(yàn)處理得到的研究區(qū)年均形變速率、累積形變量如圖2、3所示，圖中紅色區(qū)域表示下沉，綠色區(qū)域表示抬升。該區(qū)域年均形變速率范圍為-180.03～42.03 mm/a，累積形變量范圍為-167.20～63.31 mm，形變量較大。

圖2 SBAS年均形變速率

3 研究結(jié)果分析

3.1 形變特征分析

圖3 SBAS時(shí)序累積形變量

由時(shí)序累積形變量可知，隨著時(shí)間的推移，礦區(qū)地表形變量逐漸增大，形變范圍也逐漸擴(kuò)大；研究區(qū)包括A、B、C、D、E五個(gè)明顯的大型形變區(qū)域，且呈由邊緣向中心形變逐漸增大的漏斗狀，這與礦區(qū)開(kāi)采工作面的分布位置有關(guān)。

為更加直觀地對(duì)礦區(qū)的形變特征進(jìn)行分析，以研究區(qū)的B工作面為例，沿該工作面走向和傾向分別作剖線L1、L2，然后在兩條剖線上分別選取一系列點(diǎn)，提取每個(gè)點(diǎn)的形變值，最后繪制每個(gè)剖線的動(dòng)態(tài)累積曲線進(jìn)行時(shí)序分析（圖4）。

由圖4a可知，隨著時(shí)間的推移，形變區(qū)整體呈下沉趨勢(shì)，每個(gè)點(diǎn)的沉降值隨時(shí)間的推移不斷增加，大的沉降值從開(kāi)始的1.7 mm增加到94.4 mm，形成了比較明顯的下沉盆地，且有一個(gè)快速下沉區(qū)域，而其他區(qū)域點(diǎn)均以緩慢速度下沉；小的沉降值從開(kāi)始的1.4 mm增加到2.8 mm。圖4b同樣表現(xiàn)出下沉盆地特征，每期累積沉降值大致表現(xiàn)為先增大后減小的過(guò)程；下沉盆地中心沉降值從開(kāi)始的1.1 mm增加到75.5 mm；整體符合礦區(qū)開(kāi)采沉陷規(guī)律。

本文在走向線和傾向線上隨機(jī)選取4個(gè)點(diǎn)，利用擬合模型擬合得到下沉值與時(shí)間的關(guān)系，結(jié)果如圖5所示，可以看出，下沉值與時(shí)間大致呈線性關(guān)系，即隨著時(shí)間的推移，沉降量不斷增大。

圖5 時(shí)序沉降特征

3.2 降水對(duì)形變的影響分析

礦區(qū)發(fā)生形變的主要原因是礦區(qū)的開(kāi)采，但其他原因也會(huì)在一定程度上影響其形變結(jié)果。研究區(qū)以風(fēng)沙草灘、黃土丘陵為主，土壤含沙量多，土質(zhì)疏松，抗腐蝕能力差，易遭受流水侵蝕[15]，因此探究降水對(duì)礦區(qū)形變的影響很有必要。本文引入降水量參數(shù)，利用§3.1中用于擬合的4個(gè)點(diǎn)提取各點(diǎn)每月的累積形變量，并通過(guò)與在中國(guó)氣象數(shù)據(jù)網(wǎng)下載的研究區(qū)同時(shí)段降水量信息進(jìn)行對(duì)比，分析礦區(qū)形變量與降水量的關(guān)系，結(jié)果如圖6所示，可以看出，降水量與形變量并無(wú)明顯的幾何特征關(guān)系，2019年1月—2019年3月該地區(qū)幾乎沒(méi)有降水，但4個(gè)點(diǎn)不斷發(fā)生形變，且形變量不斷增大；月累積降水量在6—8月逐漸增加，且在8月達(dá)到最大值，而1、3、4號(hào)點(diǎn)的形變量均呈先增后減的趨勢(shì)，2號(hào)點(diǎn)的形變量呈先減后增的趨勢(shì)；10—12月的月累積降水量逐漸減小，4個(gè)點(diǎn)的形變量均呈先增后減的趨勢(shì)，因此降水量與形變量沒(méi)有顯著關(guān)系。這是由于研究區(qū)內(nèi)地層土體滲透系數(shù)小，降水對(duì)礦區(qū)形變量影響很小，即相對(duì)于礦區(qū)大量級(jí)形變而言，其降水的影響可忽略不計(jì)。

圖6 降水量與形變量的對(duì)比性分析

3.3 精度驗(yàn)證

本文選取研究區(qū)B工作面沿走向線和傾向線上的12個(gè)水準(zhǔn)點(diǎn)來(lái)驗(yàn)證SBAS-InSAR獲取的累積形變量。由于SAR影像獲取時(shí)間與水準(zhǔn)數(shù)據(jù)觀測(cè)時(shí)間不一致，因此首先采用三次樣條插值法獲取與SAR影像獲取時(shí)間一致的水準(zhǔn)數(shù)據(jù)，然后采用最鄰近距離法提取水準(zhǔn)點(diǎn)附近的SBAS點(diǎn)，最后對(duì)二者的差值進(jìn)行對(duì)比驗(yàn)證，如圖7所示。經(jīng)過(guò)對(duì)比發(fā)現(xiàn)，SBAS-InSAR與水準(zhǔn)數(shù)據(jù)在趨勢(shì)上基本是一致的；計(jì)算得出二者的中誤差為4.75 mm，相關(guān)系數(shù)為0.981，相關(guān)系數(shù)較高，最大差值為9.42 mm，最小差值為1.39 mm，滿足監(jiān)測(cè)精度要求。由此可以得出，SBAS-InSAR技術(shù)在礦區(qū)形變監(jiān)測(cè)中具有可靠性。

圖7 水準(zhǔn)數(shù)據(jù)和SBAS-InSAR結(jié)果的對(duì)比驗(yàn)證

4 結(jié)語(yǔ)

本文采用SBAS方法，利用29景Sentinel-1A影像對(duì)陜西省某礦區(qū)進(jìn)行了形變監(jiān)測(cè)，并從形變特征和降水對(duì)形變的影響等方面進(jìn)行分析，得到的結(jié)論為：

1）研究區(qū)存在5個(gè)較大的形變區(qū)域，通過(guò)研究提取的工作面走向和傾向線上的一系列點(diǎn)發(fā)現(xiàn)，其形變特征呈下沉盆地的漏斗狀。在走向線和傾向線上隨機(jī)選取4個(gè)點(diǎn)，利用線性擬合模型擬合發(fā)現(xiàn)下沉值和時(shí)間呈線性關(guān)系。

2）引入降水量參數(shù)，提取上述4個(gè)點(diǎn)每月的累積形變量，并通過(guò)獲取研究區(qū)同時(shí)段降水量信息，分析了礦區(qū)形變量與降水量的關(guān)系，得出降水量與形變量并無(wú)顯著關(guān)系的結(jié)論。

3）將SBAS-InSAR獲取的累積形變量值與水準(zhǔn)數(shù)據(jù)進(jìn)行對(duì)比發(fā)現(xiàn)，二者的中誤差為4.75 mm，相關(guān)系數(shù)為0.981，表明SBAS-InSAR技術(shù)在礦區(qū)形變監(jiān)測(cè)中具有可靠性。

利用SBAS-InSAR技術(shù)對(duì)礦區(qū)地表形變進(jìn)行監(jiān)測(cè)，一方面可為礦區(qū)的安全開(kāi)采和煤礦塌陷區(qū)的綜合性治理提供科學(xué)依據(jù)；另一方面對(duì)礦區(qū)人民的生命安全具有重大意義。