薛 蓮 孫建寶 沈正康

1)北京大學(xué)地球物理系，北京 100871

2)中國地震局地質(zhì)研究所，地震動力學(xué)國家重點實驗室，北京 100029

2010年1月12日海地MW7.0地震InSAR同震形變觀測及同震滑動分布反演

薛 蓮1，2)孫建寶2)沈正康1，2)

1)北京大學(xué)地球物理系，北京 100871

2)中國地震局地質(zhì)研究所，地震動力學(xué)國家重點實驗室，北京 100029

2010年1月12日GMT時間21時53分，在海地境內(nèi)(72.57°W，18.44°N)發(fā)生了MW7.0地震。文中利用干涉合成孔徑雷達(InSAR)方法獲得了覆蓋整個震區(qū)的高精度形變觀測資料，用以研究該地震的發(fā)震機理。采用ALOSPALSAR數(shù)據(jù)，分析了軌道、大氣等誤差源對干涉信號的影響，最終獲得了雷達視線向(LOS)的同震形變場。基于誤差矯正后的InSAR同震數(shù)據(jù)，反演得到了發(fā)震斷層的幾何參數(shù)及斷層面上的同震滑移分布。結(jié)果顯示斷層為N傾37°，小于USGS給出的斷層傾角70°。為了測試斷層傾角的最佳估計，另外建立了傾角為70°的平面斷層模型，以及一種鏟形斷層模型。通過將不同模型預(yù)測的InSAR干涉圖和InSAR同震觀測數(shù)據(jù)對比，認(rèn)為傾角為37°的斷層模型更為合理。斷層同震滑動主要分布在4～16km深度范圍內(nèi)，最大同震滑移量達2.8m，深度為7.2km?；瑒又饕憩F(xiàn)為逆沖兼具左旋走滑2種方式，顯示了此次地震非常復(fù)雜的運動特點。據(jù)此模型得到該地震釋放的地震矩為5.64×1019Nm，對應(yīng)矩震級為MW7.1。

海地地震 InSAR同震觀測 同震滑動分布反演

0 引言

2010 年1月12 日，海地發(fā)生了MW7.0地震。震中位置為72.57°W，18.44°N(USGS，2010)，距離首都太子港(Port-au-Prince)25km(圖1)。此次地震的規(guī)模和強度是海地21世紀(jì)以來最大的1次。該地震發(fā)生于人口密集地區(qū)，造成的死亡人數(shù)達222570人(USGS，2010)。

海地位于加勒比海北部伊斯帕尼奧拉(Hispaniola)島的西半部，地處加勒比板塊和北美板塊間的俯沖帶(圖1)。加勒比板塊向N70°E方向，以18～20mm/a的速度向北美板塊移動(Mann et al.，1995)。北美板塊和加勒比板塊邊界的運動分別以板塊邊界的走滑和匯聚為主。伊斯帕尼奧拉群島北部的Puerto Rico trench區(qū)域和南部的Muertos區(qū)域，是以俯沖為主的斷層活動帶，代表了北美板塊和加勒比板塊邊界間SW向的匯聚。伊斯帕尼奧拉島的北部，是EW走向的Septentional左旋走滑斷層(SF)，南部是EW走向的Enriquillo左旋走滑斷層(EF)。這2條走滑斷層代表了北美板塊和加勒比板塊間以走滑為主的特征(Mann et al.，1995)。根據(jù)一些地質(zhì)調(diào)查結(jié)果及GPS觀測結(jié)果，SF的走滑速率約為9mm/a(Mann et al.，2002;Prentice et al.，2003);EF斷層缺少地質(zhì)調(diào)查的滑動速率，GPS觀測數(shù)據(jù)顯示其走滑速率約為10mm/a(Mann et al.，2002)。

圖1 海地地質(zhì)構(gòu)造背景圖Fig.1 Tectonic setting of Haiti.

發(fā)震斷層EF在最近的幾十年中并沒有發(fā)生較大的地震，有關(guān)的地震記錄表明1751年、1770年和1860年發(fā)生的地震和EF斷層有關(guān)，但是沒有相關(guān)的地質(zhì)調(diào)查記錄(USGS，2010)。本次海地7.0級地震發(fā)生前，Manaker等由GPS觀測得到EF的走滑速率以及地震活動性，推知該斷層的震間應(yīng)力積累可以孕育約MW7.2地震(Manaker et al.，2008)。

震后全球科學(xué)家動用各種手段對此次地震進行了快速分析，其中InSAR影像的分析起到了重要的作用。本文利用干涉合成雷達技術(shù)(InSAR)獲得了本次地震的同震形變場，為深刻理解此次地震的形變分布和發(fā)震機制，以及未來地震災(zāi)害的演化趨勢提供了重要的觀測資料。由于發(fā)震斷層EF臨近海岸，斷層上盤局部處于海面以下，所以InSAR影像中缺失了EF斷層北側(cè)的信號。但是陸內(nèi)形變場都有比較好的觀測結(jié)果。本文基于InSAR同震形變觀測結(jié)果，采用貝葉斯概率全空間反演理論，反演獲得了這次地震的斷層幾何模型和同震滑動分布模型。

1 數(shù)據(jù)及方法

1.1 干涉雷達(InSAR)觀測方法

InSAR成像具有單方向的特點，其觀測到的相位變化是地面各個方向的形變在雷達視線向(Line-of-Sight，LOS)上的投影。由于不同方向的形變在LOS方向上的投影大小不同，存在干涉雷達的LOS向模糊問題。雷達干涉中，垂直地表位移dup、北向位移dn、東向位移de在LOS向位移dLOS上的投影(圖2)可以近似用下面的公式表示：

式(1)中dLOS以遠(yuǎn)離衛(wèi)星的方向為正，θ是射線的入射角，αh是方向角。

圖2 干涉雷達3維成像幾何(升軌)(a)與干涉雷達水平投影(升軌)(b)Fig.2 3D imaging geometry of the interferometric radar(ascending)(a)and horizontal projection of the interferometric radar(ascending)(b).

由于InSAR衛(wèi)星入射角一般＜40°，可以看出，InSAR對于垂向形變的敏感性最高，對于SN向和EW向形變的敏感性次之，且與方向角αh有關(guān)。本研究采用的ALOSPALSAR入射角約為38°，比歐洲航天局衛(wèi)星(ERS，Envisat等)的常用入射角(23°)大。根據(jù)式(1)可知，ALOS PALSAR對水平方向的靈敏度比歐洲航天局?jǐn)?shù)據(jù)源稍好一些。

1.2 InSAR數(shù)據(jù)源和InSAR數(shù)據(jù)處理

本研究主要利用JAXA的ALOSPALSAR數(shù)據(jù)(L波段)獲取本次地震的同震形變場。干涉數(shù)據(jù)源的選取需要考慮較短的垂直基線，以降低空間去相干效應(yīng)，同時增大高度模糊，使得剩余地形的影響較小(孫建寶等，2007a)。本文采用的震后數(shù)據(jù)源中，最早獲取時間為震后2d，最晚獲取時間為震后16d。同震干涉信息中包含一定的震后形變信息，但是本次地震數(shù)據(jù)獲取離發(fā)震時刻都比較近，因此震后形變信息可以忽略，對同震形變觀測的影響不大。本文選取了3對升軌干涉數(shù)據(jù)源和1對降軌干涉數(shù)據(jù)源(表1)，基本上覆蓋了整個發(fā)震區(qū)域。其中2009年9月12日的數(shù)據(jù)為FBD模式，其它數(shù)據(jù)均為FBS模式。

本研究使用JPL開發(fā)的Repeat Orbit Interferometry Package(ROI_PAC)處理軟件，利用2路差分干涉方法處理得到了同震干涉形變場 (圖3)。所用2010年數(shù)據(jù)沒有精確軌道，所以在纏繞的干涉圖中可以看到清晰的軌道剩余相位。在后期的誤差校正中必須去掉這部分軌道剩余相位，避免過大的誤差對模型反演造成較大的影響。

表1 干涉使用的ALOSPALSAR衛(wèi)星數(shù)據(jù)Table 1 The ALOSPALSAR data used to generate interferograms

圖3 海地MW 7.0地震纏繞的InSAR同震形變場Fig.3 The wrapped InSAR coseismic deformation field of the Haiti MW 7.0 earthquake.

為了增強數(shù)據(jù)信噪比，我們對纏繞的干涉數(shù)據(jù)做多視處理，以抑制噪音的影響。多視處理后，每個像元的大小約為153.8m×153.8m。用干涉相關(guān)性數(shù)據(jù)做掩膜，去除相干性較低的區(qū)域。采用全局優(yōu)化方法進行相位解纏(Chen et al.，2002)，并自動矯正大多數(shù)的相位跳躍，使整個圖像的相位跳躍最少;最后分析解纏誤差，進一步手工矯正相位跳躍(孫建寶等，2007a)，最終得到解纏的同震形變場(圖5)。由于發(fā)震斷層EF臨海，所以海岸線以外的信號丟失，而EF斷層南側(cè)有比較完整的形變觀測結(jié)果。

1.3 InSAR相位誤差矯正

干涉雷達的信號中通常包括式(2)的幾項內(nèi)容：

式(2)中φdef是形變相位，即地表形變在LOS向引起的相位變化;φatm是大氣噪音，即大氣延遲作用引起的相位誤差;Δφorb是由軌道誤差引起的剩余相位;Δφε是高程剩余相位，即由高程模型DEM中的誤差引起的相位;φn是系統(tǒng)噪音，是處理過程中引起的誤差相位。干涉圖的相位由上述幾項相位構(gòu)成，其中φdef形變相位是有用相位，其他相位的貢獻都被視為形變測量中的誤差源。

圖4 大氣垂直分層信號的特征Fig.4 Characteristics of the vertical stratification signal in the atmosphere.

其中形變相位φdef、大氣延遲相位φatm、軌道誤差相位Δφorb都具有空間相關(guān)性，但是各自的空間相關(guān)尺度不一樣，可以利用其空間相關(guān)尺度的不同將這幾種信號區(qū)分開(Hanssen，2001;孫建寶等，2007a)。軌道信號具有長波特征，其尺度大于區(qū)域性的形變尺度。同震形變場主要集中在近場，因此容易區(qū)分出近場和遠(yuǎn)場信號，并利用遠(yuǎn)場信息擬合軌道誤差信號和大氣誤差信號。高程剩余相位Δφε與垂直基線和高程模型有關(guān)。本文采用精度為1″的ASTER DEM的高程模型，引入的誤差相對較小。而且在同震信號中，高程剩余誤差的信號遠(yuǎn)小于同震信號，所以不影響同震形變場的分布。系統(tǒng)噪音φn是信號噪音、插值誤差和配準(zhǔn)誤差等引入的相位誤差(Hanssen，2001)，在相干性較好的數(shù)據(jù)中，系統(tǒng)噪音對干涉相位的貢獻非常小，一般可以忽略不計。本文主要分析由軌道誤差和大氣延遲中垂直大氣分層帶來的相位誤差。

1.3.1 軌道誤差的矯正

由于本文采用干涉源的基線較長，而且2010年后的數(shù)據(jù)軌道信息不精確，所以原始干涉數(shù)據(jù)中會有明顯的剩余軌道信號 (圖3)。衛(wèi)星軌道二階軌道系數(shù)和一階軌道系數(shù)的比值約為10-6(Hanssen，2001)，而本文的觀測區(qū)域約為200km×200km，空間跨越尺度不大，所以本文采用一階線性模型矯正軌道剩余，忽略非線性項的影響。線性軌道項可以表示為

式(3)中x，y分別是方位向(azimuth)和距離向(rang)的像元坐標(biāo)。

圖5 原始干涉解纏圖和誤差矯正圖Fig.5 The original unwrapped interferograms and the error corrected images.

選定觀測區(qū)域中的遠(yuǎn)場或無形變區(qū)，假設(shè)其觀測相位主要由軌道信號、大氣信號和隨機噪聲等誤差信號構(gòu)成。大氣信號在空間分布上不存在線性關(guān)系，所以在擬合軌道信號時，大氣信號視為噪音。在大氣信號噪音不顯著的時候，對剩余軌道的估計影響不大。利用最小二乘方法，可以擬合得到線性軌道模型參數(shù)a，b，c。

在形變觀測中，若形變特征是區(qū)域性的，形變信號和遠(yuǎn)場或無形變區(qū)能夠明顯區(qū)分時，這樣的剩余軌道矯正不會干擾真實的形變場。但是在形變區(qū)的信號沒有明顯的遠(yuǎn)場或無形變區(qū)的特征時，采用這樣的方法消除剩余軌道信號會使形變信號發(fā)生畸變，偏離真實的形變信號(Biggs et al.，2007)。在同震觀測中形變信號在遠(yuǎn)離發(fā)震區(qū)明顯衰減，所以存在區(qū)分于發(fā)震區(qū)的遠(yuǎn)場信息，利用這樣的剩余軌道矯正，不會使同震形變場發(fā)生畸變。

1.3.2 大氣垂直分層誤差的矯正

由大氣延遲引起的相位誤差在干涉圖中比較難估計和矯正。大氣延遲作用主要由2種物理過程引起，分別為大氣湍流混合(turbulentmixing)與垂直分層(vertical stratification)(Hanssen，2001)。其中大氣湍流混合使得大氣反射率各向異性，引起的大氣延遲具有任意性。如果沒有比較密集的GPS觀測臺站和其它光學(xué)遙感數(shù)據(jù)的配合，這部分信號很難用數(shù)學(xué)模型精確模擬。海地沒有密集的GPS觀測臺站以及相關(guān)的大氣模型，所以本文不做由大氣湍流混合引起的大氣相位矯正。由大氣垂直分層引起的大氣延遲相位和地形相關(guān)，在多山地的區(qū)域，或者大氣濕度比較大的區(qū)域，這部分信號的影響較顯著(Cavaliéet al.，2008;Hanssen，2001)。海地同震形變觀測中，相位和地形有很好的線性相關(guān)性 (圖4)，所以本文采用線性模型消除大氣垂直分層信號的影響。大氣垂直分層信號可以表示為

式(4)中φ'atm為大氣垂直分層信號，z為觀測區(qū)域內(nèi)地形的相對高度;k為大氣垂直分層信號與地形的相關(guān)系數(shù)，整幅干涉圖中假設(shè)都為同1個常數(shù)。

這部分信號雖然也具有空間相關(guān)性，但是它和軌道及形變的空間相關(guān)尺度不一樣，所以可以分離出這部分信號。本文采用的數(shù)據(jù)，大氣信號和地形高程有較好的線性相關(guān)性，所以采用線性模型消除大氣垂直信號對形變場沒有影響。本文采用1″Aster DEM作為地形模型，選觀測區(qū)域的地形最高點作為參考點，求出觀測區(qū)域內(nèi)相對最高點的地形變化，再采用最小二乘法，和式(4)給出的線性模型消除大氣垂直分層信號(圖3)。

從誤差矯正的結(jié)果中可以看出，軌道信號和大氣垂直分層信號的影響比較顯著，尤其是軌道信號的影響，這與本文采用的軌道數(shù)據(jù)有關(guān)。這些誤差源完全淹沒了形變場信號的特征(圖5，6)，從原始的干涉圖像中較難識別同震形變場的特征，因此必須對數(shù)據(jù)進行誤差矯正后才能對形變場特征進行分析。

2 觀測結(jié)果及定性分析

本文得到此次地震纏繞的同震形變場(圖3)，是最原始的雷達數(shù)據(jù)處理結(jié)果，沒有經(jīng)過誤差矯正，不能真實地反映同震形變場。但是可以根據(jù)干涉條紋的特點，定性判斷該地震的同震形變場特征。

升軌數(shù)據(jù)覆蓋了沿EF斷層約200km長的范圍，基本覆蓋了整個發(fā)震區(qū)域。降軌數(shù)據(jù)沒有升軌軌道覆蓋范圍廣，但基本包含了震中區(qū)域的形變場。升軌InSAR觀測圖中，圖3b軌道的條紋比其他幾幅干涉圖的條紋密集，這主要是剩余軌道信號，與形變無關(guān)。干涉條紋中明顯與整體條紋走向不一致的弧形條紋，是由發(fā)震區(qū)域內(nèi)的地表形變引起的，為此次地震的近場范圍。

圖6 經(jīng)過軌道誤差、大氣延遲誤差矯正之后的解纏InSAR同震形變場Fig.6 Unwrapped InSAR coseismic deformation field after orbital errors and atmospheric errors removed.

本文采用了ALOSPALSAR數(shù)據(jù)，其中3對為升軌數(shù)據(jù)，1對為降軌數(shù)據(jù)。雷達的中心入射角均為38°，升軌的方向角為-12.27°，降軌的方向角為-167.73°。根據(jù)雷達成像的幾何特點，分別得到了地表位移對雷達視線LOS向位移的貢獻。

由式(5)、(6)可知，無論是升軌還是降軌，地面垂直位移和EW向位移對LOS向位移的貢獻大小較為接近，存在相互抵消的情況。

矯正誤差后的解纏同震形變場(圖6)給出了基本的同震形變特征。觀測區(qū)域的最西邊a軌道的干涉數(shù)據(jù)覆蓋了震中以西的區(qū)域;從圖中可以看出Miragoane區(qū)域以西地方的形變幾乎為零。觀測區(qū)域最東邊的c軌道的干涉圖覆蓋了太子港以東的區(qū)域。該觀測區(qū)域中除了太子港附近外，其他地方幾乎沒有形變。由此推知同震形變場的分布范圍，西側(cè)不超過Miragoane區(qū)域，東側(cè)的影響范圍不超過太子港。這與余震的分布范圍一致(USGS，2010)。

發(fā)震斷層是向N傾的左旋逆沖斷層(USGS，2010)。發(fā)震斷層上盤的走滑分量和逆沖分量對升軌LOS向位移的貢獻都為負(fù)，下盤的走滑分量和逆沖分量對升軌LOS向的貢獻都為正，這與我們的觀測結(jié)果是一致的(圖6a，b，c)。而發(fā)震斷層上盤的走滑分量對降軌LOS向位移的貢獻為正，逆沖分量對降軌LOS向位移的貢獻為負(fù);下盤的走滑分量對降軌LOS向位移的貢獻為負(fù)，逆沖分量對降軌LOS向位移的貢獻為正。因此，發(fā)震斷層的走滑分量和逆沖分量對降軌LOS向位移的貢獻在一定程度上是相互抵消的。城市Léogane和Petit Goave之間區(qū)域的位移，在升軌、降軌的LOS向上的投影都為負(fù)，說明該地區(qū)的運動以逆沖為主。而在城市Léogane附近的地表運動在升軌LOS向的位移為負(fù)，在降軌LOS向的位移為正 (圖6b，d)，說明該區(qū)域運動以W向走滑為主，這和野外地質(zhì)調(diào)查結(jié)果一致(Prentice et al.，2010)。

3 反演同震模型

3.1 反演方法

本文利用獲得的InSAR觀測的同震形變場，基于彈性半空間理論，利用貝葉斯概率全空間反演方法(fully Bayesian inversion)反演斷層幾何參數(shù)及同震滑移分布(Fukuda et al.，2008)。先將斷層劃分成若干個子斷層，假設(shè)每個子斷層是均勻滑動的，不同子斷層具有非均勻滑動，彈性介質(zhì)的泊松比為0.25。斷層幾何參數(shù)、每個子斷層的滑動量和觀測數(shù)據(jù)之間的關(guān)系為

式(7)中d是觀測的同震形變場;G是將觀測值和斷層幾何參數(shù)、斷層滑移量聯(lián)系起來的格林函數(shù);m是斷層幾何參數(shù)，它和地表形變之間是非線性關(guān)系;s是每個子斷層的滑動分量;ε是觀測的誤差矩陣，假設(shè)其滿足高斯正態(tài)分布ε～N(0，Σd)，其中Σd是觀測數(shù)據(jù)協(xié)方差矩陣，本文采用對角矩陣作為協(xié)方差矩陣。

斷層劃分得越細(xì)，待求的未知量就會越多，反演就會越不穩(wěn)定。為了增強反演解的穩(wěn)定性，加入先驗約束方程，引入目標(biāo)函數(shù)

式(8)中對滑動量s采用離散的拉普拉斯平滑(Jonsson et al.，2002)，L是離散拉普拉斯算子;β2是光滑因子，決定滑動量s的平滑約束程度。

常采用的一種反演方法是選定合理的光滑因子β2，利用非線性反演法求使目標(biāo)函數(shù)Ф(式(8))最小時的斷層幾何參數(shù)和斷層滑動分布(孫建寶等，2007c;Sun et al.，2008;萬永革等，2008)。本文采用貝葉斯概率全空間反演法，并用Monte Carlo算法對解空間進行樣本采樣，求解斷層幾何參量、滑動分布和平滑因子β2及數(shù)據(jù)權(quán)重σ2的后驗概率分布 P(s，m，β2，σ2|d)(Fukuda et al.，2008)。P(s，m，β2，σ2|d) 表示在已知觀測形變量 d 的情況下，未知量 s，m，β2，σ2的概率分布。從未知量的概率分布可以計算出未知量s，m，β2，σ2的平均值和標(biāo)準(zhǔn)差。

3.2 形變場離散采樣

干涉雷達獲得的形變場是連續(xù)的，一副干涉圖像中包含了大約106個數(shù)據(jù)點，并且像元之間具有較強的相關(guān)性，若利用全部的觀測數(shù)據(jù)點進行反演，會造成計算的繁冗。本文采用四叉樹離散采樣法(Jonsson et al.，2002)對觀測數(shù)據(jù)進行離散采樣。由于觀測數(shù)據(jù)中，c軌道干涉圖(圖6c)中的大氣信號非常顯著，覆蓋發(fā)震區(qū)域較少，所以反演時不采用這對觀測數(shù)據(jù)。采樣后的觀測數(shù)據(jù)，代表了整個形變場的特征(圖7)，但是數(shù)據(jù)點的個數(shù)降為9,858個，有效地提高了計算效率，并降低了噪音對形變場的影響(孫建寶等，2007b)。為了消除同一條軌道入射角變化帶來的誤差，利用LOS形變場的采樣方式，對相應(yīng)的入射角取平均，得到了數(shù)據(jù)點對應(yīng)的觀測矢量 (圖7)。

3.3 反演模型及結(jié)果

3.3.1 模型1：貝葉斯概率全空間反演模型

根據(jù)同震形變的分布范圍，我們設(shè)定1個長76km，寬33km的斷層模型。將斷層沿走滑方向劃分為19列，沿傾向劃分為11層(圖8)。共有209個子斷層，每個子斷層的尺度是4km×3km，每個子斷層皆均勻滑動。利用Okada矩形元位錯模型(Okada，1985)生成格林函數(shù)矩陣G。反演過程中對斷層滑動量加以非負(fù)約束，使得每個子斷層的滑動分布滿足左旋走滑加逆沖的滑動特點。固定斷層中點的經(jīng)度為-72.6°，利用貝葉斯概率全空間反演觀測數(shù)據(jù)權(quán)重σ2，平滑因子β2，斷層中點的緯度，斷層傾角、走向以及每個子斷層的走滑分量和傾滑分量。

利用Monte Carlo算法對未知參數(shù)解進行空間采樣，計算出未知參數(shù)的后驗概率P(s，m，β2，σ2|d)分布。迭代過程中前150次迭代沒有收斂到最佳解區(qū)域，所以舍棄前150個樣本點。利用剩下的400個樣本來估計未知參量s，m，β2，σ2的概率分布 (圖8)?？梢钥闯雒總€未知量的樣本空間都有很好的采樣。反演得到的斷層中點的緯度為18.6°，傾角為37°，走向為260.0°，觀測數(shù)據(jù)的權(quán)重σ2為0.059，平滑因子β2為20.02。得到的斷層滑動分布、傾向分量標(biāo)準(zhǔn)差及走滑分量標(biāo)準(zhǔn)差如圖9所示。

參考USGS給出的震源機制解，發(fā)現(xiàn)反演得到的斷層傾角37°遠(yuǎn)小于USGS給出的斷層傾角70°。為了判斷傾角的合理性，我們另外建立了2個模型。

3.3.2 模型2：固定斷層傾角模型

為了判斷地震學(xué)方法給出的斷層傾角的合理性，我們建立了斷層傾角為70°的斷層模型。根據(jù)不同模型與觀測數(shù)據(jù)的擬合程度判斷模型的合理性。基于貝葉斯反演得到的結(jié)果，固定數(shù)據(jù)權(quán)重σ2為0.059，平滑因子β2為20.02。將斷層傾角固定為70°，其他的斷層幾何參數(shù)、斷層空間位置以及劃分方法和模型1一致。為了使滑動分布滿足左旋逆沖的特點，采用快速非負(fù)最小二乘算法(Bro et al.，1997)求算斷層的滑動分布 (圖10)。

3.3.3 模型3：鏟形斷層模型

因為模型1得到的斷層傾角37°遠(yuǎn)小于地震學(xué)給出的斷層傾角，推測發(fā)震斷層可能是鏟形斷層。鏟形斷層接近地表的傾角比較陡峭，而下面的傾角則比較平緩，如果使用1個斷層面擬合，會使得斷層傾角過小。因此，我們建立了上下連接的2個不同傾角的斷層面模擬鏟形斷層。斷層在空間中的位置及走向和前面2個模型一致。設(shè)定上下斷層的長度都為76km，將上部斷層的傾角固定為70°，反演求解上部斷層的寬度及下部斷層的傾角。下部斷層的寬度固定為24km，上、下斷層都分別沿走向劃分為19列，沿傾向劃分為6層，共得到228個子斷層。下部斷層的每個子斷層尺度為4km×4km，上部斷層的子斷層的長度為4km，寬度和上部斷層的總寬度有關(guān)。對上部斷層的不同寬度和下部斷層的不同傾角求取反演模型和觀測值之間的殘差。當(dāng)上部斷層寬度為14km，下部斷層傾角為0°時，擬合殘差最小。然后求解上部斷層寬度為14km，下部斷層傾角為0°時的滑動分布(圖10)。

圖7 干涉軌道形變場及其入射角離散采樣Fig.7 InSAR deformation field and its incident angle sampling.

圖8 貝葉斯反演的未知參量概率分布圖Fig.8 Probability distributions of unknown parameters estimated by Bayesian inversion.

3.4 討論

基于上述3種斷層模型及對應(yīng)的滑動分布，分別正演得到了地表形變在LOS向的位移(圖11)。分別與原始觀測數(shù)據(jù)對比，得到擬合殘差圖(圖12)。從擬合殘差圖中可以看出，3種模型與升軌a，b的觀測數(shù)據(jù)擬合都比較好，殘差都沒有系統(tǒng)性的偏差。與降軌d的觀測數(shù)據(jù)擬合，遠(yuǎn)場處擬合較好，而在斷層附近的區(qū)域存在比較明顯的剩余。但傾角為37°的斷層模型，在斷層附近的擬合殘差比其它2個模型的擬合殘差小。

通過以上3種模型的對比，我們認(rèn)為傾角為37°的斷層模型更為合理，這與利用地震波形數(shù)據(jù)反演得到的70°的斷層傾角(Sladen，2010)存在較大的差別。但是根據(jù)Eric Calais等給出的此次地震的GPS同震形變場(Calais et al.，2010)，發(fā)震斷層的遠(yuǎn)場GPS位移以垂直斷層方向匯聚為主。發(fā)震斷層近場除了Léogane附近的GPS點外，其它地方的GPS位移場有非常顯著的逆沖分量。其中靠近斷層處的最大位移達到了約50cm，且垂直斷層的逆沖分量明顯大于平行斷層的走滑分量，說明發(fā)震斷層的滑動應(yīng)是以逆沖為主。在斷層傾角為70°的情況下，通常很難產(chǎn)生大的逆沖滑動。根據(jù)USGS給出的幾個震源機制解(表2)，其中中心矩張量矩給出的斷層為S傾(表2中的USGS1)，和實際斷層傾向不一致，所以地震波數(shù)據(jù)得到的斷層面解有一定的不穩(wěn)定性。綜合幾方面的分析，我們認(rèn)為利用InSAR觀測數(shù)據(jù)反演得到的37°斷層傾角更為可信。

圖9 傾角為37°的斷層模型滑動及其標(biāo)準(zhǔn)差隨經(jīng)度、緯度、深度的分布Fig.9 Slip distribution and its standard deviationswith longitude，latitude，and depth of themodel(dipping angle 37°).

圖10 傾角為70°的斷層模型和鏟形斷層模型的滑動分布Fig.10 Slip distributions of themodelwith a dipping angle of70°and the listricmodel.

斷層傾角為37°的斷層模型反演得到的地震矩為5.64×1019Nm，相當(dāng)于矩震級MW7.1，比USGS給出4.5×1019Nm略大，這可能是因為震后形變對InSAR同震形變的影響。從其反演得到的滑動分布，可以看到有2處的滑動量較大，東邊的最大滑動處以左旋走滑為主;而在西邊最大滑移處接近于純逆沖。在深度7.2km處滑動量最大，約為2.8m，與地震波資料反演得到的同震模型(Sladen，2010)的最大滑動量的深度8km非常接近。斷層最上層接近地表處的滑動非常小，顯示破裂沒有出露地表，這與地質(zhì)考察沒有發(fā)現(xiàn)地表破裂一致(Prentice et al.，2010)。儀器震中附近(圖6中的紅色五角星處)和Petit Goave之間的區(qū)域斷層滑動為逆沖兼具走滑性質(zhì);Petit Goave以西區(qū)域的滑動接近于純逆沖。Léogane以北區(qū)域的西向運動特征 (圖6d中Léogane附近的藍(lán)色區(qū)域)受到斷層深部走滑特征的影響，而Léogane和Petit Goave區(qū)域的抬升運動主要受斷層淺層的逆沖滑動的影響。降軌LOS向擬合結(jié)果顯示近場處有明顯的殘差(圖11)，這可能是因為模型沒有完全解析出近場形變的復(fù)雜性，也可能因為這個區(qū)域的形變不完全是由斷層滑動引起的。有關(guān)的地質(zhì)調(diào)查表明這個區(qū)域有土壤液化發(fā)生(Calais et al.，2010)，這需要進一步的形變調(diào)查數(shù)據(jù)來驗證。

圖11 同震形變場的InSAR觀測結(jié)果和不同斷層模型正演結(jié)果的對比Fig.11 Comparisons of InSAR observations and the displacements predicted by the three slip models.

表2 USGS震源機制解Table 2 USGS focalmechanism solutions

4 結(jié)論

本文采用軌道誤差矯正和垂直大氣分層誤差矯正后的InSAR同震形變場，定量反演得到了斷層幾何參數(shù)以及斷層同震滑動分布。反演得到的斷層傾角為37°，小于地震波反演給出的70°斷層傾角。通過與斷層傾角為70°的斷層模型以及鏟形斷層模型對比，發(fā)現(xiàn)斷層傾角在37°時觀測數(shù)據(jù)擬合最好。同震GPS位移場顯示斷層近場以及遠(yuǎn)場都有很明顯的逆沖分量，而在較大斷層傾角的情況下，很難產(chǎn)生這么大的逆沖分量。此次地震中InSAR衛(wèi)星觀測的應(yīng)用，又一次展示了衛(wèi)星大地測量手段對于地震科學(xué)研究的重要意義。

致謝：本研究所用ALOS PALSAR數(shù)據(jù)由日本JAXA提供，版權(quán)歸日本METI和JAXA所有。

孫建寶，徐錫偉，石耀霖，等.2007a.東昆侖斷裂瑪尼段震間形變場的InSAR觀測及斷層滑動率初步估計［J］.自然科學(xué)進展，17(10)：1361—1370.

SUN Jian-bao，XU Xi-wei，SHIYao-lin，et al.2007a.The InSAR observation of interseismic deformation of the Mani section of the Eastern Kunlun Fault and preliminary estimate of the slip rate［J］.Progress in Nature Science，17(10)：1361—1370(in Chinese).

孫建寶，徐錫偉，沈正康，等.2007b.基于線彈性位錯模型及干涉雷達同震形變場反演1997年瑪尼MW7.5級地震參數(shù)-Ⅰ：均勻滑動反演［J］.地球物理學(xué)報，50(4)：1097—1110.

SUN Jian-bao，XU Xi-wei，SHEN Zheng-kang，et al.2007b.Parameter inversion of the 1997 Mani earthquake from In-SAR co-seismic deformation field based on linear elastic dislocationmodel-Ⅰ：Uniform slip inversion［J］.Chinese JGeophys，50(4)：1097—1110(in Chinese).

孫建寶，徐錫偉，沈正康，等.2007c.基于線彈性位錯模型及干涉雷達同震形變場反演1997年瑪尼MW7.5級地震參數(shù)-Ⅱ：滑動分布反演［J］.地球物理學(xué)報，50(5)：1390—1397.

SUN Jian-bao，XU Xi-wei，SHEN Zheng-kang，et al.2007c.Parameter inversion of the 1997 Mani earthquake from In-SAR co-seismic deformation field based on linear elastic dislocation model-Ⅱ：Slip distribution inversion［J］.Chinese JGeophys，50(5)：1390—1397(in Chinese).

萬永革，沈正康，王敏，等.2008.根據(jù)GPS和InSAR數(shù)據(jù)反演2001年昆侖山口西地震同震破裂分布［J］.地球物理學(xué)報，51(4)：1074—1084.

WAN Yong-ge，SHEN Zheng-kang，WANG Min，et al.2008.Coseismic slip distribution of the 2001 Kunlun mountain pass west earthquake constrained using GPSand InSAR data ［J］.Chinese JGeophys，51(4)：1074-1084(in chinese).

Biggs J，Wright T，Lu Z，et al.2007.Multi-interferogram method formeasuring interseismic deformation：Denali Fault，Alaska［J］.Geophysical Journal International，170：1165—1179.

Bro R，De Jong S.1997.A fast non-negativity-constrained least squares algorithm ［J］.Journal of Chemometrics，11：393—401.

Calais E，F(xiàn)reed A，Mattioli G，et al.2010.Transpressional rupture of an unmapped fault during the 2010 Haiti earthquake［J］.Nature Geoscience，3：794—799.

CavaliéO，Lasserre C，Doin M，et al.2008.Measurement of interseismic strain across the Haiyuan Fault(Gansu，China)，by InSAR［J］.Earth and Planetary Science Letters，275：246—257.

Chen C，Zebker H.2002.Phase unwrapping for large SAR interferograms：Statistical segmentation and generalized network models［J］.IEEE Transactions on Geoscience and Remote Sensing，40：1709—1719.

Fukuda J，Johnson K.2008.A fully Bayesian inversion for spatial distribution of fault slip with objective smoothing［J］.Bulletin of the Seismological Society of America，98：1128.

Hanssen R.2001.Radar interferometry：Data Interpretation and Error Analysis［M］.Kluwer Academic Pub，Netherland.

Jonsson S，Zebker H，Segall P，et al.2002.Fault slip distribution of the 1999 MW7.1 Hector Mine，California，earthquake，estimated from satellite radar and GPSmeasurements［J］.Bulletin of the Seismological Society of America，92：1377.

Manaker D，Calais E，F(xiàn)reed A，et al.2008.Interseismic plate coupling and strain partitioning in the northeastern Caribbean ［J］.Geophysical Journal International，174：889—903.

Mann P，Calais E，Ruegg J，etal.2002.Oblique collision in the northeastern Caribbean from GPSmeasurements and geological observations［J］.Tectonics，21：1057.

Mann P，Taylor F，Edwards R，et al.1995.Actively evolvingmicroplate formation by oblique collision and sidewaysmotion along strike-slip faults：An example from the northeastern Caribbean platemargin ［J］.Tectonophysics，246：1—69.

Okada Y.1985.Surface deformation due to shear and tensile faults in a half-space［J］.Bulletin of the Seismological Society of America，75：1135.

Prentice C，Mann P，Pena L，etal.2003.Slip rate and earthquake recurrence along the central Septentrional Fault，North American-Caribbean plate boundary，Dominican Republic［J］.Journal of Geophysical Research，108：2149.

Prentice C，Mann P，Crone A，et al.2010.Seismic hazard of the Enriquillo-Plantain Garden Fault in Haiti inferred from palaeoseismology［J］.Nature Geoscience，3(11)：789—793.

Sladen A.2010.Preliminary result01/12/2010(MW7.0)，Haiti［R］.http：∥www.tectonics.caltech.edu/slip_history/2010_haiti/index.html.

Sun J，Shen Z，Xu X，et al.2008.Synthetic normal faulting of the 9 January 2008 Nima(Tibet)earthquake from conventional and along-track SAR interferometry［J］.Geophysical Research Letters，35：L22308.

USGS.2010.Magnitude 7.0-Haiti region ［R］. http：//earthquake.usgs.gov/earthquakes/eqinthenews12010/us2010rja6.

INSAR COSEISM IC DEFORMATION OBSERVATION OF THE JAN 12TH，2010 HAITIEARTHQUAKE AND ITSCOESEISM IC SLIP DISTRIBUTION INVERSION

XUE Lian1，2)SUN Jian-bao2)SHEN Zheng-kang1，2)
1)Dept of Geophysics，School of Earth and Space Sciences，Peking University，Beijing 100871，China
2)State Key Laboratory of Earthquake Dynamics，Institute of Geology，China Earthquake Administration，Beijing 100029，China

At21：53(GMT)on Jan 12th，2010，an MW7.1 earthquake struck Haitiwith epicenter located at(72.57°W，18.44°N).We derive the high precision Interferometric Synthetic Aperture Radar(In-SAR)deformation field covering the whole rupture zone and use the observations to study the seismogenic kinematics of the quake.We analyze the influences of InSAR data errors，such as the orbiterrors and atmospheric errors，on the ALOSPALSAR data and finally obtain the coseismic deformation field in the line of sight(LOS)direction.We invert the LOS displacement for the fault geometry of the seismogenic fault and its corresponding slip distribution.We find that the fault dips to the north at 37°，lower than the 70°north-dipping solution from USGS.To test other possible dipping angles of the fault plane and find the bestestimate of the faultgeometry，we constructan other two faultmodels.One ismodeled with a 70°north-dipping fault plane and the other one is a listric faultmodelwith the hinge depth inverted.Comparing the predicted InSAR interferogram of the three faultmodels with the observed InSAR coseismic data，we argue that themodel with 37°north dipping angle ismore rational than the other two models.Our preferred model shows that the coseismic slip concentratesmainly at 4 ～16km depth，and themaximum slip is 2.8m，appearing at 7.2km depth.The coseismic displacement exhibits both reversal faulting and sinistral strike-slip，suggesting complex faultmotion in this earthquake.Themoment release of this earthquake is5.64 ×1019Nm，equivalent to an MW7.1 earthquake.

Haiti Earthquake，coseismic InSAR observation，coseismic slip distribution inversion

P315.2

A

0253-4967(2011)01-0157-18

10.3969/j.issn.0253-4967.2011.01.016

2010-06-07收稿，2011-03-08改回。

地震行業(yè)科研專項(200708002)和國家高技術(shù)研究發(fā)展計劃項目(2009AA12Z1464)共同資助。

薛蓮，女，1985年出生，2010年畢業(yè)于北京大學(xué)地球物理系，獲碩士學(xué)位，現(xiàn)就讀于加利福尼亞聯(lián)合大學(xué)圣克魯茲分校，主要研究方向為地殼形變，E-mail：xuelian.icy@gmail.com。