張志勇，謝尚平，李 曼,2，馬 鑫，李 萬

(1.東華理工大學 地球物理與測控技術(shù)學院，江西 南昌 330013；2.同濟大學 海洋與地球科學學院，上海 200092)

磁法勘探中，磁化率(κ=μr-1)變化是引起磁異常的重要原因(其中μr=μ/μ0，為無量綱參數(shù))；對磁測數(shù)據(jù)進行磁化率反演的研究工作在生產(chǎn)實踐中得到了廣泛應用(Abedi et al., 2013; Li et al., 1996; Pilkington, 1997; Portniaguine et al., 2000)。然而，大地電磁勘探一般只考慮介質(zhì)電阻率變化，并假設(shè)地下磁導率(μ)為常數(shù)，等于真空中的磁導率(μ0=4π×10-7H·m-1)。但在實際勘探環(huán)境中，巖石的磁導率隨鐵磁性礦物成分的增加而增大(Clark et al., 1991)，富含鐵磁礦物的巖石磁導率明顯高于真空中磁導率。近年來，為了提高數(shù)據(jù)處理質(zhì)量充分挖掘信息，磁化率對電磁勘探的影響得到了關(guān)注(Farquharson et al., 2003; Cao, 2005; Sasaki et al., 2010; Müller et al., 2012; Van et al., 2013)。大地電磁同時反演電阻率與磁導率的研究表明，同時反演兩個參數(shù)可以改善電阻率反演效果，同時巖石磁導率的獲取也豐富了解釋參數(shù)，提高了火成巖的成像效果(Cao et al., 2005)；磁偶源同時反演磁導率與電阻率的研究表明，對于包含磁導率異常的數(shù)據(jù)只反演電阻率會產(chǎn)生錯誤模型，同時反演兩參數(shù)可以得到可信的電阻率模型和有用的磁導率信息(Farquharson et al., 2003)；但是，磁導率對頻率變化依賴性小，單一高度的航空電磁數(shù)據(jù)同時反演電阻率與磁導率的效果較差，同時反演不同高度、不同線圈排列形式的數(shù)據(jù)可以得到更精確的結(jié)果(Sasaki et al., 2010)；海洋可控源電磁數(shù)據(jù)測量理論研究與實踐表明，布設(shè)在海底的可控源電磁裝置可有效進行淺海海底沉積物的電阻率與磁導率成像(Müller et al., 2012)；高頻電磁測深(100 kHz～100 MHz)磁場垂直分量對介質(zhì)電導率、磁導率、介電常數(shù)的靈敏度研究表明，該頻段覆蓋傳導電流與位移電流的過渡帶，對電阻率、介電常數(shù)、磁導率三個參數(shù)均較敏感(Müller et al., 2012)；對于更高頻段，鐵磁性介質(zhì)對探地雷達波(n×10 MHz～ nGHz)的波速、反射及衰減均產(chǎn)生影響(Van et al., 2013)。

鑒于此，開展了同時考慮電阻率與磁導率的二維大地電磁正則化反演研究。通過正演模擬分析了磁導率變化對卡尼亞視電阻率及其相位的影響；建立了同時反演電阻率與磁導率的正則化反演目標函數(shù)，并通過高斯-牛頓法求解；采用函數(shù)映射將電阻率與磁導率映射到反演空間，將電阻率與磁導率限定在具有實際物理意義的變化區(qū)間，通過雙共軛梯度穩(wěn)定算法求解高斯-牛頓方程，確保反演穩(wěn)定性。開展了合成數(shù)據(jù)與實測數(shù)據(jù)的反演工作，表明從大地電磁數(shù)據(jù)可以恢復地下介質(zhì)磁導率，同時反演磁導率與電阻率提高了反演效果，豐富了大地電磁數(shù)據(jù)解釋的依據(jù)。

1 二維大地電磁電阻率與磁導率影響模擬

(1)

(2)

式(1)代表E型波(也稱TE模式)，式(2)代表 H型波(也稱TM模式)。二維大地電磁問題可以通過有限單元方法求解(Wannamaker et al., 1986; De Lug?o et al., 1997; Sheen et al., 2000; Key et al., 2006; 劉云等2010; 張志勇等, 2013, 2015)。

為確定磁導率對卡尼亞視電阻率及其相位的影響，進行圖1所示模型試算(張志勇等，2013)：在ρ0=100 Ω·m，μ0=4π×10-7H·m-1，ε0=8.854 2×10-12F·m-1的均勻半空間內(nèi)設(shè)置三個截面為正方形(邊長300 m)的二度體，二度體上頂埋深220 m，間距550 m。不考慮介電常數(shù)異常，M1異常體ρ1=50 Ω·m，μ1=2μ0，M2異常體，ρ2=100 Ω·m、μ2=2μ0，M3異常體ρ3=200 Ω·m、μ3=2μ0。在地表采用50 m點距進行大地電磁測量，觀測頻率范圍為2-12～212Hz。

圖1 試算模型Fig.1 Test model

同時考慮電阻率、磁導率變化條件下，圖1試算型的計算結(jié)果如圖2。由圖2a可見TE模式視電阻率均為高阻異常，M1、M2引起的異常向低頻發(fā)展不圈閉，形態(tài)類似于只考慮電阻率變化時高阻體TM模式響應；M3引起的異常成半圈閉并向低頻發(fā)展開口。TE模式視相位(圖2b)，形成上下兩個圈閉異常且下部圈閉異常明顯；M2、M3形成圈閉異常。TM模式視電阻率(圖2c)，形成低阻向低頻發(fā)展不圈閉異常，M3形成高阻向低頻發(fā)展不圈閉異常；M2形成高阻圈閉異常；TM模式視相位(圖2d)，M1、M2形成上下兩個圈閉異常且下部圈閉異常比上部明顯；形成圈閉異常。

圖2 試算模型計算結(jié)果Fig.2 Contours of apparent resistivity and phase for test modela.TE模式視電阻率；b.TE模式視相位；c.TM模式視電阻率；d.TM模式視相位

地下電導率與磁導率的變化對視電阻率及相位均產(chǎn)生明顯影響；受磁導率影響，TE模式異常形態(tài)向無磁導率影響TM模式形態(tài)轉(zhuǎn)變，可能將低阻異常變?yōu)楦咦璁惓?；試算模型，磁導率對TM模式異常形態(tài)影響沒有TE明顯，低阻異常仍為低阻；只有磁導率差異的異常體會產(chǎn)生明顯的異常，其TE模式的異常形態(tài)類似于無磁導率影響的TM模式的視電阻率異常。

2 同時反演電阻率與磁導率的正則化反演

數(shù)值模擬表明地下介質(zhì)的磁導率差異可引起卡尼亞視電阻率異常，下面探討采用正則化方法進行電阻率與磁導率的聯(lián)合反演算法。

只考慮電阻率差異的大地電磁從一維到二維反演已發(fā)展成熟，三維正逐步進入實用階段。一維反問題可采用奇異值分解方法進行迭代反演(Jupp et al., 1975)，Jupp等(1977)擴展了這一方法開展了二維大地電磁反演。為避免反問題的多解性，Constale等(1987)提出了著名的Occam反演方法，采用模型平滑約束進行一維大電磁與直流電阻率數(shù)據(jù)正則化反演。Occam反演在大地電磁數(shù)據(jù)反演中得到了廣泛應用，在一維基礎(chǔ)上發(fā)展了二維大地磁Occam反演(Degroot-Hedlin et al., 1990)；開展了數(shù)據(jù)空間DASOCC與縮減的數(shù)據(jù)空間REBOCC反演研究(Siripunvaraporn et al., 2000, 2009;Kordy et al., 2016)，利用二維與三維條件下大地電磁數(shù)據(jù)量一般遠小于模型參數(shù)的情況來減少反演運算量。求取和存儲靈敏度矩陣是正則化反演的關(guān)鍵問題，為減少計算量和內(nèi)存需求開展了大量的研究工作：如假定電導率的水平變化遠小于垂向變化，近似計算靈敏度矩陣的快速松弛反演(RRI)(Smith et al., 1991)；互換定理是減少靈敏度矩陣計算工作量的另一有效途徑(Rodi, 1976; Jupp et al., 1977; Sasaki, 1989; De Lug?o et al., 1996, 1997)；采用基于伴隨方程計算靈敏度矩陣(McGillivray et al., 1994)，其計算量與結(jié)合了互換定理的靈敏度計算量相當，但采用均勻半空間解析解作為伴隨方程解進行計算可大大減少計算量(Farquharson et al., 1996)；非線性共軛梯度法(NLCG)因為只需要計算靈敏度矩陣及其轉(zhuǎn)置與向量的乘積，可避免直接存儲靈敏度矩陣從而減少存儲量，并最終將靈敏度矩陣及其轉(zhuǎn)置與向量的乘積轉(zhuǎn)換為正演計算(Newman et al., 2000,2003; Rodi et al., 2001; Zhdanov et al., 2004)。高斯牛頓法是另一種求解非線性反問題的有效方法，與非線性共軛梯度相比具有目標函數(shù)下降快、所需迭代次數(shù)少的優(yōu)點；該方法的主要困難，一是靈敏度矩陣的計算、滿陣乘法運算的計算量大，二是正則化反演得到的海森矩陣條件數(shù)大，求解困難(Sasaki, 1989, 2001; Usui, 2015)；通過迭代解法求解高斯-牛頓方程，可以采用NLCG類似的方法避免顯式求解靈敏度矩陣，提高計算效率(Mackie et al., 1993; Grayver et al., 2013; Jahandari et al., 2017)。

研究工作構(gòu)建了同時反演電阻率與磁導率的正則化反演目標函數(shù)，通過高斯-牛頓法求解反問題，實現(xiàn)了有限內(nèi)存條件下大地電磁電阻率與磁導率同時反演。

2.1 同時考慮電阻率與磁導率的正則化反演目標函數(shù)

同時考慮電阻率與磁導率條件下，正則化問題目標函數(shù)定義為：

φ=φd+λρφp+λμφμ

(3)

其中數(shù)據(jù)誤差，

φd=(G(ρ,μ)-dobs)TWdTWd(G(ρ,μ)-dobs)

(4)

式中，G(ρ,μ)為正演函數(shù)，dobs為觀測數(shù)據(jù)，Wd為數(shù)據(jù)方差矩陣；φρ、φμ分別為電阻率、磁導率模型誤差項，可統(tǒng)一表達為(Zhdanov, 2002)：

φρ=(ρ-ρapr)TWρTWρ(ρ-ρapr)

(5)

φμ=(μ-μapr)TWμTWμ(μ-μapr)

(6)

式中，ρ、μ分別為介質(zhì)電阻率與磁導率，ρapr、μapr為參考模型的電阻率與磁導率，Wρ、Wμ為統(tǒng)一表示的電阻率與磁導率模型正則化矩陣；λρ=γλ、λμ=(1-γ)λ,λ為正則化因子,γ是平衡電阻率模型與磁導率模型正則化項的比例系數(shù)。

正則化因子選擇可采用廣義交叉驗證準則(GCV)(Farquharson et al., 2003; Haber et al., 2000)、“L-Curve”標準(Hansen, 1992; Li et al., 2000)，本文工作采用修正的“L-Curve”算法(李曼等, 2014)。平衡系數(shù)對磁導率的反演影響很大(Farquharson et al., 2003)，Routh等(1998)進行頻譜激電的Cole-Cole四參數(shù)反演時采用各參數(shù)正則化項的比值進行平衡，Loke等(2006)則通過靈敏度矩陣模的比值進行平衡。本文采用類似Loke的方法，通過下式計算平衡系數(shù)：

(7)

式中，Gρ、Gμ為觀測數(shù)據(jù)對電阻率與磁導率的靈敏度，可采用互換定理快速計算(Rodi, 1976; Sasaki, 1989; De Lug?o et al., 1996)。

為確保電阻率與磁導率在正常的巖石與礦物取值范圍內(nèi)，反演在變換域進行(Habashy et al., 2004; Commer et al., 2008; Key, 2009, 2016)。將電阻率、磁導率進行函數(shù)變換，設(shè)反演空間變量為x，m為模型空間變量且m∈(α,β)，則可采用如下變換關(guān)系，

(8)

(9)

(10)

采用迭代法求解正則化目標函數(shù)(3)，對第k+1次正演進行Taylor極數(shù)展開，可得

G(ρk+1,μk+1) =G(ρk,μk)+GρΔρ+GμΔμ+ο(Δρ,Δμ)

(11)

忽略二次以上高次項，代入式可得高斯-牛頓方程：

HkΔm=-φ

(12)

其中，

(13)

(14)

(15)

2.2 基于雙共軛梯度穩(wěn)定算法的有限內(nèi)存優(yōu)化

地球物理反問題得到的高斯-牛頓方程式，由于海森矩陣性質(zhì)不好，如果采用直接解法求解該方程很難求得穩(wěn)定解，很多學者采用正交分解或迭代解法求不完全解，如正交分解方法(Sasaki, 1989; Grayver et al., 2013)、共軛梯度法(CG)(Mackie et al., 1993)、廣義最小余量法(GMRES)(Jahandari et al., 2017)。本次工作采用雙共軛梯度穩(wěn)定算法(BICGSTB)(Saad, 1996)迭代求解式。計算過程如下：

(2)p0:=r0.

(3)Forj=0,1,…,untilconvergene Do:

(5)sj:=rj-ajHkpj

(6)wj:=(HKsj,sj)/(Hksj,Hksj)

(7)Δmj+1:=Δmj+αjpj+wjsj

(8)rj+1:=sj-wjHKsj

(10)pj+1:=rj+1+βj(pj-wjHKpj)

(11)End Do

算法中有關(guān)海森矩陣與向量的乘積，可展開為靈敏度矩陣及其轉(zhuǎn)置與向量的乘積，最終轉(zhuǎn)換為正演計算，從而避免直接存儲海森矩陣與靈敏度度矩陣，實現(xiàn)有限內(nèi)存條下的反演計算(Mackie et al., 1993; Rodi et al., 2001)。

3 合成數(shù)據(jù)反演

對圖1所示模型計算所得數(shù)據(jù)，加入3%的隨機噪聲，采用第二部分算法，進行試算，反演結(jié)果如圖3b和圖3c，并與只反演電阻率的反演結(jié)果(圖3a)進行對比分析。

圖3 合成數(shù)據(jù)反演結(jié)果Fig.3 The inversion result of synthetic data

由圖3可見，反演算法成功得到了異常體的磁導率(圖3b)和電阻率(圖3c)。由于異常體磁導率的影響，單獨反演電阻率(圖3a)效果不理想，低阻、高磁導率異常體M1位置反演得到了一個低阻體，但其兩側(cè)出了高阻假異常；M2位置出現(xiàn)了略高于背景電阻率的異常；M3高阻、高磁導率位置，出現(xiàn)了大塊高阻，高阻區(qū)遠大于M3實際的空間大小，并受M2高磁導率異常的影響向其延伸。另外，由正則化因子、數(shù)據(jù)均方根誤差、穩(wěn)定因子隨迭代次數(shù)變化曲線(圖4)，可見研究工作采用的算法，數(shù)據(jù)均方根誤差隨迭代次數(shù)的增加穩(wěn)定下降，同時反演電阻率與磁導率圖4b可以得到更小的數(shù)據(jù)誤差。研究作采用了漸進網(wǎng)格反演策略，網(wǎng)格細化過程重新設(shè)定正則化因子，所以正則化因子與模型誤差變化曲線是分段平滑的，限于篇幅該算法另行討論。

圖4 正則化因子、數(shù)據(jù)均方根誤差、穩(wěn)定因子隨迭代次數(shù)變化曲線Fig.4 Optimum Multiplier, Root mean square error, stabilizer variation curve with iteration timesa.直接反演；b.同時反演

4 某花崗巖區(qū)聯(lián)合反演電阻率與磁導率實例

采用本文算法對某花崗巖區(qū)音頻大地電磁數(shù)據(jù)進行了聯(lián)合反演。反演剖面長度7 km，設(shè)計點距50 m，實際采集測點138個?？碧焦ぷ鞯闹饕康氖翘綔y花崗巖內(nèi)部的斷裂構(gòu)造帶與輝綠巖脈。標本測試表明，該區(qū)花崗巖電阻率變化范圍為4 200～8 000 Ω·m，輝綠巖電阻率一般高于12 000 Ω·m；花崗巖磁化率變化不明顯，但輝綠巖磁化率明顯大于花崗巖，相對磁導率可達約2.7；研究區(qū)濕潤多雨，斷裂構(gòu)造與高阻圍巖相比存在明顯的電阻率差異。物性分析表明，該區(qū)存在電磁勘探的物性基礎(chǔ)。

反演采用了所有138個測點數(shù)據(jù)，頻率從6.9 Hz到10.4 kHz共42個頻點。分別進行了單獨反演電阻率(圖5a)和聯(lián)合反演電阻率與磁導率(圖5b、c)的反演工作。總體上看，單獨反演電阻率斷面相比聯(lián)合反演所得到的電阻率斷面電阻率值要高，這與數(shù)值模擬計算結(jié)果相對應，高磁導率有類似于高阻的影響；單獨反演得到的電阻率斷面，巖體內(nèi)部出現(xiàn)近水平的低阻條帶狀異常，很大部分與地質(zhì)條件不符；聯(lián)合反演得到的電阻率斷面，淺地表低阻異常明顯，對比地質(zhì)斷面低阻異常體與出露的斷裂構(gòu)造對應，同時揭示了多條隱伏構(gòu)造。聯(lián)合反演出的磁導率斷面，以低值為主，并出現(xiàn)了幾條由地表向深部發(fā)展的高磁化率局部異常，與輝綠巖脈對應；大于526 500 m的區(qū)域出現(xiàn)了大面積磁導率中高值區(qū)，推測相較斷面其它區(qū)花崗巖完整性更好或者為不同期次侵入的花崗巖體；測線起點端，存在一長透鏡狀高磁化率區(qū)，推測與巖體內(nèi)部鐵磁性物質(zhì)升高有關(guān)。

圖5 某花崗巖區(qū)反演結(jié)果Fig.5 The inversion result of a granite area MT data

5 結(jié)論與建議

通過數(shù)值模擬，論證了磁導率異常對大地電磁勘探的影響不能忽略；采用正則化反演技術(shù)實現(xiàn)了同時考慮介質(zhì)電阻率與磁導率的大地電磁二維正則化反演。主要取得如下認識：

(1)大于真空中磁導率的局部異常體對大地電磁觀測存在明顯的影響，二倍電阻率差異條件下、二倍于真空中磁導率的異常體足以改變卡尼亞視電阻率與相位形態(tài)；

(2)通過正則化反演可以實現(xiàn)電阻率與磁導率的反演，同時反演兩個參數(shù)有利于提高反演效果；

(3)模型空間到反演空間的函數(shù)映射，將電阻率與磁導率限定在具有實際物理義的變化區(qū)間，提高反演的穩(wěn)定性；

(4)高斯-牛頓法可實現(xiàn)同時考慮電阻率與磁導率的正則化目標函數(shù)優(yōu)化，雙共軛梯度穩(wěn)定算法求解高斯-牛頓系統(tǒng)，求解過程穩(wěn)定、計算高效。