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

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

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

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

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

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

1 同震形變場獲取

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

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

紅色線表示研究區(qū)范圍內(nèi)斷層分布；紅色五角星表示主震位置；綠色圓圈表示余震分布；藍(lán)色虛線框為獲取的ASAR數(shù)據(jù)分布情況；Gowk和Mayband為斷層名稱圖1 2003年12月26日巴姆Mw6.6地震區(qū)域構(gòu)造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處理平臺,使用航天飛機雷達(dá)地形測繪使命(shuttle radar topography mission,SRTM)獲取的分辨率為90 m的數(shù)字地形高程模型(digital elevation model,DEM)對干涉圖進(jìn)行去地形處理,再通過濾波、相位解纏以及軌道誤差校正,最后得到地理編碼后的同震形變結(jié)果(圖2)。從圖2(a)可以看出,降軌獲取的形變場結(jié)果左下角存在一定的地形影響，而且由于地震震級較大,極震區(qū)破壞嚴(yán)重,導(dǎo)致升軌和降軌獲取的同震形變場中極震區(qū)都存在一定的失相干現(xiàn)象,失相干區(qū)域的存在會影響后續(xù)三維形變場解算的準(zhǔn)確度。

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

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

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)方法進(jìn)行失相干恢復(fù),首先對SAR數(shù)據(jù)獲取的形變場進(jìn)行分類,對形變場中非0區(qū)域進(jìn)行數(shù)據(jù)預(yù)處理,建立訓(xùn)練集和測試集；對形變場中0區(qū)域進(jìn)行一定的數(shù)據(jù)預(yù)處理；然后結(jié)合訓(xùn)練集、測試集以及0結(jié)果進(jìn)行核函數(shù)的選擇,進(jìn)行SVR模型訓(xùn)練和模型預(yù)測,最終獲取完整的形變場。

圖3、圖4為經(jīng)過失相干恢復(fù)和地形誤差去除后的升降軌形變場。圖3為升軌數(shù)據(jù)獲取的失相干恢復(fù)后的巴姆地震同震形變場,可以發(fā)現(xiàn)失相干恢復(fù)后的結(jié)果保證了形變場的完整性。圖4(a)為失相干恢復(fù)后的降軌形變場,可以看出,失相干恢復(fù)后,保持了形變場的完整性,但是形變場中存在明顯的地形相關(guān)誤差；圖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ū)域進(jìn)行處理后相位差的變化。最終使用進(jìn)行失相干恢復(fù)并去除地形誤差后的降軌形變場[圖4(c)]及升軌數(shù)據(jù)形變場(圖3)進(jìn)行后續(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向形變場結(jié)果[圖3、圖4(c)]可以看出,地震造成的形變中心呈蝴蝶狀,定義LOS向形變場中靠近衛(wèi)星飛行方向為正值,遠(yuǎn)離衛(wèi)星飛行方向為負(fù)值。從升軌同震形變場(圖3)可以看出,斷層左側(cè)有兩個明顯的形變中心,左上形變中心表現(xiàn)為正值,左下表現(xiàn)為負(fù)值,相比之下斷層右側(cè)形變中心量級較小,但依然能看出兩個特征不同的形變中心；從降軌同震形變場圖4(c)可以看出,升降軌同震形變場特征類似,整體表現(xiàn)為蝴蝶狀,但是受升降軌特征影響,降軌中斷層右側(cè)形變特征表現(xiàn)的更明顯,最大形變量近0.3 m。

1.2 MAI技術(shù)獲取方位向形變

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

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

圖5 MAI技術(shù)獲取的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向形變結(jié)果和方位向形變結(jié)果與地表真實東西向、南北向、垂直向形變結(jié)果的關(guān)系表達(dá)式為

(2)

(3)

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

(4)

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

X=(BTPB)-1BTPL

(5)

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

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

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

3 結(jié)論

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

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

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