汶川地震同震形變場(chǎng)的GPS和InSAR鄰軌平滑校正與斷層滑移精化反演

楊瑩輝，陳強(qiáng)＊，劉國(guó)祥，程海琴，2，劉麗瑤，胡植慶

1西南交通大學(xué)遙感信息工程系，成都 610031

2華東交通大學(xué)土木建筑學(xué)院，南昌 330013

3臺(tái)灣大學(xué)地質(zhì)科學(xué)系，臺(tái)灣臺(tái)北

1 引言

2008年5月12日，以汶川縣映秀鎮(zhèn)為震中發(fā)生了MS8.0級(jí)特大地震，地震發(fā)生機(jī)制為龍門(mén)山中央斷裂的右旋走滑逆沖，其動(dòng)力誘因?yàn)橛《劝鍓K持續(xù)向北運(yùn)動(dòng)擠壓歐亞大陸板塊，使得青藏高原東緣地塊向東運(yùn)動(dòng)，而與之相鄰的四川盆地地塊較為堅(jiān)固（徐錫偉等，2005；單新建等，2009），導(dǎo)致龍門(mén)山斷層應(yīng)力持續(xù)積累并趨于臨界狀態(tài)，最終積蓄能量從映秀—北川主斷裂帶釋放導(dǎo)致汶川地震的發(fā)生.

汶川地震引發(fā)龍門(mén)山地區(qū)強(qiáng)烈的地表形變，同震形變場(chǎng)的準(zhǔn)確獲取對(duì)于分析龍門(mén)山地區(qū)地表破裂的空間分布特征、震源參數(shù)計(jì)算以及地球物理反演具有重要的科學(xué)意義.汶川地震后，中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)項(xiàng)目組及時(shí)對(duì)龍門(mén)山地區(qū)的GPS觀測(cè)站進(jìn)行了復(fù)測(cè)，并公布了觀測(cè)結(jié)果（張培震，2008）.GPS形變觀測(cè)顯示，以映秀—北川斷裂帶為中心，兩側(cè)斷層上下盤(pán)存在顯著的相向運(yùn)動(dòng)和水平縮短現(xiàn)象，同時(shí)斷層下盤(pán)以下降為主，上盤(pán)在地表破裂線周邊呈現(xiàn)上升趨勢(shì)，但在遠(yuǎn)離破裂線后很快又轉(zhuǎn)為下降運(yùn)動(dòng)，GPS形變的這些特征揭示了斷層運(yùn)動(dòng)以逆沖為主的特性（張培震，2008）.但由于GPS稀疏分布的特性，難以獲得發(fā)震區(qū)完整的高分辨率地表形變場(chǎng)，不利于斷層面破裂的精細(xì)反演計(jì)算.汶川地震前后，ALOS衛(wèi)星搭載的PALSAR傳感器獲取了該區(qū)域大范圍SAR影像數(shù)據(jù)，InSAR技術(shù)可為龍門(mén)山地區(qū)提供大范圍、高分辨率的地表形變觀測(cè)資料（單新建等，2009）.

汶川地震InSAR同震形變場(chǎng)的觀測(cè)研究表明，由于雷達(dá)衛(wèi)星數(shù)據(jù)受軌道誤差、大氣信號(hào)延遲等負(fù)面因素的影響，從InSAR干涉圖中提取精細(xì)化的地表形變場(chǎng)干擾較大，尤其是汶川地震波及范圍廣，鄰軌PALSAR干涉圖之間存在嚴(yán)重的形變場(chǎng)不連續(xù)問(wèn)題（Zhang etal.，2008；Hashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010），積極探索鄰軌InSAR同震位移的修正方法，有助于準(zhǔn)確獲取龍門(mén)山地區(qū)精細(xì)的同震形變場(chǎng)與斷層滑移參數(shù)的精化反演.

為解決InSAR軌道誤差殘留相位與鄰軌形變條紋的不連續(xù)問(wèn)題，本文提出聯(lián)合GPS形變觀測(cè)與鄰軌相位修正的數(shù)據(jù)處理策略，利用稀疏分布的GPS位移觀測(cè)去除PALSAR軌道誤差引入的殘留平地相位，采用鄰軌約束平滑消除干涉圖之間的不連續(xù)問(wèn)題，從而獲得龍門(mén)山地區(qū)大范圍精細(xì)的同震形變場(chǎng)，在此基礎(chǔ)上，采用彈性半空間位錯(cuò)模型反演地殼發(fā)震斷層的滑移模型，并由此計(jì)算汶川強(qiáng)震的地震矩?cái)?shù)值.

2 ALOS／PALSAR影像干涉處理

使用ALOS衛(wèi)星PALSAR傳感器獲取的地震前后雷達(dá)影像開(kāi)展干涉處理，以其中連續(xù)6個(gè)升軌條帶共46幅影像作為同震形變觀測(cè)數(shù)據(jù)，各軌道影像的基本信息如表1所示.

雷達(dá)影像各條帶的場(chǎng)景覆蓋范圍如圖1紅色矩形方框所示，圖中的綠線為野外調(diào)查獲得的汶川地震地表破裂線，主要包括中央主斷裂（映秀—北川）、小魚(yú)洞斷裂和前山斷裂（灌縣—彭縣），小魚(yú)洞斷裂帶位于中央斷裂帶和前山斷裂帶之間，前山斷裂帶長(zhǎng)60余公里，在地下約10km處并入中央主斷裂帶（李勇等，2009）.圖1中的藍(lán)色與棕色箭頭表示龍門(mén)山地區(qū)的GPS同震形變觀測(cè)向量（張培震，2008），其中棕色箭頭代表單位長(zhǎng)度形變?yōu)?0cm，藍(lán)色箭頭單位長(zhǎng)度形變?yōu)?0cm.圖中的紅色小圓為MS＞4.0級(jí)的余震震中位置（Zhang etal.，2008；Hashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010）.

表1 汶川地震PALSAR影像的干涉對(duì)信息Table 1 PALSAR interferometric pairs used in the study of Wenchuan earthquake

對(duì)上述PALSAR影像進(jìn)行干涉處理，使用ASTER DEM作為外部DEM以去除地形相位分量，獲得雷達(dá)差分干涉相位圖（如圖2）.從圖2可以看出，相鄰軌道間的干涉條紋存在顯著的不連續(xù)性以及軌道殘留相位的干擾，在475—476軌道鄰接邊界（左上角的黑色方框）存在明顯的干涉條紋不連續(xù)，而圖幅中間黑色方框所標(biāo)識(shí)的473條帶則存在顯著的軌道殘留相位，導(dǎo)致遠(yuǎn)場(chǎng)存在明顯的干涉條紋.此外，在右下角的471—472條帶連接處，則存在嚴(yán)重的大氣相位，其形變條紋與473—474存在顯著的不一致性，造成形變場(chǎng)的不連續(xù).為克服上述干涉條紋形變場(chǎng)的不連續(xù)問(wèn)題，本文提出采用GPS形變位移對(duì)PALSAR軌道誤差殘留平地相位進(jìn)行建模修正，利用鄰軌平滑約束條件校正鄰軌干涉條紋.

3 軌道殘留相位建模與校正

在InSAR數(shù)據(jù)處理過(guò)程中，平地相位（也稱(chēng)為參考相位）的計(jì)算是基于SAR軌道狀態(tài)矢量（包括位置矢量與速度矢量），而ALOS衛(wèi)星軌道數(shù)據(jù)的不確定性導(dǎo)致干涉圖中存在殘留平地相位，在圖2中473條帶的干涉相位表現(xiàn)尤為顯著，在473條帶的黑色矩形框內(nèi)，地震遠(yuǎn)場(chǎng)呈現(xiàn)出顯著的周期性條紋，從形變場(chǎng)空間分布特征表明，該信號(hào)為殘留的平地相位.

圖1 汶川地震PALSAR影像場(chǎng)景覆蓋范圍Fig.1 Ground coverage of PALSAR images for the Wenchuan earthquake

圖2 InSAR同震形變干涉相位Fig.2 InSAR derived interferogram of coseismic deformation field

圖3 軌道殘留相位校正前后的干涉相位圖（a）校正前的干涉相位；（b）校正后的干涉相位.Fig.3 Local interferogram before（a）and after（b）correction from residual orbit phase modeling

為消除軌道數(shù)據(jù)誤差引入的殘留參考相位，本文提出以GPS精密形變測(cè)量值作為約束，對(duì)同一區(qū)域的InSAR與GPS形變觀測(cè)進(jìn)行差分建模，從而消除或減弱由于衛(wèi)星軌道誤差引入的殘留相位值.

首先將GPS觀測(cè)得到的三維形變轉(zhuǎn)換至雷達(dá)視線方向上的LOS形變（陳強(qiáng)等，2012a），定量的數(shù)學(xué)轉(zhuǎn)換模型如下：

式中，gLos代表GPS三維形變投影到衛(wèi)星LOS視線方向的形變，dUp、dN、dE分別為GPS測(cè)量得到的豎直向、南北向和東西向形變量，θ代表衛(wèi)星側(cè)視角，α是衛(wèi)星飛行方向與北方向的夾角，θ和α可依據(jù)衛(wèi)星參數(shù)予以計(jì)算.

從圖2中未校正的473條帶可以看出，軌道殘留相位在空間尺度上呈現(xiàn)出顯著的線性分布（張磊等，2007），因此構(gòu)建如下的軌道殘留誤差修正函數(shù)模型：

式中，gLos為GPS形變轉(zhuǎn)換到LOS方向的形變，可通過(guò)式（1）計(jì)算得到，dLos為InSAR觀測(cè)得到的形變，l、p分別表示GPS測(cè)站在雷達(dá)干涉圖中的方位向（Azimuth）和距離向（Range）坐標(biāo)，a0、a1、a2為待求模型系數(shù)，采用3個(gè)或3個(gè)以上的GPS形變觀測(cè)量，即可求解出上述系數(shù)，從而建立軌道殘留相位修正模型.

從圖2中的三維GPS形變測(cè)站分布來(lái)看，PALSAR條帶的472、473、474、475四個(gè)條帶中，GPS點(diǎn)數(shù)量均大于3個(gè)，可用于軌道殘留相位的建模，而472條帶中的H010號(hào)GPS點(diǎn)位于余震分布密集的區(qū)域，其InSAR相干性較差，形變觀測(cè)數(shù)據(jù)質(zhì)量不高，不用于軌道誤差的校正.

采用上述模型和GPS觀測(cè)數(shù)據(jù)分別對(duì)473、474和475條帶中的殘留相位予以建模，獲得校正后的差分干涉圖，如圖3所示（對(duì)應(yīng)于圖2中的矩形框區(qū)域），對(duì)比校正前后干涉相位，可以看出由于軌道不確定性所導(dǎo)致的殘留相位在校正后相位圖中已得到較好的消除，在473—474條帶的邊界已經(jīng)沒(méi)有明顯的多余條紋.

4 鄰軌平滑校正

對(duì)于PALSAR條帶內(nèi)無(wú)足夠數(shù)量GPS校正點(diǎn)的情況下，如475—476的條帶連接處，本文提出了基于鄰軌間干涉形變差異建模的相對(duì)校正方法.由于主震震級(jí)造成的地表形變遠(yuǎn)大于余震和其他因素引起的形變（如震后余滑和孔隙回彈效應(yīng)在短時(shí)間內(nèi)造成的形變量較小，與主震形變量級(jí)相比不顯著），因此可以認(rèn)為在上述PALSAR影像的時(shí)間基線跨度內(nèi)，相鄰條帶重疊數(shù)據(jù)具有平滑一致的形變量.以GPS校正后的干涉圖為基準(zhǔn)，對(duì)鄰軌干涉圖與基準(zhǔn)干涉圖之間的形變差異進(jìn)行建模，從兩幅鄰軌重疊影像中提取高相關(guān)性的目標(biāo)點(diǎn)，解算形變差異模型系數(shù)，進(jìn)而實(shí)現(xiàn)對(duì)未校正PALSAR條帶形變場(chǎng)與干涉相位的修正.該方法充分利用鄰軌重疊區(qū)域已校正的干涉相位數(shù)據(jù)，對(duì)無(wú)足夠可用數(shù)量GPS形變的干涉條帶進(jìn)行校正，可進(jìn)一步減弱未校正條帶內(nèi)的軌道相位殘差和大氣殘留相位.

圖4為鄰軌條帶的振幅影像，藍(lán)色虛線方框?yàn)橄噜彈l帶間的重疊部分，顧及條帶重疊區(qū)的同震形變應(yīng)保持一致性的原則，以其中一幅經(jīng)GPS校正后的干涉數(shù)據(jù)作為基準(zhǔn)，對(duì)另一幅未校正干涉條帶相位值進(jìn)行建模修正.例如，475條帶已經(jīng)過(guò)GPS校正，且與474和473形變條紋保持較好的一致性，則以475為基準(zhǔn)對(duì)鄰軌的476干涉條紋進(jìn)行建模修正，具體算法如下.

首先在鄰軌PALSAR影像的重疊區(qū)進(jìn)行同名點(diǎn)識(shí)別，基于相關(guān)系數(shù)最大值匹配方法開(kāi)展重疊區(qū)同名點(diǎn)搜索，結(jié)果如圖4所示，圖中的藍(lán)色和白色點(diǎn)代表所識(shí)別的同名點(diǎn).為保證后續(xù)建模的可靠性，提取相關(guān)系數(shù)大于0.4的同名點(diǎn)（圖中白色點(diǎn)）用于數(shù)學(xué)建模，計(jì)算高相關(guān)性同名點(diǎn)在兩個(gè)條帶中的形變差，依據(jù)其在未校正條帶中的方位向和距離向坐標(biāo)，構(gòu)建如下式的形變校正模型：式中，Δd表示兩條帶同名點(diǎn)形變量之差，l、p分別為同名點(diǎn)在被校正條帶中的方位向和距離向坐標(biāo)，bi（i＝0，1，2，3，4，5）為待求模型參數(shù).當(dāng)高質(zhì)量的同名點(diǎn)數(shù)量大于6，使用最小二乘方法解算出模型系數(shù)，然后使用該模型對(duì)鄰軌條紋相位進(jìn)行數(shù)據(jù)校正，獲得校正后的差分干涉圖，圖5為條帶475與476經(jīng)鄰軌校正后的干涉相位圖，從圖5中可以看出，條帶之間的不連續(xù)性問(wèn)題已得到較好解決，其鄰軌干涉條紋連續(xù)性較好.

對(duì)上述校正后的差分干涉圖進(jìn)行相鄰條帶相位數(shù)據(jù)拼接，獲得完整的地表形變干涉相位圖，如圖6所示，使用Snaphu軟件對(duì)拼接后的干涉相位進(jìn)行全局相位解纏和絕對(duì)形變解算（Chen and Zebker，2001；Chen and Zebker，2002；陳強(qiáng)等，2012b），獲得地表沿衛(wèi)星雷達(dá)視線方向（LOS）的形變場(chǎng)，顧及圖2中的JYAN GPS基站所在的InSAR像元具有較高的相干性，且遠(yuǎn)離斷層破裂線，選取該GPS觀測(cè)點(diǎn)作為相位解纏的形變基準(zhǔn)點(diǎn)，最終獲得PALSAR影像區(qū)域完整的同震形變場(chǎng)，如圖7所示.

圖4 鄰軌PALSAR振幅影像重疊區(qū)與聯(lián)系點(diǎn)Fig.4 Spatial overlapping of adjacent PALSAR amplitude maps and detected tie points

圖5 鄰軌約束校正前后的干涉相位（a）校正前的干涉圖；（b）校正后的干涉圖.Fig.5 Local interferogram before（a）and after（b）correction using adjacent geometric constraints

圖6 —7所示的校正后InSAR同震形變場(chǎng)與已有研究成果相比（單新建等，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010），整體趨勢(shì)基本一致，例如，觀測(cè)到最大的地表形變場(chǎng)位于北川區(qū)域，LOS向形變量級(jí)在1.2m左右，汶川—映秀段斷層上盤(pán)存在量級(jí)0.6～0.8m的顯著形變區(qū)，成都平原存在厘米量級(jí)的地表形變（Hashimoto etal.，2009；Liu etal.，2010；Shen etal.，2009；Furuya etal.，2010；Tong etal.，2010）.

然而，校正后的InSAR形變場(chǎng)在局部細(xì)節(jié)區(qū)域呈現(xiàn)出顯著的不同，首先，已有研究中條帶473存在的嚴(yán)重軌道誤差在本文中通過(guò)GPS校正方法得以很好去除，并在此基礎(chǔ)上，使用473條帶的觀測(cè)數(shù)據(jù)對(duì)條帶472和471進(jìn)行了鄰軌平滑校正，顯著減弱了條帶471—472中存在的大氣、軌道等系統(tǒng)誤差.此外，依據(jù)鄰軌校正方法，條帶476也基于475進(jìn)行了誤差校正.校正后的干涉相位較好地解決了已有研究成果中顯著存在的鄰軌拼接不連續(xù)現(xiàn)象，使得最終提取的形變場(chǎng)具有較好的鄰軌平滑性，有利于后續(xù)斷層滑動(dòng)分布模型的精確反演.

圖6 校正后的InSAR同震形變干涉條紋Fig.6 Corrected InSAR interferogram related to coseismic deformation

圖7 校正后的汶川地震InSAR同震形變場(chǎng)Fig.7 Corrected InSAR surface deformation field of Wenchuan earthquake

5 同震形變場(chǎng)精度評(píng)定

5.1 絕對(duì)精度分析

為定量分析上述校正后InSAR同震形變場(chǎng)的精度與可靠性，以高精度GPS觀測(cè)值作為參考，對(duì)InSAR形變與同點(diǎn)位GPS形變進(jìn)行差異比較，統(tǒng)計(jì)分析校正后InSAR同震形變場(chǎng)的實(shí)際絕對(duì)精度.

將PALSAR干涉圖中的GPS點(diǎn)位分為兩類(lèi)，干涉相干性在0.3以上的點(diǎn)位劃分為高相干點(diǎn)，低于該相干閾值則為低相干點(diǎn)，將GPS觀測(cè)得到形變位移投影轉(zhuǎn)換到雷達(dá)LOS方向上的位移，定量數(shù)據(jù)比較如表2所示.

從表2中可以看出，經(jīng)軌道殘留相位校正和鄰軌平滑校正后的同震形變場(chǎng)，高相干點(diǎn)的形變量殘差由8.95cm減小至3.59cm，精度提高了60%，低相干點(diǎn)平均殘差由25.21cm減小至15.00cm，精度提高達(dá)40%，表明本文提出的軌道殘留相位建模和鄰軌平滑方法具有較高的可靠性，同震形變場(chǎng)的精度得以顯著提高.

5.2 鄰軌形變的相對(duì)精度分析

對(duì)于地震斷層錯(cuò)動(dòng)引發(fā)的地表形變，在斷層上下盤(pán)錯(cuò)動(dòng)的邊界處通常存在顯著的形變突變，如逆斷層上盤(pán)呈抬升運(yùn)動(dòng)，而下盤(pán)則表現(xiàn)為下降，對(duì)于其余方位的形變場(chǎng)則較為平滑.因此，在InSAR同一個(gè)條帶內(nèi)的同震形變場(chǎng)，除了斷層上下盤(pán)分界處的地表破裂帶區(qū)域形變不連續(xù)之外，其余地方則呈現(xiàn)出較好的平滑度和連續(xù)性.然而，在相鄰兩個(gè)不同PALSAR條帶之間的邊界區(qū)域，則可能出現(xiàn)前述的由于軌道誤差、大氣相位等因素引起的干涉相位或形變場(chǎng)不連續(xù)現(xiàn)象，因此采用條帶邊界形變的平滑程度指標(biāo)，可對(duì)同震形變場(chǎng)的可靠性與平滑度進(jìn)行相對(duì)精度評(píng)價(jià).

本文引入拉普拉斯平滑參數(shù)作為同震形變場(chǎng)連續(xù)性的定量評(píng)價(jià)指標(biāo)，以中心像元S（i，j）為計(jì)算單元，通過(guò)鄰近像元的差分計(jì)算，獲得干涉形變場(chǎng)的拉普拉斯二階導(dǎo)數(shù)，如下式所示：（4）式中，i，j為待計(jì)算像素位置，S（i，j）為像元位置處的形變量，Δl1和 Δl2為SAR分辨單元尺寸，ρi，j為平滑因子，｜ρi，j｜的數(shù)值越小，表示形變場(chǎng)的平滑度越高（Jonsson etal.，2002）.

表2 修正前后的InSAR形變與GPS形變的差異比較Table 2 Displacement comparison between InSAR and GPS derived deformation field

對(duì)校正前后InSAR形變場(chǎng)各相鄰條帶邊界處的相位，使用上述公式（4）計(jì)算其平滑度指標(biāo)，得到如圖8所示的鄰軌邊界平滑指數(shù)分布.

從圖8可以看出，經(jīng)軌道誤差和鄰軌約束校正后，條帶邊界干涉形變場(chǎng)的非零｜ρi，j｜像元數(shù)量顯著減少，邊界平滑度得到明顯改善，對(duì)所有鄰軌邊界像元平滑因子ρi，j進(jìn)行數(shù)據(jù)統(tǒng)計(jì)，直方圖如圖9所示，統(tǒng)計(jì)數(shù)據(jù)表明，校正后形變場(chǎng)平滑因子多趨于零，標(biāo)準(zhǔn)偏差減小了33%.對(duì)平滑因子進(jìn)行分區(qū)間統(tǒng)計(jì)，在平滑因子絕對(duì)值較小的（-40，40）區(qū)間內(nèi)，校正后形變場(chǎng)平滑因子比校正前提高了7.7%，對(duì)于平滑因子絕對(duì)值較大的（-200，-40）、（40，200）區(qū)間，校正后形變場(chǎng)的平滑因子減少了6.8%.綜合以上分析，經(jīng)過(guò)軌道誤差校正和鄰軌約束校正后的形變場(chǎng)，鄰軌邊界具有更好的平滑度，形變條紋的不連續(xù)性得以較好修正.

圖8 同震形變場(chǎng)校正前后的鄰軌邊界平滑度（a）校正前的邊界平滑度；（b）校正后的邊界平滑度.Fig.8 The smoothness of coseismic interferogram respectively before（a）and after（b）correction

圖9 形變場(chǎng)的鄰軌邊界平滑因子統(tǒng)計(jì)（a）校正前鄰軌邊界平滑因子；（b）校正后鄰軌邊界平滑因子.Fig.9 Statistics of adjacent track smoothing factor of derived deformation field

6 斷層滑移模型反演

在精確獲得同震地表形變場(chǎng)基礎(chǔ)上，采用Okada彈性半空間位錯(cuò)模型反演發(fā)震斷層表面的滑移分布（Okada，1985）.通過(guò)對(duì)野外地質(zhì)調(diào)查數(shù)據(jù)分析和綜合已有研究成果（Zhang etal.，2008；Hashimoto etal.，2009；Liu etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010），本文構(gòu)建了如圖10中所示的F1—F6六段斷層分段模型，首先假設(shè)斷層表面為均勻滑動(dòng)，使用非線性遺傳算法對(duì)斷層參數(shù)進(jìn)行優(yōu)化反演，為保證反演結(jié)果的可靠性，我們賦予反演模型不同的初始值，進(jìn)行多次反演搜優(yōu)計(jì)算，并對(duì)反演參數(shù)集進(jìn)行統(tǒng)計(jì)分析，最終獲得了各子斷層的幾何參數(shù)及其對(duì)應(yīng)精度，如表3所示.

表3 汶川地震斷層幾何Table 3 Fault geometry of Wenchuan earthquake

對(duì)于中央斷裂帶（映秀—北川）長(zhǎng)大斷層來(lái)說(shuō)，斷層面的實(shí)際滑動(dòng)往往呈現(xiàn)為非均勻分布，為更精細(xì)獲得局部斷層滑動(dòng)的空間特征，對(duì)表3中各子斷層進(jìn)行4km×4km的格網(wǎng)離散化處理，共獲得1051個(gè)子斷層，以其作為基本單元精細(xì)反演斷層滑動(dòng)的空間分布，采用拉普拉斯平滑約束斷層滑動(dòng)模型，如式（5）所示：

式中，Gi為根據(jù)離散子斷層單位走滑量和逆沖量計(jì)算得到的InSAR數(shù)據(jù)反演格林函數(shù)，Hi為GPS數(shù)據(jù)反演格林函數(shù)，α、β分別對(duì)應(yīng)其權(quán)重因子，Los、def為InSAR和GPS觀測(cè)得到的地表形變量，Li為拉普拉斯約束格林函數(shù)，κ代表平滑因子分別代表斷層的走滑量和逆沖量.

由于GPS形變觀測(cè)在垂直方向的精度較弱，且垂直向觀測(cè)資料較少，因此在反演過(guò)程中僅使用了GPS平面觀測(cè)數(shù)據(jù).InSAR形變觀測(cè)數(shù)據(jù)的條帶471和472在校正后仍可能存在電離層及大氣負(fù)面影響，為避免其對(duì)反演結(jié)果造成不利影響，471—472條帶中形變場(chǎng)連續(xù)性較差的數(shù)據(jù)不參與反演計(jì)算.此外，為盡可能減弱不可靠的形變觀測(cè)結(jié)果對(duì)反演的影響，使用InSAR相干性大于0.3的形變觀測(cè)數(shù)據(jù)，并對(duì)其進(jìn)行四叉樹(shù)降采樣工作，該采樣方法在高形變梯度區(qū)域采樣量大，在低形變梯度區(qū)域采樣數(shù)據(jù)量小，能最大限度地保留形變場(chǎng)的空間特征.降采樣過(guò)程一方面可避免地震遠(yuǎn)場(chǎng)InSAR非構(gòu)造形變?cè)肼晫?duì)反演結(jié)果的影響，另一方面也可有效避免因觀測(cè)數(shù)據(jù)過(guò)多引起的計(jì)算負(fù)荷過(guò)重和計(jì)算效率低下的問(wèn)題，最終的降采樣形變場(chǎng)如圖10所示.

圖10 四叉樹(shù)降采樣InSAR形變場(chǎng)Fig.10 Distribution of InSAR sampled points based on quad-tree method

對(duì)于InSAR和GPS形變觀測(cè)數(shù)據(jù)用于斷層滑移反演的權(quán)重分配，本文將依據(jù)其觀測(cè)精度進(jìn)行確定.根據(jù)前述的形變場(chǎng)精度評(píng)定可知，高質(zhì)量InSAR觀測(cè)數(shù)據(jù)的殘差為0.036m，在實(shí)際反演中使用相干性較高的觀測(cè)點(diǎn)，此處適當(dāng)放寬其精度取為0.04m，對(duì)于GPS形變觀測(cè)數(shù)據(jù)的精度，根據(jù)相關(guān)研究可知其平面精度約為2～5mm，取3mm作為精度指標(biāo)（張培震，2008）.因此，InSAR和GPS點(diǎn)的誤差之比為13.3∶1.0，反比例計(jì)算可得其相應(yīng)權(quán)重分別為α＝0.075、β＝1.0.

在確定了數(shù)據(jù)權(quán)重之后，賦予反演模型連續(xù)變化的平滑因子κ，并依據(jù)反演結(jié)果計(jì)算獲得相應(yīng)斷層面粗糙度與模型殘差的對(duì)應(yīng)關(guān)系，如圖11所示，從圖中可看出，隨著粗糙度的增大，模型殘差急劇減小，顧及模型殘差和斷層粗糙度的折中關(guān)系，我們選取了如圖11中箭頭標(biāo)明處的粗糙度為7.5cm／km處的平滑因子為最優(yōu)的κ取值.

圖11 模型殘差與斷層面粗糙度關(guān)系曲線Fig.11 The trade-off curve with misfit function of InSAR and GPS plotted as function of average solution roughness

圖12 斷層反演滑動(dòng)及其精度分布（a）斷層滑動(dòng)三維分布；（b）斷層滑動(dòng)誤差分布；（c）斷層滑動(dòng)平面分布.Fig.12 Fault slip distribution and its error（a）Three-dimensional distribution of coseismic slip；（b）Errors of coseismic slip；（c）Coseismic slip distribution on fault surface.

確定了各項(xiàng)模型參數(shù)后，在非負(fù)最小二乘原則下進(jìn)行斷層滑動(dòng)量的反演計(jì)算，獲得汶川地震斷層滑動(dòng)的空間分布，如圖12所示，其中圖12a為反演斷層滑動(dòng)的三維分布.為了驗(yàn)證滑動(dòng)反演結(jié)果的可靠性，我們基于蒙特卡洛算法生成100組帶隨機(jī)擾動(dòng)的觀測(cè)數(shù)據(jù)集，并基于這些數(shù)據(jù)反演相應(yīng)的斷層滑動(dòng)分布，對(duì)所有計(jì)算得到的斷層滑動(dòng)分布集進(jìn)行統(tǒng)計(jì)分析，從而估算出反演模型的精度分布，如圖12b所示.統(tǒng)計(jì)表明反演模型的平均精度為0.014m，表明圖12a的反演結(jié)果是可靠的.同時(shí)需要指出的是，最大誤差分布位于青川地區(qū)16km深度處，分析此現(xiàn)象的可能原因，一方面是該區(qū)域的InSAR觀測(cè)數(shù)據(jù)由于存在大氣影響而未用于參與滑動(dòng)反演，導(dǎo)致反演模型對(duì)此部分?jǐn)鄬踊瑒?dòng)的約束較弱；另一方面是汶川主震后存在密集余震，InSAR觀測(cè)的滯后性使得相鄰條帶形變場(chǎng)中包含了不同余震的形變信息，對(duì)反演精度造成了一定影響.

為清晰顯示斷層面的滑動(dòng)分布狀態(tài)，我們繪制了如圖12c所示的斷層滑動(dòng)平面分布圖，從中可看出，斷層滑動(dòng)整體表現(xiàn)為逆沖為主兼具右旋走滑特性，平面滑動(dòng)量主要分布于映秀、北川和青川區(qū)域，斷層最大滑動(dòng)量為9.0m，位于北川地表附近，導(dǎo)致北川地區(qū)地表破裂最為嚴(yán)重.在地殼深度方向，地表以下0～16km范圍內(nèi)孕育了較大的滑動(dòng)成分，而16km以下滑動(dòng)量分布并不顯著.

由于斷層滑動(dòng)反演所使用的地表InSAR形變場(chǎng)約束以及斷層幾何模型等的不同，本文的反演結(jié)果與已有研究成果存在一定差異.首先，斷層北川段的滑動(dòng)量在沿?cái)鄬幼呦蚍植几鼜V，但在深度方向上相對(duì)較淺，而映秀—汶川段的平均滑動(dòng)量級(jí)較已有研究成果偏?。℉ashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010）；此外，青川段滑動(dòng)量分布更為深入地下，這可能與青川段震后密集的余震分布相關(guān)（朱艾斕等，2008）.

根據(jù)上述反演出的斷層滑移模型與斷層幾何參數(shù)，正演計(jì)算汶川地震的InSAR地表形變場(chǎng)，如圖13中的左下插圖（a），與圖6中校正后InSAR同震形變場(chǎng)相比，兩者具有較高的一致性，尤其是在映秀和北川等形變嚴(yán)重的區(qū)域，正演數(shù)據(jù)與InSAR實(shí)際觀測(cè)值吻合的較好.但需要說(shuō)明的是，在主斷層北端的青川區(qū)域，由于InSAR觀測(cè)數(shù)據(jù)存在嚴(yán)重大氣干擾負(fù)面影響，形變觀測(cè)值可靠性相對(duì)較弱，未參與到滑動(dòng)分布模型的反演，使該部分?jǐn)?shù)據(jù)與正演形變場(chǎng)存在顯著差異.

進(jìn)一步分析同震形變場(chǎng)的模型反演殘差與精度，圖13為GPS的反演殘差向量（各矢量箭頭）與InSAR的反演殘差數(shù)值分布.圖（b）與（c）為GPS觀測(cè)值與模型計(jì)算值的對(duì)比圖，GPS在北東方向的反演殘差分別為4.3cm和6.2cm.此外，統(tǒng)計(jì)表明約有75%的GPS平面殘差小于2cm，且僅在近斷層附近存在殘差較大的現(xiàn)象.GPS的整體形變殘差為5.5cm，與同類(lèi)已有研究基本相當(dāng).然而，本文InSAR反演形變殘差有顯著改善，從圖13中的InSAR殘差分布中可看出，大部分區(qū)域殘差小于10cm，形變殘差統(tǒng)計(jì)值達(dá)到9.2cm，與已有同類(lèi)研究成果（平均殘差約15cm）對(duì)比（Zhang etal.，2008；Hashimoto etal.，2009；Liu etal.，2009；Shen etal.，2009；Furuya etal.，2010；Tong etal.，2010），InSAR形變反演精度提高約30%，驗(yàn)證了InSAR同震形變場(chǎng)的高精度特性與可靠性，同時(shí)也表明前述InSAR形變場(chǎng)校正方法的有效性與準(zhǔn)確性.

需要指出的是，斷層滑動(dòng)模型的反演殘差一方面來(lái)源于地表形變場(chǎng)約束的觀測(cè)誤差（包括InSAR形變場(chǎng)中殘留的大氣、軌道誤差等），另一方面，主要受限于不確定的復(fù)雜斷層幾何、位錯(cuò)理論模型的系統(tǒng)誤差等.本文斷層滑動(dòng)反演的統(tǒng)計(jì)殘差為9.2cm，仍大于InSAR高質(zhì)量點(diǎn)3.6cm的觀測(cè)校正殘差，這主要是由于斷層幾何參數(shù)的精度、以及Okada位錯(cuò)模型未考慮地殼不均勻特性所帶來(lái)的系統(tǒng)誤差.其中，斷層幾何對(duì)反演結(jié)果影響最為顯著，已有研究均假定斷層幾何在地殼不同深度具有相同的斷層傾角，這與實(shí)際斷層地下延伸模式可能存在差異.因而，即使在獲得高精度地表形變場(chǎng)的情況下，對(duì)于龍門(mén)山斷裂帶復(fù)雜的斷層幾何來(lái)說(shuō)，仍可能出現(xiàn)相對(duì)較大的反演殘差.但相對(duì)于已有研究平均15cm的反演殘差，本文反演精度（9.2cm）有較大提高，這也正是InSAR形變場(chǎng)校正對(duì)于改善滑動(dòng)模型反演精度的重要貢獻(xiàn)，所獲得的斷層滑動(dòng)精細(xì)分布有助于進(jìn)一步推進(jìn)對(duì)汶川地震斷層破裂機(jī)制的認(rèn)識(shí)和深入理解.

依據(jù)上述反演所獲得的斷層幾何參數(shù)和滑動(dòng)分布模型，計(jì)算得到本次地震的地震矩為8.0469×1020N·m，對(duì)應(yīng)的矩震級(jí)Mw7.9，與已有反演結(jié)果基本相當(dāng)（Tong etal.，2010），但本文計(jì)算得到的地震矩稍顯偏大，一方面是因?yàn)楸疚氖褂昧私?jīng)鄰軌校正后的精細(xì)化地表同震形變場(chǎng)，另一方面因?yàn)楸疚乃玫臄鄬幽Ｐ蛶缀蝺A角與已有方法存在一定差異.

圖13 InSAR和GPS形變的聯(lián)合反演殘差（a）模擬InSAR干涉圖；（b）GPS北方向形變與模型正演對(duì)比；（c）GPS東方向形變與模型正演對(duì)比.Fig.13 Residuals of InSAR and GPS joint inversion（a）Simulated InSAR interferogram；（b）Model prediction versus GPS observation in north direction；（c）Model prediction versus GPS observation in east direction.

7 結(jié)論

本文采用軌道誤差校正方法與鄰軌平滑技術(shù)對(duì)汶川地震形變場(chǎng)進(jìn)行修正處理，獲得汶川地震精細(xì)的同震形變場(chǎng).GPS形變觀測(cè)與InSAR同震形變的殘差分析表明，校正后形變場(chǎng)高質(zhì)量數(shù)據(jù)殘差為3.6cm，精度提高約60%，低質(zhì)量數(shù)據(jù)精度提高約40%；鄰軌邊界形變場(chǎng)平滑度分析表明，校正后形變場(chǎng)邊界平滑因子的標(biāo)準(zhǔn)差減小33%，驗(yàn)證了本文提出的同震形變場(chǎng)校正方法的有效性和可靠性.

基于Okada彈性半空間位錯(cuò)模型的斷層滑動(dòng)反演結(jié)果表明，斷層滑動(dòng)量在沿?cái)鄬幼呦蛏现饕植加谟承?、北川和青川區(qū)域，在深度方向上主要集中在地殼淺部的0～16km范圍內(nèi)，最大滑動(dòng)量為9.0m，位于北川縣城地表附近，InSAR反演后的模型殘差達(dá)9.2cm，反演精度約有30%的提高，由此計(jì)算出的地震矩?cái)?shù)值為8.0469×1020N·m.

Chen C W，Zebker H A.2001.Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization.Journal of the Optical Society of America A，18（2）：338-351，doi：10.1364／JOSAA.18.000338.

Chen C W，Zebker H A.2002.Phase unwrapping for large SAR interferograms：statistical segmentation and generalized network models.IEEE Transactions on Geoscience and Remote Sensing，40（8）：1709-1719，doi：10.1109／TGRS.2002.802453.

Chen Q，Liu G X，Hu Z Q，etal.2012a.Mapping ground 3-D displacement with GPS and PS-InSAR networking in the Pingtung area，southwestern Taiwan，China.Chinese Journal of Geophysics（in Chinese），55（10）：3248-3258，doi：10.6038／j.issn.0001-5733.2012.10.007.

Chen Q，Yang Y H，Liu G X，etal.2012b.InSAR phase unwrapping using least squares method with integer ambiguity resolution and edge detection.Acta Geodaetica et Cartographica Sinica（in Chinese），41（3）：441-448.

Furuya M，Kobayashi T，Takada Y，etal.2010.Fault source modeling of the 2008Wenchuan Earthquake based on ALOS／PALSAR data.Bulletin of the Seismological Society of America，100（5B）：2750-2766，doi：10.1785／0120090242.

Hashimoto M，Enomoto M，F(xiàn)ukushima Y，etal.2009.Coseismic deformation from the 2008Wenchuan，China，Earthquake derived from ALOS／PALSAR images.Tectonophysics，491（1-4）：59-71，doi：10.1016／j.tecto.2009.08.034.

Jonsson S，Zebker H A，Segall P，etal.2002.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，92（4）：1377-1389，doi：10.1785／0120000922.

Li Y，Huang R Q，Densmore A L，etal.2009.Basic features and research progresses of Wenchuan Ms8.0earthquake.Journal of Sichuan University（in Chinese），41（3）：7-25.

Liu G X，Li J，Xu Z，etal.2010.Surface deformation associated with the 2008 Ms8.0Wenchuan earthquake from ALOS L-band SAR interferometry.International Journal of Applied Earth Observation and Geoinformation，12（6）：496-505，doi：10.1016／j.jag.2010.05.005.

Liu Z J，Zhang Z，Wen L，etal.2009.Co-seismic ruptures of the 12May 2008，Ms8.0Wenchuan earthquake，Sichuan：Eastwest crustal shortening on oblique，parallel thrusts along the eastern edge of Tibet.Earth Planet.Sci.Lett.，286（3-4）：355-370，doi：10.1016／j.epsl.2009.07.017.

Okada Y.1985.Surface deformation due to shear and tensile faults in a half-space.Bull.Seismol.Soc.Am，75（4）：1135-1154.

Shan X J，Qu C Y，Song X G.2009.Coseismic surface deformation caused by the Wenchuan Ms8.0earthquake from InSAR data analysis.Chinese Journal of Geophysics（in Chinese），52（2）：496-504.

Shen Z K，Sun J B，Zhang P Z，etal.2009.Slip maxima at fault junctions and rupturing of barriers during the 2008Wenchuan earthquake.Nature，10：718-724，doi：10.1038／ngeo636.

Tong X P，Sandwell D T，F(xiàn)ialko Y，etal.2010.Coseismic slip model of the 2008Wenchuan earthquake derived from joint inversion of interferometric synthetic aperture radar，GPS，and field data.Journal of Geophysical Research，115（B4）：1-19，doi：10.1029／2009JB006625.

Xu X W，Zhang P Z，Wen X Z，etal.2005.Features of active tectonics and recurrence behaviors of strong earthquakes in the western Sichuan province and its adjacent regions.Seismology and Geology（in Chinese），27（3）：446-461.

Zhang L，Wu J C，Chen Y L，etal.2007.Determination of coseismic displacement field of Taiwan Chi-Chi earthquake based on Doris.Journal of Geodesy and Geodynamics（in Chinese），27（3）：18-24.

Zhang P Z.2008.Co-seismic deformation field of 2008 Ms8.0 Wenchuan earthquake revealed by GPS data.Science in China Press（in Chinese），38（10）：1195-1206.

Zhang Y H，Gong W Y，Zhang J X，etal.2008.Measuring coseismic deformation of the sichuan earthquake by satellite differential INSAR.International Conference on Earth Observation Data Processing and Analysis，7285：4F-1：4F-8，doi：10.1117／12.816004.

Zhu A L，Xu X W，Diao G L，etal.2008.Relocation of the Ms8.0 Wenchuan earthquake sequence in part：preliminary seismotectonic analysis.Seismology and Geology（in Chinese），30（3）：759-767，doi：10.3969／j.issn.0253-4967.2008.03.014.

附中文參考文獻(xiàn)

朱艾斕，徐錫偉，刁桂苓等.2008.汶川Ms8.0地震部分余震重新定位及地震構(gòu)造初步分析.地震地質(zhì)，30（3）：759-767，doi：10.3969／j.issn.0253-4967.2008.03.014.

張培震.2008.GPS測(cè)定的2008年汶川Ms8.0級(jí)地震的同震位移場(chǎng).中國(guó)科學(xué)，38（10）：1195-1206.

張磊，伍吉倉(cāng)，陳艷玲等.2007.基于Doris的臺(tái)灣集集地震同震位移場(chǎng)計(jì)算.大地測(cè)量與地球動(dòng)力學(xué)，27（3）：18-24.

徐錫偉，張培震，聞學(xué)澤等.2005.川西及其鄰近地區(qū)活動(dòng)構(gòu)造基本特征與強(qiáng)震復(fù)發(fā)模型.地震地質(zhì)，27（3）：446-461.

單新建，屈春燕，宋小剛等.2009.汶川Ms8.0級(jí)地震InSAR同震形變場(chǎng)觀測(cè)與研究.地球物理學(xué)報(bào)，52（2）：496-504.

李勇，黃潤(rùn)秋，Densmore A L等.2009.汶川8.0級(jí)地震的基本特征及其研究進(jìn)展.四川大學(xué)學(xué)報(bào).41（3）：7-25.

陳強(qiáng)，楊瑩輝，劉國(guó)祥等.2012（b）.基于邊界探測(cè)的InSAR最小二乘整周相位解纏方法.測(cè)繪學(xué)報(bào)，41（3）：441-448.

陳強(qiáng)，劉國(guó)祥，胡植慶等.2012（a）.GPS與PS-InSAR聯(lián)網(wǎng)監(jiān)測(cè)的臺(tái)灣屏東地區(qū)三維地表形變場(chǎng).地球物理學(xué)報(bào)，55（10）：3248-3258，doi：10.6038／j.issn.0001-5733.2012.10.007.