薛 蓮 孫建寶 沈正康

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

2)中國(guó)地震局地質(zhì)研究所，地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2010年1月12日海地MW7.0地震InSAR同震形變觀測(cè)及同震滑動(dòng)分布反演

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

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

2)中國(guó)地震局地質(zhì)研究所，地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

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

海地地震 InSAR同震觀測(cè) 同震滑動(dòng)分布反演

0 引言

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

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

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

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

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

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

1.1 干涉雷達(dá)(InSAR)觀測(cè)方法

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

式(1)中dLOS以遠(yuǎn)離衛(wèi)星的方向?yàn)檎仁巧渚€(xiàn)的入射角，αh是方向角。

圖2 干涉雷達(dá)3維成像幾何(升軌)(a)與干涉雷達(dá)水平投影(升軌)(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對(duì)于垂向形變的敏感性最高，對(duì)于SN向和EW向形變的敏感性次之，且與方向角αh有關(guān)。本研究采用的ALOSPALSAR入射角約為38°，比歐洲航天局衛(wèi)星(ERS，Envisat等)的常用入射角(23°)大。根據(jù)式(1)可知，ALOS PALSAR對(duì)水平方向的靈敏度比歐洲航天局?jǐn)?shù)據(jù)源稍好一些。

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

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

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

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

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

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

1.3 InSAR相位誤差矯正

干涉雷達(dá)的信號(hào)中通常包括式(2)的幾項(xiàng)內(nèi)容：

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

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

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

1.3.1 軌道誤差的矯正

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 反演同震模型

3.1 反演方法

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

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

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

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

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

3.2 形變場(chǎng)離散采樣

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

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

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

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

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

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

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

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

3.3.3 模型3：鏟形斷層模型

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

圖7 干涉軌道形變場(chǎng)及其入射角離散采樣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種斷層模型及對(duì)應(yīng)的滑動(dòng)分布，分別正演得到了地表形變?cè)贚OS向的位移(圖11)。分別與原始觀測(cè)數(shù)據(jù)對(duì)比，得到擬合殘差圖(圖12)。從擬合殘差圖中可以看出，3種模型與升軌a，b的觀測(cè)數(shù)據(jù)擬合都比較好，殘差都沒(méi)有系統(tǒng)性的偏差。與降軌d的觀測(cè)數(shù)據(jù)擬合，遠(yuǎn)場(chǎng)處擬合較好，而在斷層附近的區(qū)域存在比較明顯的剩余。但傾角為37°的斷層模型，在斷層附近的擬合殘差比其它2個(gè)模型的擬合殘差小。

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

圖9 傾角為37°的斷層模型滑動(dòng)及其標(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°的斷層模型和鏟形斷層模型的滑動(dòng)分布Fig.10 Slip distributions of themodelwith a dipping angle of70°and the listricmodel.

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

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

表2 USGS震源機(jī)制解Table 2 USGS focalmechanism solutions

4 結(jié)論

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

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

孫建寶，徐錫偉，石耀霖，等.2007a.東昆侖斷裂瑪尼段震間形變場(chǎng)的InSAR觀測(cè)及斷層滑動(dòng)率初步估計(jì)［J］.自然科學(xué)進(jìn)展，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.基于線(xiàn)彈性位錯(cuò)模型及干涉雷達(dá)同震形變場(chǎng)反演1997年瑪尼MW7.5級(jí)地震參數(shù)-Ⅰ：均勻滑動(dòng)反演［J］.地球物理學(xué)報(bào)，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.基于線(xiàn)彈性位錯(cuò)模型及干涉雷達(dá)同震形變場(chǎng)反演1997年瑪尼MW7.5級(jí)地震參數(shù)-Ⅱ：滑動(dòng)分布反演［J］.地球物理學(xué)報(bào)，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).

萬(wàn)永革，沈正康，王敏，等.2008.根據(jù)GPS和InSAR數(shù)據(jù)反演2001年昆侖山口西地震同震破裂分布［J］.地球物理學(xué)報(bào)，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è)科研專(zhuān)項(xiàng)(200708002)和國(guó)家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2009AA12Z1464)共同資助。

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