大地電磁與接收函數(shù)結(jié)構(gòu)約束聯(lián)合反演

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

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

0 引言

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

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

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

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

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

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

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

(1)

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

(2)

(3)

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

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

(5)

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

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

(6)

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

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

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

(7)

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

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

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

2 理論模型測試

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

(8)

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

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

圖1 模型I含低阻層時接收函數(shù)與大地電磁聯(lián)合反演結(jié)果(a) 單獨反演； (b) 聯(lián)合反演，聯(lián)合反演參數(shù)β=3.0,=1.5 .中間子圖為每個部分的數(shù)據(jù)擬合結(jié)果， 上圖為接收函數(shù)的數(shù)據(jù)擬合，下圖為大地電磁視電阻率的數(shù)據(jù)擬合；右圖為兩個模型的歸一化后的空間梯度平方分布.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) 單獨反演； (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) 含低阻異常的單獨反演與聯(lián)合反演結(jié)果； (c)(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) 單獨反演； (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ū)臺站分布圖，紅色圓圈為大地電磁臺站，藍色三角為地震臺站. (b)聯(lián)合反演目標函數(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 地震臺站A504和大地電磁臺站1141n5TM模式聯(lián)合反演結(jié)果，右圖為歸一化后的空間梯度平方的分布圖(a) 單獨反演； (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 地震臺站A504 和大地電磁臺站1141n5數(shù)據(jù)擬合結(jié)果 在第一排中，灰線為疊加之前的接收函數(shù).(a) 單獨反演； (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.

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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