孫 軍，張 錦

（太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024）

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

1 研究區(qū)域概況

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

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

2 SAR數(shù)據(jù)與研究方法

2.1 SAR數(shù)據(jù)

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

2.2 研究方法

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

數(shù)據(jù)的處理和結(jié)果的表達(dá)均采用開源軟件，主要有ISCE2、StaMPS/MTI和QGIS3，主要步驟包括：

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

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

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

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

3 結(jié)果與分析

3.1 StaMPS方法參數(shù)分析

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

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

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

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

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

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

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

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

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

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

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

由圖4可得，5個參數(shù)值計算結(jié)果分別為：-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。計算結(jié)果表明時間窗口對解纏相位值影響較小。考慮到露天煤礦區(qū)短期內(nèi)形變量就達(dá)到較大值，2期影像采集時間間隔至少為12 d，因此設(shè)置15 d或者30 d的時間窗口均可。

3.2 基于StaMPS的SBAS-InSAR干涉結(jié)果分析

升降軌數(shù)據(jù)以起始日期為參考得到LOS向時序形變，升軌數(shù)據(jù)計算的露天煤礦區(qū)時序形變?nèi)鐖D5，降軌數(shù)據(jù)計算的露天煤礦區(qū)時序形變?nèi)鐖D6。

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

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

圖5中LOS向形變量范圍為-101～22 mm，由于正負(fù)數(shù)值差異較大，因此設(shè)置色標(biāo)范圍為-50～40 mm，隨著時間的推移LOS向形變逐漸增大且為負(fù)值，表明露天煤礦區(qū)PS目標(biāo)點遠(yuǎn)離SAR衛(wèi)星，主要存在沉降，沉降區(qū)域位于ATB礦和AJL礦的西側(cè)，礦坑內(nèi)東側(cè)區(qū)域失相干嚴(yán)重?zé)o法獲取形變量。

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

計算的露天煤礦區(qū)平均形變速率為-82～-176 mm/a，少數(shù)區(qū)域形變量達(dá)到-258 mm/a；開采區(qū)和剝離現(xiàn)場僅有少量或者完全沒有PS點，進(jìn)一步以ATB礦為研究區(qū)域，將ATB礦升降軌LOS向形變量分解為垂直向和東西向形變，ATB露天煤礦二維累積形變?nèi)鐖D7，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，沉降量最大達(dá)到-134 mm，對應(yīng)的圖8沿剖線AA′的沉降曲線，剖線處的沉降量最大達(dá)到-79 mm，呈現(xiàn)為“漏斗”狀；同時排土場內(nèi)部也存在2個失相關(guān)區(qū)域，主要是由于排土場存在新堆置土石初期的壓實過程，沉降速度較快，引起地面目標(biāo)時間失相干。ATB露天煤礦開采和剝離區(qū)失相干面積較大，頻繁的剝離和開采作業(yè)導(dǎo)致作業(yè)區(qū)內(nèi)部發(fā)生大梯度形變甚至地形變化。圖7中（b）為東西向的累積形變量，東西向累積形變量較小，形變量范圍為-29～33 mm，東西向形變量為正表示向東，負(fù)值表示形變方向向西。

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

ATB露天煤礦的排土場、剝離區(qū)和開采區(qū)均存在失相干區(qū)域，為獲得失相干區(qū)域的形變信息使用偏移追蹤技術(shù)對4景升降軌數(shù)據(jù)進(jìn)行了處理，利用斜距向偏移量解算了ATB露天煤礦二維方向的形變，偏移追蹤解算的ATB礦大梯度形變?nèi)鐖D9。

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

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

進(jìn)一步使用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露天煤礦開采區(qū)和剝離區(qū)地形變化較大，最大高程差達(dá)到-534 m，表示高程下降，通過文獻(xiàn)[10]中的高程模糊度公式可以估計干涉相位對地形變化的敏感性，研究區(qū)域的斜距為908 km，干涉對中最大的垂直基線為144 m，入射角為42.18°，C波段波長為5.6 cm，則計算的高程模糊度最小值為119 m，TanDEM-X DEM數(shù)據(jù)采集時間是2010年至2015年，而本文SAR數(shù)據(jù)采集時間在2019年，因此地形變化對時序差分干涉結(jié)果的影響是有限的，但是高強度、頻繁作業(yè)引起地面發(fā)生大梯度形變甚至地形變化，導(dǎo)致地面目標(biāo)后向散射強度變化較大從而產(chǎn)生誤匹配，使得局部區(qū)域偏移結(jié)果存在誤差。

4 結(jié) 語

以露天煤礦為研究區(qū)域，利用SBAS-InSAR和偏移追蹤技術(shù)監(jiān)測露天煤礦不同量級的形變，可識別潛在的災(zāi)害風(fēng)險，為露天煤礦安全生產(chǎn)提供服務(wù)和支持。

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

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

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