馬 麗，羅建洪，王 波，李蓮娥，趙彥文

（1.師宗縣自然資源局，云南 師宗 655700；2.云南省測繪地理信息科技發(fā)展有限公司，云南 昆明650034）

地面沉降通常是指由地下松散地層的固結和壓縮引起的地面逐漸沉降或迅速下沉，是造成建筑物、道路、橋梁基礎設施嚴重損壞的重大區(qū)域地質災害之一。通過工程技術手段準確有效地監(jiān)測地表沉降對于地質災害的預防和治理具有重要意義。以水準測量、GNSS測量為代表的傳統(tǒng)監(jiān)測手段難以獲取大面積的沉降監(jiān)測數(shù)據(jù)，且這些方法均需投入大量的人力物力[1]。為了能高效準確地對地表形變進行監(jiān)測，同時降低監(jiān)測成本，急需采用新方法來彌補傳統(tǒng)測量方法的不足。隨著新一代SAR衛(wèi)星Sentinel-1A的升空，大大降低了SAR影像的獲取成本，為解決上述問題提供了新的思路[2]。差分干涉合成孔徑雷達（DInSAR）技術作為一種重要的、先進的測量地表變形的大地測量技術，能有效彌補傳統(tǒng)測量方法的不足，獲取大范圍的地表形變數(shù)據(jù)；但仍受時間去相關、空間去相關、大氣延遲等誤差的影響，導致測量結果的誤差。針對這些缺點，F(xiàn)erretti A[3]等提出了永久散射體干涉測量技術（PSInSAR），能有效避免時間去相關、空間去相關、信號大氣延遲等因素引起的測量誤差，獲取大范圍高精度的地表形變數(shù)據(jù)，為防治災害提供重要依據(jù)。本文以個舊市為研究區(qū)，利用PS-InSAR技術處理2018-03-22—2019-01-16的23景Sentinel-1A衛(wèi)星影像，從而監(jiān)測研究區(qū)域的地表形變，并對形變結果進行分析。

1 研究區(qū)概況與數(shù)據(jù)來源

個舊市位于云南省東南部，云貴高原的南端，坐標范圍為23e21′～23.36′N、102e54′～103e25′E。個舊市是世界主要的錫礦生產基地，其錫礦開采歷史可追溯至2 000 a前的漢朝[4]。由于常年開采，個舊市部分地區(qū)的地表產生了非常嚴重的沉降，給道路、房屋、橋梁、水工等建筑物帶來了潛在的風險，因此及時發(fā)現(xiàn)并持續(xù)監(jiān)測地表沉降嚴重區(qū)域對于保障該地區(qū)的工業(yè)生產安全和人民群眾的生命財產安全具有重要意義。本文選取23景Sentinel-1A衛(wèi)星影像進行處理，Sentinel-1A衛(wèi)星的重返周期為12 d，工作頻率為C波段（中心頻率為 5.405 GHz）[5]。研究區(qū)域范圍如圖1中黃色方框所示，本文選用的Sentinel-1A衛(wèi)星影像完整覆蓋了研究區(qū)域。

圖1 研究區(qū)范圍

2 實驗方法與步驟

2.1 最佳主影像選取

影響影像干涉質量的主要因素包括干涉影像的時間基線、空間垂直基線、多普勒質心頻率差和熱噪聲失相關。主影像的選擇應使干涉影像的垂直基線的差異度盡可能地降低，從而提高干涉影像的總體質量。為了選擇最佳主影像，本文采用式（1）、式（2）來評估干涉影像的總體相關性[6]。

式中，為影像k和影像m的空間垂直基線；Tk,m為影像k和影像m的時間基線；fk,m為影像k和影像m的多普勒質心頻率差；Bc、Tc、fc為上述3個變量的臨界值，在實驗中分別設置為1 200 m、5 a、1 380 Hz。

通過歐洲空間局提供的SANP軟件[7]快速計算得到所有干涉影像的時間基線、空間垂直基線和多普勒質心頻率差，當γm達到最大值時，此時選取的主影像為最佳主影像。通過計算，本文以2018-07-20獲取的影像為最佳主影像，其相關參數(shù)信息如表1所示。

表1 影像相關參數(shù)信息表

2.2 差分干涉圖生成

以2018-07-20獲取的影像為最佳主影像，采用增強譜分集方法將其余22景影像（從影像）配準至最佳主影像上。該方法能有效提升配準精度，使配準誤差小于1/1 000個像元[8]。配準后將主影像與從影像兩兩進行差分干涉處理，得到22景差分干涉影像，處理時使用的DEM為美國地質調查局提供的SRTM數(shù)據(jù)。干涉影像對的時空基線圖如圖2所示，藍色實線連接主影像和從影像，可以看出，實驗采用的影像空間基線在100 m以內。

圖2 干涉影像對的時空基線圖

2.3 永久散射（PS）點選取

本文采用振幅離差法識別PS點。振幅離差法通過評估影像像元的信噪比來確定PS點，相位的信噪比可用振幅離差指數(shù)來衡量[3]，即

式中，σA、μA分別為SAR影像的振幅標準差和振幅平均值。實驗中設定DA的閾值為0.4，即當該位置振幅離差指數(shù)小于0.4時，則認定為PS點；反之，則為非PS點。

2.4 PS點的時序形變求解

ψx,i為第i張干涉影像上的第x個PS點的干涉相位，主要由5個部分組成[6]，即

式中，W{}為相位纏繞算子；φD,x,i為由地面形變引起的相位；φA,x,i為大氣延遲誤差相位；ΔφS,x,i為衛(wèi)星軌道誤差相位；Δφθ,x,i為殘余地形相位；ΔφN,x,i為由配準誤差、熱噪聲失相關等引起的噪聲相位。

通過DEM相關相位誤差校正、三維相位解纏、去除大氣誤差、去除空間相關相位誤差等步驟，可獲得第i張干涉影像獲取時刻的地表形變相位φD,x,i；再通過距離相位轉換公式即可得到各PS點的時間形變序列；最后將相位轉換為形變速率，得到年平均形變速率圖。

3 實驗結果分析

PS-InSAR方法在D-InSAR技術處理得到的22景差分干涉圖的基礎上獲得了PS點的形變速率圖和時間形變序列圖。研究區(qū)的年平均形變速率圖如圖3所示，獲取的形變速率方向為衛(wèi)星的視線方向（LOS），可以看出，在城鎮(zhèn)區(qū)域、地表裸露無植被覆蓋的礦山開采區(qū)存在大量的PS點，研究區(qū)的最大沉降形變速率為-49.89 mm/a，整個區(qū)域的沉降速率不均勻，錫礦開采區(qū)域地表沉降嚴重，城鎮(zhèn)區(qū)域地表形變量較小。選取A區(qū)域（位于個舊市城區(qū)）、B區(qū)域（位于個舊市老廠鎮(zhèn)錫礦礦區(qū)）做進一步分析發(fā)現(xiàn)，A區(qū)域地質構造穩(wěn)定，地表沉降量很小；B區(qū)域由于錫礦開采，地表沉降嚴重，該區(qū)域主要以開采金屬礦產為主，其形變量級小于煤礦開采所導致的地表沉降，在這類區(qū)域很難出現(xiàn)形變量級過大而導致的干涉影像失相干。

A區(qū)域的平均形變速率圖如圖4所示，可以看出，該區(qū)域的地表沉降較小，大部分PS點的年均形變速率為-15～15 mm/a。以T1、T2、T3點為例，本文分別繪制了各點的時間形變序列圖，如圖5所示，可以看出，位于城區(qū)的T1、T2、T3點的地表形變量變化趨勢較平緩，T1、T2、T3點的LOS向形變速率分別為1.5 mm/a、 2.58 mm/a和1.16 mm/a，說明A區(qū)域的地質條件穩(wěn)固。

圖4 A區(qū)域的平均形變速率圖

圖5 A區(qū)域T1、T2、T3點時間形變序列圖

B區(qū)域的平均形變速率圖如圖6所示，可以看出，受到長期礦山采掘的影響，該區(qū)域的地表沉降較嚴重，部分區(qū)域的地表沉降速率超過了-40 mm/a。以P1、P2、P3點為例進行形變序列分析，本文分別繪制了各點 的時間形變序列圖，如圖7所示，可以看出，由于 錫礦開采，這3點的地表沉降趨勢明顯，均產生了較嚴重的沉降；P1、P2、P3點的LOS向形變速率分別為 -48.85 mm/a、-46.86 mm/a和-44.79 mm/a；P1點在2018-12-23達到了沉降量最大值，P2、P3點在 2019-01-04達到了沉降量最大值，這3點在2018-08以后形變速率呈增加趨勢，需進一步加強觀測。

圖6 B區(qū)域的平均形變速率圖

圖7 B區(qū)域P1、P2、P3時間形變序列圖

4 結 語

本文利用PS-InSAR方法對覆蓋個舊市的23景 Sentinel-1A衛(wèi)星數(shù)據(jù)進行了處理，獲得了該地區(qū) 2018-03-22－2019-01-16的地表形變信息。研究結果表明，該地區(qū)的沉降速率在-49.89～25.52 mm/a之間，沉降嚴重的區(qū)域與錫礦開采區(qū)基本一致，而個舊市城區(qū)地質條件穩(wěn)固、形變量較小。PS-InSAR方法能對該地區(qū)的地表進行監(jiān)測，有效減少了人力物力成本；但InSAR技術也存在一些缺陷，如植被茂密地區(qū)干涉信號失相干嚴重。今后還需結合GNSS、水準測量等多種技術手段，合理利用各種測量方法的優(yōu)勢，對該地區(qū)地表沉降進行有效監(jiān)測，以便能更好地為該地區(qū)的災害防治提供技術支撐。