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

季靈運，朱良玉，劉傳金，張文婷，邱江濤，徐曉雪

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

0 引 言

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

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

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

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

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

1.2 InSAR技術(shù)獲取同震形變場的缺點

InSAR技術(shù)受到自身條件的限制，在獲取同震形變場的實際應(yīng)用中也存在一些缺點。

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

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

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

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

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

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

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

(1)

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

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

1.3.1 多視角SAR影像獲取三維同震形變場

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

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

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

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

ΦMAI=Φf-Φb

(2)

(3)

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

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

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

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

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

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

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

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

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

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

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

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

2.1 彈性半空間斷層模型

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

圖件引自文獻[80]圖4 斷層位錯模型Fig.4 Dislocation Model of Fault

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

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

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

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

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

2.3 反演方法與先驗約束

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

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

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

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

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

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

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

3 結(jié) 語

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

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

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

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

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