賀黎明, 裴攀科, 吳立新, 張香凝

(1. 東北大學 資源與土木工程學院， 遼寧 沈陽 110819; 2. 中南大學 地球科學與信息物理學院， 湖南 長沙 410083)

近年來,隨露天礦產資源的大量開發(fā)利用，礦區(qū)發(fā)生的滑坡災害越來越多，造成巨大的破壞和經濟損失[1].分析滑坡演化過程在風險預防中起關鍵作用[2]，其中對滑坡前地表運動特征進行分析顯得尤為重要，該特征可以幫助預測邊坡變形的發(fā)展趨勢[3]，對確保露天礦山安全生產具有重要意義.

傳統(tǒng)的滑坡形變監(jiān)測方法一般需知滑坡位置等先驗知識，主要方法有實地測量法、北斗/GNSS監(jiān)測法、分布式光纖法等.這些監(jiān)測方法都是將局部特征點作為數據采集目標，從局部映射到整體進行分析，無法全面且精準反映整個研究區(qū)域內滑坡形變的情況，不能應用于人員難以進入且處于高風險的易滑坡區(qū)，容易造成“測區(qū)未滑，滑區(qū)未測”的現象.隨遙感技術的發(fā)展，大量新的滑坡變形監(jiān)測手段涌現并得到廣泛應用，如雷達差分干涉測量[4]、無人機航空攝影測量[5]、三維激光掃描[6]等，但每一種方法都有其自身的局限性和適用范圍[7]，如何精準預測滑坡并降低滑坡災害所造成的經濟財產損失及人員傷亡仍是一個亟待解決的問題.利用時序InSAR獲取礦區(qū)地表變形場信息[8-9]，對礦區(qū)滑坡前地表運動特征進行分析，結合先驗知識對滑坡進行早期預測、預警是一種重要且有效的方法[10].

實驗以鞍山市鞍千啞巴嶺露天采場邊坡為研究對象，基于44景Sentinel-1雷達影像，對鞍千啞巴嶺露天采場2019年11月25日發(fā)生滑坡的西南側邊坡進行了遙感監(jiān)測.利用多軌SAR數據提取邊坡形變信息，對不同軌道SAR數據集反演結果進行可靠性分析，提出評價不同軌道SAR數據集監(jiān)測結果可靠性的判別依據.在此基礎上，發(fā)現鞍千啞巴嶺露天采場西南側邊坡滑坡前地表運動異常加速的變化特征，并對誘發(fā)滑坡的因素進行分析.研究結果表明，多軌時序InSAR技術可為礦區(qū)地質災害預測預警提供可靠科學依據[11]，對于保障礦山安全生產具有重要意義.

1 研究區(qū)域

研究區(qū)域為隸屬于鞍山市鞍千礦業(yè)有限責任公司的啞巴嶺露天采場，采場位于遼寧省鞍山市東部，距離鞍山市中心的直線距離為12 km，生產規(guī)模為1×107t/a.該采場的范圍為：北緯41°05′50″～41°06′30″，東經123°07′50″～123°08′40″，礦體長約1 900 m，平均水平厚度約120 m，海拔標高約為100～200 m，礦體傾角一般大于75°，傾向為NE或SW，如圖1所示.由于長期持續(xù)開采，啞巴嶺采場對比2008年投產以來發(fā)生了很大變化，其大部分區(qū)域已經變成深凹露天開采區(qū)，形成很多易滑坡危險區(qū).

圖1 鞍千啞巴嶺露天采場地理位置圖Fig.1 Location map of the Anqian Yabaling open pit

2 數據集與方法

2.1 SAR數據集

研究使用哨兵1A/B(Sentinel-1 A/B)C波段(波長約為5.6 cm)SAR影像來監(jiān)測獲取鞍千啞巴嶺采場邊坡的時間序列形變圖.針對覆蓋整個鞍千啞巴嶺采場的研究區(qū)，在歐空局網站收集了2019年6月7日至2019年11月22日共29景干涉寬幅模式(IW)的升軌(T25，T98)單視復數(single look complex，SLC)影像和2019年6月6日至11月21日共15景干涉寬幅模式(IW)的降軌(T105)單視復數影像，影像時間與空間分布分別如圖2和圖3所示.

圖2 SAR影像時間分布圖Fig.2 SAR images time distribution

圖3 實驗數據覆蓋范圍Fig.3 Coverage range of the experimental data

2.2 數據處理方法

由于SAR衛(wèi)星影像數據的限制，我國將InSAR技術應用于滑坡預測領域中的時間較晚，近些年隨Sentinel衛(wèi)星影像數據免費公開使用及國產雷達衛(wèi)星投入使用，InSAR在該領域的研究進展迅速[12].將干涉技術與合成孔徑雷達(SAR)技術結合形成的時序InSAR技術為地質災害所引起的大規(guī)模位移監(jiān)測提供了新手段，可通過不同時間從SAR傳感器收集的影像相位信息來監(jiān)測細微的地表形變，通過多個不同軌道提供的多視角結果可更全面地評價監(jiān)測結果的可靠性.

數據處理分為三大步驟，分別為基礎資料收集、SBAS-InSAR數據處理、滑坡識別與滑坡前地表運動特征分析，技術路線如圖4所示.

圖4 技術路線流程圖Fig.4 Technical flow chart

SBAS-InSAR技術即小基線集技術，是基于早期小基線差分干涉圖集反演地表形變的方法所演變而來的一項技術.SBAS-InSAR技術的基本原理是基于一定數量同一研究區(qū)域的SAR影像，通過處理大量短時空基線干涉像對，應用奇異值分解(singular value decomposition, SVD)得到相干點的位移速率和位移時間序列.

SBAS-InSAR算法基于N+1幅單視復數SAR影像，通過限定最大時間和空間基線的閾值組合成大量短基線像對以減弱時空去相關現象.一方面，由于像對空間基線較小從而避免了部分地形因素對殘余相位的影響；另一方面，在時間基線域中將干涉圖分離成幾個獨立的子集合，只要干涉圖子集合在時間域上重疊，就可以正確分離出對應相位值[13]，從而使用最小二乘法得到每一個子集合的形變時間序列.在M個干涉圖中選擇一組共同的相干像對，通過相位展開從干涉相位中分離出絕對相位值，從而將所有像素都映射到一個具有高相干性和先驗信息的參考基準點.基于所選擇的高相干點目標建立觀測方程，再將所有的小基線子集合通過奇異值分解法進行解算，從而得到所需要的時間形變序列結果.SBAS-InSAR技術將每一組D-InSAR所得到的形變結果作為SBAS-InSAR處理的觀測值，這大大增加了時間采樣率，同時通過約束基線的方法以達到時空去相關的目的，有效提高了處理效率.

3 基于時序InSAR的礦區(qū)滑坡前形變監(jiān)測

通過處理實驗數據得到了三組不同軌道影像數據集所獲取的啞巴嶺露天采場視線(line of sight，LOS)的時間序列形變結果，如圖5～7所示.其中值為負值(紅色)表示遠離衛(wèi)星的地面運動，而正值(藍色)則表示朝向衛(wèi)星的地面運動.三組不同軌道數據集所獲取的InSAR監(jiān)測結果表明啞巴嶺采場西南側邊坡存在異常變形，變形范圍不僅局限于邊邦陡坡，下面將依次對每一幅時間序列圖進行詳細分析.

圖5 啞巴嶺露天采場升軌1(T25)數據集LOS向形變 時間序列圖Fig.5 Time series graphs of deformation in LOS direction of the Yabaling open-pit slope ascending 1(T25)dataset

啞巴嶺露天采場升軌1(T25)數據集所獲取的LOS向形變時間序列圖如圖5所示，整個監(jiān)測期間采場西南側邊坡中心區(qū)域出現了較大量級和范圍的沉降.根據意大利IFFI工程所改進的滑坡活動性分級體系[14]，充分證明該區(qū)域為發(fā)生滑坡的潛在危險區(qū)，可初步根據監(jiān)測結果快速識別滑坡危險區(qū)，為進一步分析礦區(qū)滑坡前地表運動特征提供研究范圍.

啞巴嶺露天采場升軌2(T98)數據集LOS向形變時間序列如圖6所示，可知在整個監(jiān)測期間，啞巴嶺采場西南側邊坡中心區(qū)域隨時間逐漸出現了較大量級和范圍的沉降，最終在2019年11月3日至11月15日的監(jiān)測結果中達到最大形變值.值得注意的是，最后一幅影像的拍攝時間距離滑坡事件發(fā)生還有10天，按照這種沉降模式的發(fā)展趨勢，后續(xù)的10天中其最大沉降值將會繼續(xù)增加.

圖6 啞巴嶺露天采場升軌2(T98)數據集LOS向形變 時間序列圖Fig.6 Time series graphs of deformation in LOS direction of Yabaling open-pit slope ascending 2(T98)dataset

啞巴嶺露天采場降軌(T105)數據集所獲取的LOS向形變時間序列如圖7所示，在整個監(jiān)測期間，啞巴嶺采場西南側邊坡中心區(qū)域的LOS向形變值變化范圍及量級均較小.采場西南側邊坡中心區(qū)域整體變化趨勢與升軌數據集的監(jiān)測結果相反，形變值多為正值(藍色)，表現為朝向衛(wèi)星的地面運動，無法通過該數據集的監(jiān)測結果直接快速判定易滑坡危險區(qū)的范圍，需要對該數據集的監(jiān)測結果再進行詳細分析.

圖7 啞巴嶺露天采場降軌(T105)數據集LOS向形變 時間序列圖Fig.7 Time series graph of deformation in LOS direction of Yabaling open-pit slope descending(T105)dataset

綜上所述，升軌1(T25)與升軌2(T98)數據集所獲取的鞍千啞巴嶺露天采場的形變場基本一致，在采場的西南側邊坡中心區(qū)域出現了較大范圍和量級的地表形變位移.同樣，在兩組升軌(T25，T98)數據集的監(jiān)測結果與降軌(T05)數據集所獲取的監(jiān)測結果中也發(fā)現了部分區(qū)域在形變范圍和量級上存在較大差異.為了定量分析該采場滑坡前的地表形變時空演化規(guī)律，需要對這三個不同軌道數據集時序InSAR處理后的結果再進行分析.

提取跨越整個異常形變區(qū)域的三條剖線AA′，BB′及CC′如圖8所示，分別在每一條剖線上選取5個特征點(共計15個特征點)對其進行時間序列形變分析.

圖8 剖面線及特征點地理位置分布圖Fig.8 Geographical location of the profile lines and feature points

1) 縱剖面AA′：升軌1(T25)、升軌2(T98)及降軌(T105)縱剖面AA′特征點位移時間序列圖如圖9所示.

升軌1(T25)和升軌2(T98)數據集的結果顯示所有特征點均表現出了持續(xù)下沉的趨勢，且該斷面的中心區(qū)域P04，P07兩個特征點的累計形變值最大，如圖9a,9b所示.從圖9c的降軌(T105)

圖9 升軌1(T25)、升軌2(T98)及降軌(T105) 縱剖面AA′特征點位移時間序列圖Fig.9 Ascending1(T25), ascending 2(T98)and descending(T105)displacement time series diagram of the characteristic pointsalong the longitudinal section AA′(a)—升軌1(T25); (b)—升軌2(T98); (c)—降軌(T105).

數據集所獲取的形變結果可知，該斷面所有特征點的形變范圍較小，整體呈線性變化，且變化趨勢很小.采礦所引發(fā)的地表形變與時間常常呈非線性關系.從圖9b升軌2(T98)中的LOS向形變時間序列中不難發(fā)現該斷面所有的特征點在2019年6月12日至2019年10月10日(紅色虛線處)均處于一個緩慢線性下沉的趨勢，但在2019年10月10日(紅色虛線處)至2019年11月15日忽然出現了加速沉降的現象，這一形變演化過程對應于礦區(qū)單點沉降中的初始沉降、加速沉降兩個過程.說明在發(fā)生滑坡前該斷面所有的點均存在一個加速下沉的特征，這符合采礦引發(fā)地表形變的規(guī)律，可以作為分析礦區(qū)滑坡前地表運動特征的參考依據.

2) 縱剖面BB′：升軌1(T25)、升軌2(T98)及降軌(T105)縱剖面BB′特征點位移時間序列如圖10所示.

圖10 升軌1(T25)、升軌2(T98)及降軌(T105) 縱剖面BB′特征點位移時間序列圖Fig.10 Ascending1(T25), ascending2(T98)and descending(T105)displacement time series diagram of characteristic points along the longitudinal section BB′(a)—升軌1(T25); (b)—升軌2(T98); (c)—降軌(T105).

升軌1(T25)和升軌2(T98)數據集的結果顯示所有特征點均表現出持續(xù)下沉的趨勢，但位于縱剖面兩端的兩點P02和P04的累計形變值較小，中心區(qū)域P05，P08，P11的累計形變值較大，說明形變的中心區(qū)域位于特征點P05，P08，P11之間.同樣的，在圖10b升軌2(T98)中不難發(fā)現P05，P08，P11在2019年6月12日至2019年10月10日(紅色虛線處)期間，所有特征點均處于一個緩慢線性沉降的趨勢，但在2019年10月10日(紅色虛線處)至2019年11月15日，P05，P08，P11這三個特征點忽然出現了加速沉降的現象，說明在發(fā)生滑坡前該斷面中心區(qū)域的點均存在一個加速下沉的特征.從圖10c降軌(T105)數據集所獲取的形變結果可知，在整個監(jiān)測期間，該斷面所有特征點的形變范圍相比縱剖面AA′的范圍有所擴大，整體仍呈線性變化，且變化趨勢很小.

3) 縱剖面CC′：升軌1(T25)、升軌2(T98)及降軌(T105)縱剖面CC′特征點位移時間序列如圖11所示.

升軌1(T25)數據集的結果顯示，所有特征點的累計沉降量均小于縱剖面AA′與BB′，如圖11a所示.升軌2(T98)數據集的結果顯示，位于縱剖面CC′上所有的特征點在整個監(jiān)測期間均以線性形式進行沉降，但沉降速率及沉降范圍相比縱剖面AA′與CC′都要小，未能發(fā)現存在加速下沉的特征點，如圖11b所示.降軌(T105)數據集的結果顯示，在整個監(jiān)測期間，該斷面所有特征點仍呈線性變化，且變化趨勢很小，如圖11c所示.以上結果表明，縱剖面CC′已經是形變范圍的邊界區(qū)域，屬于線性穩(wěn)定沉降區(qū)域.

圖11 升軌1(T25)、升軌2(T98)及降軌(T105) 縱剖面CC′特征點位移時間序列圖Fig.11 Ascending1(T25), ascending2(T98) and descending(T105)displacement time series diagram of characteristic points along the longitudinal section CC′(a)—升軌1(T25); (b)—升軌2(T98); (c)—降軌(T105).

通過對縱剖面AA′、BB′與CC′的分析以及從升軌2(T98)數據集所獲取的所有特征點的形變時間序列圖來看，不難發(fā)現P01，P04，P05，P07，P08，P10，P11，P13這8個特征點在整個監(jiān)測期間存在滑坡前的加速時期(2019年10月10日至2019年11月15日)，如圖12所示.圖中，紅色點表示這些出現異常加速現象的點，白色點表示在整個監(jiān)測期間表現為線性勻速沉降的點.

圖12 不同沉降模式特征點分布圖Fig.12 Distribution map of characteristic points of different settlement modes

以線性行為為特征區(qū)域，是極限平衡狀態(tài)下的不穩(wěn)定性指標.沒有跡象表明這種行為可以持續(xù)多長時間，或可能會演變?yōu)榧铀傩袨榛蚍€(wěn)定行為.在這種情況下，如果這些穩(wěn)定區(qū)域圍繞加速區(qū)域，就極有可能受到加速區(qū)域移動的影響[15].從2019年10月10日至2019年11月25日(滑坡事件發(fā)生時)在啞巴嶺礦區(qū)滑坡位置頂部后緣位置存在一個扇形加速沉降區(qū)域，如圖13所示.這成為直接誘導滑坡事件發(fā)生的主要因素，根據這一現象可以圈定出一個可能發(fā)生滑坡的范圍，即易滑坡區(qū).

圖13 易滑坡區(qū)示意圖Fig.13 Sketch map of the landslide area

4 InSAR多軌道觀測值之間的可靠性檢驗

為了解釋來自三組不同觀測軌道InSAR測量結果之間的可靠性，比較分析了從三組不同的哨兵1A / B(Sentinel-1 A / B)數據集獲取的鞍千啞巴嶺采場西南側邊坡形變區(qū)域中所有像素的InSAR平均形變測量值，整個過程分為以下兩個步驟：

1) 將時序InSAR反演所獲得的LOS向形變場重采樣為100 m×200 m的網格，從而達到最小化地理定位誤差(解決測量之間的空間不匹配問題)的目的；2)異常形變區(qū)域(藍色)光學影像如圖14所示，選擇了三個軌道共有的局部形變區(qū)域中所有形變監(jiān)測數據不為空值的像素共計324個.異常形變區(qū)域所有像素平均形變時間序列圖如圖15所示，統(tǒng)計了三條不同軌道數據集在該形變區(qū)域內所有像素LOS向平均形變值的時間序列，并分別進行線性擬合.根據線性擬合的結果可以發(fā)現：升軌1(T25)與升軌2(T98)的判定系數R2分別為0.89與0.98，基本表現出相同的沉降趨勢，而降軌(T105)與升軌1(T25)和升軌2(T98)趨勢線呈相反關系.這些現象表明從不同數據集獲得的不同沉降模式彼此相關.升軌的兩個數據集監(jiān)測所得到的結果在2019年6月至2019年11月之間表現出相似的線性運動，使這兩個數據集所得到的監(jiān)測結果具有可信度.

圖14 異常形變區(qū)域光學影像Fig.14 The optical image of the abnormal deformation area

圖15 異常形變區(qū)域所有像素平均形變時間序列圖Fig.15 Time series graph of the average deformation of all pixels in the abnormal deformation area

同時不難發(fā)現，只有升軌2(T98)數據集的監(jiān)測結果顯示了在滑坡前存在一個突然加速的區(qū)域使得整個滑坡事件得以發(fā)生.對整個形變區(qū)域平均形變時間序列線性擬合結果進行分析，升軌2(T98)數據集的判定系數R2為0.98，這說明整個大的形變區(qū)域是以線性模式進行均勻沉降，這符合滑坡事件發(fā)生的基本沉降規(guī)律，所以針對鞍千啞巴嶺露天采場，以整個大范圍的平均形變時間序列的線性擬合效果的差異來判斷不同軌道數據集之間監(jiān)測結果的可靠性是合理的.

比較相同或不同軌道之間時序InSAR反演結果是定量檢查測量一致性的常見做法[16].本次實驗所觀察到的三組數據集之間的差異可能源于采樣不足、地理參考不確定、季節(jié)性偏差、長時間間隙、不同的影像獲取時間及不同的軌道參數等因素.在本次實驗數據集中觀察到的升軌和降軌所獲取的結果之間產生相反的變化趨勢.原因主要包括：三組不同軌道數據集擁有不同的入射角及航向角，其對研究區(qū)域真實形變的敏感度不同，由于采礦所引發(fā)的地表沉降導致其LOS向形變中可能會包含水平形變的影響，這些因素共同造成了差異性.

5 礦區(qū)滑坡因素分析

影響巖質邊坡穩(wěn)定性的因素可分為兩個方面，即內部與外部因素.內部因素為啞巴嶺采場西南側巖質邊坡的一些固有屬性，而外部因素通常會引起邊坡失穩(wěn)從而發(fā)生滑坡，如工程爆破、降雨、地震等.盡管目前沒有特別有力的證據表明鞍千啞巴嶺露天采場在2019年11月25日的滑坡是由于采場西南側邊坡的不穩(wěn)定因素而引起的，但考慮到該礦區(qū)的特殊性質，強烈和頻繁的開采活動及礦石爆破工作有可能構成了邊坡失穩(wěn)的主要外部因素.

孔隙壓力是滑坡的主要驅動因素[17]，季節(jié)性滑動速率的變化與降水驅動的瞬時孔隙壓力的增加有關.鞍山地區(qū)位于遼河平原的東南部邊緣，暖溫帶大陸性季風普遍影響當地氣候.這使得該地區(qū)四季分明、干冷同期、雨熱同季.鞍山地區(qū)2019年平均降水量約為686 mm，其中85%的降雨集中在5月至9月，但該年11月份的降水量卻異常增多，占據非雨季降水量的34%.2019年11月25日位于該地區(qū)的啞巴嶺露天采場邊坡發(fā)生一次大型滑坡.遼寧鞍山53339號區(qū)站2019年11月的日降水量數據和鞍山地區(qū)11月份的平均溫氣與天氣數據如圖16所示.

圖16 鞍山地區(qū)2019年11月部分氣象數據Fig.16 Part of the meteorological data of Anshan area in November 2019

由圖可知：11月12日開始到滑坡前陸陸續(xù)續(xù)一直在下雨或下雪，11月19日氣溫升高，在之前積雪融化的基礎上又開始降雨.由于前期降水持續(xù)時間較長，降水量也很大，尤其是在發(fā)生滑坡的前兩天，直接導致啞巴嶺西南側邊坡裂隙中進入大量的雨水，長時間大范圍的降水導致軟弱層(滑面)形成，誘導本來處于線性穩(wěn)定沉降的坡面發(fā)生滑坡事件.

6 結 論

1) 時序InSAR監(jiān)測結果顯示，鞍千啞巴嶺露天采場西南側邊坡存在異常形變，可通過該監(jiān)測結果初步定位滑坡區(qū)域.

2) 啞巴嶺采場西南側邊坡頂部后緣位置在發(fā)生滑坡前的一段時間范圍內(約45 d)出現異常強烈的加速現象，進而導致原本處于線性形變的相對穩(wěn)定區(qū)域出現了滑坡.

3) 不同軌道數據集針對同一研究區(qū)域所獲取的結果具有差異性，形變區(qū)域中所有像素的平均形變時間序列的線性擬合回歸效果越好，對應軌道數據集的監(jiān)測結果越具有可靠性.

4) 啞巴嶺采場滑坡的主要內部影響因素是西南側邊坡巖質的不穩(wěn)定結構，長時間大范圍的降水導致軟弱層(滑面)形成，因此持續(xù)的大規(guī)模降水是啞巴嶺礦區(qū)發(fā)生滑坡的外部重要誘因.