大地電磁與接收函數結構約束聯合反演

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

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

0 引言

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

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

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

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

1 聯合反演方法

接收函數和大地電磁的聯合反演通過最優化共同的目標函數來實現：

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

(1)

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

(2)

(3)

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

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

(5)

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

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

(6)

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

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

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

(7)

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

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

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

2 理論模型測試

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

(8)

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

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

圖1 模型I含低阻層時接收函數與大地電磁聯合反演結果(a) 單獨反演； (b) 聯合反演，聯合反演參數β=3.0,=1.5 .中間子圖為每個部分的數據擬合結果， 上圖為接收函數的數據擬合，下圖為大地電磁視電阻率的數據擬合；右圖為兩個模型的歸一化后的空間梯度平方分布.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不含低阻層的接收函數與大地電磁聯合反演結果(a) 單獨反演； (b) 聯合反演，聯合反演參數β=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接收函數與大地電磁聯合反演結果 (a)(b) 含低阻異常的單獨反演與聯合反演結果； (c)(d) 含高阻異常的單獨反演與聯合反演結果，聯合反演參數β=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接收函數與大地電磁聯合反演結果(a) 單獨反演； (b) 聯合反演，聯合反演參數β=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) 工區臺站分布圖，紅色圓圈為大地電磁臺站，藍色三角為地震臺站. (b)聯合反演目標函數隨β演化曲線， 紅色虛線表示被選中的β值，φ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模式聯合反演結果，右圖為歸一化后的空間梯度平方的分布圖(a) 單獨反演； (b) 聯合反演，聯合反演參數β=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數據擬合結果 在第一排中，灰線為疊加之前的接收函數.(a) 單獨反演； (b) 聯合反演.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中的模型進行修改，保留速度模型的低速異常，去掉了電阻率模型的低阻異常.可以看出，聯合反演的結果與真實模型對應良好，沒有出現假異常，但是由于速度模型的約束，在對應的地方出現了一定的跳躍變化，這也對應著模型附近的空間梯度異常.這說明聯合反演對于速度模型和電阻率模型不耦合的情況有一定的辨識能力.

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

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

3 實際數據應用

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

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

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

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

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

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

4 結論

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

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

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