基于軌道精煉控制點精選的極艱險區(qū)域時序InSAR地表形變監(jiān)測

潘建平，鄧福江，徐正宣，向淇文，涂文麗，付占寶

（1.重慶交通大學土木工程學院，重慶 400074；2.中鐵二院工程集團有限責任公司地勘院，四川 成都 610031）

0 引言

合成孔徑雷達干涉測量（Interferometry Synthetic Aperture Radar，InSAR）作為近年來新發(fā)展起來的一種地表監(jiān)測技術(shù)，具有全天時、全天候、監(jiān)測精度高、監(jiān)測范圍廣等優(yōu)點[1]，由其發(fā)展而來的合成孔徑雷達差分干涉測量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）通過對同一研究區(qū)的兩景SAR影像進行差分干涉，獲取精度為厘米級甚至毫米級的地表形變信息，極大得提升了InSAR技術(shù)在地表形變監(jiān)測領(lǐng)域的應用前景[2]。而以D-InSAR為基礎(chǔ)發(fā)展起來的時序InSAR技術(shù)又解決了長期監(jiān)測中時間失相干、空間失相干、大氣延遲相位等因素的影響，成為了一種有效探測長時間緩慢地表形變的監(jiān)測技術(shù)[3?5]。

常見的時序InSAR處理技術(shù)有永久散射體合成孔徑雷達干涉測量（Persistent Scatters InSAR,PS-InSAR）技術(shù)[6]和小基線集合成孔徑雷達干涉測量（Small Baseline Subset InSAR，SBAS-InSAR）技術(shù)[7]。PS-InSAR技術(shù)利用時間序列影像中的永久性散射體（Persistent Scatters，PS）點來進行時序分析，可以有效地減少時間失相干、空間失相干以及大氣延遲相位對地表形變監(jiān)測結(jié)果的影響。2017年白澤朝等[8]利用PS-InSAR技術(shù)成功獲取了天津市2015年至2016年地面沉降速率。SBASInSAR技術(shù)通過設(shè)置時空基線閾值，從而利用短時間內(nèi)時序影像的空間相干性，有效減小了時間失相干、空間失相干以及大氣延遲相位對地表形變監(jiān)測結(jié)果的影響，同時還增加了時間采樣率[9?10]。2019年劉曉杰等[11]利用SBAS-InSAR技術(shù)和Sentinel-1A數(shù)據(jù)得到了廈門新機場海洋填海區(qū)地面沉降特征。潘光永等[12]利用SBAS-InSAR技術(shù)對濟南井田礦區(qū)進行了地面沉降監(jiān)測。PS-InSAR技術(shù)和SBAS-InSAR技術(shù)各有優(yōu)點，因此兩種技術(shù)在不同研究區(qū)域下的結(jié)合與交叉應用是當前的研究熱點。2019年王舜瑤等[13]利用一種顧及永久性散射體的SBAS-InSAR方法獲取了鄭州市區(qū)的地表形變信息，發(fā)現(xiàn)形變結(jié)果與水準測量吻合較好，說明了該方法的可靠性。

文章以折多山區(qū)域為研究對象。極艱險區(qū)域永久性散射體分布稀疏，PS-InSAR技術(shù)無法獲得較好的地表形變結(jié)果，同時該區(qū)域高差極大，地形復雜，不利于SBAS-InSAR技術(shù)選取穩(wěn)定的地面控制點（GCP），從而降低了地表形變監(jiān)測的精度。為此，采用了一種改進的SBAS-InSAR技術(shù)，通過篩選合適的PS點作為GCP點引入到SBAS-InSAR技術(shù)處理流程，從而有效提高該區(qū)域地表形變監(jiān)測精度。

1 技術(shù)方法

文章采用SBAS-InSAR技術(shù)處理流程為主體，引入顧及研究區(qū)實際概況的GCP點篩選研究，提高SBASInSAR技術(shù)地表形變監(jiān)測精度，技術(shù)流程如圖1所示。改進的SBAS-InSAR技術(shù)原理重點涉及GCP點篩選和SBAS-InSAR技術(shù)原理兩個部分。

圖1 改進的SBAS-InSAR技術(shù)的基本流程圖Fig.1 Basic flow chart of improved SBAS-InSAR Technology

1.1 GCP點篩選

GCP點被引入SBAS-InSAR軌道精煉步驟來進行殘余相位估計和重去平處理，從而提高地表形變結(jié)果的精度。GCP點應當處于地形平緩，沒有相位躍變且遠離形變區(qū)域的位置，根據(jù)GCP點選取特點，可以采用由PS-InSAR技術(shù)處理后提取的PS點作為GCP候選點。

PS點是長時間序列影像中具有較高后向散射特性以及穩(wěn)定性的點，其在不同時期的影像干涉對中均表現(xiàn)出高一致性與高相干性。振幅離差指數(shù)法[14]、相干系數(shù)閾值法[15]均是常見的PS點選取方法。

（1）振幅離差指數(shù)法

PS點較高的后向散射特性和穩(wěn)定性體現(xiàn)在其回波的相位信息的長時間序列上具有一定的統(tǒng)計特性。根據(jù)這一特性，可以用振幅離差指數(shù)DA來 定量表示，具體表達示如下：

式中：DA——振幅離差指數(shù)；

σA——振幅標準差；

μA——振幅均值。

當某像素點的振幅離差指數(shù)DA的大小在0.25～0.4時，該點可以被視為PS點，但為了保證PS點的質(zhì)量和最終形變監(jiān)測結(jié)果的可靠性，至少需要25幅SAR影像參與計算。

（2）相干系數(shù)閾值法

相干系數(shù)是衡量干涉影像像對的干涉質(zhì)量的重要指標，它主要用于描述干涉對中主、副影像同一區(qū)域的相似程度。相干系數(shù)數(shù)值分布區(qū)間為[0,1]，0表示完全不相干，1表示完全一致。根據(jù)PS點的特性，在相干性好的區(qū)域?qū)南喔上禂?shù)也相對較高[16]，因此可以采用相干系數(shù)為PS點選取的判斷閾值，其表達式為：

式中：m、n—主、輔影像內(nèi)需要計算相干性的數(shù)據(jù)塊大?。?/p>

i、j——數(shù)據(jù)塊內(nèi)的行列號；

M(i,j)、S(i,j)— 主、輔影像數(shù)據(jù)在 (i,j)位置處的復數(shù)值；

*——復數(shù)的共軛算子；

|～|2——數(shù)據(jù)的二階范數(shù)。

通過以任一像元為中心的i×j窗口大小進行計算，即可得到該像元的相關(guān)系數(shù)值γ 。γ值越高，像元越穩(wěn)定，受噪聲影響較小，干涉相位的質(zhì)量越高，高于設(shè)定的相干系數(shù)閾值的像元點，即可判定為PS點。

實際處理中，聯(lián)合多種PS點選取方法，通過多層閾值過濾的方式，減少各PS點選取方法考慮特性單一帶來的影響，能夠有助于提高PS點判斷的準確性。選取出的PS點作為GCP候選點引入后續(xù)得處理中。

高山區(qū)域地形起伏較大，部分裸露巖石位于非平坦區(qū)域，此類因雷達散射特性穩(wěn)定的永久性散射體，依然有被提取并作為GCP候選點的可能。然而此類GCP候選點不符合GCP點的要求，若作為GCP點參與后續(xù)流程解算，會降低形變結(jié)果精度。因此可以根據(jù)光學影像來去除此類GCP候選點，最終得到GCP點并引入SBASInSAR技術(shù)處理流程。

1.2 SBAS-InSAR技術(shù)原理

覆蓋研究區(qū)的所有時間序列的SAR影像，假設(shè)有N+1幅SAR影像，其時間序列為：

通過選取其中一幅SAR影像作為超級主影像，將其他N幅SAR影像都配準到該超級主影像的雷達坐標系上來，像按照一定的時間、空間基線閾值來組合SAR影像，得到M幅干涉對，并對其進行差分干涉處理。假設(shè)第i景差分干涉圖由t1，t2(t1

式中：i∈[1,2,···,M]；

dt— —t時刻相對于t0時刻的雷達視線向累積形變量；

式中AV為[M×N]的矩陣，每一行對應每個干涉對。當干涉對屬于同一子基線集時，有r(A)=N，由最小二乘法可得：

式中δφ為解纏后干涉相位。然而，受時空基線的限制并非任意干涉對同屬于某一小基線集，A為秩虧陣，導致最小二乘法的解不唯一。因此，采用奇異值分解（singular value decomposition，SVD）來求解A的廣義逆矩陣，因而可得到最小二乘范數(shù)解：

式（7）中U為[M×N]的 正交矩陣且的對角元素為奇異值 σi(i=1,···,N)，將求相位解轉(zhuǎn)換為求相位變換速率，則可得：

式（8）中B為[M×N]的 矩陣，對B進行奇異值分解，可以解出各時間段內(nèi)相位變化速率v，最后計算得到最終平均形變速率以及形變速率時間序列。

2 實驗過程

2.1 研究區(qū)概況

研究區(qū)位于四川省甘孜自治州康定市西面的折多山區(qū)域，面積約為1 000 km2。研究區(qū)域地形高差大，溝壑密布，山嶺縱橫，活動斷裂帶構(gòu)造大量發(fā)育，屬于極艱險區(qū)域。研究區(qū)地理位置如圖2所示。該區(qū)域復雜的地理環(huán)境，使得常規(guī)工程測量施測難以開展，InSAR技術(shù)所具備的全天時、全天候、大范圍、高精度監(jiān)測的優(yōu)勢能在該區(qū)域得到較好的體現(xiàn)。

圖2 研究區(qū)地理位置圖Fig.2 Geographical location of the study area

2.2 數(shù)據(jù)源

本文采用由歐空局提供的Sentinel-1A雷達衛(wèi)星影像升軌數(shù)據(jù)，成像模式為干涉寬視場（IW）模式，極化方式為VV同極化，重訪周期為12天，時間跨度為2017年1月至2019年7月，共計76景。衛(wèi)星精密軌道數(shù)據(jù)采用歐空局提供的POD Precise Obrit Ephemerides（POD精密定軌星歷數(shù)據(jù)），共計76個。DEM數(shù)據(jù)采用美國航空航天局（NASA）和國防部國家測繪局（NIMA）聯(lián)合測量的SRTM1，分辨率為30 m。

2.3 數(shù)據(jù)處理

實驗數(shù)據(jù)處理采用了SARscape軟件，按照改進后的SBAS-InSAR處理流程，得到研究區(qū)時間序列形變結(jié)果和地表平均形變速率，并采用GMT軟件繪制成年均地表形變速率圖。

（1）數(shù)據(jù)預處理。為減少計算冗余和處理時間，需要按研究區(qū)范圍進行預先的裁剪。

（2）連接圖生成。對預處理后的研究區(qū)影像進行基線估算，設(shè)置時空基線閾值，并生成連接圖。本文空間基線設(shè)置為臨界基線的3%，時間基線設(shè)置為48天，最終形成了289對干涉對，時空基線分布如圖3、圖4所示。

圖3 時間基線連接圖Fig.3 Time baseline connection diagram

圖4 空間基線連接圖Fig.4 Spatial baseline connection diagram

（3）干涉處理。對所有的干涉像對進行干涉處理，并通過查看差分干涉結(jié)果圖，剔除了49對相干性差、解纏不理想的影像干涉對，為軌道精煉和重去平，以及SBAS反演估算做好數(shù)據(jù)準備。

（4）GCP點篩選。首先通過PS-InSAR技術(shù)獲取PS點來作為GCP候選點，根據(jù)GCP點的特性，綜合考慮振幅離差，相干性以及形變速率，來獲取PS點信息。PS點振幅離差閾值在0.25至0.4之間；相干性閾值不能設(shè)置過高，否則處于高山山地的極艱險區(qū)域很難獲取一定數(shù)量的PS點，需要反復試驗和調(diào)整，使獲取的PS點數(shù)量滿足GCP點數(shù)量要求；PS點必須是穩(wěn)定點，形變速率需要小于1 mm/a。因此，本文設(shè)置振幅離差閾值Dv為3.5，相干性閾值Tγ為 0.87，形變速率閾值Tv為±1 mm/a，共獲取了51個PS點，也就是51個GCP候選點。然后結(jié)合光學影像查看51個GCP候選點的位置分布，發(fā)現(xiàn)其中9個點沒有位于地勢較緩的平坦區(qū)域，例如圖5所示GCP點處于高山半坡，將此類GCP候選點剔除后最終得到42個穩(wěn)定的GCP點，其結(jié)果與人工選取GCP點的結(jié)果對比如圖6所示。

圖5 結(jié)合Google earth光學影像篩選剔除的GCP點Fig.5 GCP points filtered out by Google Earth optical image

圖6 GCP點選取結(jié)果對比（左為人工選取，右為本文方法選?。〧ig.6 Comparison of GCP point selection results (left is manual selection, right is method selection in this paper)

（5）軌道精煉與重去平。將篩選得到的GCP點轉(zhuǎn)換到SAR影像雷達斜坐標系中，采用3次軌道精煉多項式估算軌道誤差和相位偏移量，消除斜坡相位，最后基于控制點對數(shù)據(jù)進行重去平。

（6）SBAS反演估算。先是通過第一次反演估算形變速率和殘余地形，并設(shè)置相應的小波分解等級去除殘余地形的影響；再是通過第二次反演，進行時間域高通濾波和空間域低通濾波計算，去除大氣相位影響。最終得到時間序列形變結(jié)果。

（7）地理編碼。將SBAS反演估算得到的形變結(jié)果轉(zhuǎn)換到參考DEM坐標系下，最后通過GTM軟件繪制成年均地表形變速率圖(圖7)。

圖7 折多山區(qū)域年均地表形變速率圖Fig.7 Annual surface deformation rate of zheduoshan area

3 結(jié)果與分析

對比PS-InSAR技術(shù)、SBAS-InSAR技術(shù)和改進的SBAS-InSAR技術(shù)的地表形變監(jiān)測結(jié)果，可以發(fā)現(xiàn)研究區(qū)域基于SBAS-InSAR技術(shù)提取的地表形變點密度明顯高于PS-InSAR技術(shù)獲取的地表形變點密度，SBASInSAR技術(shù)監(jiān)測結(jié)果中該區(qū)域中疑似滑移區(qū)域體現(xiàn)得更加明顯和清晰。這是由于位于高山山地地帶的極艱險區(qū)域永久性散體較為稀少，PS-InSAR獲取的PS點信息也相對較少。區(qū)域內(nèi)山體海拔較高，多存在冰雪覆蓋區(qū)域，受氣溫影響，地表土壤呈現(xiàn)出較為明顯的季節(jié)性反復凍融，地表形變量往往較大，在采用單一主影像的PS-InSAR技術(shù)的長時間基線干涉圖中，這些區(qū)域常出現(xiàn)失相干現(xiàn)象，難以獲得一定數(shù)量的PS點。通過對比分析可知，SBAS-InSAR在該研究區(qū)域的監(jiān)測效果更好。改進的SBAS-InSAR技術(shù)的地表形變結(jié)果很好的保留了SBAS-InSAR技術(shù)的優(yōu)勢，得到了同樣密度較高的地表形變點，整體的形變速率結(jié)果與SBAS-InSAR技術(shù)處理結(jié)果基本一致（圖8）。

圖8 不同方法的折多山區(qū)域年均地表形變速率圖Fig.8 Annual surface deformation rate of zheduoshan area with different methods

進一步的分析改進SBAS-InSAR技術(shù)與SBAS-InSAR技術(shù)的處理結(jié)果，對比兩者的形變速率、形變速率精度，如表1、表2所示。

表1 兩種時序InSAR技術(shù)的形變速率統(tǒng)計Table 1 Deformation rate statistics of two time series InSAR techniques

表2 兩種時序InSAR技術(shù)的形變速率精度統(tǒng)計Table 2 Deformation rate accuracy statistics of two time series InSAR techniques

通過對比改進前后的監(jiān)測結(jié)果，該區(qū)域的平均形變速率由?1.208 922變?yōu)?1.099 689 1，標準差從9.224 748變?yōu)?.172 168。該區(qū)域形變結(jié)果表明，人工選取的GCP點處于輕微沉降區(qū)域，而GCP點代表穩(wěn)定的位置，這使得使用人工選取GCP點的監(jiān)測結(jié)果中形變點沉降速率相對偏大；而使用精選后GCP點的監(jiān)測結(jié)果更接近形變點的真實沉降速率。并且，在形變速率精度量級上，改進后的形變速率精度平均值從改進前的3.851 148變?yōu)?.091 360，標準差從1.501 882變?yōu)?.279 466。因此結(jié)果表明，改進后的SBAS-InSAR技術(shù)的監(jiān)測結(jié)果在形變速率、形變速率精度上均有較好的改善。高山山地地帶的地理環(huán)境條件使得傳統(tǒng)SBAS-InSAR技術(shù)中人工選取GCP點較為困難，而利用篩選得到的永久性散體點作為GCP點來代替人工選取GCP點，有助于提高地表形變監(jiān)測精度。這說明基于軌道精煉控制點精選的SBAS-InSAR技術(shù)在極艱險區(qū)域的地表形變監(jiān)測中具有更好的適用性。

4 結(jié)論

本文對比了PS-InSAR技術(shù)和SBAS-InSAR技術(shù)在折多山區(qū)域的地表形變監(jiān)測結(jié)果，發(fā)現(xiàn)SBAS-InSAR技術(shù)在該研究區(qū)域獲得的地表形變點密度更高，整體效果更好；又考慮到極艱險區(qū)域的地理條件，設(shè)計了GCP點精選方案來代替GCP人工選點，對SBAS-InSAR技術(shù)進行了改進。通過分析，本文得到以下結(jié)論：

（1）在極艱險區(qū)域的地表形變監(jiān)測中，SBAS-InSAR技術(shù)能夠獲取較高地表形變點密度的監(jiān)測結(jié)果，而PSInSAR技術(shù)的應用效果不佳。SBAS-InSAR技術(shù)在極艱險區(qū)域比PS-InSAR技術(shù)更具有應用優(yōu)勢。

（2）本文提出的基于軌道精煉控制點精選的SBASInSAR技術(shù)在極艱險區(qū)域的應用中，保留了SBASInSAR技術(shù)在極艱險區(qū)域地表形變監(jiān)測中的優(yōu)勢，取得較高密度地表形變監(jiān)測結(jié)果的同時，避免了GCP人工選點引起的誤差，提高了地表形變監(jiān)測精度，具有良好的應用價值。