甘露, 吳慶舉, 黃清華, 張慧茜， 唐榮江

1 中國(guó)地震局地球物理研究所， 北京 100081 2 北京大學(xué)地球與空間科學(xué)學(xué)院， 北京 100871

0 引言

隨著地球物理技術(shù)的發(fā)展，對(duì)地下目標(biāo)詳細(xì)的、多屬性參數(shù)探測(cè)的需求日益增加，這在很大程度上促進(jìn)了自動(dòng)化綜合地球物理新方法的發(fā)展，特別是聯(lián)合反演方法的更新(Julià et al., 2000; Gallardo and Meju, 2003; Moorkamp et al., 2007; Moorkamp, 2007; Haber and Holtzman Gazit, 2013).本文嘗試了對(duì)不同地球物理參數(shù)敏感的大地電磁和接收函數(shù)數(shù)據(jù)的聯(lián)合反演.大地電磁數(shù)據(jù)對(duì)地下地質(zhì)體的電性結(jié)構(gòu)敏感，而接收函數(shù)對(duì)速度的突變界面敏感.雖然在地殼和地幔中，這兩種參數(shù)的變化并不是嚴(yán)格空間相關(guān)的，通常來說，對(duì)于大地電磁，電性結(jié)構(gòu)的異常通常存在于巖體的次要成分，比如石墨、硫化物，或者滲入巖石內(nèi)部的流體，這些通常對(duì)地震方法來說探測(cè)比較困難；但是，同樣也存在許多物理參數(shù)將會(huì)同時(shí)影響電導(dǎo)率和速度，例如溫度和孔隙度等.此外，在一些主要的巖性邊界，電導(dǎo)率和速度這兩種參數(shù)都會(huì)發(fā)生改變 (Marquis et al., 1995; Lahti et al., 2005; Tournerie and Chouteau, 2005)，即使影響電導(dǎo)率的是石墨，也可以認(rèn)為是構(gòu)造活動(dòng)決定了石墨體的存在和規(guī)模，該構(gòu)造活動(dòng)也同樣決定了地震波速度(Carcione et al., 2007).此外，兩個(gè)主要的邊界即莫霍面和巖石圈邊界也明顯與電導(dǎo)率的變化相關(guān)(Jones, 1999; Jones and Ferguson, 2001; Gatzemeier and Moorkamp, 2005).因此，可以預(yù)料到，接收函數(shù)與大地電磁聯(lián)合反演得到的兩種模型，存在一定的結(jié)構(gòu)相似性.

大地電磁和接收函數(shù)聯(lián)合反演的目標(biāo)是重建地殼乃至上地幔結(jié)構(gòu)，一個(gè)主要的挑戰(zhàn)是：在實(shí)際情況中，高速(低速)異常和高阻(低阻)異常不一定完全對(duì)應(yīng) (Jones et al., 2009).在地殼內(nèi)部，對(duì)于速度結(jié)構(gòu)來說，隨著深度的增加，由于巖石壓力增加，導(dǎo)致巖石更加致密，密度和波速也增大.即使在地幔熱物質(zhì)存在的區(qū)域，波速會(huì)下降但不會(huì)影響速度的整體變化趨勢(shì).電阻率的大小主要受溫度、壓力以及含水量的綜合影響，電阻率與壓力成正相關(guān)，與溫度和含水量成負(fù)相關(guān).隨著深度的增加，溫度和壓力都會(huì)增加，含水量也會(huì)發(fā)生變化，電阻率受到三者綜合影響總體上呈下降趨勢(shì)，局部出現(xiàn)反彈.地震波速度和電阻率在深部的變化規(guī)律使得它們很難通過經(jīng)驗(yàn)公式聯(lián)系起來，通常情況下經(jīng)驗(yàn)公式來源于測(cè)井?dāng)?shù)據(jù)，而測(cè)井?dāng)?shù)據(jù)通常較淺，一般不超過8 km；另一方面，如果使用共享部分屬性(例如層厚等)信息實(shí)現(xiàn)聯(lián)合(Moorkamp, 2007)，需要設(shè)定較少的層數(shù)實(shí)現(xiàn)約束，其實(shí)質(zhì)是將約束條件下的欠定問題變?yōu)檫m定或者超定問題，然后尋找共同界面下的最優(yōu)解.目前主流的大地電磁與接收函數(shù)聯(lián)合反演是通過共享界面來達(dá)到聯(lián)合反演的目的，而未加入模型耦合項(xiàng)，假設(shè)不同觀測(cè)數(shù)據(jù)對(duì)相同的異常都具有敏感性 (Moorkamp et al., 2007, 2010; Moorkamp, 2007; Zevallos et al., 2009; 彭淼等, 2012).然而這種方法限制了模型反演的靈活性.在該方法中，使用多目標(biāo)函數(shù)遺傳算法NSGA-2(Deb et al., 2002)尋找的帕累托最優(yōu)解是基于少量的模型層數(shù)，難以找到約束條件下兩個(gè)目標(biāo)函數(shù)的全局最優(yōu)解.當(dāng)層數(shù)較多時(shí)，兩個(gè)模型是解耦的，分別都能找到各自的最優(yōu)解.因此，聯(lián)合反演的一個(gè)最大挑戰(zhàn)在于找到既可以提供模型足夠的自由度去擬合所有觀測(cè)特征，又可以耦合兩種不同屬性的模型(Moorkamp, 2007).為了在維持模型一定自由度情況下，提高模型的耦合程度，本文提出新的約束算子用于速度和電阻率模型相似度約束.

接收函數(shù)與大地電磁的反演不但針對(duì)不同的物性，它們各自本身屬性也具有較大差異.接收函數(shù)屬于地震信號(hào)序列，控制方程為波動(dòng)方程，信號(hào)的強(qiáng)弱受波阻抗界面的控制；電磁場(chǎng)的控制方程在低頻段的地下介質(zhì)中滿足傳導(dǎo)方程，信號(hào)強(qiáng)弱受區(qū)域內(nèi)地球內(nèi)部電阻率的綜合影響.這些差異使得兩種反演方法具有完全不同的特性.接收函數(shù)對(duì)于地層分界面敏感，對(duì)絕對(duì)速度的約束較弱，這導(dǎo)致它雖然具有較高的分辨率，但是嚴(yán)重受到初始模型的影響(如，Ammon et al., 1990; Kind et al., 1995; Tomlinson et al., 2006)；大地電磁對(duì)于電性界面沒有很好的分辨能力，但是對(duì)于絕對(duì)電阻率屬性約束較好，且受初始模型的影響相對(duì)接收函數(shù)較弱.綜上可以推測(cè)，大地電磁和接收函數(shù)的聯(lián)合反演具有一定的互補(bǔ)性.本文將通過理論模型和實(shí)際資料的測(cè)試，探索聯(lián)合反演能否降低二者單獨(dú)反演中的非唯一性.

本文采用擬牛頓法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，尋找最優(yōu)解.擬牛頓法是梯度類反演中最具代表性的，也是大型非線性問題中最優(yōu)秀的方法之一.擬牛頓法的本質(zhì)是牛頓法的近似，它通過近似迭代求解牛頓法中計(jì)算量最大的Hessian矩陣，保證了模型的二階(曲率)下降.與線性化的方法和全局搜索類算法相比，梯度類的算法通?；赥ikhonov 正則化方法，可以很方便地建立不同目標(biāo)函數(shù)以及約束項(xiàng)求和形式.為了設(shè)定足夠的反演參數(shù)，同時(shí)促進(jìn)兩個(gè)模型找到最優(yōu)耦合形式，一個(gè)合理的方法是建立起電阻率模型和速度模型的聯(lián)合算子項(xiàng)，作為目標(biāo)函數(shù)需要優(yōu)化的多項(xiàng)式之一，這種方案是目前結(jié)構(gòu)相似算法普遍采用的方案，例如Gallardo和Meju(2003)提出的交叉梯度聯(lián)合反演.受到交叉梯度方法的啟發(fā)，本文將設(shè)計(jì)速度模型和電阻率模型的一維梯度耦合算子，使用擬牛頓法同時(shí)優(yōu)化大地電磁、接收函數(shù)以及聯(lián)合算子項(xiàng)，實(shí)現(xiàn)大地電磁與接收函數(shù)的聯(lián)合反演.

1 聯(lián)合反演方法

接收函數(shù)和大地電磁的聯(lián)合反演通過最優(yōu)化共同的目標(biāo)函數(shù)來實(shí)現(xiàn)：

φ(mse,mem)=φse(mse)+φem(mem)+φjoint(mse,mem),

(1)

其中，φse為接收函數(shù)的目標(biāo)函數(shù)，φem為大地電磁的目標(biāo)函數(shù)，φjoint為電阻率和地震模型的聯(lián)合算子，具體展開形式如下：

(2)

(3)

此外，Vse和Vem分別對(duì)應(yīng)兩種數(shù)據(jù)的協(xié)方差矩陣的對(duì)角線元素(Egbert and Kelbert, 2012)，以Vse為例，Vem的形式與Vse一致，其表達(dá)式為：

Vse=diag{E[(X-E[X])(X-E[X])T]}

(5)

這里X為由觀測(cè)數(shù)據(jù)組成的向量.在實(shí)測(cè)數(shù)據(jù)中，由于每個(gè)元素都有大量的采樣，因此可以對(duì)每個(gè)元素求取方差.該矩陣可以讓方差較大的觀測(cè)數(shù)據(jù)在參與目標(biāo)函數(shù)的計(jì)算中占有較小的比重，以保證當(dāng)觀測(cè)數(shù)據(jù)質(zhì)量較差時(shí)反演能夠穩(wěn)定進(jìn)行.在合成數(shù)據(jù)的測(cè)試中，由于數(shù)據(jù)是沒有誤差的，Vse和Vem取單位矩陣.

在反演中，速度模型和電阻率模型共享層厚和權(quán)重，假設(shè)反演層數(shù)為k，目標(biāo)函數(shù)的梯度可以寫作：

(6)

由于假設(shè)模型為一維層狀介質(zhì)，這對(duì)于接收函數(shù)和大地電磁的反演來說計(jì)算量較小，因此，本文總是同時(shí)反演三項(xiàng)目標(biāo)函數(shù)，速度模型和電阻率模型在反演中同步更新.

另一個(gè)值得注意的問題是，在聯(lián)合反演中，目標(biāo)函數(shù)中的每一項(xiàng)應(yīng)保持在同一個(gè)數(shù)量級(jí)上，以防止目標(biāo)函數(shù)僅以擬合較大項(xiàng)為目的而忽略了較小項(xiàng)的擬合.接收函數(shù)的值一般位于[-0.3,0.8]之間；對(duì)于大地電磁，假設(shè)視電阻率的取值范圍為[1,100000]Ωm，那么對(duì)其取對(duì)數(shù)后的取值范圍為[0,5].這個(gè)范圍大于接收函數(shù)的觀測(cè)值，因此這里用二者觀測(cè)數(shù)據(jù)的最大值對(duì)各自的觀測(cè)數(shù)據(jù)進(jìn)行歸一化，同時(shí)消除采樣點(diǎn)數(shù)的影響.此外，協(xié)方差矩陣也可能造成目標(biāo)函數(shù)項(xiàng)之間的量級(jí)差，特別是值較小的方差可能造成目標(biāo)函數(shù)的值突然增大，為了減小方差影響，對(duì)原數(shù)據(jù)加權(quán)矩陣作以下處理:

V′=γ2N(V+E),

(7)

其中γ為觀測(cè)數(shù)據(jù)的最大值；E為單位矩陣，N為各自的采樣點(diǎn)數(shù).γ2N消除了接收函數(shù)和視電阻率數(shù)據(jù)量級(jí)的差異，以及采樣點(diǎn)數(shù)的差異，E消除了小方差的影響.

從公式(4)可以看出，φjoint通過極小化電阻率模型和速度模型的縱向梯度絕對(duì)值之間的差異來實(shí)現(xiàn)二者在空間上的耦合.其中，方括號(hào)里的平方是為了同化梯度的方向，即不論模型逐漸增大還是逐漸減小，只要梯度的絕對(duì)值一致，都認(rèn)為它們具有相似性.這更加符合地殼尺度的反演，因?yàn)樵诘貧?nèi)部速度隨著深度的增加總體上逐漸上升，而電阻率總體呈下降趨勢(shì).方括號(hào)外的平方是為了保持目標(biāo)函數(shù)為凸函數(shù)，使其在反演中總能尋找到極小值.在反演中可以設(shè)計(jì)足夠?qū)訑?shù)的模型，使得算法找到擬合觀測(cè)數(shù)據(jù)的模型參數(shù)；同時(shí)聯(lián)合算子保證了兩種模型的相似性，即既能夠維持足夠的模型自由度，又同時(shí)耦合了兩種模型.

本文采用擬牛頓法來極小化目標(biāo)函數(shù)，擬牛頓法的基本思想是根據(jù)在每次迭代中儲(chǔ)存前m次的模型以及梯度修正量，以盡可能小的計(jì)算代價(jià)構(gòu)建近似Hessian矩陣，進(jìn)而求得近似的牛頓下降方向，保證模型的二階下降屬性(Byrd et al.,1994,1995).除了計(jì)算下降方向以外，還需要一個(gè)合適的步長(zhǎng)以保證模型的充分下降.通?；赪olf準(zhǔn)則，采用二次插值算法，回溯迭代來搜索一個(gè)合適的步長(zhǎng).本文對(duì)有限內(nèi)存擬牛頓法(L-BFGS)開源代碼(Byrd et al., 1994, 1995; Zhu et al., 1997; https:∥github.com/stephenbeckr/L.BFGS.B.C)進(jìn)行修改后完成接收函數(shù)與大地電磁聯(lián)合反演.梯度類反演方法總是從一個(gè)初始模型開始，不斷迭代更新以逐漸逼近真實(shí)解.接收函數(shù)由于對(duì)絕對(duì)速度不敏感，導(dǎo)致反演結(jié)果高度依賴初始模型(Ammon et al., 1990; Kind et al., 1995; Tomlinson et al., 2006)，且大地電磁反演結(jié)果也一定程度上受到初始模型的影響，因此選擇合適的初始模型對(duì)反演具有重要意義.對(duì)于電阻率模型，通常采用均勻半空間為初始模型，以適用于地殼尺度的大地電磁反演(Zhang et al., 2016；葉濤等，2021)；對(duì)于速度模型，初始模型修改自全球速度模型ak135(Kennett et al., 1995).

2 理論模型測(cè)試

為了驗(yàn)證算法的有效性，本節(jié)設(shè)計(jì)了耦合的理論速度模型和電阻率模型進(jìn)行測(cè)試，電阻率模型和速度模型總是在異常處耦合，但是異常的正負(fù)不完全對(duì)應(yīng)，這更加符合真實(shí)地殼結(jié)構(gòu).每個(gè)模型的反演層數(shù)設(shè)定為80層，反演深度70 km.接收函數(shù)正演時(shí)窗為-5～46.2 s，采樣間隔0.1 s，共512個(gè)采樣點(diǎn)數(shù)，以方便進(jìn)行快速傅里葉變換.反演中僅需要反演部分時(shí)窗，即-5～25 s，采樣點(diǎn)數(shù)為300.接收函數(shù)的射線參數(shù)設(shè)定為0.055，高斯濾波系數(shù)設(shè)定為2.5.大地電磁數(shù)據(jù)的頻率范圍為 0.001～100 Hz，以對(duì)數(shù)間隔共設(shè)定40個(gè)采樣點(diǎn).權(quán)重參數(shù)α1=1.0，α2=1.0，正則化因子λ1=λ2=0.1.目標(biāo)函數(shù)的終止受控于梯度的最大分量和目標(biāo)函數(shù)的變化率，根據(jù)Byrd等(1995)，當(dāng)滿足以下條件之一時(shí)，目標(biāo)函數(shù)終止：

(8)

這里fk為第k次迭代的目標(biāo)函數(shù)，g[i]為梯度的第i個(gè)分量，l為模型的層數(shù)；pgtol為人為設(shè)定參數(shù)，在本文中設(shè)定為10-5；epsmch為機(jī)器精度(machine precision), 64位機(jī)通常為2-53，factr設(shè)定為1010.可以看出，第一個(gè)條件的物理意義為梯度足夠小，第二個(gè)為目標(biāo)函數(shù)無法下降或者反彈.

圖1設(shè)計(jì)了一組耦合的低速-低阻模型，可以看出電阻率模型和速度模型在形態(tài)上具有較大差異，然而在某些界面處，特別是低速-低阻處模型結(jié)構(gòu)是耦合的，這在地殼結(jié)構(gòu)中是非常常見的.由于25 km處低阻層較薄，單獨(dú)的大地電磁反演無法有效約束出低阻層的具體位置，而聯(lián)合反演中25 km處低阻異常明顯，大地電磁在深部的分辨率明顯提高.從空間梯度分布圖可以看出，電阻率模型梯度的平方隨深度變化，其極值對(duì)應(yīng)著梯度變化的極大值，也就是模型劇烈變化的區(qū)域.與單獨(dú)反演相比，聯(lián)合反演后速度模型和電阻率模型的空間梯度更加吻合，這說明在反演中聯(lián)合反演算子使得兩種不同屬性的模型朝著相似的方向變化.

圖1 模型I含低阻層時(shí)接收函數(shù)與大地電磁聯(lián)合反演結(jié)果(a) 單獨(dú)反演； (b) 聯(lián)合反演，聯(lián)合反演參數(shù)β=3.0,=1.5 .中間子圖為每個(gè)部分的數(shù)據(jù)擬合結(jié)果， 上圖為接收函數(shù)的數(shù)據(jù)擬合，下圖為大地電磁視電阻率的數(shù)據(jù)擬合；右圖為兩個(gè)模型的歸一化后的空間梯度平方分布.Fig.1 Joint inversion of receiver function and magnetotelluric for model I with low resistivity layer(a) Separate inversion; (b) Joint inversion with β=3.0,=1.5. The middle sub-figure is the data fitting result of each part, the upperFigure is the data fitting of the receiver function, the lowerFigure is the data fitting of the magnetotelluric apparent resistivity; the rightFigure is the distribution of normalized spatial gradient square from the two models.

圖2 模型I不含低阻層的接收函數(shù)與大地電磁聯(lián)合反演結(jié)果(a) 單獨(dú)反演； (b) 聯(lián)合反演，聯(lián)合反演參數(shù)β=3.0,=1.0.各子圖的解釋與圖1相同.Fig.2 Joint inversion of receiver function and magnetotelluric for model I without low resistivity layer(a) Separate inversion; (b) Joint inversion with β=3.0,=1.0. The description of each subgraph is the same as that in Fig.1.

圖3 模型II接收函數(shù)與大地電磁聯(lián)合反演結(jié)果 (a)(b) 含低阻異常的單獨(dú)反演與聯(lián)合反演結(jié)果； (c)(d) 含高阻異常的單獨(dú)反演與聯(lián)合反演結(jié)果，聯(lián)合反演參數(shù)β=1.0,=1.5.Fig.3 Joint inversion of receiver function and magnetotelluric for model II (a)(b) Separate and joint inversion results with low resistivity anomalies; (c)(d) Separate and joint inversion results with high resistivity anomalies,β=1.0,=1.5.

圖4 模型III接收函數(shù)與大地電磁聯(lián)合反演結(jié)果(a) 單獨(dú)反演； (b) 聯(lián)合反演，聯(lián)合反演參數(shù)β=1.0,=1.5.各子圖的解釋與圖1相同.Fig.4 Joint inversion of the receiver function and magnetotelluric for model III(a) Separate inversion; (b) Joint inversion with β=1.0,=1.5. The description of each subgraph is the same as that in Fig.1.

圖5 (a) 工區(qū)臺(tái)站分布圖，紅色圓圈為大地電磁臺(tái)站，藍(lán)色三角為地震臺(tái)站. (b)聯(lián)合反演目標(biāo)函數(shù)隨β演化曲線， 紅色虛線表示被選中的β值，φjoint下降98%Fig.5 (a) Distribution of stations in the work area. The red circle shows magnetotelluric stations, and the blue triangle shows the seismic stations. (b) β evolution curve for joint inversion objective function terms. Red dotted line indicates the selected β, which corresponds to φjoint decreased by 98%

圖6 地震臺(tái)站A504和大地電磁臺(tái)站1141n5TM模式聯(lián)合反演結(jié)果，右圖為歸一化后的空間梯度平方的分布圖(a) 單獨(dú)反演； (b) 聯(lián)合反演，聯(lián)合反演參數(shù)β=2.5,=1.0.Fig.6 Joint inversion of seismic station A504 and magnetotelluric data of station 1141n5 for TM mode. The rightFigure shows the distribution of normalized spatial gradient square(a) Separate inversion; (b) Joint inversion with β=2.5, =1.0.

圖7 地震臺(tái)站A504 和大地電磁臺(tái)站1141n5數(shù)據(jù)擬合結(jié)果 在第一排中，灰線為疊加之前的接收函數(shù).(a) 單獨(dú)反演； (b) 聯(lián)合反演.Fig.7 Data fitting of seismic station A504 and magnetotelluric station 1141n5 In the first row, the gray line is the receiver function before superposition.(a) Separate inversion; (b) Joint inversion.

為了進(jìn)一步測(cè)試模型不耦合的情況，在圖2中，對(duì)圖1中的模型進(jìn)行修改，保留速度模型的低速異常，去掉了電阻率模型的低阻異常.可以看出，聯(lián)合反演的結(jié)果與真實(shí)模型對(duì)應(yīng)良好，沒有出現(xiàn)假異常，但是由于速度模型的約束，在對(duì)應(yīng)的地方出現(xiàn)了一定的跳躍變化，這也對(duì)應(yīng)著模型附近的空間梯度異常.這說明聯(lián)合反演對(duì)于速度模型和電阻率模型不耦合的情況有一定的辨識(shí)能力.

此外，在實(shí)際情況中，電阻率和速度不一定成正相關(guān).為了驗(yàn)證聯(lián)合反演算法是否適用于不同情況，分別設(shè)計(jì)了高阻和低阻異常進(jìn)行測(cè)試(圖3)，兩種電阻率模型都對(duì)應(yīng)著低速異常.可以看出，當(dāng)?shù)妥璁惓４嬖跁r(shí)，單獨(dú)反演可以約束出異常的位置，而聯(lián)合反演使得異常進(jìn)一步突出，且在45 km的界面處出現(xiàn)明顯跳變，聯(lián)合總體上優(yōu)于單獨(dú)反演.然而，從圖3c,d可以看出，聯(lián)合反演對(duì)于高阻異常的約束能力不如低阻異常，聯(lián)合反演沒有出現(xiàn)明顯的高阻異常，僅在異常處呈現(xiàn)一個(gè)跳變，這種跳變是速度模型約束的結(jié)果.盡管如此，聯(lián)合反演結(jié)果總體上優(yōu)于單獨(dú)反演，在單獨(dú)反演中的電阻率模型平緩，沒有任何異常；而聯(lián)合反演能更好的識(shí)別出主要電性界面的位置.同樣，在圖4中，聯(lián)合反演相比單獨(dú)反演能更好的約束出電性界面的位置，這種界面往往是速度和電阻率模型共同的界面，可以從空間梯度平方的分布圖上明顯識(shí)別出來.在圖1b,2b,4b中，聯(lián)合反演視電阻率的擬合效果均比單獨(dú)反演稍差，這是因?yàn)樵诼?lián)合反演中不僅需要極小化觀測(cè)數(shù)據(jù)與正演數(shù)據(jù)之差，同時(shí)需要極小化聯(lián)合算子，因此，聯(lián)合反演的數(shù)據(jù)擬合程度一般會(huì)略低于單獨(dú)反演.

從以上討論可以看出，在聯(lián)合反演中，速度模型對(duì)電阻率模型有較好的約束能力，而速度模型受電阻率模型影響較小，其主要原因在于大地電磁數(shù)據(jù)對(duì)深部結(jié)構(gòu)的約束能力較弱，電磁場(chǎng)的觀測(cè)資料相對(duì)于接收函數(shù)對(duì)深部局部異常變化的敏感性較低.本文聯(lián)合反演實(shí)質(zhì)是通過空間梯度上的一致性，達(dá)到兩種模型相互約束的目的.大地電磁反演中的多解性主要體現(xiàn)在局部的模型變化，例如薄層、互層、異常體等可以對(duì)應(yīng)著相同的觀測(cè)數(shù)據(jù)，而接收函數(shù)的多解性主要體現(xiàn)在由于絕對(duì)速度的偏差導(dǎo)致模型整體的平移，例如對(duì)于簡(jiǎn)單兩層地殼模型，一個(gè)相對(duì)低速莫霍面較淺的模型和另一個(gè)相對(duì)高速、但是莫霍面較深的模型對(duì)應(yīng)著同一個(gè)觀測(cè)數(shù)據(jù)(Ammon et al., 1990).兩種方法多解性問題的體現(xiàn)是不一樣的，而聯(lián)合算子使用梯度并不能對(duì)整體速度結(jié)構(gòu)進(jìn)行約束，這解釋了為什么電阻率模型對(duì)速度模型的約束不明顯.

3 實(shí)際數(shù)據(jù)應(yīng)用

本節(jié)對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行接收函數(shù)與大地電磁聯(lián)合反演，以進(jìn)一步測(cè)試算法的適用性.在華北地區(qū)選擇了一組相鄰大地電磁和地震臺(tái)站，其具體位置如圖5所示，兩種數(shù)據(jù)都具有較高的信噪比.

接收函數(shù)數(shù)據(jù)來源于中國(guó)地震局地球物理研究所布設(shè)于華北的A504臺(tái)站.首先篩選出來自該臺(tái)站的震級(jí)在5.5級(jí)以上，震中距在30°到95°之間的事件，并對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，主要包括事件提取、去傾斜、去趨勢(shì)、去均值、濾波、降采樣等；然后采用迭代時(shí)間反褶積(Ligorría and Ammon, 1999)計(jì)算接收函數(shù)，高斯濾波系數(shù)采用的是2.5，然后截取P波初至的前5 s和后25 s作為接收函數(shù)的有效時(shí)窗，得到采樣頻率為10 Hz、共300個(gè)采樣點(diǎn)數(shù)的接收函數(shù).

對(duì)于大地電磁實(shí)測(cè)數(shù)據(jù)，首先對(duì)電磁場(chǎng)各個(gè)分量的時(shí)間序列進(jìn)行傅里葉變換(Thomson, 1982)，得到頻率域的電磁場(chǎng)數(shù)據(jù)，然后使用平均交叉譜(Sims et al., 1971)估計(jì)得到阻抗張量，并使用遠(yuǎn)參考點(diǎn)法(Gamble et al., 1979)以進(jìn)一步減少由非平面波場(chǎng)或儀器噪聲等原因造成的計(jì)算偏差.

聯(lián)合算子權(quán)重β的取值對(duì)反演結(jié)果有重要影響，隨著β的不斷增大，目標(biāo)函數(shù)可能陷入局部極小，數(shù)據(jù)擬合不充分，然而過小的β值將導(dǎo)致聯(lián)合算子不能很好的約束兩個(gè)模型朝著相似的方向演化.因此，本文的調(diào)整β值的準(zhǔn)則為：一方面需要保證模型可以很好的解釋觀測(cè)數(shù)據(jù)；另一方面，希望聯(lián)合算子在反演中約束兩種不同屬性的模型朝相似方向演化.因此，這里計(jì)算了不同β情況下，不同目標(biāo)函數(shù)項(xiàng)的RMS曲線，盡可能選擇一個(gè)合理的β值使得聯(lián)合算子項(xiàng)既占有足夠的比重來對(duì)模型進(jìn)行約束，又確保大地電磁和地震數(shù)據(jù)充分下降(圖5b).最終選擇β=2.5，此時(shí)φjoint從1.08下降到0.027，下降程度98%.

圖6為接收函數(shù)和大地電磁聯(lián)合反演結(jié)果，可以看出，與單獨(dú)反演相比，聯(lián)合反演的電阻率模型在莫霍面附近出現(xiàn)明顯的跳變，與地震模型匹配良好.在單獨(dú)反演中，兩者梯度變化的極值完全不相關(guān)，說明模型完全解耦；而在聯(lián)合反演中，兩種模型的梯度匹配良好，在地表，40 km莫霍面，以及65 km處梯度同時(shí)出現(xiàn)極大值，表明速度模型和電阻率模型在該地區(qū)發(fā)生明顯變化，而其余地區(qū)兩者梯度變化平緩.由于速度模型在地表以及莫霍面的約束，電阻率模型在這兩處相對(duì)于單獨(dú)反演出現(xiàn)更加劇烈的變化；而在深部65 km處，電阻率模型變化率較大，引起這種變化主要的原因是觀測(cè)數(shù)據(jù)對(duì)最后一層電阻率相比倒數(shù)幾層更為敏感.電阻率的梯度對(duì)速度模型形成約束，聯(lián)合反演后速度模型出現(xiàn)低速異常，然而大地電磁的異常在該深度上存在不確定性，因此不能判斷速度的異常為真實(shí)異常.此外，由圖7可以看出，單獨(dú)反演和聯(lián)合反演都能較好的擬合觀測(cè)數(shù)據(jù)，且聯(lián)合反演中φjoint出現(xiàn)明顯下降，這說明聯(lián)合反演找到了耦合的電阻率模型和速度模型.

實(shí)際數(shù)據(jù)的擬牛頓聯(lián)合反演測(cè)試表明，該聯(lián)合反演算法一定程度上提升了大地電磁深部的分辨率，特別是在莫霍面附近電阻率模型被約束出明顯的異常，這對(duì)于大地電磁的解釋來說至關(guān)重要.在深部，大地電磁對(duì)于速度模型的約束相對(duì)微弱，主要原因在于大地電磁對(duì)于深部的局部異常例如薄層或突變界面等敏感度較弱，模型的局部變化對(duì)于視電阻率的影響是微弱的，這使得目標(biāo)函數(shù)更容易約束電阻率模型以實(shí)現(xiàn)聯(lián)合算子下降.

4 結(jié)論

合成數(shù)據(jù)和實(shí)際觀測(cè)資料證實(shí)了擬牛頓算法聯(lián)合反演大地電磁和接收函數(shù)的有效性，大地電磁模型分辨率得到一定提升.本文的算法不僅能夠維持模型本身足夠的自由度，又能夠最大程度地獲取耦合的速度模型和電阻率模型.

聯(lián)合反演算法非常適用于識(shí)別高阻體中的低阻夾層，這種異?？赡鼙粐鷰r的高阻異常所淹沒，很難通過單獨(dú)反演得到低阻層的具體位置，而聯(lián)合反演可以利用速度模型有效識(shí)別；對(duì)于高阻夾層，聯(lián)合反演效果不如對(duì)低阻異常約束明顯，但總體上優(yōu)于單獨(dú)反演結(jié)果；聯(lián)合反演可以有效地約束出異常體的邊界，相比于單獨(dú)反演分辨率有一定的提升.總的來說，借助地震模型的約束，聯(lián)合反演可以一定程度上提高大地電磁在深部的分辨率.

在華北地區(qū)，選取了一個(gè)地震和大地電磁相鄰臺(tái)站的數(shù)據(jù)，進(jìn)行了大地電磁與接收函數(shù)的聯(lián)合反演.反演結(jié)果表明，聯(lián)合反演得到的速度模型和電阻率模型相比于單獨(dú)反演呈現(xiàn)出更加耦合的特征，特別是對(duì)于莫霍面的約束，這有利于深部構(gòu)造的解釋，且可以借助RMS-β曲線選取最優(yōu)的β參數(shù)值，它同時(shí)也是速度模型和電阻率模型是否耦合的重要指標(biāo).此外，通過對(duì)聯(lián)合算子演化過程以及兩種模型空間梯度分布的綜合分析，可以判斷電阻率模型和速度模型的耦合程度，這有利于地球物理反演結(jié)果的綜合解釋.