耿 丹，涂 寬，諶 華，鄭 健

（二十一世紀空間技術(shù)應(yīng)用股份有限公司，北京 100096）

引言

礦產(chǎn)開采沉陷所引發(fā)的地表大面積沉降會給礦區(qū)及周邊環(huán)境帶來嚴重的負面影響。對礦區(qū)地表形變進行長時間序列監(jiān)測，分析其時空演化特征，可為礦區(qū)地面沉降、礦洞坍塌、泥石流、滑坡等地質(zhì)災(zāi)害防治提供有效支撐，減少對礦區(qū)周邊人民生命財產(chǎn)及工程設(shè)施的威脅，提高礦區(qū)土地的使用價值，促進礦區(qū)的可持續(xù)發(fā)展［1-2］。

高分辨率光學(xué)及雷達遙感技術(shù)已廣泛應(yīng)用于地表形變監(jiān)測、地質(zhì)災(zāi)害監(jiān)測、信息提取等方面［3-6］，其中合成孔徑雷達差分干涉測量（Differentialinterfer Interferometric Synthetic Aperture Radar，D?InSAR）技術(shù)能夠快速獲取大面積的整體地表形變信息，且監(jiān)測成本低，能克服各種天氣進行全天時監(jiān)測，在地表形變監(jiān)測方面已開展廣泛應(yīng)用［7-9］，但是該方法受時間、空間失相干和大氣延時相位的影響，難以在長時間序列的地表緩慢變形監(jiān)測中得到理想結(jié)果［10-14］。小基線集差分干涉（Small Baseline Subsets Interferometric Synthetic Ap?erture Radar，SBAS?InSAR）同時具有D?InSAR技術(shù)的優(yōu)點，而且受垂直基線的影響較小，能有效減小空間失相干和大氣相位影響，監(jiān)測長時間大范圍地表緩慢形變，監(jiān)測有效性及精度較高，同時能夠獲得研究區(qū)的形變時間演化特征，在地表形變監(jiān)測方面應(yīng)用廣泛［15-16］。孫曉鵬等［17］利用SBAS?InSAR方法和ENVISAT ASAR數(shù)據(jù)，監(jiān)測成都平原2008年到2010年的地表形變；彭米米等［18］采用SBAS?InSAR方法和哨兵1A數(shù)據(jù)監(jiān)測大西安地區(qū)2015年到2017年的地表形變并對時序特征進行研究；華怡穎等［12］以大寧礦區(qū)為例，利用小基線集技術(shù)對研究區(qū)9景PALSAR數(shù)據(jù)進行時序處理，得到研究區(qū)沉降速率圖、時序累計沉降圖。上述學(xué)者都利用SBAS?InSAR技術(shù)對地表沉降進行監(jiān)測和時序分析，但卻不能有效的對地表沉降的空間演化特征進行分析。

標準差橢圓（Standard Deviational Ellipse，SDE）是分析地理要素空間分布特征的經(jīng)典方法之一，已廣泛應(yīng)用于經(jīng)濟、人口、農(nóng)業(yè)、城市發(fā)展等各行各業(yè)。近年來，已有學(xué)者將標準差橢圓方法應(yīng)用于自然災(zāi)害的空間分布特征研究。如：熊俊楠等［19-20］利用標準差橢圓等方法對重慶市1950?2015年、西藏自治區(qū)1983?2015年山洪災(zāi)害進行時空分布特征分析；周超凡等［21］利用標準差橢圓方法對2003?2010年北京市地面沉降進行空間演化特征研究。

綜上所示，為研究礦區(qū)地表沉降時空演化特征，本文選取阜新某礦區(qū)為試驗區(qū)，采用SBAS?InSAR技術(shù)提取礦區(qū)長時間序列地表形變信息，基于形變剖面提取典型特征點分析礦區(qū)沉降時間演化特征，利用標準差橢圓方法分析礦區(qū)形變空間分布演化特征，結(jié)合高分辨率光學(xué)遙感影像對監(jiān)測結(jié)果進行驗證，為礦區(qū)的地質(zhì)災(zāi)害防治提供參考依據(jù)。

1 試驗區(qū)及數(shù)據(jù)源介紹

1.1 試驗區(qū)介紹

本文選取阜新市某礦區(qū)作為試驗區(qū)（圖1），該試驗區(qū)地勢由東南向西北傾斜，海拔100～621 m，中部地勢較低。礦區(qū)地層由三部分構(gòu)成：上部主要為砂巖、砂礫巖；中部是含有少量砂巖、礫巖和砂質(zhì)頁巖等的煤層；下部基本是細砂巖［22］。根據(jù)已有調(diào)查成果，礦區(qū)有多處矸石山和采坑，矸石山與采坑相互連接，整體沿北東-南西向呈條帶狀展布［23］。

圖1 試驗區(qū)范圍Fig.1 Scope of test area

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

本文采用由歐空局免費提供的Sentinel-1B星Level-1級的單視復(fù)數(shù)圖像（Single Look Complex，SLC）作為數(shù)據(jù)源［24-26］，共獲取2019年10月14日-2020年3月30日期間的15景影像數(shù)據(jù)。同時采用2020年12月獲取的北京二號0.8 m高分辨率光學(xué)遙感數(shù)據(jù)產(chǎn)品進行專家判識。雷達及光學(xué)衛(wèi)星相關(guān)參數(shù)見表1。為提高地面形變結(jié)果精度，選用SRTM 30m分辨率DEM數(shù)據(jù)來減少地形因素對SBAS?InSAR監(jiān)測結(jié)果的影響。

表1 影像相關(guān)參數(shù)Table 1 Relevant parameters of image

2 研究方法

本文以Sentinel-1B SAR影像和北京二號高分辨率光學(xué)遙感影像為數(shù)據(jù)源，基于SBAS?InSAR技術(shù)和標準差橢圓方法開展礦區(qū)地表形變時空演化特征研究。首先提取礦區(qū)形變時間序列和礦區(qū)年均形變速率，并基于GIS的空間分析功能對礦區(qū)地表形變的時間和空間演化特征進行分析，最后利用北京二號光學(xué)影像對監(jiān)測結(jié)果進行驗證。研究技術(shù)路線見圖2。

圖2 礦區(qū)地表形變時空演化特征研究技術(shù)路線Fig.2 Technical route of research on temporal and spatial evolution characteristics of land subsidence in mining area

2.1 SBAS-InSAR技術(shù)

SBAS?InSAR小基線集差分干涉技術(shù)將所有雷達影像數(shù)據(jù)組合成多個集合，集合內(nèi)數(shù)據(jù)之間的時間和空間基線距盡量小，再解算每個小集合的地表形變時間序列［27?38］。主要步驟包括：（1）按照時間順序排列研究區(qū)的雷達單復(fù)數(shù)影像數(shù)據(jù)，設(shè)定時空基線閾值，進行差分干涉；（2）選取一定數(shù)量的高相干點進行軌道精煉和重去平；（3）在穩(wěn)定區(qū)域選擇GCP控制點，對差分干涉影像進行相位解纏；（4）利用SVD方法解算解纏相位，求解形變速率和高程系數(shù)，去除大氣相位和地形殘余相位；（5）地理編碼得到WGS84坐標系下的礦區(qū)年平均形變速率、形變時間序列信息。

2.2 標準差橢圓方法

1926年Lefever提出了標準差橢圓的方法，其基本思想是以橢圓的中心、面積、長軸、短軸、方位角來分析礦區(qū)形變的中心、范圍、形狀、方向等演化過程和趨勢［29?31］。以礦區(qū)地表形變量為權(quán)重，求解標準差橢圓的中心、長軸、短軸和方位角的計算公式如下［32］：

（1）標準差橢圓的中心

式中，（SX,S Y）即為標準差橢圓的中心，（x i，y i）是時序點i的坐標是所有時序點平均中心，n是時序點的數(shù)量。

（2）標準差橢圓的方位角

式中，θ即為標準差橢圓的方位角和是時序點坐標（x i，y i）與所有時序點平均中心的偏差。

（3）標準差橢圓長軸和短軸的標準差

式中，σx和σy分別為沿長軸和短軸的標準差。

3 結(jié)果分析與驗證

3.1 SBAS-InSAR處理結(jié)果

本文利用SBAS?InSAR技術(shù)對15景Sentinel-1B數(shù)據(jù)處理，對影像進行多視、濾波、相位解纏、干涉對挑選等，最終獲取86個干涉對，像對的時空基線見圖3，部分差分干涉像對見圖4。

圖3 干涉像對時空基線連接圖Fig.3 Time space baseline connection diagram of interference image pair

圖4 差分干涉像對Fig.4 Differential interference image pair

礦區(qū)2019年10月-2020年3月的地表形變年平均速率見圖5。整個監(jiān)測期內(nèi)，礦區(qū)整體存在不同程度的形變，地表不均勻形變嚴重，形成了多個形變中心。礦區(qū)整體形變沿北東-南西方向呈條帶狀分布，與開采坑口空間展布方位一致。礦區(qū)內(nèi)形變嚴重區(qū)域形變速率達到50～60 mm/a，主要分布在形變的中心位置。

3.2 礦區(qū)形變時序演化特征分析

礦區(qū)整體存在不同程度的形變，部分區(qū)域形變量較大，造成在形變中心及周邊區(qū)域失相干現(xiàn)象嚴重，在圖5中表現(xiàn)為無值區(qū)，無法分析及確定形變中心位置。為分析礦區(qū)整體形變時序演化，利用Kriging插值方法對2019年10月-2020年3月的累計形變量進行空間插值分析。Kriging插值是一種能夠準確產(chǎn)生預(yù)測表面的地統(tǒng)計方法，在土壤科學(xué)和地質(zhì)學(xué)中尤為適用。

圖5 礦區(qū)年平均速率分布圖Fig.5 Distribution of annual average rate in mining area

基于礦區(qū)地表累計形變分布圖及插值圖，根據(jù)形變中心位置（圖5中無值區(qū)）沿礦區(qū)的北東-南西方向作一條縱向貫穿整個礦區(qū)的剖面線AA’，盡可能多的通過分布在該方向上的形變中心；依據(jù)以上原則，再沿北西-南東方向作四條垂直于剖面線AA’的剖面線BB’、CC’、DD’、EE’，如圖6黑色線所示（CC’和DD’剖面線之間區(qū)域由于形變中心聚集在AA’剖面線附近，因此不再繪制剖面線進行分析），并得到五條剖面線的形變剖面圖（圖6），進而對整個礦區(qū)的地表形變時序形變特征進行分析。

圖6 礦區(qū)地表累計形變量插值圖Fig.6 Interpolation map of accumulated surface settlement in mining area

縱向貫穿礦區(qū)的剖面線AA’形成了7個沉降漏斗，見圖7（a），分別對應(yīng)圖6中的特征形變點A1-A7點。分析這7個特征形變點的地表形變時間序列變化（圖7（b））發(fā)現(xiàn)，特征點的地表形變量在監(jiān)測期內(nèi)呈現(xiàn)出逐漸增大的趨勢，且形變量基本都到達150 mm以上，最大累計形變量為217 mm，位于A5號特征點。隨著各特征點形變范圍的不斷擴大，相鄰沉降漏斗很可能會合并，如A3和A4特征點。貫穿礦區(qū)的剖面線BB’形成了2個沉降漏斗，見圖7（c），分別對應(yīng)圖5中的特征形變點B1、B2。由這2個特征形變點的形變時間序列變化圖（圖7（d））發(fā)現(xiàn)，監(jiān)測期內(nèi)B1和B2都呈現(xiàn)出逐漸增大的趨勢，且B2點從2020年1月開始形變速度加快，超過B1點，最終形變量達到197 mm。剖面線CC’形成了3個沉降漏斗，見圖7（e），分別對應(yīng)圖5中的特征形變點C1、C2、C3。由這3個特征形變點的形變時間序列變化圖（圖7（f））發(fā)現(xiàn)，監(jiān)測期內(nèi)3個特征點都呈現(xiàn)出逐漸增大的趨勢，且C2點從2019年12月開始形變速度明顯加快，最終形變量達到212 mm。剖面線DD’形成了5個沉降漏斗，見圖7（g），分別對應(yīng)圖5中的特征形變點D1-D5，由這5個特征形變點的形變時間序列變化圖（圖7（h））發(fā)現(xiàn)，監(jiān)測期內(nèi)特征點都呈現(xiàn)出逐漸增大的趨勢，且D1點從2020年3月18日形變量超過D5特征點，達到205 mm。剖面線EE’形成了2個沉降漏斗見圖7（i），分別對應(yīng)圖5中的特征形變點E1、E2，由圖7（j）可以看出監(jiān)測期內(nèi)特征點都呈現(xiàn)出逐漸增大的趨勢，且形變速率基本一致，形變量均超過200 mm。

圖7 礦區(qū)剖面圖及剖面上特征點累計形變量Fig.7 Mining area profile and cumulative settlement of characteristic points on the profile

3.3 礦區(qū)形變空間演化特征分析

采用標準差橢圓方法進一步分析礦區(qū)形變的空間分布及形變演化趨勢，以地表形變量為權(quán)重，計算標準差橢圓來分析地表形變隨時間推移的空間演變趨勢。

以2019年10月-2020年3月各月份的地表累計形變圖制作礦區(qū)各月份地表形變的標準差橢圓，見圖8。從2019年10月到2020年3月，形變量在持續(xù)增大，且在礦區(qū)不均勻分布。本文從標準差橢圓的圓心變化、面積、長短軸變化以及方位角變化等幾個角度來定量分析2019年10月-2020年3月期間礦區(qū)地表形變的空間分布特征和發(fā)展趨勢。2019年10月-2020年3月標準差橢圓參數(shù)數(shù)值見表2，各參數(shù)變化趨勢見圖9。

圖8 2019年10月-2020年3月礦區(qū)地表形變空間演化圖Fig.8 Spatial evolution of land subsidence in mining area from October 2019 to March 2020

表2 2019年10月-2020年3月標準差橢圓參數(shù)列表Table 2 Standard deviation ellipse parameter list from October 2019 to March 2020

圖9 標準差橢圓參數(shù)變化圖Fig.9 Variation of standard deviation ellipse parameters

（1）礦區(qū)形變中心變化分析

從礦區(qū)2019年10月-2020年3月年各月份的累計形變標準差橢圓的中心經(jīng)緯度以及遷移軌跡（圖10），2019年10月至2020年3月，標準差橢圓的中心逐漸向南西方向移動，礦區(qū)形變中心在這6個月內(nèi)整體向南西方向移動了505.62 m。

圖10 標準差橢圓中心位置變化圖Fig.10 Variation of center position of standard deviation ellipse

（2）礦區(qū)形變空間范圍變化分析

標準差橢圓的面積及長軸變化反應(yīng)礦區(qū)形變空間范圍的變化。礦區(qū)2019年10月-2020年3月累計形變量的標準差橢圓面積統(tǒng)計結(jié)果及長軸長度見表2，橢圓面積變化圖及長軸變化見圖9（a）（b）。2019年10月-11月，標準差橢圓的面積減少，對應(yīng)的長軸長度也相應(yīng)減少，表明在此期間礦區(qū)地表形變的范圍在減少；但從2019年11月-2020年3月，標準差橢圓的面積由8.06 km2增大到9.32 km2，長軸長度由2.21 km增大到2.64 km，表明在此期間礦區(qū)形變面積在逐漸擴大，且長軸方向的形變也在逐漸增大。

（3）礦區(qū)形變空間形狀變化分析

標準差橢圓的短軸比長軸值的變化可以表示礦區(qū)形變空間分布形狀的變化。比值越接近于1，表明礦區(qū)地表形變在各個方向演化較為均勻；而比值越接近于0，表明礦區(qū)地表形變的方向性分布特征越明顯。礦區(qū)2019年10月-2020年3月的地表累計形變標準差橢圓短軸與長軸的比值見圖9（c）。2019年10月到同年11月，短軸比長軸值變大，說明礦區(qū)地表形變在北東-南西方向形變減緩，而在北西-南東方向形變加劇；而在2019年11月-2020年3月，短軸與長軸比值逐步下降，礦區(qū)地表形變空間分布形狀表現(xiàn)出扁化趨勢，表明地表形變在長軸方向（北東-南西）上逐步發(fā)展，而在短軸方向上（北西-南東）發(fā)展相對減緩。

（4）礦區(qū)形變空間方向變化分析

標準差橢圓的方位角變化反應(yīng)礦區(qū)地表形變空間發(fā)展的主方向變化。礦區(qū)2019年10月-2020年3月的地表累計形變標準差橢圓的方位角變化見圖9（d）。2019年10月-2020年3月，礦區(qū)地表形變空間分布方位角呈逐步增大趨勢，方位角由2019年10月的16.7°增大到2020年3月的21.2°。為更直觀地表達地表形變空間分布方位角的變化趨勢，在圖11中以圖的形式來示意各期次標準差橢圓角度的增量變化，可以看出各期次標準差橢圓在空間上表現(xiàn)為逆時針旋轉(zhuǎn)，2019年10月到2020年3月方位角增大了4.5°。

圖11 準差橢圓方位角變化趨勢示意圖Fig.11 Variation trend of azimuth angle of standard deviation ellipse

3.4 結(jié)果驗證

利用2019年12月獲取的北京二號0.8 m高分辨率光學(xué)遙感影像，疊加礦區(qū)地表形變速率見圖12，可以看出：礦區(qū)整體存在不同程度的形變，部分區(qū)域形變量較大，造成在形變中心及周邊區(qū)域失相干現(xiàn)象嚴重，在圖11中表現(xiàn)為無值區(qū)。對圖中9個地表形變中心分析（圖13），發(fā)現(xiàn)形變嚴重區(qū)域在高分辨率光學(xué)影像上均表現(xiàn)為已開采礦坑或正在進行開采區(qū)域，地表形變中心位置與礦區(qū)開采位置一致。

圖12 礦區(qū)地表形變中心分布圖Fig.12 Mine location map

圖13 SABS?InSAR結(jié)果與光學(xué)影像對比圖Fig.13 Comparison of SABS InSAR results and optical images

另外，以沉降中心位置1為例，利用Google earth上與監(jiān)測期相近的3期光學(xué)影像（2019年8月、2019年9月和2020年6月）進行礦區(qū)開采活動驗證，如圖14。從3期影像的紅線范圍內(nèi)可以看出，該位置開采程度隨時間推移逐漸加劇。

圖14 開采位置1在Google earth中不同時相影像對比圖Fig.14 Comparison of different time phase images of mining location 1 in Google Earth

本次試驗，利用0.8 m高分辨率的北京二號光學(xué)遙感影像和多期次Google earth影像同時驗證了利用SBAS-InSAR進行礦區(qū)地表形變監(jiān)測的準確性和技術(shù)的可行性。

4 結(jié)論

本文利用15景Sentinel-1B干涉寬幅雷達降軌數(shù)據(jù)及0.8 m高分辨率光學(xué)影像，基于SBAS?InSAR獲取阜新某礦區(qū)2019年10月14日-2020年3月30日地表形變信息，并對礦區(qū)地表形變的時空演化特征進行分析，利用同期北京二號高分辨率光學(xué)影像對SBAS?InSAR監(jiān)測結(jié)果進行驗證，得到以下結(jié)論：

（1）SBAS?InSAR監(jiān)測結(jié)果表明，礦區(qū)整體存在明顯的不均勻地表形變現(xiàn)象，最大地表形變速率達到50-60 mm/a之間，主要分布在礦區(qū)形變的中心位置。

（2）對貫穿礦區(qū)的5條剖面及剖面線上形變特征點的時序形變分析，共識別出礦區(qū)19個沉降漏斗，位置分散，最大形變量達到217 mm。隨著特征點形變范圍的不斷擴大，相鄰特征點對應(yīng)的沉降漏斗在后期很可能會合并。從剖面圖也可以看出，隨著剖面線方向的推進，形變量在漏斗邊緣地區(qū)較低，而在距離形變中心較近的區(qū)域形變量在不斷增大，這與礦區(qū)開采的沉降規(guī)律一致。

（3）基于標準差橢圓方法對試驗區(qū)累計形變標準差橢圓的中心、面積、形狀和方向4個方面的變化分析發(fā)現(xiàn)，2019年10月-2020年3月，試驗區(qū)地表形變的形變中心由北東向南西方向遷移，形變的空間范圍也在逐步擴大；且礦區(qū)地面形變在長軸（北東-南西）逐步發(fā)展，短軸（北西-南東）發(fā)展相對減緩，導(dǎo)致橢圓形狀出現(xiàn)扁化趨勢。

（4）結(jié)合北京二號0.8 m高分辨率光學(xué)影像專家判識，發(fā)現(xiàn)礦區(qū)形變嚴重區(qū)域為已開采或者正在開采區(qū)域，驗證了利用SBAS?InSAR進行礦區(qū)地表形變監(jiān)測的準確性和技術(shù)的可行性。

綜上所述，本文利用SBAS?InSAR技術(shù)對試驗區(qū)15景Sentinel?1B降軌數(shù)據(jù)進行處理，獲取由煤礦開采活動引起的地表形變信息和典型形變特征點，定量分析了礦區(qū)地表形變的時空演化特征，結(jié)合同期高分辨率光學(xué)遙感影像對監(jiān)測結(jié)果進行驗證，為礦區(qū)的地質(zhì)災(zāi)害防治提供了參考依據(jù)。

由于本次研究缺乏監(jiān)測期內(nèi)的地面實測數(shù)據(jù)，不能有效檢驗SBAS?InSAR監(jiān)測結(jié)果的精度，后續(xù)將重點收集試驗區(qū)地面實測數(shù)據(jù)，進一步完善研究成果。