升降軌InSAR數(shù)據(jù)約束下的2007年阿里地震反演分析

溫揚茂，許才軍，劉 洋，馮萬鵬，李志才

1.武漢大學(xué)測繪學(xué)院，湖北武漢430079；2.武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，湖北 武漢430079；3.格拉斯哥大學(xué)地理與地球科學(xué)學(xué)院，英國 格拉斯哥G12 8QQ；4.國家基礎(chǔ)地理信息中心大地測量部，北京100048

1 引 言

北京時間2007年5月5日16時51分，在西藏阿里地區(qū)日土縣與改則縣交界處發(fā)生了Mw 6.1級地震（本文稱為阿里地震），地震震中位于81.99°E，34.26°N，震源深度為7.0km（來自美國地質(zhì)調(diào)查局（United States Geological Survey，USGS）。由于此次地震所在地的大部分區(qū)域為無人區(qū)，沒有造成人員傷亡。該地震是自2005年以來發(fā)生在青藏高原一系列正斷層事件（主要包括2005年仲巴Mw6.2級地震、2008年改則Mw6.4級地震與Mw5.9級地震、2008年于田Mw7.1級地震、2008年仲巴Mw6.7級地震和2008年當(dāng)雄Mw6.3級地震等）中震級較小的一個，可能進一步佐證了青藏高原板塊內(nèi)部在重力作用下東西向擴張的事實［1］。

2007年阿里地震震區(qū)位于班公錯—怒江斷裂帶以北活動構(gòu)造帶上，該構(gòu)造帶主要由發(fā)育較好的NE和NW向活動斷裂組成，如托和平錯斷裂、美馬錯—德克瑪洛斷裂和魯瑪江冬錯斷裂等（圖1）。這一區(qū)域的隆起帶隆升較強烈，呈單斜式隆升，并且斷陷帶也較發(fā)育。其中，第四紀沖、洪積扇分布廣泛并有較大湖泊分布［2］。在過去的30年中，2007年阿里地震發(fā)生前，該地區(qū)的地震活動并不強烈。據(jù)USGS地震目錄記載，在1975年至2005年的30年間，主震震中1°范圍內(nèi)發(fā)生的最大地震為1985年的Mw5.0級地震。因此，研究這樣一個少有的中強地震的同震位移場與斷層參數(shù)對認識該地區(qū)的構(gòu)造運動特征具有重要意義。

圖1 2007年阿里地震震中區(qū)域地質(zhì)構(gòu)造背景Fig.1 Tectonic setting of epicenter region of the 2007Ali earthquake

此外，由于阿里地震發(fā)生在藏北無人區(qū)，野外地質(zhì)調(diào)查和地球物理數(shù)據(jù)的采集都非常困難，對它的研究起初只能依賴于全球地震臺網(wǎng)（GSN）記錄的寬頻帶遠震體波數(shù)字資料，據(jù)此給出粗略的矩張量解和最佳雙力偶解。而作為新興的空間大地測量技術(shù)的InSAR技術(shù)，能夠獲取高精度、高分辨率的地表形變信息，且不需要提供地面控制點，非常適合偏遠地區(qū)的地震形變監(jiān)測［3－5］。需要特別指出的是，InSAR僅能觀測到一維形變量，它獲取到的是雷達波視線向（line－of－sight，LOS）的形變，是地表三維形變（通常指東西向、南北向和垂向形變）在視線向上的投影（圖2）。在絕大多數(shù)情況下，這種一維形變不能完全反映地表的真實形變，有時候甚至?xí)鹌?。不同成像幾何條件下的干涉數(shù)據(jù)（如具有不同成像幾何特征的升、降軌數(shù)據(jù)）能夠提供更為豐富的形變信息，可以有效地降低視線向的模糊問題［6］，從而對地震震源機制和破裂過程的反演等提供更可靠的數(shù)據(jù)約束。本文利用EnviSat衛(wèi)星的升、降軌數(shù)據(jù)來獲取2007年阿里地震的高質(zhì)量同震地表形變場，然后采用線性和非線性反演方法來確定地震的斷層參數(shù)和滑動分布，反演結(jié)果將深化對阿里地震的成因和機理的認識。

圖2 右視SAR的升、降軌成像幾何Fig.2 The geometry of the ascending and descending right－looking SAR acquisitions

2 InSAR數(shù)據(jù)分析

地震發(fā)生前后，歐洲航天局（European Space Agency，ESA）的EnviSat衛(wèi)星C波段ASAR傳感器對阿里地震進行了觀測。為了研究本次地震的幾何結(jié)構(gòu)和運動學(xué)特征，筆者通過中國科技部和歐洲航天局聯(lián)合資助的“龍計劃”三期國際合作項目收集了覆蓋此次地震發(fā)生區(qū)域的雷達影像，為了盡可能獲取高相干數(shù)據(jù)和減小非同震形變信號的影響，最終選取了兩組像對（其中升、降軌各一組）來進行干涉處理（表1）。

表1 研究中所用EnviSat衛(wèi)星數(shù)據(jù)的基本信息Tab.1 Details of EnviSat satellite images used in this study

干涉處理使用的平臺為瑞士GAMMA軟件［7］?；谠撥浖?，采用二通法［8］對獲取的衛(wèi)星SAR影像進行差分干涉處理來獲取阿里地震的同震地表形變場。在數(shù)據(jù)處理過程中，使用ESA提供的DOR精密軌道和NASA提供的90m分辨率的SRTM DEM數(shù)據(jù)［9］來去除地形相位的影響。由于干涉像對的空間基線較長，因此盡管使用了較為精確的DOR軌道數(shù)據(jù)，但是在得到的干涉圖中仍然存在由于殘余軌道誤差所引起的殘余相位影響。對于這個影響，本文采用多項式模型來擬合去除。同時為了降低干涉相位的噪聲水平，提高干涉圖的信號質(zhì)量，采用基于能量譜的局部自適應(yīng)濾波對干涉圖進行濾波。最終，采用枝切法來解纏得到差分干涉相位。阿里地區(qū)地處高寒山區(qū)，平均海拔為5000m，氣候干燥，大氣影響相對較小。盡管如此，為了進一步削弱干涉相位中的大氣相位貢獻，還根據(jù)地形與大氣相位的相關(guān)性擬合去除了與地形相關(guān)的大氣相位［10］。最終得到了地理編碼后的升、降軌阿里地震的高精度同震地表形變場（圖3）。

圖3 2007年阿里地震的升軌T155A（a）和降軌T205D（b）同震地表形變場Fig.3 Coseismic interferograms associated with the 2007Ali earthquake from（a）ascending track T155Aand（b）descending track T205D

從同震形變場（圖3）中可以看出，該地震引起的地表形變總體表現(xiàn)為橢圓狀分布，其NNW向長約40km，NEE向長約20km。除了部分被冰雪和湖泊所覆蓋的區(qū)域外，整個干涉形變相位連續(xù)，條紋光滑清晰，特征明顯。在T155A的形變圖（圖3（a））中，視線向（LOS）的最大位移達到5.8cm；而在T205D（圖3（b））中，LOS的最大位移為6.7cm。由斷層兩盤的相對運動模式估計的斷層運動模式與震源機制解相一致。

在InSAR同震地表形變場中，主要包含有軌道誤差、大氣誤差、DEM誤差、熱噪聲和數(shù)據(jù)處理中的配準、重采樣誤差等多種誤差，這些誤差對形變相位的影響方式和量級各不相同。為了估計InSAR同震形變場的精度水平，本文采用1D方差－協(xié)方差函數(shù)［11］：來描述形變相位的誤差特征，其中Ci，j為相距r的像素i和j之間的協(xié)方差函數(shù)，σ2是整個形變場的方差，α是誤差衰減距離。據(jù)此給出的升、降軌InSAR同震形變場的中誤差分別為3.8mm和4.3mm，方差－協(xié)方差衰減距離分別為2.3km和2.6km（表1）。該結(jié)果略優(yōu)于2008年青海大柴旦地震的同震地表形變場的精度水平［4］，接近InSAR自身的精度水平，這可能也意味著在所得到的升、降軌InSAR地表形變場（圖3）中主要包含的是地震同震形變信息。

3 地震反演模型

為了提高反演效率，本文采用基于格林函數(shù)分辨率的降采樣法［12］來對同震形變場（圖3）進行降采樣，并按照得到的采樣點位置來計算實際的衛(wèi)星入射角及其軌道方位角，最終得到了1060個觀測數(shù)據(jù)來反演震源參數(shù)，其中包括503個升軌InSAR觀測數(shù)據(jù)和557個降軌InSAR觀測數(shù)據(jù)。

為了獲取地震的震源參數(shù)，本文采用兩步法策略［13－14］來進行反演分析。首先是基于矩形位錯模型［15］采用非線性反演方法確定斷層的幾何結(jié)構(gòu)；然后通過線性反演方法來估計斷層面上的精細滑動分布。這里特別需要指出的是，對于線性反演而言，基于非線性反演給出的斷層幾何結(jié)構(gòu)（尤其是傾角）也許并不是最優(yōu)［16］，它需要在線性反演中重新進行估計。

本文采用 PSOKINV 程序［17－18］來進行反演，該程序采用多峰值顆粒群（MPSO）優(yōu)化算法來確定斷層幾何參數(shù)，如斷層位置、長度、深度、走向、傾向、滑動角和滑動量等參數(shù)，具有收斂效率高、控制參數(shù)少和程序簡潔易行等特點。最佳擬合的均勻斷層模型參數(shù)見表2，表明發(fā)震斷層是一個以正傾滑為主兼有少量右旋走滑的斷層，其走向近 NW－SE方向，傾角為39.2°±5.0°。斷層的深度為6.0km，意味著該發(fā)震斷層并沒有出露于地表。反演得到的地震矩為1.05×1018Nm，小于GCMT給出的1.54×1018Nm，但大于USGS給出的6.46×1017Nm。

表2 InSAR和地震學(xué)給出的阿里地震發(fā)震斷層參數(shù)Tab.2 Fault parameters for the Ali earthquake from InSAR and seismic resource

在線彈性位錯模型中，當(dāng)斷層的幾何結(jié)構(gòu)確定后，斷層面上的滑動量（走滑和傾滑分量）與地表形變之間呈線性關(guān)系。在確定斷層精細滑動分布模型過程中，本文采用之前均勻滑動模型給出的斷層位置和走向等參數(shù)，但對其中的斷層傾角進行重新估計。將斷層長度延長到30km，斷層面寬度延長到25km，然后將斷層面離散成1km×1km大小的750個斷層片。這些斷層片上的滑動量與觀測值之間的關(guān)系為

式中，G為格林函數(shù)；κ2為光滑因子；L為拉普拉斯二階平滑算子［19］；D為InSAR觀測值；S為待求的滑動量。

由于均勻滑動模型給出的斷層傾角并非是線性反演模型所對應(yīng)的最優(yōu)傾角［16－17］，這里采用綜合考慮模型粗糙度和模型擬合殘差的方法［18］確定最優(yōu)斷層傾角和光滑因子（κ2）。在該方法中，定義目標(biāo)函數(shù)lg（ε＋ψ），其中ε為模型擬合殘差，ψ為滑動分布模型的粗糙度。圖4為斷層傾角、光滑因子（κ2）和目標(biāo)函數(shù)lg（ε＋ψ）的關(guān)系圖，文中選取目標(biāo)函數(shù)最小位置處的斷層傾角（43°）和光滑因子（κ2＝2）作為線性反演模型最優(yōu)斷層傾角和光滑因子。

線性反演給出的阿里地震同震滑動分布模型見圖5（a）。從圖5（a）中可以看出，同震滑動分布主要集中在7～12km深度位置（對應(yīng)10～18km寬度位置），其平均滑動量為0.14m，平均滑動角為－131.7°；最大滑動量位于9km深度位置，為0.31m，滑動角為－111.2°。地震釋放的能量為1.24×1018Nm，相當(dāng)于矩震級 Mw 6.03。為了估計線性反演模型所給出滑動分布的精度水平，文中采用蒙特卡洛誤差傳遞方法［20］來估計模型的精度，得到的同震滑動的誤差分布見圖5（b）。從圖5（b）中可以看到，線性反演的同震滑動的誤差分布比較均勻，其平均誤差為2cm，最大誤差為5.6cm。

圖4 斷層傾角、光滑因子與函數(shù)lg（ε＋ψ）之間的關(guān)系（紅色五角星為全局最優(yōu)位置）Fig.4 Contour map of lg（ε＋ψ）with variations of dips and smoothing factors（The red star indicates the point of global minimum）

圖5 2007年阿里地震同震滑動分布（a）及其誤差（b）Fig.5 Coseismic slip distribution（a）and its errors（b）of the 2007Ali earthquake

圖6顯示的是采用線性反演給出的同震滑動分布模型模擬的升、降軌干涉圖以及殘差分布，擬合結(jié)果表明該滑動分布模型能較好地解釋阿里地震的InSAR同震地表形變場，模型生成的干涉圖非常清晰，殘差結(jié)果也很小，觀測值與擬合值之間的相關(guān)性為96.4%。T155A和T205D的殘差均方誤差分別為2.7mm和3.0mm，與InSAR同震地表形變場的精度相當(dāng)。

圖6 分布滑動模型擬合的升、降軌形變場（（a）和（b））和殘差形變場（（c）和（d））Fig.6 Modeled（upper left，upper right）and residual（bottom left，bottom right）ascending and descending displacement maps

4 結(jié) 論

2007年5月5日阿里Mw 6.1級地震發(fā)生在藏北無人區(qū)，對人類的生產(chǎn)和生活造成的影響非常小，但是它對理解藏北地區(qū)的運動特征有著重要的意義。盡管該地區(qū)的惡劣自然條件使野外地質(zhì)考察難以完全展開，但是本研究同時利用升、降軌道EnviSat衛(wèi)星SAR數(shù)據(jù)獲取了該地震的高精度同震地表形變場。在此基礎(chǔ)上，分別采用MPSO非線性和最小二乘線性方法反演給出了地震的最佳斷層幾何模型和均勻滑動分布，為研究地震破裂對周邊斷層的應(yīng)力狀態(tài)的影響、震后形變機理分析和災(zāi)害評估等研究提供重要數(shù)據(jù)支撐。

地震發(fā)生后，盡管采用全球地震臺網(wǎng)（GSN）記錄的寬頻帶遠震體波數(shù)字資料可以快速地給出粗略的矩張量解和最佳雙力偶解，但是由于在亞洲地區(qū)地震臺站稀少以及青藏高原地區(qū)地殼的不均勻性使得基于地震波資料給出的地震震中位置有著較大的不確定性。本文采用升、降軌InSAR資料給出了2007年阿里地震的精確震中位置，發(fā)現(xiàn)震中位于USGS和NEIC基于地震波資料給出的震中位置的西北5km處，這個結(jié)果類似于2008年于田Mw7.1級地震的結(jié)果［1］。

為了分析不同類型觀測數(shù)據(jù)對反演結(jié)果的影響，本文還分別單獨采用升軌和降軌數(shù)據(jù)來進行同震滑動分布反演，這些不同數(shù)據(jù)反演得到的殘差均方誤差以及滑動量見表3。從表3中可以看到，單獨采用升軌數(shù)據(jù)進行的反演可以很好地擬合升軌自身的觀測值，卻不能較好地擬合降軌觀測值，而基于降軌數(shù)據(jù)的反演結(jié)果也是如此。但是聯(lián)合升、降軌數(shù)據(jù)反演給出的同震滑動分布模型雖然對于單類數(shù)據(jù)的擬合效果不如采用單類數(shù)據(jù)的反演結(jié)果，但是其總體擬合效果要更好，其殘差均方誤差比采用單類數(shù)據(jù)的結(jié)果降低了0.2～0.3mm（6%～10%）。此外，不同數(shù)據(jù)反演給出的最大滑動量和平均滑動量大小基本一致。

表3 不同數(shù)據(jù)反演給出的殘差和滑動量Tab.3 The residuals and slips from inversions with different observations

InSAR地表形變觀測表明2007年阿里地震形成了40km×20km大小的形變區(qū)，引起約4cm的下沉。反演結(jié)果表明發(fā)震斷層是美馬錯—德克瑪洛斷裂的次生斷層，走向158°，傾角43°，傾向西南向，以正斷層為主，兼有少量右旋走滑分量。分布式滑動分布模型顯示同震滑動分布主要集中在7～12km深度位置，而在斷層的上部沒有顯著的滑動分布，存在著明顯的類似于走滑型地震的“淺部滑動缺失”現(xiàn)象［6］，這可能意味著該發(fā)震斷層的上部區(qū)域在將來有著較大的地震可能性。

致謝：感謝中國科技部和歐洲航天局（ESA）聯(lián)合資助項目（龍計劃三期）提供的Envisat衛(wèi)星SAR數(shù)據(jù)（ID：10607）。

［1］ELLIOTT J R，WALTERS R J，ENGLANDP C，et al.Extension on the Tibetan Plateau：Recent Normal Faulting Measured by InSAR and Body Wave Seismology［J］.Geophysics Journal International，2010，183（2）：503－535.

［2］HAN Tonglin.Active Tectonic of Xizang［M］.Beijing：The Geological Publishing House，1987.（韓同林.西藏活動構(gòu)造［M］.北京：地質(zhì)出版社，1987.）

［3］WEN Yangmao，XU Caijun.Ms 7.9Mani Earthquake Slip Distribution Inversion by a Sensitivity－based Iterative Fitting Method［J］.Geomatics and Information Science of Wuhan University，2009，34（6）：732－735.（溫揚茂，許才軍.基于敏感度的迭代擬合法反演瑪尼Ms 7.9級地震滑動分布［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2009，34（6）：732－735.）

［4］WEN Yangmao，XU Caijun，LIU Yang，et al.Source Parameters of 2008Qinghai Dachaidan Mw6.3Earthquake from InSAR Inversion and Automated Fault Discretization Method［J］.Geomatics and Information Science of Wuhan University，2012，37（4）：458－462.（溫揚茂，許才軍，劉洋，等.利用斷層自動剖分技術(shù)的2008年青海大柴旦Mw6.3級地震InSAR反演研究［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2012，37（4）：458－462.）

［5］WANG Leyang，XU Caijun，WEN Yangmao.Fault Parameters of 2008Qinghai Dacaidan Mw 6.3Earthquake from STLN Inversion and InSAR Data［J］.Acta Geodaetica et Cartographica Sinica，2013，42（2）：168－176.（王樂洋，許才軍，溫揚茂.利用STLN和InSAR數(shù)據(jù)反演2008年青海大柴旦Mw 6.3級地震斷層參數(shù)［J］.測繪學(xué)報，2013，42（2）：168－176.）

［6］FIALKO Y，SANDWELL D，SIMONS M，et al.Threedimensional Deformation Caused by the Bam，Iran，Earthquake and the Origin of Shallow Slip Deficit［J］.Nature，2005，435（7040）：295－299.

［7］WERNER C，WEGMULLER，STROZZI T，et al.GAMMA SAR and Interferometric Processing Software［C］∥Proceedings of the ERS－ENVISAT Symposium.Gothenburg：［s.n.］，2000.

［8］MASSONNET D，ROSSI M，CARMONA C，et al.The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry［J］.Nature，364（6433）：138－142.

［9］FARRT G，ROSENP A，CARO E，et al.The Shuttle Radar Topography Mission［J］.Reviews of Geophysics，2007，45（2）：20－24.

［10］WEN Y，XU C，LIU Y，et al.Coseismic Slip in the 2010 Yushu Earthquake（China），Constrained by Wide－swath and Strip－map InSAR［J］.Natural Hazards and Earth System Sciences，2013，13：35－44.

［11］HANSSEN R F.Radar Interferometry：Data Interpretation and Error Analysis［M］.New York：Kluwer Academic Publishers，2001.

［12］LOHMAN R B，SIMONS M.Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation：Noise Structure and Data Down Sampling［J］.Geochemistry，Geophysics，Geosystems，2005，6（1），doi：10.1029/2004GC000841.

［13］WRIGHT T J，LU Z，WICKS C.Source Model for the Mw 6.7，23October 2002，Nenana Mountain Earthquake（Alaska）from InSAR［J］.Geophysical Research Letters，2003，30（18），doi：10.1029/2003GL018014.

［14］WEN Yangmao，HE Ping，XU Caijun，et al.Source Parameters of the 2009L’Aquila Earthquake，Italy from ENVISAT and ALOS Satellite SAR Images［J］.Chinese Journal of Geophysics，2012，55（1）：53－65.（溫揚茂，何平，許才軍，等.聯(lián)合ENVISAT和ALOS衛(wèi)星影像確定L’Aquila地震震源機制［J］.地球物理學(xué)報，2012，55（1）：53－65.）

［15］OKADA Y.Surface Deformation due to Shear and Tensile Faults in a Half－space［J］.Bulletin of the Seismological Society of America，1985，75（4）：1135－1154.

［16］BüRGMANN R，AYHAN M E，F(xiàn)IELDING E J，et al.Deformation during the 12November 1999Düzce，Turkey，Earthquake，from GPS and InSAR Data［J］.Bulletin of the Seismological Society of America，2002，92（1）：161－171.

［17］FENG W P，LI Z H，ELLIOTT J R，et al.The 2011Mw 6.8Burma Earthquake：Fault Constraints Provided by Multiple SAR Techniques［J］.Geophysics Journal International，2013，195（1）：650－660.

［18］FENG Wanpeng，LI Zhenhong.A Novel Hybrid PSO/Simplex Algorithm for Determining Earthquake Source Parameters Using InSAR Data［J］.Progress in Geophysics，2010，25（4）：1189－1196.（馮萬鵬，李振洪.InSAR資料約束下震源參數(shù)的PSO混合算法反演策略［J］.地球物理學(xué)進展，2010，25（4）：1189－1196.）

［19］JONSSON S，ZEBKER H，SEGALL P，et al.Fault Slip Distribution of the 1999Mw7.1Hector Mine，California Earthquake，Estimated from Satellite Radar and GPS Measurements［J］.Bulletin of the Seismological Society of America，2002，92（4）：1377－1389.

［20］PARSONS B，WRIGHT T，ROWE P，et al.The 1994 Sefidabeh（Eastern Iran）Earthquakes Revisited：New Evidence from Satellite Radar Interferometry and Carbonate Dating about the Growth of an Active Fold above a Blind Thrust Fault［J］.Geophysics Journal International，2006，164（1）：202－217.