基于InSAR同震形變觀測反演2010年新西蘭南島Mw7.1Darfield地震同震破裂分布

羅旭巍，孫建寶，沈正康，＊，Yo Fukushima

1 北京大學地球與空間科學學院，北京 100871

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

3 Disaster Prevention Research Institute，Kyoto University，Uji，Kyoto，Japan

1 引 言

2010年9月4日新西蘭南島Canterbury平原發(fā)生了Mw7.1級地震，震源深度約為10km，震中位于Darfield東南10km處，地震造成2人受傷［1］.但它的余震，2011年2月21日Christchurch Mw6.3級地震造成了180余人死亡［2］.

新西蘭島嶼橫跨太平洋板塊和澳洲板塊（圖1），地殼運動十分活躍.新西蘭南島和北島之間，太平洋板塊以約43mm／a的速度向澳洲板塊俯沖；而在南島的西南端，澳洲板塊以35mm／a的速度插入太平洋板塊下方［3］（圖1）.由于上述兩大板塊匯聚在南島附近且存在斜向運動，新西蘭南島南端形成了一條東西走向的阿爾卑斯走滑斷裂帶，GPS測量數(shù)據(jù)顯示其存在約40mm／a的右旋走滑速率［4－5］.除了阿爾卑斯走滑斷裂帶吸收的震間滑動外，還有約10mm／a的跨邊界帶縮短分布在阿爾卑斯山走滑斷裂帶東南部的山地內(nèi)［5］.因而新西蘭成為世界上地震頻繁的國家之一（每年超過15000個ML≥2的地震）［7］，歷史上發(fā)生的大震包括1855年的Mw8.2 Wairarapa地震和1931年的 Mw7.8Hawke′s Bay地震.

2010年的Mw7.1Darfield地震發(fā)生在地震活動相對寂靜的地區(qū)，發(fā)震斷層震前并未為人所知.地震之后一些地質(zhì)學家和地球物理學家迅速進行了相關研究，調(diào)查了斷裂的地表破裂軌跡［8］，給出了地震的震源機制解［7，9］.同時，利用地震波形以及GPS和InSAR位移場反演了該地震的同震破裂分布［1，4，10］.本文中，我們利用一套時空基線更好的InSAR衛(wèi)星觀測數(shù)據(jù)，結(jié)合地質(zhì)調(diào)查，通過最大后驗概率反演方法同時確定斷層的幾何參數(shù)，光滑參數(shù)，不同InSAR數(shù)據(jù)集的權(quán)重，以及這次地震的同震滑動分布.最后，將形變反演結(jié)果與野外觀測結(jié)果進行系統(tǒng)比較，并研究該地震對周圍地區(qū)應力場的加載作用.

圖1 新西蘭地區(qū)構(gòu)造形變主要特征（箭頭表示板塊相對運動方向）Fig.1 Main features of plate tectonics in New Zealand（Arrows show relative plate motion directions）

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

本研究采用了日本空間探測局（JAXA）的ALOS衛(wèi)星上PALSAR傳感器獲得的雷達數(shù)據(jù)，在覆蓋地震發(fā)震時間范圍內(nèi)，選取了三對垂直基線較短且完整覆蓋地震形變區(qū)的數(shù)據(jù)源（表1），包括兩對升軌數(shù)據(jù)和一對降軌數(shù)據(jù).

表1 衛(wèi)星數(shù)據(jù)參數(shù)Table 1 Parameters of Raw Data

選擇的兩條升軌數(shù)據(jù)都滿足垂直基線較短，震后時間較短，主輔圖間隔時間較短三個要素，有助于保持震前震后雷達信號的干涉相干性，且同震形變場中混入的震后形變較小，從而獲得高信噪比和高精度的同震形變場.與前人所采用的雷達數(shù)據(jù)不同的是，本研究所用的336軌道數(shù)據(jù)垂直基線為302m，成像時間間隔為6個月，而前人采用的數(shù)據(jù)雷達數(shù)據(jù)成像時間間隔達到了2年6個月［10］.另外，在數(shù)據(jù)處理過程中，我們采用了二階雷達圖像配準算法，使得雷達干涉信號的相干性有所增強.特別是本研究所用的336軌道數(shù)據(jù)，只有在使用二階配準算法的情況下才能得到完整的形變場，而前人使用的一階配準模型無法得到有效的形變信號.因此，本研究中的數(shù)據(jù)選擇和處理方面比前人有明顯改善，而這些顯著改善對于更好約束斷層的幾何特征和滑動分布非常關鍵［11－12］.

本研究使用美國JPL實驗室開發(fā)的Repeat Orbit Interferometry Package［13］（ROI－PAC）軟 件包以及Gamma InSAR處理軟件包，處理獲得同震形變干涉條紋（圖2），同時用雷達幅值信息互相關的方法得到雷達斜距向和方位向的水平位移數(shù)據(jù)（圖3（c，d））.InSAR相位通過全局最優(yōu)算法解纏，獲得了覆蓋地震形變區(qū)的包含三分量位移信息的雷達視線向圖2（Line－of－Sight，LOS）位移［14］.

差分干涉雷達相位信息中包括若干分量：φ＝φdef＋Δφε＋φatm＋Δφorb＋φn.其中φdef是形變相位，為地表形變在LOS向引起的相位變化；φatm是大氣噪音，由大氣延遲作用引起的相位誤差；Δφorb是由軌道誤差引起的剩余相位；Δφε是高程剩余相位，由高程模型DEM中的誤差引起的相位，φn是系統(tǒng)噪音，為處理過程中引起的誤差相位.φdef形變相位是有用相位，其它的相位的貢獻都被視為形變測量中的誤差源.高程剩余相位Δφε與垂直基線和高程模型有關.本文采用精度為3″的SRTM高程模型［15］，并將其插值到1″用于消除干涉相位中的地形貢獻，其引入的誤差相對較小，對同震形變場的影響可以忽略.系統(tǒng)噪音φn是信號噪音、插值誤差和配準誤差等引入的相位誤差［13］，在相干性較好的數(shù)據(jù)中，系統(tǒng)噪音對干涉相位的貢獻非常小，一般可以忽略不計.因此本文處理的數(shù)據(jù)中，大氣誤差φatm和 軌道誤差φorb是主要的誤差來源.

觀測信號中的形變相位φdef，大氣延遲相位φatm，軌道誤差相位Δφorb都具有空間相關性，但是各自的空間相關尺度不一樣，可以利用其空間相關尺度的不同將這幾種信號區(qū)分開［16－17］.軌道誤差信號具有長波特征，其尺度大于區(qū)域性形變特征尺度.同震形變主要集中在近場，選擇一個足夠遠的區(qū)域使得同震形變近似為零，利用最小二乘方法將此遠場區(qū)域誤差用一個線性或者二階的剩余軌道面進行擬合，得到剩余軌道相位，并從干涉形變場中消除，從而消除軌道誤差的影響.另外，大氣噪聲中有一部分與地形起伏線性相關［13］，因此可以將上述改正的遠場信號再用最小二乘方法擬合其與地形的線性相關性，將得到的大氣分量從總的形變場中進一步消除掉.

3 同震破裂模型反演

對于線性、非線性混合反演的問題，目前有兩種比較常見的概率反演算法，分別是完全貝葉斯概率反演方法（Fully Bayesian Inversion，以下簡稱FBI）［18］和線性－非線性混合蒙特卡洛反演方法（Mixed linear－nonlinear Monte－Carlo method，以下簡稱MCI）［19］.前者將所有參數(shù)，包括斷層幾何參數(shù)和滑動參數(shù)，數(shù)據(jù)相對權(quán)重和光滑因子等都作為非線性參數(shù)反演；而后者將滑動參數(shù)作為線性參數(shù)（基于最小二乘方法）、其它參數(shù)作為非線性參數(shù)（基于蒙特卡洛方法）反演.FBI方法對所有的未知量都采用蒙特卡洛隨機采樣進行試探，以期獲得最終的收斂解，其求解大量未知量的地球物理問題時，存在耗時長、收斂慢的問題.也正是因為這一不足，F(xiàn)ukuda和Johnson又將FBI方法改良為了 MCI.在 MCI中，未知量被分成線性未知量和非線性未知量兩類.MCI算法大大降低了未知力量的個數(shù)，需要通過蒙特卡洛方法反演的未知參量的個數(shù)，因而能夠有效減少迭代次數(shù).兩種方法使用相同的參數(shù)采樣方法，即馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法.它們有一個共同的問題，即參數(shù)采樣步長選擇不當可能使得反演不收斂.如果采樣步長設定過大，反演參數(shù)采樣可能趨向于小概率區(qū)域，得到一個滑動分布過于粗糙但數(shù)據(jù)擬合很好的解.如果采樣步長過小，每一步的采樣缺乏效率，需要大量的嘗試才能得到收斂解.實際操作中，往往需要多次嘗試才能最終確定合理步長.對于存在較大誤差的InSAR觀測數(shù)據(jù)，很難在較短的時間內(nèi)獲得合理的反演結(jié)果.鑒于此，我們最終采用了最大后驗概率反演方法（three－step Maximum－A－Posteriori Probability inversion）.

3.1 最大后驗概率反演簡介

形變場的最大后驗概率反演［21］過程分為三步.

第一步，采用自適應模擬退火法（Adaptive Simulated Annealing，以下簡稱 ASA）［22］搜索形變模型的最大后驗概率，將滿足最大后驗概率的斷層滑動模型作為形變場反演結(jié)果，反演中模型的后驗概率公式沿用了Fukuda和Johnson［18］在MCI方法中的定義.每一步采樣中，斷層傾角、數(shù)據(jù)權(quán)重和平滑因子作為非線性量隨機采樣，而斷層錯動量則通過最小二乘法求得.不斷重復這一過程，直到得到穩(wěn)定的收斂結(jié)果.反演第一步是全局最優(yōu)化反演，其效率遠遠高于FBI和MCI中的MCMC采樣方法.ASA反演中自動確定參數(shù)步長和各個參數(shù)在反演中的敏感性，因此該步反演能夠很快收斂到最大后驗概率模型附近，同時局部后驗概率極大值和孤立的參數(shù)采樣區(qū)域都被逐一測試和跳過.

第二步，將第一步得到的模型作為初始值，采用類似于FBI的方法對各個未知量進行了帶約束的非線性反演，將斷層幾何參數(shù)和滑動參數(shù)均作為非線性參數(shù)反演.第一步中為了反演快速收斂到合理的區(qū)域，斷層滑動參數(shù)在反演中用最小二乘方法一次反演得到，但最小二乘法無法對滑動參數(shù)施加約束，可能會在一些數(shù)據(jù)約束不夠強的區(qū)域（如斷層較深部位），或者觀測數(shù)據(jù)局部噪聲過大時，出現(xiàn)不合理的解，如局部的滑動方向與周圍區(qū)域的滑動方向相反，需要進一步的修正；因此第二步中將斷層上每一個矩形元的滑動參數(shù)均作為非線性參數(shù)在給定的范圍內(nèi)逐一單獨采樣反演，便于對滑動參數(shù)的值進行進一步的修正，反演進一步趨于系統(tǒng)解的最大后驗概率.

第三步，將第二步中得到的斷層幾何模型固定下來，只反演錯動量及數(shù)據(jù)權(quán)重、平滑因子等非線性參數(shù).該步反演可以避免斷層幾何參數(shù)和滑動參數(shù)之間的相互影響（trade－off），同時進一步使得反演趨于最大后驗概率模型.

3.2 斷層分段依據(jù)和斷層反演模型設置

對于主要發(fā)震斷層Greendale斷層，前人做了很多卓有成效的調(diào)查工作，Elliot等［10，23］已經(jīng)準確地描述出了本次地震在地表的破裂軌跡（圖3a）.本研究中的兩條升軌干涉圖信號具有很高的信噪比，形變干涉條紋連續(xù)清晰.337軌道的干涉條紋圖中（圖2b）可以看到，形變條紋中存在一條近東西走向的分界線，分界線兩側(cè)具有不同的條紋色彩變化模式，顯示它們具有不同的LOS移動方向.該分界線與野外觀測發(fā)現(xiàn)的Greendale斷層一致.Greendale斷層以南的條紋呈現(xiàn)以斷層為界的橢圓形，且條紋非常連續(xù)，表明這個區(qū)域的形變主要由Greendale斷層的滑動產(chǎn)生.斷層北側(cè)，可以明顯看到形變干涉條紋的不規(guī)則扭曲，表明這個區(qū)域內(nèi)可能存在多條斷層，并產(chǎn)生了破裂.

從337軌道的信號相干性（圖3b）圖上，可以看到6個低相干性條帶，其中1—4是貫通的，與地質(zhì)調(diào)查的Greendale斷層基本一致（圖3a），5號低相干帶位于Greendale斷層的西南端，可以解釋在337、631（圖2b，圖2c）兩條軌道干涉圖中在（171.9°E，43.57°S）附近出現(xiàn)的單獨閉合條紋帶.6號低相干帶接近儀器震中，從336、337兩幅干涉圖（圖2a，圖2b）中可以明顯看到幾條與周圍顏色反向的條紋，說明此處存在一條較小斷層.以上幾條斷層在337軌道斜距向位移場中也能夠看到（圖3c）.從337軌道方位方向上（圖3d）還能看到一條位移突變分界線，相同位置上336、337的干涉圖條紋向南彎曲，顯示此處存在一條微弱滑動的斷層，在對應位置添加了一條斷層單元7.根據(jù)余震分布［24］，將斷層兩端稍加延伸，使得整個斷層模型東西總長約45km.破裂主體的Greendale斷層，根據(jù)走向變化由西向東分成四段，分別編號為1至4，三個隱伏斷層由西向東分別編號為5至7.經(jīng)過多次測試之后，發(fā)現(xiàn)斷層上的錯動量基本集中在10km深度以內(nèi)，同時，考慮到InSAR數(shù)據(jù)對斷層錯動的分辨率［11］，最終將斷層深度固定為12km，相關參數(shù)如表2所示.斷層面劃分成約1km×1km的矩形位錯元，傾角作為未知量參與非線性反演.

表2 斷層單元相關參數(shù)Table 2 Parameters of fault segments

3.3 位移場離散采樣

干涉雷達獲得的形變場是連續(xù)的，一幅干涉圖像中包含了大約106個數(shù)據(jù)點，像元之間具有較強的相關性，若利用全部的觀測數(shù)據(jù)點進行反演，會造成計算困難.本文采用四叉樹離散采樣法［25］對觀測數(shù)據(jù)離散采樣.該方法利用形變量變化梯度值為圖像采樣依據(jù)，采樣后的觀測數(shù)據(jù)，既能有效提高計算效率，同時不會在近場劇烈形變區(qū)域丟失信息，可以代表整個形變場的特征，降低了噪音對形變場的影響［26］.經(jīng)過測試，選用形變梯度0.06為閾值，即當形變區(qū)域梯度大于0.06時將區(qū)域分割為4個小塊，直至形變梯度小于0.06為止，同時忽略了干涉圖中相關系數(shù)小于0.25區(qū)域內(nèi)的數(shù)據(jù)點.最終從336、337和631三個軌道的干涉圖里分別提取出了1152，2513，733個數(shù)據(jù)點加入反演.

4 反演結(jié)果及討論

4.1 斷層幾何形態(tài)及同震破裂分布

通過三步反演，最終得到了令最大后驗概率趨于最大值的解.每一步反演的相關結(jié)果如表3所示.

表3 三步反演的結(jié)果Table 3 Results of the three－step method

本次事件總共釋放的地震矩達到了5.0×1019N·m，對應的矩震級約為7.1.模型的滑動分布如圖4a、表4所示.

表4 各斷層整體錯動特征Table 4 Overall Fault Slip Parameters

圖4 （a）各斷層面上位錯分布圖；（b）模型中Greendale斷層地表錯動與其他結(jié)果對比圖Fig.4 （a）Slip distribution on each fault segment；（b）Comparison of surface rupture calculated in our model with other data

對于Darfield地震的同震破裂模型，已經(jīng)有兩篇研究成果發(fā)表［4，10］.其中 Beavan等采用了 GPS數(shù)據(jù)（尚未公開）以及 ALOS 337軌道獲得的InSAR數(shù)據(jù)進行聯(lián)合反演，Elliot等采用了ALOS三個軌道的InSAR數(shù)據(jù)以及P波和SH波形數(shù)據(jù)分別反演.本文同樣采用了三個軌道的InSAR數(shù)據(jù)，但部分數(shù)據(jù)不同于Elliot等的模型中所用數(shù)據(jù)［10］，本研究中的336軌道震前、震后數(shù)據(jù)離發(fā)震時間更近，因此獲得的數(shù)據(jù)具有更高的性噪比，對于研究斷層東段的破裂有很大幫助.

本研究的反演結(jié)果顯示，Greendale斷層總長度為45km，1號子斷層位于它的最西端，傾角為東北傾70°，最東邊有一個最大值約3.5m的破裂集中區(qū)，往西逐漸減少.最西端（第1段）的3km范圍滑動量接近于0，與地表觀測吻合較好.Greendale斷層中央的2號和3號子斷層是主要破裂段，傾角分別為南傾87°和85.5°.2號子斷層面上的破裂較集中，滑動發(fā)生在深度5km以上區(qū)域，最大滑動約為7.0m.3號子斷層面上的滑動最為顯著，出現(xiàn)了整個斷層的最大滑動，約8.5m，該段部分區(qū)域的破裂在深度9.0km處仍達到了3.1m.Greendale斷層東段，即4號子斷層，是一個南傾69.5°的斷層，最大破裂約為4.2m，發(fā)生在與3號斷層相連的地方.在它的末端，出現(xiàn)了一部分的逆沖分量，最大約為1.0m.

斷層段6長度為6km，傾角為59°（傾向東南）.反演顯示該段以逆沖為主，最大逆沖分量約為5.3m，整個斷層面上都有比較大的滑動，并且一直延伸到10km深度，該斷層段是此次地震中初始破裂所在位置，因此區(qū)域地震波形反演得到了以逆沖為主的初始破裂震源機制解［1，9］.該斷層的西南端與Greendale斷層相連.

5號斷層段也是一個逆沖斷層，傾角為50°（傾向西北）.它的破裂主要出現(xiàn)在中部5km以上的較小區(qū)域內(nèi)，最大的滑動為3.8m，相當于一個Mw＝6.2級地震，最大滑移量比Beavan等［4］的模型值略大，與Elliott等［10］的模型值比較一致，能夠很好地解釋干涉條紋中5號斷層段附近的獨立封閉形變條紋.

7號斷層段在InSAR干涉圖上沒有明顯的顯示，因此并沒有加入到Beaven等［4］的模型中，計算顯示的破裂也很小，最大破裂不超過2.0m，以走滑破裂為主，對整個形變場貢獻較小.

4.2 LOS位移場殘差

圖5是正演得到 的三條軌道上預測的InSAR形變場以及殘差.336和337軌道，信號較好，獲得了大范圍清晰形變場，模型正演產(chǎn)生的條紋整體形態(tài)與觀測比較吻合.與前人的結(jié)果［4，10］比較，本文的模型能夠更好地解釋形變觀測.特別是，前人的破裂模型無法有效消除6號斷層段東側(cè)的剩余殘差，而從本文的模型殘差圖中可以看到，無論是336軌道還是337軌道，該區(qū)域近場殘差最大僅為20.0cm左右.

631軌道基線較長，且震前震后雷達數(shù)據(jù)獲取時間間隔較長，因而相干性不高，只得到了不連續(xù)的形變條紋，但能看到大體形變趨勢.從模型計算預測的形變條紋看，其基本趨勢與觀測數(shù)據(jù)是一致的，特別是5、6號逆沖斷層段引起的單獨條紋帶比較明顯.Greendale斷層東西兩端存在一部分殘差，說明形變場還存在一部分形變無法用現(xiàn)有模型有效解釋，該殘差可能是由于形變場中包含的非彈性形變信息，如沙土液化等因素所致.沙土液化現(xiàn)象在本次地震中大量存在，也造成了一定的地面形變［4］.

4.3 與前人結(jié)果的比較

斷層幾何參數(shù)方面，本研究所用的4條斷層位置與Elliott等［10］所用的斷層模型基本一致，而與Beavan等［4］所用模型僅相差一條整體破裂較小的7號斷層段.從斷層幾何形態(tài)來看，本文的模型與前人研究給出的Greendale斷層傾角變化一致，6號斷層的傾角在本文的模型和Elliott等的反演結(jié)果［10］中均為約50°，而Beavan等［4］計算給出的傾角為75°.考慮到該段以逆沖斷層活動為主，較緩的傾角可能更加合理.

從錯動分布來看，我們計算得到靠近震源的斷層，即6號斷層，整個斷層面上都有較大的錯動，在約3km深度的地方存在3.0m左右的滑動.從殘差上也能看到，這種滑動分布能比較合理地解釋6號斷層東側(cè)的位移.Elliott等［10］的結(jié)果顯示在這條斷層上只在4km深度以下的部分區(qū)域才存在較明顯的錯動，可能低估了這條斷層段的滑動量.對于Greendale斷層，Beavan等［4］，Elliott等［10］以及我們的結(jié)果都顯示Greendale斷層的最大破裂出現(xiàn)在中部偏東的區(qū)域.Beavan等［4］的模型最大滑動量（約5.0m）遠小于Elliott等模型的約8.2m以及我們的結(jié)果約8.5m.這可能是Beavan等［4］所采用的數(shù)據(jù)缺乏足夠的近場信息造成的.在4號子斷層的中部3～4km深度，本文模型顯示存在超過4m的破裂，超過了其它兩篇文獻中模型的滑動量，該差異可能來源于本文采用的336軌道數(shù)據(jù)的形變信號優(yōu)于Elliott等［10］采用的形變數(shù)據(jù)，因而本文模型能提供4號子斷層附近更多的破裂信息.

我們選取Greendale斷層0～1km深度斷層面上的位移近似表征它的表面滑動，與地質(zhì)調(diào)查結(jié)果［8］、航拍數(shù)據(jù)等遙感數(shù)據(jù)［10］以及 Elliott等Greendale斷層模型中0～1km深度上的位移［10］進行對比（圖4b）.結(jié)果表明，觀測數(shù)據(jù)之間匹配較好，地表的最大破裂出現(xiàn)在2號和3號子斷層交界的區(qū)域，最大值在5.5～6m之間，之后向兩邊減小.在我們的模型中，地表的最大錯動出現(xiàn)在3號斷層的最西端，滑動量為6.2m，與觀測結(jié)果基本吻合，而在向兩邊減小的趨勢也與觀測結(jié)果大致吻合.Elliott等［10］的模型中可能由于過強的滑動約束使斷層滑動過于光滑，該模型中的地表錯動在1號、2號、3號子斷層都小于地表野外觀測數(shù)據(jù)，而在4號子斷層上又高于地表觀測，峰值約4.0m出現(xiàn)在3、4號子斷層交界的地方，比觀測結(jié)果偏東7km左右.因此，本文的模型結(jié)果比Elliott等［10］的反演結(jié)果更合理，能更好解釋地表位移數(shù)據(jù).在2號子斷層中部偏東部分，本文模型中給出的地表位移高于地表觀測結(jié)果2.0m左右，這可能與2號斷層和6號斷層在此處相鄰，它們的滑移存在較高相關造成模型誤差增大有關.

4.4 Darfield地震對周圍區(qū)域應力場的影響

Darfield地震是該地區(qū)應力場高度集聚后的調(diào)整，同時會對其附近應力場中余震的發(fā)生產(chǎn)生一定影響［27］.Darfield地震的余震既有逆沖型地震，也有走滑型地震，反映了該地區(qū)比較復雜的區(qū)域應力狀態(tài).我們選取Greendale斷層的震源機制（走滑為主）作為接收斷層上地震的震源機制，計算Darfield地震對于周圍區(qū)域應力場的影響［28－29］.接收斷層的參數(shù)選如表5.

圖5 位移條紋和模型殘差（由纏繞條紋表示，一個周期：－11.81－11.81cm）Fig.5 Deformation fringes and residues of each path（Showed by wrapped fringes，value in one cycle is from －11.81cm to 11.81cm）

表5 選定的接收斷層參數(shù)Table 5 Parameter of receiving fault

圖6a中顯示了這一系列斷裂對6.0km深度范圍應力場的加卸載作用.與余震的分布對照可以看到，約有65%的余震發(fā)生距在發(fā)震斷層3km距離范圍內(nèi)，這是由于主震發(fā)生之后斷層附近積累的應力沒有完全釋放，而發(fā)震后的斷層面容易發(fā)生破裂，因此在斷層面上發(fā)生了后續(xù)破裂.如Greendale斷層西側(cè)5號斷層段附近，雖然主震減小了該區(qū)域的應力，但仍能明顯地看到斷層面上發(fā)生了密集余震.遠離斷層的余震帶主要是以43.55°S，172.05°E為中心的一個矩形，這恰好與庫侖應力的一個增加區(qū)域（紅區(qū)）相吻合，因此我們推測，這個區(qū)域內(nèi)可能存在一些隱伏斷層，震后釋放了Darfield地震施加的部分應力.Greendale斷層東南端的藍區(qū)也出現(xiàn)了較多的余震，這個余震區(qū)域呈帶狀分布，向東北延伸到2011年Christchurch地震的發(fā)震斷層，我們認為這部分小震預示著區(qū)域應力釋放的方向從Greendale斷層轉(zhuǎn)變到一條西南—東北走向的斷層上，Christchurch地震的發(fā)生可能與此有關.

2011年Christchurch地震區(qū)域正好處于一個庫侖應力增大的區(qū)域.選定截面為Christchurch地震的斷層面，可以看到在斷層面上庫侖應力增加值在0.5～1.2bar之間，在震源處［7］增大量約為1.0 bar（圖6b）.可見僅僅從庫侖應力的變化來考慮，Darfield地震的發(fā)生對于Christchurch地震起到了一定的觸發(fā)作用.

4.5 區(qū)域應力場分配與震源機制

圖6 （a）斷層破裂造成周圍區(qū)域庫侖應力場變化（圓點為震后一個月內(nèi)余震分布）；（b）Darfield地震引起的Christchurch地震斷層面上的庫倫應力改變Fig.6 （a）The change of the Coulomb stress around the faults（The dots are the after－shocks in one month）；（b）The change of the Coulomb stress on the fault of the Christchurch earthquake caused by the Darfield earthquake

澳洲板塊和太平洋板塊的邊界帶貫穿新西蘭南島.根據(jù)GPS數(shù)據(jù)計算發(fā)現(xiàn)（圖1），在此地區(qū)太平洋板塊相對澳洲板塊存在約42mm／a方向為南偏西76°的運動［3，29］.邊界帶內(nèi)的阿爾卑斯斷裂帶的走向為北偏東60°［7］，存在約40mm／a的右旋走滑［3］.因此邊界帶內(nèi)阿爾卑斯斷裂帶以外仍有約10mm／a垂直于邊界帶的相對縮短，即北偏西約60°的壓應變在新西蘭南島內(nèi)積累.Darfield地震破裂主要發(fā)生在走向33°的6號逆沖斷層段，以及近東西走向的Greendale右旋走滑斷層上，它是板塊邊界帶區(qū)域應力的一種表現(xiàn).區(qū)域震源機制解顯示此次地震的初始破裂是發(fā)生在6段斷層段上的逆沖破裂［1，9］，其最大主壓應力的方向垂直于斷層的走向，因此6號斷層段的震源機制對應于北偏西57°的壓應力.其后發(fā)生的破裂主要為近東西走向的Greendale斷層上的右旋滑移，其震源機制對應于北偏西約45°的壓應力.因此這二個主破裂段釋放的應力場就是板塊邊界帶應變分配在該地區(qū)殘留的構(gòu)造應力場，Darfield地震的復雜性體現(xiàn)了區(qū)域構(gòu)造應力場的特點.

5 結(jié) 論

通過完整覆蓋發(fā)震區(qū)域且信噪比較高的SAR數(shù)據(jù)對，獲得了高空間分辨率的同震形變場.根據(jù)形變信息，可以對Darfield地震的斷層形態(tài)和同震破裂分布進行研究.本文中斷層錯動模型在地表的滑動，相比前人結(jié)果，與地質(zhì)調(diào)查和遙感觀測數(shù)據(jù)更為吻合.根據(jù)錯動模型計算得到的庫倫應力變化場跟遠離發(fā)震斷層的余震有較好的對應.2011年2月21號的Christchurch地震震源處于Darfield地震同震應力增加區(qū)域內(nèi)，庫侖應力增加約1.0bar，證明Darfield地震對于余震Christchurch地震存在應力加載效應.

（References）

［1］Gledhill K，Ristau J，Reyners M，et al.The darfield（Canterbury，New Zealand）Mw7.1earthquake of september 2010：apreliminary seismological report.Seismological Research Letters，2011，82（3）：378－386，doi：10.1785／gssrl.82.6.378.

［2］Reyners M. Lessons from the destructive Mw6.3 Christchurch，New Zealand，earthquake.Seismol.Res.Lett.，2011，82（3）：371－372.

［3］Norris R J，Cooper A F.Late quaternary slip rates and slip partitioning on the alpine fault，New Zealand.J.Structural Geol.，2001，23（2－3）：507－520.

［4］Beavan J，Samsonov S，Motagh M，et al.The MW7.1 Darfield（CANTERBURY）earthquake：Geodetic observations and preliminary source model.Bull.N.Z.Soc.Earthq.Eng，2010，43：228－235.

［5］Wallace L M，Beavan J.A large slow slip event on the central Hikurangi subduction interface beneath the Manawatu region，North Island，New Zealand.Geophysical Research Letters，2006，33（11），doi：10.1029／2006GL026009.

［6］Wallace L M，Beavan J，McCaffrey R，et al.Balancing the plate motion budget in the South Island，New Zealand using GPS，geological and seismological data.Geophysical Journal International，2007，168（1）：332－352.

［7］USGS http：／／earthquake. usgs. gov／earthquakes／eqinthenews／2010／us2010atbj／

［8］Quigley M，Dissen R V，Villamor P，et al.Surface rupture of the greendale fault during the darfield （canterbury）earthquake，New Zealand：initial findings.Bulletin of the New Zealand Society for Earthquake Engineering，2010，43（4）：236－242.

［9］Holden C，Beavan J，F(xiàn)ry B，et al.Preliminary source model of the Mw7.1Darfield earthquake from geological，geodetic and seismic data ／／ Proceedings of the Ninth Pacific Conference on Earthquake Engineering Building an Earthquake－Resilient Society，Auckland，New Zealand，2011.

［10］Elliott J R，Nissen E K，England P C，et al.Slip in the 2010—2011Canterbury earthquakes，New Zealand.Journal of Geophysical Research，2012，117，doi：10.1029／2011JB008868.

［11］Delouis B，Giardini D，Lundgren P，et al.Joint inversion of InSAR，GPS，Teleseismic，and Strong－Motion data for the spatial and temporal distribution of earthquake slip：Application to the 1999izmit mainshock.Bulletin of the Seismological Society of America，2002，92（1）：278－299.

［12］Delouis B，Nocquet J M，Vallée M，et al.Slip distribution of the February 27，2010 Mw＝8.8maule earthquake，central chile，from static and high－rate GPS，InSAR，and broadband teleseismic data.Geophysical Research Letters，2010，37（17），L17305，doi：10.1029／2010GL043899.

［13］Rosen P A，Hensley S，Peltzer G，et al.Updated repeat orbit interferometry package released.Eos，Transactions American Geophysical Union，2004，85（5）：47.

［14］Fialko Y，Simons M，Agnew D.The complete（3－D）surface displacement field in the epicentral area of the 1999 Mw7.1 Hector Mine Earthquake，California，from space geodetic observations.Geophysical Research Letters，2001，28（16）：3063－3066.

［15］Farr T G，Rosen P A，Caro E，et al.The shuttle radar topography mission.Rev.Geophys.，2007，45（2），doi：10.1029／2005RG000183.

［16］Hanssen R F.Radar Interferometry－Data Interpretation and Error Analysis（V.2）.Dordrecht：Kluwer Academic Pub.，2001.

［17］孫建寶，徐錫偉，石耀霖等.東昆侖斷裂瑪尼段震間形變場的InSAR觀測及斷層滑動率初步估計.自然科學進展，2007，17（10）：1361－1370.Sun J B，Xu X W，Shi Y L，et al.InSAR observation of the inter－seismic deformation of Mani segment in the East Kunlun fault and estimation of the slip rate.Progress in Natural Science （in Chinese），2007，17（10）：1361－1370.

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

［19］Fukuda J，Johnson K M.Mixed linear－non－linear inversion of crustal deformation data：Bayesian inference of model，weighting and regularization parameters. Geophysical Journal International，2010，181（3）：1441－1458.

［20］Okada Y.Internal deformation due to shear and tensile faults in a half－space，Bull.Seism.Soc.Am.，1992，82（2）：1018－1040.

［21］Sun J B，Shen Z K，Bürgmann R，et al.A three－step Maximum－A－Posteriori probability method for INSAR data inversion of coseismic rupture with application to the April 14，2010 Mw6.9Yushu，China earthquake.Journal of Geophysical Research：Solid Earth，doi：10.1002／jgrb.50244.

［22］Ingber L.Simulated annealing：Practice versus theory.Math.Comput.Model.，1993，18（11）：29－57.

［23］Quigley M，van Dissen R，Litchfield N，et al.Surface rupture during the 2010 Mw7.1Darfield （Canterbury）earthquake：Implications for fault rupture dynamics and seismic－h(huán)azard analysis.Geology，2012，40（1）：55－58.

［24］GeoNet http：／／www.geonet.org.nz／

［25］Jónsson S，Zebker H，Segall P，et al.Fault slip distribution of the 1999 Mw7.1Hector Mine，California，earthquake，estimated from satellite radar and GPS measurements.Bulletin of the Seismological Society of America，2002，92（4）：1377－1389.

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

Sun J B，Xu X W，Shen Z K，et al.Parameter inversion of the 1997Mani earthquake from INSAR co－seismic deformation field based on linear elastic dislocation model—Ⅰ.uniform slip inversion.Chinese J.Geophys.（in Chinese），2007，50（4）：1097－1110.

［27］King G C P，Stein R S，Lin J.Static stress changes and the triggering of earthquakes.Bulletin of the Seismological Society of America，1994，84（3）：935－953.

［28］Lin J，Stein R S.Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike－slip faults.Journal of Geophysical Research，2004，109（B2），doi：10.1029／2003JB002607.

［29］Sella G F，Dixon T H，Mao A L.REVEL：A model for recent plate velocities from space geodesy.Journal of Geophysical Research，2002，107（B4），doi：10.1029／2000JB000033.