基于PS-InSAR方法的個(gè)舊市地表形變監(jiān)測

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

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

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

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

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

圖1 研究區(qū)范圍

2 實(shí)驗(yàn)方法與步驟

2.1 最佳主影像選取

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

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

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

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

2.2 差分干涉圖生成

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

圖2 干涉影像對的時(shí)空基線圖

2.3 永久散射（PS）點(diǎn)選取

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

式中，σA、μA分別為SAR影像的振幅標(biāo)準(zhǔn)差和振幅平均值。實(shí)驗(yàn)中設(shè)定DA的閾值為0.4，即當(dāng)該位置振幅離差指數(shù)小于0.4時(shí)，則認(rèn)定為PS點(diǎn)；反之，則為非PS點(diǎn)。

2.4 PS點(diǎn)的時(shí)序形變求解

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

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

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

3 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

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

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

4 結(jié) 語

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