劉東烈，陸超然，郭金城，蔡杰華，董 杰，張 路，廖明生

（1. 貴州省自然資源衛(wèi)星應用中心，貴州 貴陽 550001；2. 航天宏圖信息技術股份有限公司，北京 100195；3. 武漢大學測繪遙感信息工程國家重點實驗室，湖北 武漢 430079；4. 武漢大學遙感信息工程學院，湖北 武漢 430079）

貴州省位于我國西南地區(qū)，地形地貌復雜、雨源性水系發(fā)育，加之水利工程、城市建設、礦山開采等促進經濟發(fā)展的工程項目日趨增多，在自然和人為因素的雙重影響下，地質災害頻繁發(fā)生[1]，給人民生命財產安全和社會經濟發(fā)展帶來了不可逆轉的傷害，開展大范圍地質災害隱患排查刻不容緩。合成孔徑雷達干涉（InSAR）采用主動式微波成像技術，空天作業(yè)方式實現了全天時、全天候、高精度、大范圍對地觀測[2]，可滿足貴州省“天無三日晴、地無三尺平”的復雜觀測條件。差分干涉測量技術（D-InSAR）是目前InSAR提取地表形變中的常用技術，已成功應用于地震[3]、礦區(qū)沉降[4]、滑坡[5]等地質災害監(jiān)測中；但隨著時間和空間基線的拉長，差分干涉對的相干性降低[6]，且兩期影像成像時大氣條件的差異也會對形變結果造成極大干擾[7]，影響解纏和解譯過程，從而限制了其應用范圍[8]。隨著SAR 傳感器的不斷進步和SAR數據的積累，可靠性更高的時序InSAR（TSI）技術應運而生，如永久散射體干涉測量[9]、小基線集[10]（SBAS）、干涉點目標分析[11]和StaMPS方法[12-14]等。這些時間序列方法通過對一系列差分干涉相位進行聯合分析，有效抑制了失相干、DEM誤差和大氣相位延遲等的影響[15]。近年來，眾多學者采用InSAR 技術開展了地質災害隱患早期識別研究，如Shi X G[16]等集成3個相鄰軌道的ALOS/PALSAR數據獲取了三峽地區(qū)秭歸—奉節(jié)段的斜坡形變，并探測到30 處活動性滑坡；Zhao C Y[17]等采用TSI技術在金沙江烏東德水電站附近區(qū)域發(fā)現了22處滑坡，并探究了金坪子滑坡的時空形變特征；廖明生[18]等針對滑坡形變監(jiān)測中的關鍵問題，先后提出了相干散射體TSI 方法[19]、基于點目標的像素偏移量追蹤方法[20]、垂直分層大氣延遲改正模型[21]等，為滑坡形變監(jiān)測應用中常見的重難點問題提供了有效的解決方案。在山區(qū)開展InSAR 處理分析時，需引入外部DEM 去除地形相位。常用的開放DEM 數據的空間分辨率和高程測量精度有限，且或多或少已“過時”，與真實地表形態(tài)存在差異。已有學者對比了不同DEM 數據對TSI 結果的影響[22-23]，但未使用高分辨率高精度DEM （LiDAR DEM）進行分析。

本文針對山區(qū)地質災害隱患早期識別的需求與現狀，首先在貴州省西南山區(qū)的興仁縣開展D-InSAR與TSI的對比分析，以驗證TSI方法在山區(qū)地質災害隱患探測中的優(yōu)勢和必要性；然后基于貴州省水城縣重點區(qū)域AW3D30、SRTM-3 和LiDAR DEM，系統研究了外部DEM 對山區(qū)TSI 數據處理結果的影響，為山區(qū)SAR數據處理提供參考。

1 研究方法

D-InSAR對兩景雷達影像進行相位差分，處理簡單，能快速獲取觀測時段內的地表形變信息，但失相干、大氣相位延遲、DEM誤差等因素干擾了形變信號的提取。StaMPS-SBAS方法在多對差分干涉結果的基礎上進一步開展相位分析（圖1），假設形變相位空間相關，首先采用低通自適應濾波器迭代計算相位的穩(wěn)定性，并以此為測度選取相干點目標；再估計并剔除DEM誤差；然后在平面、空間和時間3個維度上開展相位解纏，并將小基線集的多主影像相位轉到單主影像相位；最后進行時空濾波，去除大氣延遲相位，獲取形變速率和形變序列。該方法能削弱大氣延遲、DEM誤差、去相干噪聲等干擾相位的影響，可獲取精準、可靠的時序形變結果。

圖1 StaMPS-SBAS方法流程圖

在上述過程中，需引入外部DEM模擬和去除地形相位，通常是將DEM 與主影像配準，結合成像幾何，根據式（1）計算地形相位：

式中，h為地面高；λ為雷達波長；B⊥為基線垂直于視線向（LOS）的分量；R為雷達天線到地面目標的距離；θ為雷達入射角。

由此可知，DEM 直接影響了干涉相位，DEM 誤差越大，干涉相位中的殘余地形誤差相位越大；而TSI 分析則是在大量差分干涉結果的基礎上進行時空域濾波、迭代相位分析等工作，DEM誤差對其影響無法直觀、定量地表示出來，仍需探究。

2 研究區(qū)域與數據來源

2.1 研究區(qū)概況

研究區(qū)位于貴州省西南山區(qū)，分別隸屬于黔西南布依族苗族自治州興仁縣和六盤水市水城縣。興仁縣是地質災害多發(fā)縣，地形起伏劇烈，地層和地質構造復雜，地貌多樣且侵蝕剝蝕地貌十分發(fā)育；雨量充沛，年均降水量為1 286 mm，最大降水量達1 618 mm；金礦和煤礦資源富集，礦區(qū)開采活動強烈；地質災害主要包括滑坡、崩塌、泥石流、地裂縫、地面塌陷、地面沉降等類型，其中滑坡災害占大部分，空間上呈局部密集式分布特征，主要表現為人類工程活動強烈的地區(qū)成群出現，其中礦山開采導致的地質災害占比最高。水城縣地形復雜，部分山區(qū)呈現出斷崖式的劇烈起伏，坡度甚至達到45°~60°，山體巖溶地貌發(fā)育，地表穩(wěn)定性差；雨量充沛，年均降水量為1 100 mm，降雨集中在5—9 月；煤炭資源豐富，存在不同程度的開采活動。2019年7 月23 日雞場鎮(zhèn)坪地村岔溝組發(fā)生了特大型山體滑坡，造成了重大人員傷亡，該地區(qū)是貴州省地質災害防治的重點區(qū)域。

2.2 SAR數據

本文采用4 景L 波段的ALOS-2/PALSAR2 超精細模式（Ultrafine）數據[24]和117 景升軌Sentinel-1 寬幅干涉模式數據[25]（表1）。Sentinel-1 數據面向全球開放，重訪周期為12 d，軌道穩(wěn)定，充足、連續(xù)、穩(wěn)定的觀測數據保證了一定時間內的相位相干性。Sentinel-1 數據的時空基線見圖2，共構成345 對小基線對。

表1 星載SAR數據參數

圖2 Sentinel-1數據集的時空基線組合

2.3 DEM數據

本文分別利用SRTM-3、AW3D30 和LiDAR DEM數據模擬InSAR 地形相位，相關參數對比見表2，高程基準均為EGM96水準面，無需進行基準校正。由于LiDAR DEM 數據的精度遠高于AW3D30 和SRTM-3，結果分析中將LiDAR DEM 作為真值，水城縣的Li-DAR DEM、AW3D30 和SRTM-3 見圖3，隨著DEM 分辨率的提高，地形反映得更加真實細致，可輔助隱蔽性地質災害的識別[26]。

表2 LiDAR DEM、AW3D30和SRTM-3參數對比

圖3 水城縣的LiDAR DEM、AW3D30和SRTM-3

3 D-InSAR與TSI滑坡探測的對比

本文在興仁縣境內選取4 個具有代表性的區(qū)域，對D-InSAR和StaMPS-SBAS時序分析的滑坡探測與監(jiān)測結果進行對比分析，以證明TSI 在山區(qū)形變范圍圈定、大氣延遲削弱、微小形變探測和長時間序列特征分析方面的優(yōu)勢。

1）TSI 圈定形變范圍分析。P1 區(qū)域位于興仁縣潘家莊鎮(zhèn)褚皮田村。由圖4 可知，D-InSAR 探測的形變范圍與TSI 基本一致，但D-InSAR 獲取形變范圍明顯更小，只能獲取兩次成像期間的累積形變，而TSI 能準確、全面地給出長時序形變的空間范圍。

2）TSI抗誤差干擾分析。P2區(qū)域位于興仁縣潘家莊鎮(zhèn)扯尼姑村。圖5a中紅色框南北區(qū)域分別出現了沉降和抬升的形變信號，而圖5b僅探測到南部的沉降信號，北部未見明顯抬升。對比多對干涉圖發(fā)現，該信號僅出現在與2019-01-20 影像干涉的相位圖中；且根據距P2 區(qū)域18 km 的興仁縣氣象站歷史數據顯示，2019-01-20前后出現間歇性雨雪天氣，因此推斷北部出現的異常相位是由氣象條件差異引起的，而TSI 能削弱該影響。

圖5 P2區(qū)域D-InSAR與TSI結果

3）TSI 探測微小形變分析。P3 區(qū)域位于P1 和P2 之間，由圖6b 可知，該區(qū)域發(fā)生了形變速率為-17 mm/a 的微小形變，D-InSAR 方法難以有效探測出這樣的微小形變；而TSI 則能滿足微小形變的探測與分析。

圖6 P3區(qū)域D-InSAR與TSI結果

4）TSI長時間形變特征分析。P4區(qū)域位于興仁縣城南街道保駒村，由圖7b可知，該區(qū)域在2019年1月前基本未發(fā)生形變，此后開始滑動，LOS向形變速率達到了-83 mm/a，TSI 能反映隱患點的形變序列特征，為災害的穩(wěn)定性預警提供重要的指示依據。

圖7 P4區(qū)域D-InSAR與TSI結果

4 DEM對TSI分析的影響

本文以LiDAR DEM 為地表高程真值，對比分析了3 種DEM 的高程及其對相干點選取、DEM 誤差估計、形變速率、形變序列的影響。

4.1 多源DEM數據高程對比

AW3D30 和SRTM-3 的高程誤差分布和統計直方圖見圖8，紅色、藍色分別表示LiDAR DEM高于、低于其他兩種DEM，可以看出，SRTM-3 的高程誤差明顯大于AW3D30，誤差標準差達到24.71 m；在坡度較陡的區(qū)域地形誤差均較大（藍框區(qū)域），符合AW3D30、SRTM-3 在坡度大區(qū)域高程精度較低的結論[27,29]；部分區(qū)域展現出了無明顯分布規(guī)律的小范圍誤差（紅框區(qū)域），可能為地表形態(tài)發(fā)生了較大變化。

圖8 AW3D30和SRTM-3高程誤差分布和統計直方圖

4.2 DEM對相干點選取的影響

由于StaMPS 通過分析相位穩(wěn)定性來選取相干點，外部DEM可通過影響地形殘余相位的方式間接影響相干點選取的數量和質量。在相同的參數條件下，Li-DAR DEM、AW3D30 和SRTM-3 選取的相干點數量分別為268 425個、268 998個、274 419個，平均相干性均為0.31，總體差異不大。3 種DEM 選取的相干點分布的局部細節(jié)見圖9，可以看出，LiDAR DEM 與AW3D30 基本一致，而與SRTM-3 差異明顯；道路上分布有相干點，LiDAR DEM、AW3D30選取的相干點與道路重疊度較高，而SRTM-3 在地形變化區(qū)域選取的相干點與道路明顯偏移，地理編碼精度較低?？傮w而言，不同DEM對選點的數量和相干性水平沒有顯著影響，但在局部區(qū)域，低精度DEM獲取的相干點可能存在地理編碼不準確的問題。

圖9 采用3種DEM時的相干點局部分布

4.3 DEM誤差估計

采用StaMPS-SBAS 估計的DEM 誤差見圖10，對比圖8可知，估計的DEM地形誤差與實際地形差異的局部特征基本一致，如圖8 中的地形差異較大區(qū)域（藍框）和地形變化區(qū)域（紅框）在圖10 中也有相應特征；但StaMPS-SBAS 估計的DEM 誤差明顯低于實際地形誤差，尤其是低精度的SRTM-3。

圖10 采用AW3D30和SRTM-3的地形誤差估計結果

StaMPS-SBAS估計的地形誤差與真實地形誤差的關系見圖11，可以看出，StaMPS-SBAS 估計的DEM誤差與真實誤差基本呈正相關關系，但明顯低于真實地形誤差，DEM真實誤差為50 m時，AW3D30估值約為22 m，SRTM-3 估值僅為8 m，說明StaMPS-SBAS方法能去除一定的DEM誤差，但會低估其真實值，該問題在低精度DEM中表現得更嚴重。

圖11 AW3D30和SRTM-3估計的地形誤差與真值的關系

4.4 DEM對形變速率的影響

采用LiDAR DEM、AW3D30、SRTM-3進行TSI分析獲取的LOS 向形變速率圖見圖12。以LiDAR DEM獲取的平均形變速率為準，以10 m為半徑尋找3組結果的同名點并作差，得到形變速率差異圖（圖13）：AW3D30與LiDAR DEM獲取的形變速率差異較小，平均絕對誤差約為3.55 mm/a；SRTM-3獲取的形變速率差異集中在±20 mm/a；對比圖8 發(fā)現，形變速率差異與DEM 誤差的空間分布不相關，說明StaMPS-SBAS方法能在一定程度上反演出DEM誤差的整體趨勢與局部特征，但低估了DEM誤差的真實值，從而引起形變速率的估計偏差。

圖12 采用3種DEM的LOS向形變速率圖

圖13 LiDAR DEM與AW3D30、SRTM-3形變速率差異和統計直方圖

4.5 DEM對形變序列的影響

對研究區(qū)內的兩個較大的形變區(qū)和一個穩(wěn)定區(qū)域進行形變序列對比分析，3個區(qū)域的TSI結果見圖14~16。Q1點位于水城縣發(fā)耳鎮(zhèn)尖山營發(fā)耳煤礦，由于煤層開采形成采空區(qū)，坡體的應力場發(fā)生改變，巖體產生變形塌陷，光學影像上可以看出該區(qū)域存在多條地裂隙，坡體上有多條鏟溝，坡腳存在大量的堆積碎屑，滑坡面積較大。由于地表發(fā)生了劇烈變化，在主坡體上并未選取到相干點，僅在上部未滑坡區(qū)域和下部碎屑提取到部分相干點。該區(qū)域長期處于不穩(wěn)定變形狀態(tài)，2018年5 月后坡體形變速率發(fā)生多次改變，AW3D30、SRTM-3 與LiDAR DEM 的形變序列也開始出現差異。Q2位于水城縣都格鎮(zhèn)新盤村保興煤礦，形變區(qū)域位于山體坡腳，不具備滑坡、崩塌等地質災害發(fā)生的地表形態(tài)，因此推斷該形變是由采礦活動引起的。然而地下礦產開采往往導致地表向下塌陷，圖15中顯示該區(qū)域地表出現了抬升形變，經查詢，該礦場建設有瓦斯氣發(fā)電站，瓦斯壓力增大可能引起地表的抬升形變。該區(qū)域以45 mm/a 的形變速率持續(xù)發(fā)生形變，LiDAR DEM 與AW3D30、SRTM-3 的形變序列基本一致。Q3 點位于水城縣雞場鎮(zhèn)旗幟村。TSI結果顯示該點未發(fā)生形變，其累積形變序列僅有微弱的周期性變化，LiDAR DEM 與AW3D30、SRTM-3 的形變序列差異較小。

圖14 Q1點TSI結果

圖15 Q2點TSI結果

圖16 Q3點TSI結果

AW3D30、SRTM-3 與LiDAR DEM 在3 個區(qū)域形變序列差異見表3，在時間序列上SRTM-3 與LiDAR DEM的差異更大，形變速率未發(fā)生變化的Q2和Q3點的形變序列差異明顯小于發(fā)生了形變速率變化的Q1點。

表3 AW3D30、SRTM-3與LiDAR DEM的形變序列差異的平均絕對誤差/mm

綜上所述，DEM誤差對形變序列估計的影響主要集中在形變速率發(fā)生變化的時刻：當地表形變速率穩(wěn)定時，不同DEM獲取的形變序列基本一致；當地表形變速率發(fā)生變化時，不同DEM獲取的形變序列會出現階躍式差異，并維持這一形變差異直至速率再次發(fā)生變化。

5 結 語

本文在貴州省興仁縣開展了D-InSAR 和TSI 滑坡探測對比分析，證明TSI 具備抗誤差干擾能力，能探測微小形變，并提供直觀的形變序列演化趨勢，為山區(qū)地質災害隱患排查工作提供了重要的科學技術參考。同時，近年來大量雷達衛(wèi)星數據的開放，為TSI分析提供了基礎數據支持，充分滿足大范圍地質災害普查工作的需求。TSI 數據處理流程中需引入外部DEM 模擬和去除地形相位，目前常用的公開DEM 精度有限，可能對InSAR 形變探測與分析產生不利影響。本文采用LiDAR DEM、AW3D30 和SRTM-3 等外部DEM 在貴州省水城縣開展了StaMPS-SBAS 處理，通過對比分析相干點選取情況、DEM誤差估計、形變速率和形變序列發(fā)現，StaMPS-SBAS方法無法完全剔除DEM誤差，從而導致形變速率、形變序列估計存在一定偏差。因此，在開展精細化InSAR 形變分析時，應盡可能選取高精度的外部DEM數據，以保證正確估計地質災害隱患形變特征。