巴姆地震三維形變場獲取

張慶云, 張景發(fā), 李永生, 李兵權, 羅三明, 馬慶尊

(1.中國地震局第一監(jiān)測中心, 天津 300180； 2.應急管理部國家自然災害防治研究院, 北京 100085)

隨著合成孔徑雷達干涉測量技術(interferometric synthetic aperture radar, InSAR)的不斷發(fā)展,在地震發(fā)生后,基于合成孔徑雷達(synthetic aperture radar,SAR)數(shù)據(jù)快速獲取地震同震形變場,可以對地震造成的破壞程度和范圍進行直觀描述。雖然使用差分干涉測量(differential interferometric synthetic aperture radar,D-InSAR)進行同震形變場獲取的技術已經(jīng)相對成熟,但對于地形復雜地區(qū),仍然會存在地形誤差影響；而且傳統(tǒng)D-InSAR技術獲取的同震形變場是衛(wèi)星視線向(line of sight,LOS)結果,無法反映同震形變的全貌,因此三維形變場的解算對深入了解發(fā)震構造以及地震造成的形變具有重要意義。

目前結合InSAR技術獲取三維形變場的方法主要有:①不同傳感器升降軌InSAR數(shù)據(jù)分析,該方法的缺點是對數(shù)據(jù)要求較高,而且由于SAR獨特成像原理,導致其對南北向不敏感,所以結合升降軌的三維解算在南北向存在較大誤差,同時目前不同傳感器的入射角和方位角相差不大,也會導致解算存在誤差[1-2]；②InSAR技術聯(lián)合外部數(shù)據(jù)全球定位系統(tǒng)(global positioning system,GPS)數(shù)據(jù)、或者根據(jù)模型模擬得到的三維結果進行三維解算,該方法缺點是有些地區(qū)沒有GPS覆蓋,模型模擬會引入一定程度的模型誤差[3-6]；③不同InSAR技術聯(lián)合進行三維解算,D-InSAR技術獲取LOS向形變,多孔徑干涉雷達測量技術(multiple aperture interferometry,MAI)可以獲取方位向形變場。偏移跟蹤技術(offset-tracking)可以獲取地震的方位向和距離向形變結果,綜合使用升降軌數(shù)據(jù)的LOS向形變以及方位向或者距離向形變可以獲取三維形變結果[7-10]。

經(jīng)過中外學者多年努力,目前InSAR三維形變場獲取技術已日益成熟趨于完善,在同震、地面沉降和滑坡等各類地質(zhì)災害形變監(jiān)測中得到廣泛應用。He等[11]、張慶云等[12-13]利用InSAR技術獲取了2017年伊朗地震的三維形變場；Chen等[14]、姚佳明等[15]和董龍凱[16]利用InSAR技術進行了礦區(qū)三維形變場解算研究；Shi等[17]和王晨興等[18]基于InSAR技術研究了滑坡三維形變特征,并取得了可靠的分析結果；楊麗葉等[19]利用InSAR技術分析了錯郎瑪冰川的三維形變特征；Pepe[20]和馬艷鴿[21]分析了InSAR技術獲取三維形變場的方法原理。雖然已有的研究成果已經(jīng)相對成熟,但是仍存在不足之處,如MAI技術和D-InSAR技術對SAR影像相干性要求相對較高,導致地震、礦區(qū)等形變量較大區(qū)域,存在嚴重失相干；以及在地形復雜區(qū)域,地形誤差如何準確去除等。失相干以及各類誤差的去除對三維形變場解算結果的準確性都有重要影響,這些都是當前仍需要解決的主要問題。

針對InSAR技術獲取的形變場中存在失相干及地形等相關誤差去除問題,以常用的第3種三維形變場解算方法為基礎,以2003年12月26日巴姆地震為例,研究地形誤差去除方法、完整形變場獲取技術,基于D-InSAR和MAI技術進行三維形變場解算,分析地震三維形變場結果。

1 同震形變場獲取

2003年12月26日,伊朗東南部巴姆(Bam)古城發(fā)生Mw6.6地震,震中位置在29.01°N,58.26°E,地震造成四萬多人死亡,五萬多人受傷,成千上萬人無家可歸,巴姆古城內(nèi)幾乎所有的建筑物都被損毀,地震造成了嚴重的人員傷亡和財產(chǎn)損失。

地震發(fā)生后,Zare等[22]通過分析伊朗國家強震臺網(wǎng)獲取的強震記錄,認為Bam斷層長約10 km,且經(jīng)過巴姆市；Tatar等[23]分析了巴姆地震余震分布特征；Wang等[24]使用先進合成孔徑雷達(advanced synthetic aperture radar,ASAR)數(shù)據(jù)分析了巴姆地震發(fā)震特征；季靈運等[25]、王洪友[26]使用ASAR數(shù)據(jù)獲取了巴姆地震三維形變場,但是上述研究存在一定的局限性,形變場中存在一定程度的失相干現(xiàn)象和地形誤差。在已有研究基礎上,對地形誤差進行進一步去除,同時對失相干區(qū)域進行恢復,然后基于多種技術聯(lián)合的方法進行三維形變場解算。圖1為巴姆地震區(qū)域構造圖。

紅色線表示研究區(qū)范圍內(nèi)斷層分布；紅色五角星表示主震位置；綠色圓圈表示余震分布；藍色虛線框為獲取的ASAR數(shù)據(jù)分布情況；Gowk和Mayband為斷層名稱圖1 2003年12月26日巴姆Mw6.6地震區(qū)域構造Fig.1 Regional structure of Bam Mw6.6 earthquake on December 26, 2003

1.1 LOS向形變場獲取

地震發(fā)生后,歐洲航天局ENVISAT衛(wèi)星獲取了地震震前和震后SAR數(shù)據(jù),由于該區(qū)域植被稀疏且以平原為主,因此適合用SAR數(shù)據(jù)分析其形變場。選擇合適的SAR數(shù)據(jù),獲取地震的同震形變場,表1為主要數(shù)據(jù)參數(shù)。

使用兩軌差分方法,以GAMMA軟件為InSAR處理平臺,使用航天飛機雷達地形測繪使命(shuttle radar topography mission,SRTM)獲取的分辨率為90 m的數(shù)字地形高程模型(digital elevation model,DEM)對干涉圖進行去地形處理,再通過濾波、相位解纏以及軌道誤差校正,最后得到地理編碼后的同震形變結果(圖2)。從圖2(a)可以看出,降軌獲取的形變場結果左下角存在一定的地形影響，而且由于地震震級較大,極震區(qū)破壞嚴重,導致升軌和降軌獲取的同震形變場中極震區(qū)都存在一定的失相干現(xiàn)象,失相干區(qū)域的存在會影響后續(xù)三維形變場解算的準確度。

圖2 基于ASAR數(shù)據(jù)獲取的巴姆地震同震原始形變場Fig.2 Coseismic original deformation field of Bam earthquake based on ASAR data

為保證獲取可靠的形變場,針對降軌中存在地形影響的區(qū)域進行誤差去除。首先勾選出存在地形影響的區(qū)域,然后對存在地形誤差區(qū)域的參考DEM與形變場經(jīng)緯度之間建立函數(shù)模型,進行地形誤差去除。最終使用的模型為

Z=[Ddisp-0.021 237-0.002 076LlatDdem-(-0.002 617 8)LlonDdem-(-9.989×10-6)Llat2Ddem-(-1.608 2×10-5)Llon2Ddem-(-2.580 5×10-5)Llon2Ddem-(-0.106 15)Ddem]

(1)

式(1)中:Z為去除地形誤差之后點的形變量；Ddisp為未去除地形誤差時點的形變量；Llat為點的緯度；Llat為點的經(jīng)度；Ddem為點的高程。

針對極震區(qū)失相干現(xiàn)象,采用基于非線性支持向量機(support vector regression,SVR)方法進行失相干恢復,首先對SAR數(shù)據(jù)獲取的形變場進行分類,對形變場中非0區(qū)域進行數(shù)據(jù)預處理,建立訓練集和測試集；對形變場中0區(qū)域進行一定的數(shù)據(jù)預處理；然后結合訓練集、測試集以及0結果進行核函數(shù)的選擇,進行SVR模型訓練和模型預測,最終獲取完整的形變場。

圖3、圖4為經(jīng)過失相干恢復和地形誤差去除后的升降軌形變場。圖3為升軌數(shù)據(jù)獲取的失相干恢復后的巴姆地震同震形變場,可以發(fā)現(xiàn)失相干恢復后的結果保證了形變場的完整性。圖4(a)為失相干恢復后的降軌形變場,可以看出,失相干恢復后,保持了形變場的完整性,但是形變場中存在明顯的地形相關誤差；圖4(b)為存在地形誤差的降軌數(shù)據(jù)纏繞相位；圖4(c)為使用函數(shù)模型去除地形誤差后降軌數(shù)據(jù)形變場；圖4(d)為去除地形誤差后纏繞相位。從圖4(b)、圖4(d)對比可以看出,去除地形誤差前,存在地形誤差的區(qū)域與周邊區(qū)域存在很大的相位差。經(jīng)過去地形誤差處理后,整體形變場除中心主形變區(qū)外,周邊區(qū)域的形變量都趨于0。從圖4(d)可以看出，原本存在相位差的地形區(qū)域進行處理后相位差的變化。最終使用進行失相干恢復并去除地形誤差后的降軌形變場[圖4(c)]及升軌數(shù)據(jù)形變場(圖3)進行后續(xù)三維分析處理。

圖3 2003年12月26日巴姆Mw6.6地震升軌數(shù)據(jù)LOS向形變場Fig.3 Los direction deformation field of the orbit ascending data of Bam Mw6.6 earthquake on December 26, 2003

從D-InSAR獲取的LOS向形變場結果[圖3、圖4(c)]可以看出,地震造成的形變中心呈蝴蝶狀,定義LOS向形變場中靠近衛(wèi)星飛行方向為正值,遠離衛(wèi)星飛行方向為負值。從升軌同震形變場(圖3)可以看出,斷層左側有兩個明顯的形變中心,左上形變中心表現(xiàn)為正值,左下表現(xiàn)為負值,相比之下斷層右側形變中心量級較小,但依然能看出兩個特征不同的形變中心；從降軌同震形變場圖4(c)可以看出,升降軌同震形變場特征類似,整體表現(xiàn)為蝴蝶狀,但是受升降軌特征影響,降軌中斷層右側形變特征表現(xiàn)的更明顯,最大形變量近0.3 m。

1.2 MAI技術獲取方位向形變

MAI技術[27]是利用雷達波束多普勒頻移的正負,分裂波束的InSAR處理算法,從而獲得目標地物方位向形變的信息。將兩幅單視復數(shù)影像根據(jù)多普勒頻移的符號,分裂成前后視兩對主從影像,并分別進行干涉。由于前后視干涉圖的空間基線幾乎相等,又受到了同等的大氣影像,而且它們探測到的距離向的形變也是相同的,所以將前視和后視的干涉圖進行相位差分,得到方位向形變的相位信息。

圖5為MAI技術獲取的巴姆地震方位向形變場,受數(shù)據(jù)影響,升軌數(shù)據(jù)對形變場的覆蓋并不完整,降軌數(shù)據(jù)可以完整覆蓋形變場。從圖5可以看出,升降軌數(shù)據(jù)獲取的方位向形變場在斷層兩側都有明顯的特征,升軌中[圖5(a)],斷層左側表現(xiàn)為正值,斷層右側表現(xiàn)為負值,降軌數(shù)據(jù)與之相反[圖5(b)],最大方位向形變量為1 m。

圖5 MAI技術獲取的2003年12月26日巴姆Mw6.6地震方位向形變場Fig.5 Azimuth deformation field of Bam Mw6.6 earthquake on December 26, 2003 acquired by MAI Technology

2 三維形變場解算

LOS向形變結果和方位向形變結果與地表真實東西向、南北向、垂直向形變結果的關系表達式為

(2)

(3)

通過獲取升降軌LOS向和方位向形變場,根據(jù)升降軌情況,可以得到4個方程,轉(zhuǎn)化為矩陣形式可表示為

(4)

根據(jù)最小二乘原理,三維形變場的最優(yōu)解為

X=(BTPB)-1BTPL

(5)

式(5)中：P為觀測值方差組成的權陣。

在SAR形變場結果中,選擇合適的窗口,計算窗口中心點標準差,從而估計出數(shù)據(jù)的近似標準差。

根據(jù)式(4)、式(5)方程,使用多種InSAR技術聯(lián)合的三維地表形變解算方法,獲取巴姆地震的三維形變場,由于受升軌數(shù)據(jù)的影響,獲取的三維形變場結果主要覆蓋斷層右側信息,左側信息不完全。從三維解算結果(圖6)可以看出,在東西向上[圖6(a)],定義向東為正,向西為負,主要表現(xiàn)為3個形變中心,斷層左側有一個向西形變中心,斷層右側有2個形變中心,其中右上為向西形變,右下為向東形變,最大向西和向西形變量為0.25 m；在南北向上[圖6(b)],定義向北為正,向南為負,地震造成的形變沿斷層方向斷層左側為向北運動,右側向南運動,最大向南形變量為0.4 m；垂直向上[圖6(c)],定義向上為正,向下為負,主要表現(xiàn)為三個形變中心,斷層左側為向上抬升,斷層右側兩個形變中心,右上為向下沉降,右上為向上抬升,最大抬升形變量為0.24 m,最大下沉形變量為0.2 m。從三維形變場結果可以看出,巴姆地震的發(fā)震斷層主要表現(xiàn)為近南北向的右旋走滑運動,同時存在一定的東西向水平旋轉(zhuǎn)運動。獲得的三維形變場結果與洪順英[28]、季靈運等[25]和Fialko等[29]的研究成果一致。

3 結論

(1)引入非線性支持向量機回歸(SVR)對形變場中失相干區(qū)域進行恢復,保證了形變場的完整性,為后續(xù)研究分析提供基礎。

(2)構建DEM與形變場經(jīng)緯度之間函數(shù)模型,對形變場中存在的地形誤差進行去除,提高了形變場的可靠性。

(3)聯(lián)合D-InSAR獲取的LOS向同震形變場與MAI技術獲取的方位向形變場,使用多種InSAR技術聯(lián)合的方法進行三維形變場解算,對地震造成的地面形變進行更可靠的分析,通過巴姆地震的三維形變場解算結果可以看出,巴姆地震發(fā)震斷層為近南北向的右旋走滑斷層。完整三維形變場的獲取,為后續(xù)發(fā)震斷層研究、地質(zhì)構造分析提供理論基礎。