InSAR同震形變場(chǎng)及其在震源參數(shù)確定中的應(yīng)用研究進(jìn)展

季靈運(yùn)，朱良玉，劉傳金，張文婷，邱江濤，徐曉雪

(1.防災(zāi)科技學(xué)院 地球科學(xué)學(xué)院,河北 三河 065201；2.中國(guó)地震局第二監(jiān)測(cè)中心,陜西 西安 710054)

0 引 言

從大地測(cè)量的角度，大地震發(fā)生后，地震科技工作者需要回答兩個(gè)問(wèn)題，一是同震形變的空間范圍、量級(jí)和特征，二是確定發(fā)震斷層的幾何特征和運(yùn)動(dòng)性質(zhì)。要回答這兩個(gè)問(wèn)題，首先需要獲取地震引起的地表形變場(chǎng)，得到地震同震形變空間展布范圍，在此基礎(chǔ)上反演發(fā)震斷層的震源參數(shù)，確定發(fā)震斷層性質(zhì)。從活動(dòng)構(gòu)造的角度，地震震源參數(shù)能夠反映發(fā)震斷層的產(chǎn)狀和運(yùn)動(dòng)性質(zhì)，為認(rèn)識(shí)區(qū)域構(gòu)造應(yīng)力場(chǎng)和斷層運(yùn)動(dòng)特征提供最直接的證據(jù)。從解剖地震的角度，獲得精細(xì)的震源參數(shù)，對(duì)于認(rèn)識(shí)和理解地震孕育過(guò)程和發(fā)生機(jī)理具有重要意義。一次大地震的發(fā)震斷層往往涉及到一條斷層的多個(gè)斷層段或者多條斷層，大地測(cè)量數(shù)據(jù)能夠捕捉發(fā)震斷層的近場(chǎng)運(yùn)動(dòng)特征，可以識(shí)別這種復(fù)雜的破裂特征[1]。合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar，InSAR)技術(shù)能夠獲取高空間分辨率的地殼形變場(chǎng)，可為精確確定斷層破裂特征提供重要約束。本文首先梳理了利用InSAR技術(shù)獲取同震形變場(chǎng)的研究現(xiàn)狀，總結(jié)了InSAR技術(shù)獲取同震形變場(chǎng)的優(yōu)缺點(diǎn)、InSAR技術(shù)獲取三維同震形變場(chǎng)的現(xiàn)狀；回顧了地震震源參數(shù)正反演模型及方法的研究歷程，包括從簡(jiǎn)單的彈性半空間位錯(cuò)模型到接近真實(shí)地球介質(zhì)的數(shù)值模型，分析了不同反演方法和先驗(yàn)約束對(duì)結(jié)果的影響；最后，從同震形變場(chǎng)分離、自動(dòng)化InSAR數(shù)據(jù)處理與震源參數(shù)反演等方面進(jìn)行了展望。

1 InSAR技術(shù)獲取同震形變場(chǎng)

InSAR技術(shù)是20世紀(jì)70年代發(fā)展起來(lái)的一種空間大地測(cè)量技術(shù)，其利用衛(wèi)星發(fā)射電磁波在地表的反射信號(hào)進(jìn)行相位差分，計(jì)算地表微小變形[2]。1993年Massonnet等利用兩景存檔的ERS-1數(shù)據(jù)成功獲取了美國(guó)Landers地震的同震形變場(chǎng)，其結(jié)果與野外實(shí)測(cè)數(shù)據(jù)一致性較好且擁有更高的空間分辨率，豐富了觀測(cè)信息的同時(shí)極大地降低了觀測(cè)成本，開啟了InSAR技術(shù)在同震斷層運(yùn)動(dòng)領(lǐng)域應(yīng)用的先河[2-3]。在中國(guó)，InSAR技術(shù)在同震形變測(cè)量中也發(fā)揮了不可替代的作用[4-17]。

1.1 InSAR技術(shù)獲取同震形變場(chǎng)的優(yōu)點(diǎn)

InSAR技術(shù)是一種主動(dòng)遙感技術(shù)。雷達(dá)傳感器通過(guò)發(fā)射和接收電磁波信號(hào)以獲取SAR影像，不依賴可見(jiàn)光，具有全天時(shí)的優(yōu)勢(shì)；因?yàn)槔走_(dá)波波長(zhǎng)為厘米級(jí)，可以穿透雨霧云雪，所以InSAR技術(shù)具有全天候的優(yōu)勢(shì)[2,18]；并且，利用InSAR技術(shù)獲得的地表形變場(chǎng)精度較高，可達(dá)毫米級(jí)[19]。由于SAR影像空間分辨率可達(dá)米級(jí)，影像幅寬達(dá)到數(shù)百千米，所以InSAR與水準(zhǔn)、GNSS采集的離散點(diǎn)位信息相比，具有高空間分辨率、廣域覆蓋的優(yōu)勢(shì)。并且，隨著雷達(dá)衛(wèi)星的不斷發(fā)射升空，同一地區(qū)可以獲取多平臺(tái)、多軌道和多時(shí)相的SAR數(shù)據(jù)，使InSAR技術(shù)獲取同震形變場(chǎng)的時(shí)間間隔為1～12 d[20]。以歐洲航天局哨兵影像為代表的SAR影像可免費(fèi)獲取，大大降低了InSAR地震監(jiān)測(cè)的成本[10,21-22]。因此，InSAR技術(shù)獲取同震形變場(chǎng)具有時(shí)效性和便利性。綜上所述，在同震形變監(jiān)測(cè)領(lǐng)域，InSAR技術(shù)具有以下優(yōu)勢(shì)：第一，一般情況下，InSAR技術(shù)可以獲取斷層近場(chǎng)和近斷層的同震形變場(chǎng)，而近斷層形變能夠精確刻畫斷層地表破裂特征、走向和空間偏移量分布模式，在一定程度上反映了斷層的物理性質(zhì)；第二，隨著越來(lái)越多SAR衛(wèi)星的發(fā)射，InSAR技術(shù)獲取同震形變場(chǎng)的時(shí)間間隔不斷縮短，使得InSAR技術(shù)觀測(cè)的同震形變場(chǎng)中包含的震后形變分量越來(lái)越少，可以獲得更加精確的同震破裂，減少震后形變的貢獻(xiàn)；第三，由于能夠進(jìn)行空間連續(xù)面狀觀測(cè)，所以InSAR技術(shù)可以獲取相鄰斷層的活動(dòng)信息。

1.2 InSAR技術(shù)獲取同震形變場(chǎng)的缺點(diǎn)

InSAR技術(shù)受到自身?xiàng)l件的限制，在獲取同震形變場(chǎng)的實(shí)際應(yīng)用中也存在一些缺點(diǎn)。

第一，InSAR技術(shù)僅能獲取雷達(dá)視線方向(Line-of-sight，LOS)的一維形變信息。由于受到SAR衛(wèi)星極軌飛行和側(cè)視成像的制約，InSAR技術(shù)觀測(cè)到的形變信息是地表真實(shí)形變?cè)诶走_(dá)視線方向上的投影，且對(duì)南北向形變不敏感，所以僅靠單一軌道的InSAR同震形變有可能對(duì)斷層運(yùn)動(dòng)性質(zhì)誤判甚至錯(cuò)判[18,20,23-24](圖1)。例如，四川汶川8.0級(jí)地震的實(shí)例顯示，該發(fā)震斷層的逆沖運(yùn)動(dòng)與右旋走滑運(yùn)動(dòng)在InSAR圖像上表現(xiàn)為兩種相反的特征，這也是本次地震中InSAR技術(shù)不能給出確切最大位錯(cuò)的主要原因[25]。

圖1 逆沖型和走滑型同震形變?cè)贚OS向的表現(xiàn)示意圖Fig.1 Sketch Views Showing Radar’s LOS Displacement Caused by Coseismic Slip on a Thrust Fault and a Strike-slip Fault

第二，InSAR相位失相干。由于InSAR技術(shù)是利用雷達(dá)波在地面的反射波進(jìn)行相位差分，地表覆蓋物(如茂密植被)容易造成雷達(dá)波的散射，所以時(shí)間間隔稍長(zhǎng)就可能引起相位失相干而無(wú)法獲得可靠的干涉測(cè)量結(jié)果。此外，地震前、后獲取的兩幅SAR影像軌道之間距離越長(zhǎng)，會(huì)導(dǎo)致干涉相位噪聲的水平越高，嚴(yán)重的會(huì)導(dǎo)致無(wú)法形成有效的干涉測(cè)量[26-27]。2008年新疆烏恰6.8級(jí)地震[28]、2008年新疆于田地震[29-30]、2008年四川汶川8.0級(jí)地震[31]等都存在震中區(qū)域嚴(yán)重的失相干。為了減少失相干的影響，通常盡量選擇時(shí)間跨度小、空間垂直基線較短的震前、震后兩景SAR影像進(jìn)行差分干涉，并對(duì)得到的干涉圖進(jìn)行多視和濾波；對(duì)于震級(jí)較大的地震，也可以選取波長(zhǎng)更長(zhǎng)的L波段ALOS-1/2衛(wèi)星數(shù)據(jù)以減少失相干。此外，多時(shí)相InSAR技術(shù)(Multi-temporal InSAR，MT-InSAR)有時(shí)也被用于同震形變觀測(cè)，可以有效減少噪聲影響[32-35]。

第三，InSAR大氣相位延遲(Atmospheric Phase Screen,APS)。SAR衛(wèi)星信號(hào)在穿透大氣時(shí)，對(duì)流層中水汽、電離層中帶電物質(zhì)會(huì)對(duì)雷達(dá)波的傳播速度和傳播路徑產(chǎn)生影響，反映在干涉圖中為大氣相位延遲，降低同震形變監(jiān)測(cè)的精度。例如，2017年中國(guó)新疆精河6.6級(jí)地震的SAR數(shù)據(jù)受到大氣噪聲的嚴(yán)重影響，其中一個(gè)干涉圖的同震位移幾乎無(wú)法分辨[36]。2016年中國(guó)臺(tái)灣美濃6.4級(jí)地震和2015年尼泊爾8.1級(jí)地震的SAR干涉圖顯示，電離層對(duì)在低緯度的升軌數(shù)據(jù)影響較大[37-38]。對(duì)流層延遲是大氣延遲的主要組成部分，其具有隨機(jī)性，難以使用解析法直接估計(jì)。對(duì)于同震形變監(jiān)測(cè)，目前常用的對(duì)流層延遲估計(jì)方法主要有：①經(jīng)驗(yàn)改正法，基于對(duì)流層延遲與地形的強(qiáng)相關(guān)關(guān)系，建立地形相關(guān)模型估計(jì)對(duì)流層延遲中的垂直分層延遲；②基于外部數(shù)據(jù)法，包括地面站的氣象觀測(cè)數(shù)據(jù)、GPS數(shù)據(jù)、紅外輻射數(shù)據(jù)及氣象模型數(shù)據(jù)等，其中的GACOS(Generic Atmospheric Correction Online Service for InSAR)大氣產(chǎn)品目前已廣泛應(yīng)用于地震監(jiān)測(cè)中，該產(chǎn)品利用歐洲中期天氣預(yù)報(bào)中心(European Centre for Medium-range Weather Forecasts,ECMWF)大氣產(chǎn)品和 GNSS數(shù)據(jù)生成高精度大氣對(duì)流延遲相位，具有全球覆蓋、近乎實(shí)時(shí)運(yùn)作的特性。與對(duì)流層延遲不同，電離層誤差對(duì) LOS向和方位向形變觀測(cè)都有影響。目前，電離層延遲校正方法有利用電離層模型、方位向和距離向偏移量估計(jì)方法和電離層模型的法拉第效應(yīng)估計(jì)與去除等[39-42]。

第四，單次地震同震形變場(chǎng)難以有效分離。由于InSAR技術(shù)是根據(jù)地震前、后的兩景SAR影像獲取同震形變場(chǎng)，而SAR衛(wèi)星飛行具有周期性，所以同震形變場(chǎng)往往包含了影像獲取時(shí)間間隔內(nèi)發(fā)生的所有形變。對(duì)于震群型地震或有較大震級(jí)余震的地震，無(wú)法分離出主震的形變場(chǎng)，即經(jīng)常會(huì)對(duì)同震引起的形變估計(jì)過(guò)高。例如，1996年喀喇昆侖山西口7.1級(jí)地震和2008年西藏改則6.4級(jí)地震的InSAR同震形變場(chǎng)均包含了強(qiáng)余震形變[43-44]。

第五，受制于D-InSAR測(cè)量的精度(厘米級(jí))，InSAR技術(shù)一般只能用于5級(jí)以上地震的同震形變監(jiān)測(cè)，難以觀測(cè)到地表變形過(guò)小的低震級(jí)地震。

1.3 InSAR技術(shù)獲取三維同震形變場(chǎng)重建方法

假定遠(yuǎn)離SAR衛(wèi)星的運(yùn)動(dòng)在LOS向的觀測(cè)結(jié)果為負(fù)，反之為正，則根據(jù)SAR成像幾何關(guān)系，地面上某一點(diǎn)在LOS向(DLOS)和方位向(DAZ)上的觀測(cè)值與其在U、E、N三方向上的形變分量(Du、De、Dn)的關(guān)系表達(dá)式為[45]

(1)

式中：θ為雷達(dá)入射角；α為衛(wèi)星飛行的方位角。

根據(jù)式(1)，當(dāng)獲取某點(diǎn)3個(gè)及以上線性無(wú)關(guān)的LOS向或方位向觀測(cè)值，即可求解相應(yīng)的三維形變分量[46]。近年來(lái)，隨著InSAR技術(shù)的進(jìn)步，三維同震形變場(chǎng)獲取技術(shù)也發(fā)展較快，主要包括多視角SAR影像法、聯(lián)合升降軌D-InSAR形變與方位向形變法、借助外部數(shù)據(jù)或形變模型法等。這些方法各有優(yōu)缺點(diǎn)，需要視情況采取不同的方法。

1.3.1 多視角SAR影像獲取三維同震形變場(chǎng)

1.3.2 聯(lián)合升降軌D-InSAR形變與方位向形變獲取三維同震形變場(chǎng)

(1)聯(lián)合升降軌D-InSAR和像素偏移量估計(jì)法獲取三維同震形變場(chǎng)。像素偏移量估計(jì)法獲取同震形變的原理是基于SAR影像幅度信息，利用互相關(guān)(Cross-correlation)或相干追蹤(Coherence Tracking)技術(shù)計(jì)算地震前、后兩景影像的偏移量，進(jìn)而獲得方位向和距離向二維形變場(chǎng)。由于不需要相位信息，所以像素偏移量估計(jì)法不受失相干影響，測(cè)量精度約為像元分辨率的1/10[48]；受SAR影像空間分辨率的制約，像素偏移量估計(jì)法的監(jiān)測(cè)精度比D-InSAR技術(shù)低一個(gè)數(shù)量級(jí)，目前主要用于監(jiān)測(cè)形變量級(jí)較大的同震形變。Michel等最早將像素偏移量估計(jì)法用于同震形變計(jì)算，獲取了1992年美國(guó)Landers地震的方位向和距離向形變[49]。為了抑制SAR影像中斑點(diǎn)噪聲對(duì)像素偏移量估計(jì)法計(jì)算精度的影響，Wang等對(duì)震前和震后的多景SAR影像取平均，這樣在不降低空間分辨率的情況下提高SAR影像質(zhì)量，提高了像素偏移量估計(jì)法監(jiān)測(cè)形變的精度[50]。目前，聯(lián)合升降軌D-InSAR和像素偏移量估計(jì)法在獲取大地震三維同震形變場(chǎng)中得到了廣泛應(yīng)用，例如2016年新西蘭Kaikoura 7.8級(jí)地震[51]。

另一方面，由于D-InSAR技術(shù)可以獲取精度更高的距離向形變，像素偏移量估計(jì)法獲得的距離向偏移量估計(jì)結(jié)果常被丟棄；值得注意的是，當(dāng)?shù)乇砦灰七^(guò)大或形變不連續(xù)時(shí)，距離向偏移量也可以作為D-InSAR觀測(cè)結(jié)果的有效補(bǔ)充[52]。此外，像素偏移量估計(jì)法估計(jì)得到的位移還可以作為先驗(yàn)信息，為確定斷層幾何參數(shù)提供先驗(yàn)約束[53]或在相位變化梯度較大處輔助相位解纏[54]。

(2)聯(lián)合D-InSAR和多孔徑SAR干涉法獲取三維同震形變場(chǎng)。多孔徑SAR干涉法利用SAR傳感器前視、后視成像的不同姿態(tài)獲取子孔徑前視、后視的干涉信號(hào)，計(jì)算得出地面目標(biāo)的方位向位移，是繼像素偏移量估計(jì)法后又一可獲得SAR方位向形變信息的方法。多孔徑SAR干涉法最早由Bechor 等提出，用于獲取1999年美國(guó)Hector礦地震的二維形變場(chǎng)[55]。其基本原理為：首先采用方位向帶通濾波對(duì)原始SAR影像進(jìn)行分頻，獲得子孔徑前視和后視SLC影像，計(jì)算子孔徑SAR影像的前視(Φf)和后視(Φb)干涉圖，最后從前視、后視干涉圖的相位差(ΦMAI)中提取方位向形變信息(DAZ)。其表達(dá)式為

ΦMAI=Φf-Φb

(2)

(3)

式中：l為真實(shí)天線長(zhǎng)度；n為歸一化斜視指數(shù)。

近年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)多孔徑SAR干涉法不斷改進(jìn)，以減小前視、后視干涉圖基線差及電離層誤差等對(duì)觀測(cè)精度的影響[39]。由于前視、后視干涉圖的有效視數(shù)和相干性都低于常規(guī)干涉圖，所以多孔徑SAR干涉法觀測(cè)的精度低于D-InSAR技術(shù)，為亞分米級(jí)[56]。同時(shí)，多孔徑SAR干涉法利用了相位信息，在相干性較好的區(qū)域，監(jiān)測(cè)精度優(yōu)于像素偏移量估計(jì)法。

對(duì)于近年來(lái)廣泛應(yīng)用的Sentinel-1 SAR衛(wèi)星，由于其采用的大幅寬TOPS成像模式使得方位向分辨率降低至20 m，所以常規(guī)多孔徑SAR干涉法不再適用。對(duì)此，Grandin 等提出了脈沖重疊-多孔徑干涉技術(shù)，成功獲取了2015年智利Illapel 8.3級(jí)地震的三維形變場(chǎng)[56]。脈沖重疊-多孔徑干涉技術(shù)利用了方位向上相鄰脈沖間在重疊區(qū)成像時(shí)的視角差來(lái)計(jì)算方位向位移，精度可以達(dá)到多孔徑SAR干涉法的10倍[57]，但由于重疊區(qū)范圍有限，所以該方法無(wú)法獲得連續(xù)的面狀方位向位移(圖2、3)。脈沖重疊-多孔徑干涉技術(shù)也可以與像素偏移量估計(jì)法相結(jié)合以提高偏移量估計(jì)的精度[58]。此外，得益于TOPS觀測(cè)技術(shù)，Sentinel-1 SAR數(shù)據(jù)還可以實(shí)現(xiàn)距離向分頻(Range Split-bandwidth Interferometry,Range-SBI)技術(shù)，得到較高精度的距離向位移(圖2)；當(dāng)遇到形變量過(guò)大導(dǎo)致D-InSAR無(wú)法解纏的情況時(shí)，Range-SBI技術(shù)可以作為有效替代[52]。

InSAR技術(shù)獲取形變的精度方面，D-InSAR技術(shù)獲取的LOS向形變精度最高，理論上為厘米級(jí)；其次是多孔徑SAR干涉法和脈沖重疊-多孔徑干涉技術(shù)，其中多孔徑SAR干涉法獲取的方位向形變誤差至少為D-InSAR技術(shù)的2倍；像素偏移量估計(jì)法獲取的方位向和LOS向形變精度最低，約為SAR數(shù)據(jù)分辨率的1/10，但是其抵抗失相干能力更強(qiáng)。在實(shí)際應(yīng)用中需要根據(jù)研究區(qū)域特性選擇最適合的觀測(cè)技術(shù)，也可以綜合應(yīng)用多種InSAR技術(shù)實(shí)現(xiàn)優(yōu)勢(shì)互補(bǔ)[51-52,59-60]。

1.3.3 聯(lián)合InSAR和GPS數(shù)據(jù)解算三維形變

GPS 可以精確測(cè)定地表三維形變，但是受觀測(cè)條件和成本的限制，大部分地區(qū)GPS點(diǎn)位空間密度較低；若將高空間分辨率的 InSAR觀測(cè)結(jié)果與具有三維形變分量的GPS觀測(cè)結(jié)果進(jìn)行有效融合，則可以獲得高空間分辨率的三維形變場(chǎng)。關(guān)于聯(lián)合InSAR和GPS數(shù)據(jù)獲取高空間分辨率三維形變場(chǎng)，國(guó)內(nèi)外已經(jīng)開展大量的研究工作，主要分為兩類。第一類為解析法，通過(guò)構(gòu)建設(shè)計(jì)矩陣聯(lián)合InSAR和GPS求解三維形變分量，Samsonov 等利用內(nèi)插GPS和升降軌InSAR數(shù)據(jù)，結(jié)合貝葉斯估計(jì)法求解了美國(guó)南加州地區(qū)的三維形變場(chǎng)[61]。Shen等基于該方法設(shè)計(jì)了一套GPS與InSAR融合的開源軟件，并將其成果用于美國(guó)南加州地區(qū)的三維地殼形變監(jiān)測(cè)[62]。

另一類以Gudmundsson等最早開發(fā)的SISTEM系列方法為主，該方法認(rèn)為在一定范圍內(nèi)研究區(qū)域應(yīng)變是一定的，在小應(yīng)變假設(shè)下建立觀測(cè)方程，從而利用GPS和InSAR求解研究區(qū)的三維形變場(chǎng)、應(yīng)變和旋轉(zhuǎn)參數(shù)[63]。自該方法提出以來(lái)，在數(shù)據(jù)選取和反演方法上都取得了較大發(fā)展。Guglielmino等采用非內(nèi)插的GPS數(shù)據(jù)，利用加權(quán)最小二乘法求解了意大利Etna火山地區(qū)的三維形變場(chǎng)[64]。Hu 等利用方差分量估計(jì)法求解美國(guó)南加州地區(qū)三維形變[39]。Fuhrmann等利用GNSS、水準(zhǔn)測(cè)量和升降軌InSAR數(shù)據(jù)獲取了德法瑞交界處Upper Rhine地塹區(qū)的三維形變場(chǎng)[65]。Song等基于非均勻三角格網(wǎng)建立觀測(cè)方程，獲取了海原斷裂及鄰區(qū)的三維形變場(chǎng)[66]。Luo等提出一種擴(kuò)展的SIS-TEM方法(ESISTEM)，引入InSAR數(shù)據(jù)作為約束，獲取了2015年尼泊爾8.1級(jí)地震的同震三維形變場(chǎng)[67]。

除了外部形變資料，還可以結(jié)合地質(zhì)、地球物理資料等先驗(yàn)信息，以此為約束條件，計(jì)算InSAR三維同震形變場(chǎng)。孫建寶等將地表形變模型與 D-InSAR 數(shù)據(jù)相結(jié)合，獲得了2003年伊朗巴姆地震的同震三維形變場(chǎng)[68]。此外，如果已知地震引起的南北向形變量級(jí)較小時(shí)，可忽略南北向變形對(duì)InSAR觀測(cè)值的貢獻(xiàn)，只估計(jì)垂直向和東西向形變，即2.5維形變[69-71]。

圖件引自文獻(xiàn)[52]圖2 基于升軌Sentinel-1數(shù)據(jù)的2016年日本熊本地震同震形變場(chǎng)Fig.2 Coseismic Deformation Field of the 2016 Kumamoto (Japan) Earthquake Based on Ascending Sentinel-1 Data

圖件引自文獻(xiàn)[52]圖3 基于升軌Sentienl-1數(shù)據(jù)的2016 年日本熊本Kumamoto地震局部同震形變場(chǎng)Fig.3 Local Coseismic Deformation Field of the 2016 Kumamoto (Japan) Earthquake Based on Ascending Sentinel-1 Data

2 發(fā)震斷層震源參數(shù)反演模型及方法

在獲取高空間分辨率、高精度同震形變場(chǎng)的基礎(chǔ)上，反演發(fā)震斷層的震源參數(shù)(包括發(fā)震斷層幾何參數(shù)和同震位錯(cuò)分布)，為進(jìn)一步厘定發(fā)震構(gòu)造、認(rèn)識(shí)地震發(fā)生過(guò)程、揭示地震發(fā)生機(jī)理提供重要依據(jù)。本文針對(duì)利用InSAR同震形變場(chǎng)進(jìn)行震源參數(shù)反演的模型和方法進(jìn)行綜述，主要包括正演模型、反演方法、先驗(yàn)約束和多種觀測(cè)共同約束等方面。

2.1 彈性半空間斷層模型

1910年Reid根據(jù)1906年美國(guó)舊金山大地震前后的形變觀測(cè)資料，提出彈性回跳理論解釋地震發(fā)生機(jī)制，首次認(rèn)識(shí)到地震地表變形與斷層運(yùn)動(dòng)相關(guān)[72]。此后，為了建立地震地表變形與地震斷層之間的數(shù)學(xué)模型，大量學(xué)者開展了深入研究。1958年Steketee首次將晶體位錯(cuò)理論引入地震斷層領(lǐng)域[73]，該模型在均勻彈性介質(zhì)下，給出了斷層位錯(cuò)產(chǎn)生形變場(chǎng)的解析公式，此后大量學(xué)者建立了從二維到三維、從點(diǎn)源到矩形位錯(cuò)、從全空間到半空間的解析公式[74-79]。其中，日本學(xué)者Okada于1985年整理并總結(jié)了前人工作，給出了一套實(shí)用的位錯(cuò)變形計(jì)算公式[80]。目前，Okada提出的公式是應(yīng)用最廣泛的位錯(cuò)理論經(jīng)典表達(dá)式。該公式主要描述地殼中斷層面上發(fā)生位移，地表觀測(cè)點(diǎn)產(chǎn)生位移?；驹砣鐖D4所示，位于地下深處斷面發(fā)生位移，在地表任一點(diǎn)引起的位移通過(guò)假定介質(zhì)為均勻彈性半空間而求得?；诘乇碛^測(cè)，大量學(xué)者利用這一公式反演確定深部斷層的滑動(dòng)分布。

圖件引自文獻(xiàn)[80]圖4 斷層位錯(cuò)模型Fig.4 Dislocation Model of Fault

在InSAR與GPS技術(shù)發(fā)展起來(lái)之前，獲取發(fā)震斷層的靜態(tài)滑動(dòng)主要依靠地震臺(tái)站記錄的波形資料結(jié)合彈性位錯(cuò)理論進(jìn)行反演。隨著高空間分辨率InSAR同震形變場(chǎng)獲取越來(lái)越方便，眾多學(xué)者采用均勻彈性半空間位錯(cuò)模型對(duì)國(guó)內(nèi)外地震進(jìn)行了同震斷層滑動(dòng)分布反演。例如，Hamiel等利用GPS和InSAR數(shù)據(jù)對(duì)1999年土耳其Izmit 7.4級(jí)地震的斷層滑動(dòng)進(jìn)行約束[81];Toraldo Serra等對(duì)2011年新西蘭Christchurch 6.2級(jí)地震進(jìn)行斷層滑動(dòng)分布約束[1]；Wang等對(duì)美國(guó)Parkfield地震進(jìn)行滑動(dòng)反演[82]；Shen等利用不同的反演方法對(duì)2008年中國(guó)四川汶川8.0級(jí)地震的斷層滑動(dòng)分布進(jìn)行了估算[83-86]。這些研究中的斷層模型均采用了均勻彈性半空間模型。

均勻彈性半空間模型是對(duì)真實(shí)地殼介質(zhì)模型的簡(jiǎn)化，模型相對(duì)簡(jiǎn)單，運(yùn)算效率高，得到了廣泛應(yīng)用，缺點(diǎn)是沒(méi)有考慮地球介質(zhì)的分層特性。已有研究顯示，在斷層參數(shù)反演中考慮介質(zhì)分層結(jié)構(gòu)因素加以綜合分析，往往能夠更好地模擬地震同震形變場(chǎng)[87-89]。一些學(xué)者基于數(shù)值積分法[90-91]或采用收斂算法中 Lipshitz-hankel 積分的有限項(xiàng)嘗試解決該問(wèn)題[92]，但隨著層數(shù)的增加，求解波數(shù)譜的計(jì)算量呈指數(shù)增加，因此，上述兩種算法解算超過(guò)四層結(jié)構(gòu)因素是不現(xiàn)實(shí)的。Thomson等提出矩陣傳播算法，并將其應(yīng)用在分層介質(zhì)中計(jì)算傳播的地震波[93-94]，該算法具有簡(jiǎn)潔直觀的優(yōu)點(diǎn)，且理論上可應(yīng)用在層數(shù)不受限制的分層模型中，解決了計(jì)算量大的問(wèn)題，但在實(shí)際應(yīng)用中存在數(shù)值不穩(wěn)定的缺點(diǎn)[95]。Zhu等提出矩陣傳播的正交歸一法解決這一問(wèn)題，該方法簡(jiǎn)潔明了，且能夠大大提高計(jì)算效率[96-97]。

2.2 考慮三維介質(zhì)各向異性的數(shù)值模型

隨著數(shù)值計(jì)算技術(shù)的快速革新，構(gòu)建接近真實(shí)地球的介質(zhì)屬性，對(duì)于獲取精細(xì)化斷層滑動(dòng)分布逐漸成為可能。Masterlark等以1995年墨西哥Colima-Jalisc地震為例，采用Abaquas軟件對(duì)比了均勻模型與各向異性模型對(duì)同震滑動(dòng)的影響，發(fā)現(xiàn)均勻模型會(huì)高估斷層的滑動(dòng)分布量級(jí)[98]。Wald等利用三維有限元技術(shù)，結(jié)合GPS、水準(zhǔn)和地震波資料聯(lián)合反演了北橋地震的斷層滑動(dòng)分布，經(jīng)過(guò)測(cè)試發(fā)現(xiàn)各向異性的影響不可忽略[99]。Sato等利用三維有限元技術(shù)分析了彈性參數(shù)的三維各向異性對(duì)日本東北地震同震地表位移的影響，發(fā)現(xiàn)在相同斷層位移的情況下，彈性參數(shù)三維各向異性在地表位移的影響最大可達(dá)均勻彈性模型位移的40%[100]。此后，利用三維有限元技術(shù)對(duì)斷層同震滑動(dòng)進(jìn)行約束的方法逐漸被運(yùn)用到全球其他地震，如2005年印度尼西亞蘇門答臘Nias-Simeulue 8.7級(jí)地震[101]、2003～2004年中國(guó)青海德令哈震群[102]、2011年日本東北9.0級(jí)地震[103-104]、2010年智利Maule 8.8級(jí)地震[105]、2015年智利Illapel 8.3級(jí)地震[106]、2015年尼泊爾8.1級(jí)地震[107]、2016意大利Amatrice 6.2級(jí)地震[108-109]、2018年墨西哥Pinotepa 7.2級(jí)地震[110]。除了考慮三維各向異性的彈性模型進(jìn)行斷層滑動(dòng)反演外，也有學(xué)者采用彈塑性材料進(jìn)行同震形變場(chǎng)的模擬研究。Nevitt等利用三維彈塑性有限元技術(shù)對(duì)2014年美國(guó)加利福利亞South Napa地震進(jìn)行了觀測(cè)資料擬合(移動(dòng)雷達(dá)掃描和地震波數(shù)據(jù))，發(fā)現(xiàn)彈塑性模型能夠更好地?cái)M合近場(chǎng)的地表破裂，而彈性模型因?yàn)闊o(wú)法產(chǎn)生非連續(xù)變形而無(wú)法擬合地表破裂分布[111]。除了考慮更為精細(xì)的介質(zhì)屬性分布，自動(dòng)化、流程化獲取震源參數(shù)也是一個(gè)重要的發(fā)展方向。圖5為Tung等利用接近真實(shí)地球介質(zhì)的三維有限元模型快速反演了2018年墨西哥Pinotepa 7.2級(jí)地震的滑動(dòng)分布，該方法在預(yù)先建立好模型的基礎(chǔ)上，能夠自動(dòng)化、流程化地給出滑動(dòng)分布結(jié)果[110]。

與傳統(tǒng)模型相比，各向異性地殼介質(zhì)模型能夠更好地解釋觀測(cè)資料，分析發(fā)震構(gòu)造，且已經(jīng)在全球廣泛使用。因此，采用更為精細(xì)的地殼結(jié)構(gòu)和接近實(shí)際的斷層模型，是未來(lái)基于同震形變場(chǎng)快速反演震源參數(shù)的發(fā)展方向。

2.3 反演方法與先驗(yàn)約束

地球物理反演問(wèn)題一般存在多解性，其中的一個(gè)原因是缺乏足夠的約束[112]。因此，為了確定最優(yōu)解，需要對(duì)反演參數(shù)施加一定的先驗(yàn)約束，以消除系數(shù)矩陣的奇異性。在實(shí)際應(yīng)用中，反演方法和約束方法的使用對(duì)斷層參數(shù)反演結(jié)果產(chǎn)生一定影響。

地震震源參數(shù)反演包含斷層幾何參數(shù)反演和斷層滑動(dòng)分布反演[113]。在斷層幾何參數(shù)未知的情況下，斷層幾何參數(shù)反演是斷層滑動(dòng)分布反演的前提。斷層幾何參數(shù)反演一般采用非線性反演算法，包括最速下降方法[97]、蒙特卡洛貝葉斯法[114]、網(wǎng)格搜索法等；而對(duì)于線性反演，主要采用最小二乘法、約束最小二乘法。也有學(xué)者將非線性反演和線性反演共同進(jìn)行，稱作并行貝葉斯反演方法[115]，或者將斷層幾何參數(shù)進(jìn)行參數(shù)化后與其他線性參數(shù)一起反演[86]。

zPE為Pinotepa地震震源深度；zSlab1.0為板片深度;z為深度；δPE為Pinotepa地震傾角；δSlab1.0為Slab1.0模型估計(jì)的走向;δ為傾角；φPE為Pinotepa地震走向；φ為走向；Middle America Trench為中美海溝；stress-free land/seafloor surface為自由邊界/海底表面；Mexico為墨西哥；Gulf of Mexico為墨西哥灣；Pacific Ocean為太平洋；Mexico City為墨西哥城；Guatemala為危地馬拉；fine-mesh region為精細(xì)網(wǎng)格；pinned farfield boundary為遠(yuǎn)場(chǎng)邊界固定；length為長(zhǎng)度；width為寬度；strike為走向；dip為傾角；圖件引自文獻(xiàn)[110]圖5 2018年墨西哥Pinotepa7.2級(jí)地震三維有限元模型Fig.5 3D Finite Element Model of the 2018 M 7.2 Pinotepa (Mexico) Earthquake

斷層滑動(dòng)分布反演的約束方法目前常用的有平滑約束、最大位移約束、地震矩約束、斷層性質(zhì)約束、考慮斷層邊界的約束以及余震約束等。其中，平滑約束是最基本、最常用的約束，在大部分反演工作中均予以考慮，在實(shí)際應(yīng)用中，有學(xué)者發(fā)展了自適應(yīng)平滑約束方法[116]。采用平滑約束最關(guān)鍵的是確定平滑因子，比較常用的方法是L曲線法；也有學(xué)者采用方差分量估計(jì)法，即將平滑因子和不同數(shù)據(jù)權(quán)重比一起進(jìn)行優(yōu)化，以確定平滑因子[117]。通常，在利用InSAR資料開展斷層滑動(dòng)分布反演之前，地震學(xué)給出的震源機(jī)制和斷層滑動(dòng)一般可作為先驗(yàn)信息對(duì)位錯(cuò)模型進(jìn)行一定約束，比如地震矩約束[97]；同時(shí)，也可以利用地震現(xiàn)場(chǎng)調(diào)查的最大破裂位移，對(duì)斷層滑動(dòng)進(jìn)行約束[118]；另外，重新定位后的余震分布和震源機(jī)制解也是斷層滑動(dòng)分布反演的重要約束[10,119]。

2.4 多種數(shù)據(jù)聯(lián)合反演

本文主要綜述利用InSAR觀測(cè)資料進(jìn)行震源參數(shù)反演，如果區(qū)域資料比較豐富，為了獲得更為可靠的結(jié)果，還會(huì)結(jié)合GPS、水準(zhǔn)、地震波等多種資料進(jìn)行聯(lián)合反演。對(duì)此，國(guó)內(nèi)外已開展了大量相關(guān)研究，如2004年美國(guó)Parkfield地震[120]、2005年巴基斯坦地震[121]、1997年中國(guó)西藏瑪尼7.5級(jí)地震[122]、2008年中國(guó)四川汶川8.0級(jí)地震[123-124]、2011年日本東北9.0級(jí)地震[103]、2015年尼泊爾8.1級(jí)地震[107]、2010年中國(guó)青海玉樹7.1級(jí)地震[53]，2015年智利Illapel 8.3級(jí)地震等[106]。從上述研究可以看出，相對(duì)于地震波數(shù)據(jù)來(lái)講，InSAR觀測(cè)資料因?yàn)槟軌蛱峁└呖臻g分辨率的近場(chǎng)約束，所以其對(duì)斷層淺層的約束能力更強(qiáng)。因此，InSAR技術(shù)在研究近場(chǎng)地震形變分布和破裂過(guò)程中具有不可替代的作用。但是，對(duì)于震級(jí)較大的地震來(lái)講，如果震中的同震位移大于InSAR監(jiān)測(cè)形變的最大閾值，則會(huì)造成震中區(qū)域的失相干，從而無(wú)法完全捕捉近斷層的同震形變。而對(duì)于同震位移場(chǎng)來(lái)講，無(wú)法完全捕捉發(fā)震斷層淺部的滑動(dòng)，則會(huì)導(dǎo)致無(wú)法準(zhǔn)確預(yù)測(cè)余震的破裂方向，因此，采用多種數(shù)據(jù)源(光學(xué)影像等)獲取完整的近場(chǎng)同震形變場(chǎng)，是未來(lái)精細(xì)化研究地震破裂過(guò)程的一個(gè)重要方向。

值得注意的是，在采用多種觀測(cè)資料進(jìn)行聯(lián)合反演時(shí)，各種觀測(cè)資料權(quán)重的確定尤為重要。由于不同觀測(cè)資料的誤差評(píng)價(jià)標(biāo)準(zhǔn)不同，數(shù)據(jù)量和空間分布形態(tài)存在較大差異，所以需要對(duì)不同數(shù)據(jù)之間的權(quán)重進(jìn)行確定，以平衡不同類型數(shù)據(jù)對(duì)反演結(jié)果的影響。目前來(lái)說(shuō)，一般采用層次分析法進(jìn)行數(shù)據(jù)定權(quán)。層次分析法具有精度高、方便快捷等優(yōu)點(diǎn)，其通過(guò)對(duì)不同類型數(shù)據(jù)權(quán)重因子兩兩比較確定最終權(quán)重，可以避免主觀因素對(duì)整體的影響[125]。當(dāng)觀測(cè)資料種類較多時(shí)，兩兩相比的方法效率較低。為了解決該問(wèn)題，許才軍等利用方差分量估計(jì)法，將平滑因子與各種數(shù)據(jù)類型權(quán)重一起進(jìn)行迭代求解[117]。

3 結(jié) 語(yǔ)

(1)與連續(xù)GPS觀測(cè)不同，InSAR觀測(cè)的空間分辨率很高，但其時(shí)間分辨率受制于衛(wèi)星的重訪周期，因此，在監(jiān)測(cè)同震形變時(shí)，往往把強(qiáng)余震引起的變形場(chǎng)疊加到一起，從而無(wú)法準(zhǔn)確識(shí)別主震和余震產(chǎn)生的形變場(chǎng)。要解決這一問(wèn)題，除了提高InSAR觀測(cè)的時(shí)間分辨率外，也可以采用更為復(fù)雜的模型或更為多元的觀測(cè)進(jìn)行分離。具體來(lái)講，在模型方面可以將同震形變模型與震后余滑模型作為參數(shù)一起進(jìn)行疊加，然后進(jìn)行反演計(jì)算。然而這種方法一是忽略了震后松弛和孔隙彈性模型的影響，二是對(duì)大量的余震無(wú)法很好地參數(shù)化，在實(shí)際構(gòu)建模型中尚有許多尚未解決的問(wèn)題。在多種資料聯(lián)合反演方面，由于不同學(xué)科的跨度較大，處理形變和波形資料都需要專業(yè)人員完成，所以如何加強(qiáng)不同學(xué)科之間的交叉實(shí)現(xiàn)多種資料的聯(lián)合處理是未來(lái)發(fā)展的方向。

(2)目前獲取InSAR同震形變場(chǎng)的方法比較成熟，但是無(wú)人工干預(yù)、自動(dòng)化下載、數(shù)據(jù)處理和反演計(jì)算的軟件或系統(tǒng)還沒(méi)有報(bào)道，這是未來(lái)發(fā)展的趨勢(shì)。目前，在InSAR同震形變場(chǎng)自動(dòng)獲取方面的技術(shù)已經(jīng)成熟，但如何根據(jù)發(fā)震位置的構(gòu)造背景設(shè)置斷層初始參數(shù)進(jìn)行斷層參數(shù)反演還需要人工參與。目前，每個(gè)地震的發(fā)震斷層需要結(jié)合震源機(jī)制解與重新定位后的余震分布確定。因此，構(gòu)建盡可能完備的活斷層庫(kù)，或根據(jù)重新定位后的小震分布自動(dòng)確定斷層位置，是目前制約利用InSAR技術(shù)自動(dòng)化確定震源參數(shù)的關(guān)鍵，也是迫切需要解決的問(wèn)題之一。

(3)在震源參數(shù)反演時(shí)，目前國(guó)內(nèi)大部分反演模型采用均勻彈性半空間模型，但是建立區(qū)域接近真實(shí)的地殼介質(zhì)模型也是未來(lái)發(fā)展的方向。具體來(lái)講，為了獲得盡可能真實(shí)的震源參數(shù)，構(gòu)建中國(guó)大陸重點(diǎn)危險(xiǎn)區(qū)精細(xì)化的三維地殼速度結(jié)構(gòu)，并建立三維有限元模型和斷層模型，生成格林函數(shù)庫(kù)，是快速獲得震后斷層滑動(dòng)分布的重要手段之一。目前，中國(guó)大陸重點(diǎn)危險(xiǎn)區(qū)的深部結(jié)構(gòu)并沒(méi)有一個(gè)完備的數(shù)據(jù)庫(kù)，這也是未來(lái)需要完善的方向。

(4)科學(xué)研究層面，同震滑動(dòng)分布研究對(duì)于了解地震破裂過(guò)程、總結(jié)地震發(fā)生規(guī)律、分析區(qū)域構(gòu)造特征等具有重要意義。尤其是加強(qiáng)對(duì)青藏高原陸內(nèi)碰撞帶等典型構(gòu)造區(qū)大地震同震形變特征的分析，能夠提高對(duì)高原隆升及區(qū)域構(gòu)造演化過(guò)程的認(rèn)識(shí)。例如，可以通過(guò)分析InSAR技術(shù)獲取的精細(xì)斷層滑動(dòng)矢量方向，確定發(fā)震斷層的運(yùn)動(dòng)性質(zhì)，為研究發(fā)震斷層的構(gòu)造特征提供約束；可以結(jié)合震前形變場(chǎng)獲取的斷層閉鎖程度、同震滑動(dòng)分布和震后形變時(shí)間序列觀測(cè)資料確定發(fā)震斷層的摩擦參數(shù)，進(jìn)而為構(gòu)建基于速度-狀態(tài)摩擦本構(gòu)的地震預(yù)測(cè)模型提供約束；也可以利用震后高密度InSAR形變時(shí)間序列約束青藏高原不同發(fā)震斷層區(qū)域深部的流變結(jié)構(gòu)，從而將長(zhǎng)期構(gòu)造形變與震后短期形變聯(lián)系起來(lái)，為認(rèn)識(shí)青藏高原隆升和區(qū)域構(gòu)造演化提供支持依據(jù)。

在長(zhǎng)安大學(xué)七十周年華誕之際，謹(jǐn)以此文祝福母校永遠(yuǎn)輝煌，再譜絢麗華章！七十載春華秋實(shí)，風(fēng)雨兼程；七十載開拓進(jìn)取，滄桑砥礪！清晰地記得，2000年9月6日的清晨，經(jīng)過(guò)兩天汽車、火車的輾轉(zhuǎn)，我從東北一個(gè)小鄉(xiāng)村到達(dá)了古都西安，開啟了我七年的長(zhǎng)安大學(xué)求學(xué)生涯。七年的求學(xué)經(jīng)歷，讓我從一名稚氣未脫的少年蛻變成一名學(xué)有所成的青年。工作十幾年來(lái)，我始終牢記自己是一名長(zhǎng)大人，始終將母校校訓(xùn)“弘毅、明德、篤學(xué)、創(chuàng)新”融入工作和生活，深入靈魂！母校的生活讓人終生難忘，食堂懸掛的“今天我以長(zhǎng)安大學(xué)為榮，明天長(zhǎng)安大學(xué)以我為榮”標(biāo)語(yǔ)，已成為我畢生奮斗的動(dòng)力源泉，深入骨髓！張伯聲先生的塑像挺立在圖書館前，為我指明了人生努力的方向！忘不了階梯教室里，恩師授課的神采飛揚(yáng)；忘不了測(cè)圖實(shí)習(xí)時(shí)，恩師的辛勤汗水；忘不了科學(xué)研究中，恩師的指點(diǎn)迷津；忘不了日常生活里，恩師的諄諄教誨！母校給予的一切，讓我受益終身！