基于SBAS技術(shù)的九寨溝震后地表形變監(jiān)測

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

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

2017-08-08T21：19于四川九寨溝(103.82°E，33.20°N)發(fā)生Ms7.0地震，震源深度約為20 m。地震發(fā)生后，九寨溝景區(qū)遭受到較大的破壞，本次地震造成了25人死亡，520余人受傷，同時7萬余房屋受到不同程度的損壞[1-2]。而地震之后一段時間內(nèi)當(dāng)?shù)氐匦蔚淖兓蔀闉?zāi)后重建的一項(xiàng)重要監(jiān)測內(nèi)容。

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

本次實(shí)驗(yàn)選取25景Sentinel衛(wèi)星數(shù)據(jù)，使用小基線集合成孔徑雷達(dá)干涉測量技術(shù)對九寨溝地震區(qū)及其附近山區(qū)進(jìn)行測量，得到該地區(qū)震后近1年內(nèi)的時間序列形變，為探究該地區(qū)震后的沉降規(guī)律、災(zāi)后重建及地質(zhì)災(zāi)害預(yù)防提供技術(shù)和數(shù)據(jù)支持。

1 研究區(qū)域及實(shí)驗(yàn)數(shù)據(jù)

1.1 研究區(qū)域

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

圖1 研究區(qū)域覆蓋范圍圖

1.2 實(shí)驗(yàn)數(shù)據(jù)

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

表1 Sentinel-1A影像信息

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

2 SBAS技術(shù)原理

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

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

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

(1)

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

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

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

(2)

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

4)采用選取多個相對穩(wěn)定的地面點(diǎn)的方法求取干涉圖里的大氣相位以及噪聲項(xiàng)成分。在求取誤差項(xiàng)以后，使用高次多項(xiàng)式插值的方法，對研究區(qū)域中的大氣相位進(jìn)行插值，達(dá)到恢復(fù)研究區(qū)域中每一景干涉圖所包含的大氣相位及噪聲的目的。

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

(3)

(4)

Ax=Δφ.

(5)

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

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

(6)

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

Bφ=δφ.

(7)

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

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

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

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

表2 九寨溝采樣區(qū)域信息及采樣點(diǎn)分布信息

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

圖4 區(qū)域1采樣點(diǎn)時間序列形變圖

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

圖5 區(qū)域2采樣點(diǎn)時間序列形變圖

圖6 區(qū)域3采樣點(diǎn)位置圖

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

圖7 區(qū)域3采樣點(diǎn)時間序列形變圖

4 結(jié)論

本次研究基于歐空局Sentinel-1A衛(wèi)星影像，運(yùn)用SBAS技術(shù)對九寨溝地震區(qū)及其周邊地區(qū)進(jìn)行了時間序列形變的監(jiān)測，基本查明震后災(zāi)區(qū)附近的地表形變趨勢，研究結(jié)果表明：

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

2)SBAS技術(shù)能有效監(jiān)測大范圍的時序形變，其監(jiān)測精度可以達(dá)到mm級。另外，InSAR技術(shù)在監(jiān)測大范圍形變的同時，還可以監(jiān)測到小面積、大幅度的地表形變，在地表持續(xù)監(jiān)測上可以減少人力物力財(cái)力的投入，實(shí)現(xiàn)對形變異常區(qū)域更精準(zhǔn)的定位。

除此之外，可以對當(dāng)?shù)貧v史影像進(jìn)行數(shù)據(jù)處理，還能進(jìn)行后續(xù)形變監(jiān)測，從而得到更全面的形變規(guī)律，為山區(qū)形變監(jiān)測和地質(zhì)災(zāi)害預(yù)警及防治提供更強(qiáng)大的理論基礎(chǔ)和數(shù)據(jù)支持。