基于SBAS-InSAR和偏移追蹤技術的露天煤礦地面形變監測

孫 軍，張 錦

（太原理工大學 礦業工程學院，山西 太原 030024）

合成孔徑雷達干涉技術（interferometric synthetic aperture radar，InSAR）已應用于地面形變監測，具有分辨率高、監測范圍廣、多極化觀測的特點。時序差分干涉方法較多，主要包括永久散射體差分干涉技術（persistent scatterers InSAR，PSInSAR）、分布式散射體技術（SqueeSAR）、小基線集技術（small baseline subsets InSAR，SBAS-InSAR）等[1-4]。Hooper等[5]提出了斯坦福永久散射體法（stanford method for PS,StaMPS），認為小基線集組合干涉方法可以增加多余觀測；Hooper等[6]根據小基線集原理進一步提出基于StaMPS的SBAS方法，利用少量SAR影像（至少5期）就可以形成較多干涉對，將振幅差分離差作為閾值選擇慢失相干濾波相位（slowly-decorrelating filtered phase，SDFP）像素點，可以有效地增加PS點的數量，同時StaMPS方法對DEM精度要求不高，可用于露天煤礦區大梯度形變監測。在井工煤礦開采沉陷監測方面，InSAR與偏移追蹤技術相結合能夠實現高強度采區的監測[7-8]，進一步通過對升降軌數據的處理可以解算煤礦開采區的二維形變[9]；而露天煤礦形變監測中InSAR技術多用于監測露天煤礦大梯度沉降，分析邊坡穩定性和非法采礦等[10-12]，與偏移追蹤技術相結合用于形變監測的研究較少。為此，使用開源SAR數據處理軟件ISCE2和StaMPS4.1分別對覆蓋平朔礦區的18期和15期升降軌數據進行干涉處理，重點分析了安太堡露天煤礦的時序形變，根據升降軌觀測幾何解算出數據采集時間段內露天煤礦二維方向的累積形變量，通過偏移追蹤技術計算了露天煤礦內失相干區域的大梯度形變。

1 研究區域概況

平朔礦區位于山西省朔州市朔城區、平魯區境內，礦區范圍為東經112°10′～113°30′、北緯39°23′～39°37′，地貌類型為黃土低山丘陵，海拔高度1.3～1.4 km，地形起伏較大。礦區東西寬18 km,南北長21 km，總面積380 km2。目前平朔礦區露天煤礦共有3座，露天煤礦由南至北分別為安家嶺礦（AJL）、安太堡礦（ATB）和東露天煤礦（DLT）。除此之外還有地下開采的井工礦，井工一礦位于AJL西側；井工二礦位于AJL和ATB之間；井工三礦位于ATB北側及DLT西側。由于井工和露天開采并存，平朔礦區環境地質問題比較復雜，分布范圍廣泛，種類繁多，在露天煤礦區主要存在邊坡穩定性和地面沉陷問題。研究區域如圖1。

選擇安太堡露天煤礦為研究區域，ATB露天煤礦東西長5.7 km，南北長3.1 km，該露天礦的開采方向為自西向東，分別分布有排土場、開采區和剝離區。目前ATB中主要剝離方向為東北方向，排土作業區位于礦坑的西北側，各作業區之間分布有較多運輸道路，道路兩邊分布有開采平盤和邊坡。

2 SAR數據與研究方法

2.1 SAR數據

所使用的SAR影像由Sentinel-1采集，Sentinel-1是由歐洲航天局發射的對地觀測衛星，主要搭載C波段SAR傳感器。Sentinel-1是由A、B 2顆衛星組成的觀測星座，單顆衛星的回訪周期為12 d，升降軌組合衛星對地重訪周期可縮短為6 d。可以提高對地觀測頻率。地面形變監測中使用干涉寬幅（interferometric wide swath，IW）模式下采集的SAR數據，IW模式采用漸進式掃描（terrain observation by progressive Scans，TOPS）對地面觀測獲取3個子條帶。相比于傳統ScanSAR模式，TOPSAR技術可以實現相同的分辨率和覆蓋范圍，避免ScanSAR采集模式的扇貝效應，信噪比也有所提高，IW模式下SLC影像幅寬為250 km，空間分辨率約為5 m×20 m。文中升軌數據總共18期，采集時間為2019年1月8日至2019年7月31日，極化方式為雙極化即VV+VH；降軌數據共有15期，采集時間為2019年1月2日至2019年7月25日，升降軌數據分別間隔204 d和205 d，升降軌SAR數據參數見表1；此外使用精密定軌星歷數據更新原數據的軌道矢量，使用30 m分辨率的SRTM DEM模擬和去除地形相位。

2.2 研究方法

時序差分干涉方法中使用基于StaMPS的SBAS方法，該方法根據小基線集原理組成多主影像干涉對，通過對像素點的相位穩定性分析選擇PS點，主要包括4部分：①依據時間和空間基線形成干涉圖；②依據振幅差分離差指數選擇PS侯選點（persistent scatterers candidates，PSC）；③分析PSC的相位穩定性，通過估計PSC的噪聲，計算概率選出最終的PS點；④對干涉相位進行解纏，去除殘余相位和噪聲相位，計算形變相位和形變量。

數據的處理和結果的表達均采用開源軟件，主要有ISCE2、StaMPS/MTI和QGIS3，主要步驟包括：

1）利用ISCE2對SAR影像預處理，形成差分干涉對。利用精密軌道數據將主、副影像精確配準，對應像素相位共軛相乘形成初始干涉對，去除干涉相位中的參考橢球相位和地形相位，形成單主影像差分干涉對。

2）在StaMPS/MTI中形成小基線集干涉對，通過時間基線、空間基線和相干性閾值確定小基線集組合數量。Hooper等[13]設置的閾值為5年、1.1 km；Tang Wei等[10]使用哨兵數據分析露天煤礦形變時設置的閾值為50 d、100 m。本文設置的時間基線閾值為60 d，同時為避免形成多個小基線集，設置空間基線閾值為220 m，相關性閾值設為0.5。由升軌數據形成的干涉對數目為36，最大的時間基線和空間基線分別為132 d和-144 m；降軌數據形成的干涉對數目為23，干涉對的時間和空間基線最大值為108 d、89.80 m。

3）通過StaMPS方法構建的相位模型分離出形變相位，根據最小二乘原理從小基線集中計算出每個SDFP點視線向（line of sight，LOS）形變量和形變速率。在QGIS3中將升降軌結果通過空間插值方法采樣到同名點，采用Samsonov等[14]提出的補償因子方法統一數據采集時間，由于SAR衛星對南北向形變不敏感，因此根據LOS向與垂直向、東西向的幾何關系計算了露天煤礦二維累積形變量。

利用基于SAR影像強度的偏移追蹤技術計算了露天煤礦區的大梯度形變，偏移追蹤技術利用移動窗口內的強度信息計算歸一化互相關值，互相關值最大的位置即為最佳匹配位置，從而估計出斜距向和方位向偏移量。基于SNAP8.0實現偏移量計算，主要有3個步驟：①主影像和副影像的預處理：包括更新軌道數據、分離子條帶、輻射校正、提取bust；②主副影像的配準：采用DEM輔助配準法可以減小因地形起伏引起的偏移；③以主影像為參考，在影像中均勻分布像素點，選擇一定尺寸的窗口，逐像素計算偏移量。偏移追蹤技術的精度約為SAR影像分辨率的1/10～1/20，Sentinel-1影像經過多視處理后像素大小為13.9 m，因此偏移追蹤測量精度要高于1.39 m，同時文中選擇的匹配窗口大小為32像素×32像素，過采樣因子為16，偏移速率閾值為1.5 m/d，與前述類似，將由升降軌數據計算的斜距向偏移量分解到垂直向、東西向。

3 結果與分析

3.1 StaMPS方法參數分析

StaMPS方法根據像素點的殘余噪聲相位的穩定性選擇PS點，同時考慮到露天煤礦采掘和排土作業頻繁，一定時間段內局部區域地形會發生變化，因此分析了PSC點的噪聲估計中最大地形誤差對相位擬合和相干性估計的影響。在相位解纏步驟中又計算了不同大小的采樣格網和時間窗口對應的解纏結果。StaMPS方法是對已去除參考橢球和地形相位的干涉圖展開分析的，將剩余相位分為空間相關相位和空間不相關相位2部分，通過低通和自適應濾波去除空間相關相位，再利用最小二乘法去除地形誤差相位后得到殘余噪聲相位，通過式（1）統計PSC點的噪聲相位在時序中的變化。

式中：γx為殘余相位隨著時間的變化，也即反映了相位的穩定性；N為干涉對的數量；φx,i為纏繞相位，包括形變相位、大氣相位、視角誤差相位、軌道誤差相位和噪聲相位5部分；下標x，i為第i幅干涉圖中的第x個像素；為φx,i的空間相關部分；為空間不相關視角誤差相位，主要包含地形誤差相位為殘余的噪聲相位。

式中：λ為波長；B⊥,x,i為垂直基線值；為空間不相關視角誤差。

最大地形誤差用于線性擬合DEM誤差相位，設置最大地形誤差分別為5、10、20、30、40、50 m，其中20 m為默認值，分別計算DEM誤差相位的擬合值，計算結果顯示擬合值無明顯差異，主要原因是空間基線較短，對高程誤差不敏感。Hooper等[6]研究Alcedo火山的形變中空間基線范圍為-917～976 m，文中基線范圍為-99～106 m，空間不相關視角誤差的估計如圖2。

圖2 空間不相關視角誤差的估計Fig.2 Estimated spatial uncorrelated look angle error

由圖2可知，某個像素點的γx值與空間不相關視角誤差近似于線性關系，該點的纏繞相位相對穩定，分布于-1～1 rad之間，擬合的結果相對較好。參考相關文獻設置的最大地形誤差為5～20 m，本文中設置20 m的地形誤差用于估計噪聲。

StaMPS方法中采用3D相位解纏算法，相位解纏中采樣格網主要用于對PS點進行重采樣，不同尺寸采樣格網計算的解纏相位如圖3。

圖3 不同尺寸采樣格網計算的解纏相位Fig.3 Unwrapping phase in different sampling grids

圖3（a）～圖3（d）是格網尺寸為25、50、100、200 m對應的解纏結果，格網尺寸25、50 m時，解纏的相位值范圍均為-17.4～6.2 rad；格網尺寸100 m和200 m對應的結果為-13.1～5.2 rad和-10.5～5.2 rad，總體看解纏結果受采樣格網尺寸的影響較小，相位發生變化的區域主要是礦坑排土場，這些區域形變梯度較大并且存在失相干區域，解纏結果存在一定誤差，同時較大的格網具有抑制噪聲和平滑的效果，Hooper、Lazecky等[6,15]研究火山和煤礦開采區沉降時，將采樣格網均設置為100 m，本文中也采用該值。

時間窗口用于在一定時間范圍內估計PS點在每個干涉對的相位噪聲，從而在時序中對相位進行平滑。分別設置時間窗口為15、30、60、120、730 d并計算解纏相位，其中默認值為730 d，不同時間窗口對應的解纏相位如圖4。

圖4 不同時間窗口計算的解纏相位Fig.4 Unwrapping phase in different time windows

由圖4可得，5個參數值計算結果分別為：-13.9～5.9 rad、-13.5～5.7 rad、-13.9～5.7 rad、-12.9～5.2 rad、-13.1～5.2 rad。計算結果表明時間窗口對解纏相位值影響較小。考慮到露天煤礦區短期內形變量就達到較大值，2期影像采集時間間隔至少為12 d，因此設置15 d或者30 d的時間窗口均可。

3.2 基于StaMPS的SBAS-InSAR干涉結果分析

升降軌數據以起始日期為參考得到LOS向時序形變，升軌數據計算的露天煤礦區時序形變如圖5，降軌數據計算的露天煤礦區時序形變如圖6。

圖5 升軌數據計算的露天煤礦區時序形變Fig.5 Time series deformation in open pit coal mine calculated by ascending data

圖6 降軌數據計算的露天煤礦區時序形變Fig.6 Time series deformation in open-pit coal mine calculated by descending data

圖5中LOS向形變量范圍為-101～22 mm，由于正負數值差異較大，因此設置色標范圍為-50～40 mm，隨著時間的推移LOS向形變逐漸增大且為負值，表明露天煤礦區PS目標點遠離SAR衛星，主要存在沉降，沉降區域位于ATB礦和AJL礦的西側，礦坑內東側區域失相干嚴重無法獲取形變量。

由圖6可知，由降軌數據計算的LOS向時序形變量范圍為-86～27 mm，升降軌計算的形變量值的差異主要是由于衛星軌道和入射角不同，對地面形變的敏感性不同，但形變趨勢基本一致即ATB礦和AJL礦西側區域存在較大量級的沉降，礦坑內部失相干嚴重，通過衛星影像比對發現西側沉降區主要是ATB礦和AJL礦的排土場，失相干嚴重區域位于礦坑內的排土和剝離區。

計算的露天煤礦區平均形變速率為-82～-176 mm/a，少數區域形變量達到-258 mm/a；開采區和剝離現場僅有少量或者完全沒有PS點，進一步以ATB礦為研究區域，將ATB礦升降軌LOS向形變量分解為垂直向和東西向形變，ATB露天煤礦二維累積形變如圖7，ATB露天煤礦沿剖線AA′的累積沉降量如圖8。

圖7 ATB露天煤礦二維累積形變Fig.7 2D deformation of ATB Open-pit Coal Mine

圖8 ATB露天煤礦沿剖線AA′的累積沉降量Fig.8 Deformation along AA′profile of ATB Open-pit Coal Mine

圖7（a）為垂直向的累積形變量，可知排土場存在1個較大的沉降中心，沉降量范圍為-119～-60 mm，沉降量最大達到-134 mm，對應的圖8沿剖線AA′的沉降曲線，剖線處的沉降量最大達到-79 mm，呈現為“漏斗”狀；同時排土場內部也存在2個失相關區域，主要是由于排土場存在新堆置土石初期的壓實過程，沉降速度較快，引起地面目標時間失相干。ATB露天煤礦開采和剝離區失相干面積較大，頻繁的剝離和開采作業導致作業區內部發生大梯度形變甚至地形變化。圖7中（b）為東西向的累積形變量，東西向累積形變量較小，形變量范圍為-29～33 mm，東西向形變量為正表示向東，負值表示形變方向向西。

3.3 基于偏移追蹤技術的露天礦形變分析

ATB露天煤礦的排土場、剝離區和開采區均存在失相干區域，為獲得失相干區域的形變信息使用偏移追蹤技術對4景升降軌數據進行了處理，利用斜距向偏移量解算了ATB露天煤礦二維方向的形變，偏移追蹤解算的ATB礦大梯度形變如圖9。

圖9 偏移追蹤解算的ATB礦大梯度形變Fig.9 Large gradient deformation of ATB Mine by offset tracking

由圖9可得，垂直向形變量范圍為-46.89～53.29 m，東西向形變量范圍為-39.48～48.50 m，最大的垂直向和東西向形變量主要位于開采和剝離現場，同時形變分布呈現出斑塊狀，這與偏移追蹤算法通過窗口逐像素移動計算有關。盡管設置的偏移追蹤窗口從32像素到128像素，但是偏移量追蹤結果仍然受到局部噪聲的影響，同時過大的窗口會產生空值，因此本文設置的窗口尺寸為32像素×32像素。

進一步使用2期時間間隔為48 d，空間基線長為116 m的SLC影像反演了ATB露天煤礦的DEM，將反演的DEM與90 m分辨率的TanDEM-X DEM相減得到高程差圖，反演的DEM和TanDEM-X DEM的高程差如圖10。

圖10反演的DEM和TanDEM-X DEM的高程差Fig.10 Elevation difference between inverted DEM and TanDEM-X DEM

從圖10中可以看出，ATB露天煤礦開采區和剝離區地形變化較大，最大高程差達到-534 m，表示高程下降，通過文獻[10]中的高程模糊度公式可以估計干涉相位對地形變化的敏感性，研究區域的斜距為908 km，干涉對中最大的垂直基線為144 m，入射角為42.18°，C波段波長為5.6 cm，則計算的高程模糊度最小值為119 m，TanDEM-X DEM數據采集時間是2010年至2015年，而本文SAR數據采集時間在2019年，因此地形變化對時序差分干涉結果的影響是有限的，但是高強度、頻繁作業引起地面發生大梯度形變甚至地形變化，導致地面目標后向散射強度變化較大從而產生誤匹配，使得局部區域偏移結果存在誤差。

4 結 語

以露天煤礦為研究區域，利用SBAS-InSAR和偏移追蹤技術監測露天煤礦不同量級的形變，可識別潛在的災害風險，為露天煤礦安全生產提供服務和支持。

1）SBAS-InSAR技術計算的露天煤礦形變精度可達毫米級；偏移追蹤技術計算的形變量在分米級以上，2種技術在形變監測中可實現互補。

2）SBAS-InSAR中通過限制時空基線，可以減小地形相位的貢獻，降低地形誤差對形變結果的影響；露天礦邊端幫、運輸道路、邊坡等均具有較強的散射性，而偏移追蹤易受局部強散射體的影響，同時哨兵數據的分辨率較低，導致在局部地形變化區域偏移追蹤結果偏大且存在誤差。

3）與井工開采引起的地面形變不同，露天煤礦短期內形變量就達到較大值，對SAR衛星的時間分辨率和數據的空間分辨率要求較高，另外露天煤礦由不同的功能區組成，剝離區、開采區和排土場形變趨勢不同，因此后續需要結合形變模型，對不同區域的形變展開針對性研究。