郭世鵬，康偉，黃繼茂，張庭葦，李云，張王菲

(西南林業(yè)大學(xué) 地理與生態(tài)旅游學(xué)院，云南 昆明 650224)

0 引 言

SAR干涉測量技術(shù)(interferometric synthetic aperture radar，InSAR)可以監(jiān)測地球陸地表面的微小變化，精度可達(dá)毫米級，同時SAR技術(shù)具有全天時、全天候和大面積監(jiān)測等優(yōu)點。由于使用該技術(shù)不需要人進(jìn)入監(jiān)測地區(qū)，因此逐漸成為獲取地震發(fā)生區(qū)震后位移測量的首選方法[1-2]。D.Massonnet等[3]采用時隔幾個月的ERS-1衛(wèi)星干涉像對揭示了Landers地震的同震位移場，基于此獲得的形變圖與GPS監(jiān)測結(jié)果非常一致；張紅等[1]采用ERS-1/2干涉像對提取了張北地區(qū)的地震同震位移場，進(jìn)一步證實了差分干涉SAR(D-InSAR)技術(shù)在地震監(jiān)測和研究中的潛力；臧妻斌等[4]利用兩軌法對意大利拉奎拉地震形變展開研究，成功獲取了同震位移場在雷達(dá)視線方向(light of sight,LOS)的位移；張瑞等[5]應(yīng)用D-InSAR技術(shù)得到了玉樹地震區(qū)域的地表形變,并對地震的震源機制和發(fā)震機理展開分析。D-InSAR技術(shù)由于干涉像對獲取的容易性和獲取位移場算法的簡單性，在未來的地震監(jiān)測研究和監(jiān)測中具有極大潛力。這些研究證實了D-InSAR技術(shù)在地震形變場監(jiān)測中的有效性，同時也指出由于SAR成像方式的影響，即部分區(qū)域影像噪聲較高，相干性低，使得監(jiān)測結(jié)果精度降低，因此需要地表測量數(shù)據(jù)進(jìn)行相應(yīng)的校正或驗證[1-2,6-8]。由于地震產(chǎn)生的同震位移場的區(qū)域通常較大，地表合適的測量數(shù)據(jù)(如GPS、水準(zhǔn)測量等測量數(shù)據(jù))獲取困難，難以檢驗采用D-InSAR獲取的地震形變信息的不確定性和精度[2]。近年來針對升降軌數(shù)據(jù)獲取的地表形變特征及采用彈性半空間模型(Okada模型)模擬三維位移場相融合的方法為以上問題提供了解決思路[9-13]。本文基于D-InSAR獲取地震位移場的成像幾何關(guān)系和測量原理，建立Okada模型模擬地震三維位移場與D-InSAR在LOS方向的數(shù)學(xué)關(guān)系，通過升降軌相關(guān)性分析、Okada模型模擬地震位移場來間接驗證D-InSAR在地震引發(fā)的同震位移反演中的有效性。

1 研究區(qū)概況及數(shù)據(jù)獲取

1.1 研究區(qū)概況

研究區(qū)位于中國臺灣省高雄市(緯度22°09′43″N～23°21′20″N，經(jīng)度120°06′25″E～120°53′38″E)，震中位于美濃區(qū)(22°55′48″N，120°32′24″E)，如圖1所示。地震區(qū)地處左鎮(zhèn)斷層、潮州斷層及旗山斷層之間的地震密集帶，距高雄市區(qū)約30 km。此次地震達(dá)到里氏6.7級，震源深度為15 km。

圖1 研究區(qū)地理位置

1.2 研究區(qū)D-InSAR數(shù)據(jù)獲取

1.2.1 D-InSAR數(shù)據(jù)獲取

本研究獲取了升軌和降軌的Sentinel-1A衛(wèi)星InSAR數(shù)據(jù)各一景。Sentinel-1是由位于同一軌道平面內(nèi)，軌道相位差為180°的1A和1B組成的雙星星座，雙星星座重訪周期為6 d，A，B星各自重訪周期為12 d，高緯度地區(qū)最快重訪時間可達(dá)1～3 d。Sentinel-1A衛(wèi)星搭載一臺C波段(頻率為5.404 GHz)合成孔徑雷達(dá)，包含條帶(SM)、干涉寬幅(IW)、超寬幅(EW)、波浪(WV)4種成像模式，最大覆蓋幅寬可達(dá)400 km。其中IW模式為Sentinel-1A的主要模式，該模式成像由9個burst拼接而成，幅寬為375 km。[14]本研究獲取了Sentinel-1A衛(wèi)星IW模式的2對較短時間間隔的干涉SAR影像對，詳細(xì)參數(shù)信息見表1。升降軌覆蓋研究區(qū)的位置見圖1(b)。

表1 SAR影像信息

1.2.2 其他輔助數(shù)據(jù)獲取

本次研究同時獲取了輔助的DEM數(shù)據(jù)(圖1(a))，用于模擬地形相位和地理編碼，該數(shù)據(jù)為美國航空航天局(NASA)提供的SRTM DEM(分辨率為30 m)數(shù)據(jù)。為了提高軌道誤差對基線估算精度的影響，本文同時獲取了兩景影像的精密軌道數(shù)據(jù)(precise orbit data，POD)，用于軌道誤差改正，進(jìn)而提高基線估算的精度。

2 InSAR數(shù)據(jù)處理

2.1 InSAR數(shù)據(jù)處理流程

InSAR數(shù)據(jù)處理流程如圖2所示，包括兩個步驟：第一步，采用二軌差分法分別對升降軌數(shù)據(jù)進(jìn)行LOS方向地震前后二維位移量的提??；第二步，采用升降軌數(shù)據(jù)和Okada模型反演形變信息并進(jìn)行二維位移場提取結(jié)果驗證。具體處理環(huán)節(jié)如下。

圖2 InSAR數(shù)據(jù)處理流程圖

第一，首先利用ENVI軟件處理InSAR單視復(fù)數(shù)(single look complex，SLC)數(shù)據(jù)；第二，根據(jù)獲取的精密軌道數(shù)據(jù)對用于基線估算的軌道信息進(jìn)行修正，得到更精確的軌道參數(shù)；第三，通過距離向為1和方位向為5的多視處理提高相干性的統(tǒng)計精度；第四，采用Goldstein濾波對多視后的影像進(jìn)行去噪[15]，最小費用流(minimum cost network flow phase unwrapping algorithm，MCF)進(jìn)行相位解纏獲取形變相位；第五，經(jīng)過地理編碼獲得在LOS向的位移信息；第六，分別使用升降軌數(shù)據(jù)獲得的位移信息進(jìn)行相關(guān)性分析；第七，將降軌數(shù)據(jù)獲得的位移信息帶入Okada模型模擬地震同震場的信息,以此驗證采用D-InSAR反演的位移場位移結(jié)果是否具有現(xiàn)實指導(dǎo)意義[8]。

2.2 Okada模型算法

Okada通過分析各種類型的斷層錯動引起位移的解析表達(dá)式，剔除了通用點源及有限面源彈性模型，以走滑方向、傾滑方向和張性方向斷層單位位錯描述的彈性半無限空間表征地表位移公式，進(jìn)而模擬地震場的時空演化過程[9-10,12]。該算法描述的位移場ui可表示為

(1)

圖3 Okada模型坐標(biāo)系

(2)

(3)

(4)

(5)

(6)

(7)

式(5)中,

Sloc=aU1+bU2+cU3，

(8)

其中，a,b,c可以通過雷達(dá)成像幾何關(guān)系和U1,U2,U3幾何關(guān)系計算獲得。

3 結(jié)果與討論

3.1 InSAR技術(shù)提取結(jié)果

圖4為覆蓋整個地震形變區(qū)域的干涉紋圖，上盤干涉條紋密集且閉合，呈現(xiàn)4個同心圓環(huán)狀，一個同心圓代表半個雷達(dá)波長的形變量。下盤干涉條紋稀疏且條紋不太明顯，存在非閉合條紋，右下部分地區(qū)呈現(xiàn)較大失相干現(xiàn)象。從整體看，升降軌獲得的干涉圖大致相同，反演的同心圓數(shù)目也相同。但由于衛(wèi)星拍攝角度等多種因素和植被山體疊掩現(xiàn)象，升軌干涉紋圖和降軌干涉紋圖略有區(qū)別。雖然上盤和下盤均形成干涉條紋，但兩者差異明顯。上盤干涉條紋和下盤條紋的差異表明此次地震為逆沖型地震，該結(jié)果與全球多家監(jiān)測機構(gòu)給出的相關(guān)震源機制解釋一致。

圖4 地震形變區(qū)域干涉紋圖

筆者將升降軌兩組干涉紋圖經(jīng)過相位解纏、地理編碼到相同地理坐標(biāo)系下，并將相位值轉(zhuǎn)為位移量(圖5)。圖5(a)為升軌影像提取的LOS方向的位移信息，圖5(b)為降軌影像提取的位移信息。從圖5 D-InSAR反演的同震位移場可以看出，兩組D-InSAR結(jié)果較為一致，上盤主要表現(xiàn)為抬升，下盤則為下沉。升降軌位移場不僅具有相同的位移趨勢，并且上盤和下盤的位移值也基本吻合，其中升降軌監(jiān)測的位移圖中在斷層的上盤區(qū)域分別抬升14.3 cm和17.2 cm，下盤分別下降10.1 cm和15.7 cm，兩者在上盤和下盤的抬升和下降主要是由于衛(wèi)星獲取影像時間差異和升降軌入射角度差異造成的。由于SAR斜距成像，對距離的變化較敏感，因此斷層西部上盤在LOS方向遠(yuǎn)離升降軌，表現(xiàn)為較大抬升；相反，東部下盤LOS方向遠(yuǎn)離降軌，靠近升軌，使得降軌下沉值大于升軌下沉值。

圖5 InSAR同震位移場提取結(jié)果

從地質(zhì)構(gòu)造看，臺灣西南部地處復(fù)雜的破碎活動斷層，并且西南地區(qū)發(fā)生的地震基本都與這些斷層走滑相關(guān)。本次地震震中位于旗山斷層和潮州斷層的交界附近，結(jié)合升降軌LOS位移場的位置關(guān)系可以看出，同震位移場的沉降位移區(qū)位于地震帶的東南盤，而隆升位移區(qū)位于地震帶的西北盤，這一特征與逆沖斷層的性質(zhì)相符[14-16]。此外，沉降區(qū)與隆升區(qū)分界明顯，并且分界位置處于旗山斷層帶附近。因此，推測出此次地震很有可能是由多個斷層帶位移綜合引起的。由于研究區(qū)所在的臺灣島地處歐亞板塊與菲律賓板塊的相互擠壓帶,構(gòu)造復(fù)雜，在北部和東部，菲律賓海板塊沿琉球海溝向西北偏北的歐亞大陸俯沖；島嶼以南，南海(在歐亞板塊上)在馬尼拉海溝向東延伸至菲律賓海板塊之下，板塊的擠壓作用導(dǎo)致臺灣地形南北狹長。在地震發(fā)生的地方，菲律賓板塊以每年8 cm左右的速率向西北-東南方向推擠，因此推測此次地震是菲律賓海板塊與東亞大陸在臺東縱谷碰撞仰沖而強烈向西推壓影響的結(jié)果[17-18]。

3.2 升降軌互補驗證位移精度

由于缺乏當(dāng)?shù)貙嶋H的水準(zhǔn)測量數(shù)據(jù)，本次研究基于統(tǒng)計方法從沿最大抬升區(qū)和最大沉降區(qū)驗證D-InSAR提取的震后位移信息的精度。消除地形誤差后，在升降軌數(shù)據(jù)的反演結(jié)果中共選出54對典型位移點進(jìn)行相關(guān)性分析，分析結(jié)果如圖6所示。兩組數(shù)據(jù)之間相關(guān)性較高，R2=0.978 4。升降軌數(shù)據(jù)的高相關(guān)性表明采用Sentinel-1A數(shù)據(jù)的D-InSAR技術(shù)在反演地震同震位移場中LOS方向位移的有效性。

圖6 升降軌位移結(jié)果相關(guān)性分析

3.3 Okada模型反演位移場驗證精度分析

為進(jìn)一步分析地震機制與D-InSAR反演地震位移場的精確性，研究選用降軌數(shù)據(jù)反演的位移場提取用于Okada模型的反演震源參數(shù)。將D-InSAR同震位移場進(jìn)行采樣后，參考CMT(centroid moment tensor)提供的經(jīng)度、緯度、走向、傾向、滑動角、深度、斷層長度、寬度以及滑動量等求解非線性反演斷層。

依據(jù)CMT官網(wǎng)推薦的地震參數(shù)推薦區(qū)間反復(fù)迭代，不斷縮小區(qū)間范圍，直至得到合適搜索區(qū)間，并根據(jù)搜索區(qū)間范圍確定參數(shù)值進(jìn)行反演，將反演結(jié)果的參數(shù)設(shè)為最優(yōu)值。經(jīng)過Okada模型的非線性反演，得到如表2所示的參數(shù)值。表2同時比較了反演結(jié)果與CMT及USGS(united states geological survey)公布的震源參數(shù)，從表2可知，反演結(jié)果與CMT接近，與USGS公布的震源深度和傾角有一定差異，但從地震海灘球、震級、地震矩綜合分析可知，三者一致性較好。

表2 CMT、USGS和Okada形變場參數(shù)

采用Okada模型反演出的地震產(chǎn)生的破裂長度為9 890.4 m，寬度為3 976.2 m，震源深度為14.59 km，震級為MW6.4，累計釋放地震矩為5.25×1018N·m，非線性正演斷層模型反演結(jié)果如圖7所示。

由圖7斷層滑動模型結(jié)果可知，地震形成的斷層面主要位于震中西北方向，并處在潮州斷層和旗山斷層附近，斷裂發(fā)生在向西北-東南方向。向東北方向淺傾，或在向西陡傾的南北走向構(gòu)造上。此外，地震還伴有走滑現(xiàn)象，斷層的滑動角為1.26°，傾角為22°，滑動量為4.47 m。由此可以斷定，本次地震的發(fā)震機制為逆沖斷層，這與D-InSAR獲取的逆沖走滑斷層發(fā)震機制相符，進(jìn)一步說明了D-InSAR反演同震位移場的可行性。

圖7 Okada模型地震參數(shù)反演結(jié)果

由以上分析可以確定，高雄地震由該區(qū)域斷層帶引起，發(fā)震機制為逆沖走滑斷層，滑動量為4.47 m，InSAR得到上盤最大抬升量約為17 cm，下盤最大沉降量約為 12 cm。已有研究采用ALOS-2 PALSAR-2的D-InSAR的結(jié)果表明，該區(qū)域的最大抬升高度約為13 cm，與本研究的研究結(jié)果相差在5 cm內(nèi)[8]。而同樣采用2016年2月14日和2月02日的Sentinel-1A的反演結(jié)果表明研究區(qū)總體位移量為-4～12 cm，與本研究的結(jié)果基本一致[7]。

4 結(jié) 語

基于Sentinel-A升降軌數(shù)據(jù)，利用D-InSAR技術(shù)得到LOS向同震位移場，基于同震位移場及干涉相位圖，確定了研究區(qū)地震場的范圍，地震引起的地面抬升和沉降程度，結(jié)合研究區(qū)的地質(zhì)構(gòu)造特征，揭示了地震的形成機制。為了驗證D-InSAR技術(shù)得到同震位移場的有效性，采用54個典型點的升降軌形變量進(jìn)行了相關(guān)性比較，其中升降軌典型特征點形變量線性擬合R2為0.978 4，說明了D-InSAR技術(shù)在地震測量中有較大的應(yīng)用潛力。采用D-InSAR技術(shù)降維獲得的參數(shù)輸入Okada模型反演位移場特征，反演結(jié)果與已有CMT和USGS官方公布結(jié)果一致性較好，Okada模型反演參數(shù)顯示，地震位于旗山斷層和潮州斷層的斷裂帶上，深度為14.59 km，發(fā)震機制為走滑逆沖斷層，與D-InSAR反演位移場分析結(jié)果基本一致。研究表明，利用D-InSAR可以有效進(jìn)行地震同震位移場的反演，而在沒有對應(yīng)區(qū)域測量驗證數(shù)據(jù)時，可以采用升降軌技術(shù)驗證研究區(qū)同震位移場的反演結(jié)果，也可采用Okada模型的反演結(jié)果定量和定性分析D-InSAR技術(shù)獲得的位移場有效性。