基于SBAS技術的九寨溝震后地表形變監測

胡 波，吳 洋，魏德宏，楊 斌，余俊鵬，陳志謀

(1.廣東工業大學，廣東 廣州 510006；2.偉志股份公司，福建 晉江 362200)

2017-08-08T21：19于四川九寨溝(103.82°E，33.20°N)發生Ms7.0地震，震源深度約為20 m。地震發生后，九寨溝景區遭受到較大的破壞，本次地震造成了25人死亡，520余人受傷，同時7萬余房屋受到不同程度的損壞[1-2]。而地震之后一段時間內當地地形的變化成為災后重建的一項重要監測內容。

合成孔徑雷達干涉測量(Interferometry Synthetic Aperture Radar，InSAR)是一門測量地形及地表形變的技術。經過近30年的發展，InSAR技術已經成為大面積區域地面沉降監測的重要手段[3]。相比較常規測量技術而言，InSAR技術具有全天候、寬覆蓋、高精度、低成本的優勢[4]。但是，傳統的差分合成孔徑雷達干涉測量(Differential Interferometry SAR，D-InSAR)技術容易受到大氣擾動、時間失相關和空間失相關等因素的影響，其測量的精度受到一定限制，并且無法獲得監測區域地表形變的時間演化情況[5]。針對D-InSAR中存在的諸多問題，科研工作者相繼提出了時序InSAR技術。其中，具有代表性的時序InSAR技術有永久散射體干涉測量[6](Permanent Scatter Interferometry, PS)技術和小基線集[7-8](Small Baseline Subset, SBAS)技術。目前時序InSAR技術已廣泛應用于城市地下水開采[9]、地質災害[10-11](滑坡、火山、地震等)和城市地表[12-13]等形變測量當中。

本次實驗選取25景Sentinel衛星數據，使用小基線集合成孔徑雷達干涉測量技術對九寨溝地震區及其附近山區進行測量，得到該地區震后近1年內的時間序列形變，為探究該地區震后的沉降規律、災后重建及地質災害預防提供技術和數據支持。

1 研究區域及實驗數據

1.1 研究區域

本次研究的區域位于四川省阿壩州松潘縣北部，九寨溝地震震中附近。中心坐標為103.80°E，33.12°N，其覆蓋面積約為6 712 km2。研究區域包含九寨溝景區，區域內地形以山脈為主，如圖1所示。

圖1 研究區域覆蓋范圍圖

1.2 實驗數據

本次研究選用了25景Sentinel-1A的升軌數據(見表1)，采用C波段觀測。為提高Sentinel-1A影像數據處理精度，本研究使用衛星的精密軌道數據(Precise Orbit Ephemerides,POD)。實驗中所使用的DEM數據是由日本JAXA公司提供的ALOS數字表面模型“ALOS World 3D- 30 m”。干涉影像對的時間與空間基線關系如圖2所示。

表1 Sentinel-1A影像信息

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

2 SBAS技術原理

Berardina和Lanari等人于2002年首次提出了SBAS技術[7]，該技術以多幅影像作為主影像生成干涉對，并利用干涉圖中高相干點恢復研究區域的時間序列形變信息。同時，SBAS技術利用多視處理降低差分干涉圖的相位噪聲，并應用奇異值分解(Singular Value Decomposition, SVD)方法對相位進行求解得到形變相位速率，進而解得整個觀測時段的形變時間序列[14]。

SBAS能有效提高研究區域形變的時間和空間分辨率，其步驟如下：

1)針對傳感器、觀測條件、研究區域的不同設定空間基線閾值和時間基線閾值以及平均相干性閾值，可以排除空間和時間失相干的像對，同時生成M對小基線干涉對，干涉對數量應滿足：

(1)

式中：N+1為影像個數。在地形起伏較大的區域，需借助外部DEM數據去除地形起伏導致的地形相位影響，從而得到n幅差分干涉圖的形變相位；

2)對于第i景差分干涉圖中，方位向坐標A以及距離向坐標R的像元的干涉相位值可以表示為：

δφi(A,R)≈[d(TB,A,R)-d(TA,A,R)]+

(2)

3)利用得到的干涉圖進行二次篩選，將干涉圖中干涉條件差或是解纏結果中存在大面積整周跳變結果的干涉圖剔除；

4)采用選取多個相對穩定的地面點的方法求取干涉圖里的大氣相位以及噪聲項成分。在求取誤差項以后，使用高次多項式插值的方法，對研究區域中的大氣相位進行插值，達到恢復研究區域中每一景干涉圖所包含的大氣相位及噪聲的目的。

本次實驗采用三次多項式插值，為了去除DEM誤差相位和大氣延遲帶來的影響，可假設地表形變的低頻部分為：

(3)

(4)

Ax=Δφ.

(5)

5)在去掉殘余地形相位、大氣相位和相干噪聲后，可以將式(2)簡化為：

δφi(A,R)=d(TB,A,R)-d(TA,A,R).

(6)

用矩陣形式表達式(6)，可寫成：

Bφ=δφ.

(7)

式中含有M個觀測向量和N個未知數，計算結果時，可以使用最小二乘法解得結果，然而在差分干涉中，為了抑制時空基線過長帶來的去相干，常出現M幅干涉圖并不處于同一子集的情況，該情況下用最小二乘得到的結果并不唯一，此時使用SVD方法可以得到未知參數的最小范數解。

3 平均形變速率及其時間序列形變反演

圖3 九寨溝震區2017-08-011—2018-07-01平均速率圖

本研究利用SBAS技術得到的該地區整體垂直方向平均速率如圖3所示。圖3顯示該區域在2017-08-11—2018-07-01期間整體形變量較小，大部分區域形變速率保持在-50～50 mm/a。同時，研究區域的東北部地區下沉速率較大，其年形變速率達到-70～-50 mm/a。除此之外，在中部地區和南部地區出現3塊面積較小，但下沉速率較快的區域，它們的面積小于1 km2，但最快下沉速率已經超過120 mm/a。為了進一步探究研究區域整體形變隨時間的變化，選取3塊區域共11個采樣點進行分析(見圖3)，采樣區域信息及采樣點分布如表2所示。

表2 九寨溝采樣區域信息及采樣點分布信息

圖3中區域1占據了研究區域的大部分面積，該區域主要表現為較為穩定的形變，其累積形變量多在-50～50 mm之間。為探究該區域地表形變規律，本文提取4個采樣點，其中P1位于震中(103.82°E，33.20°N)附近，P2為該區域中下沉明顯的點，而P3和P4分別位于區域的西北部和西南部。從P1～P4的時間序列形變圖(見圖4)中可以看出，P3和P4的時序形變較為穩定，波動幅度較小。而P1點由于位于九寨溝地震震中附近，在地震之后仍伴有較大的起伏，其最大抬升量和最大下沉量的差值達到37 mm。而P2雖然表現為較為平穩的下沉，但其累積下沉量仍超過80 mm，達到82 mm。

圖4 區域1采樣點時間序列形變圖

區域2位于研究區東北部，九寨溝縣北部。該區域表現為整體下沉，且形變速率多在-40～-70 mm/a之間。在該區域中，本文選取3個點(P5～P7)作為采樣點進行時序分析。從圖5中可以看出，2017-08-11—2018-01-02，該地區處于快速下沉狀態，其累積下沉量分別達到41 mm、49 mm和58 mm。而2018-01-02后，該區域開始趨于穩定，其形變量波動小于±5 mm。

圖5 區域2采樣點時間序列形變圖

圖6 區域3采樣點位置圖

區域3位于研究區域南部(見圖6)，其最大下沉區域面積均不超過1 km2，但其形變速率較大，且累計形變超過-100 mm，本研究同樣在區域3和4中分別取兩個采樣點進行時序分析，時序形變圖如圖7所示，從圖7中可以看出，P8～P11始終處于下沉狀態中，且下沉量不斷累積。其中P11下沉最為明顯，其最大下沉量達到118 mm。由于該區域中兩個下沉明顯的部分沉降速率較快且保持持續下沉，故可作為重點對象進行監測和后續的考察。

圖7 區域3采樣點時間序列形變圖

4 結論

本次研究基于歐空局Sentinel-1A衛星影像，運用SBAS技術對九寨溝地震區及其周邊地區進行了時間序列形變的監測，基本查明震后災區附近的地表形變趨勢，研究結果表明：

1)研究區域中大部分區域形變速率均維持在-50～50 mm/a。研究區的東北部表現為大面積下沉，且與中部區域相比下沉幅度更大。另外研究區域中存在3處面積不足1 km2的地區(103°41′46″E，33°13′8″N；103°52′E，32°54′N；103°59′E，32°51′N)出現了速率為-80～130 mm/a的大幅度下沉，在2017-08-11—2018-07-01期間最大累積下沉量也達到118 mm。

2)SBAS技術能有效監測大范圍的時序形變，其監測精度可以達到mm級。另外，InSAR技術在監測大范圍形變的同時，還可以監測到小面積、大幅度的地表形變，在地表持續監測上可以減少人力物力財力的投入，實現對形變異常區域更精準的定位。

除此之外，可以對當地歷史影像進行數據處理，還能進行后續形變監測，從而得到更全面的形變規律，為山區形變監測和地質災害預警及防治提供更強大的理論基礎和數據支持。