彭榮華， 胡祥云， 韓波

1 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢 430074 2 不列顛哥倫比亞大學(xué)地球、海洋與大氣科學(xué)學(xué)院， 溫哥華V6T 1Z4， 加拿大 3 中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島 266100

?

基于高斯牛頓法的頻率域可控源電磁三維反演研究

彭榮華1,2， 胡祥云1*， 韓波3

1 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院， 武漢 430074 2 不列顛哥倫比亞大學(xué)地球、海洋與大氣科學(xué)學(xué)院， 溫哥華V6T 1Z4， 加拿大 3 中國海洋大學(xué)海洋地球科學(xué)學(xué)院， 青島 266100

三維反演解釋是電磁法勘探發(fā)展的重要趨勢，而如何提高三維反演的可靠性、穩(wěn)定性和計算效率是算法開發(fā)者們目前的研究重點.本文實現(xiàn)了一種頻率域可控源電磁(CSEM)三維反演算法.其中正演基于擬態(tài)有限體積法離散化，利用直接矩陣分解技術(shù)來求解大型線性系統(tǒng)方程，不僅準(zhǔn)確、穩(wěn)定，而且特別有利于含有大量發(fā)射場源位置的CSEM勘探情況；對目標(biāo)函數(shù)的最優(yōu)化采用高斯牛頓法(GN)，具有近似二次的收斂性；使用預(yù)條件共軛梯度法(PCG)求解每次GN迭代所得到的法方程，避免了顯式求解和存儲靈敏度矩陣，減小了計算量.以上這些方法的結(jié)合應(yīng)用，使得本文的三維反演算法準(zhǔn)確、穩(wěn)定且高效.通過陸地和海洋CSEM勘探場景中的典型理論模型的反演測試，驗證了本文算法的有效性.

頻率域可控源電磁法； 三維反演； 高斯牛頓法； 直接矩陣分解法； 預(yù)條件共軛梯度法

Quantitative interpretation of large-scale CSEM data in the frequency domain requires efficient and stable 3D forward modeling and inversion codes. Considerable efforts have been contributed to developing numerical algorithms concerning 3D inversion of CSEM data that can accurately and efficiently recover subsurface electrical structure over last decade. Most existing 3D inversion algorithms employ Krylov subspace iterative methods as their forward solvers. Iterative techniques usually need little memory due to only matrix-vector products storage required, and they are fast for computation of single field solution. However, there are two main issues when using iterative solvers for 3D CSEM problems: (1) The ill-conditioning of linear systems arising from discretization of Helmholtz equation can lead to poor behavior of iterative solvers and even divergence in some cases. (2) Iterative solvers are very time-consuming for multi-source problems. These difficulties become major impediments when solving large-scale multi-source CSEM problems using iterative solvers.

Given the availability of more powerful workstations or computer clusters, direct methods have been increasingly used for solving 3D CSEM problems. Compared to iterative methods, direct solvers have several distinct advantages: (1) They provide more accurate solutions. (2) They are less prone to ill-conditioning of matrix, making them more robust for highly heterogeneous models or non-uniform grids. And (3) they separate the solving of linear system into an expensive matrix factorization and comparably inexpensive forward-backward substitution steps, which make them more suitable for multi-source CSEM surveying.

In this paper, we present an efficient inversion algorithm for 3D inversion of multi-frequency and multi-source CSEM data, which is based on a direct solver for solving the forward problem. The Gauss-Newton (GN) optimization algorithm is applied considering its high convergence rate, thus limiting the number of expensive matrix factorization required when using a direct solver. A preconditioned conjugate gradient solver (PCG) is used to solve the normal equation resulted from linearization at each GN iterate, in order to avoid computing and storing sensitivity matrix explicitly. This scheme only requires matrix-vector products of Jocabian and its transpose with vectors, which are equivalent to one forward and one adjoint problem. In addition, the matrix factorization obtained when solving forward problem can be used in subsequent PCG iterates, which dramatically speeds up PCG iteration and reduces computational costs. Numerical experiments on synthetic data from land and marine CSEM surveys show that our inversion algorithm has an excellent convergence rate and only tens to several tens of iterates are needed to reach desired misfit, demonstrating its efficiency and stability.

1 引言

頻率域可控源電磁法(controlled-source electromagnetic, CSEM)由于優(yōu)良的分辨能力而廣泛應(yīng)用于陸地和海洋礦產(chǎn)和油氣資源勘探中(Constable，2010；Yang and Oldenburg, 2012； Hu et al., 2013；Grayver et al., 2014).隨著電磁勘探越來越多地深入到復(fù)雜地質(zhì)環(huán)境中，一維和二維地質(zhì)假設(shè)難以準(zhǔn)確描述勘探目標(biāo)體全方位的空間展布特征(Ledo, 2005； Siripunvaraporn et al., 2005)，電磁數(shù)據(jù)的三維反演越來越受到重視.CSEM的三維反演研究在最近十年取得了極大的進展，國內(nèi)外一些學(xué)者陸續(xù)實現(xiàn)了有效的CSEM三維反演算法.例如Haber等(2007)實現(xiàn)了時間域CSEM三維反演；Commer和Newman(2008)、Plessix和Mulder(2008)實現(xiàn)了海洋CSEM三維反演；林昌洪等(2012)、翁愛華等(2012)和趙寧等(2016)針對陸地CSAMT勘探進行了三維反演研究；劉云鶴和殷長春(2013)則完成了航空CSEM數(shù)據(jù)的三維反演算法等等.盡管已有這些成功的研究工作，但越來越大的數(shù)據(jù)規(guī)模(例如最新的海洋CSEM發(fā)射與接收一體式的拖曳式測量)以及越來越真實且復(fù)雜的地電模型假設(shè)(例如各向異性、高電性差異等)對CSEM三維反演的可靠性、穩(wěn)定性和計算效率仍然是巨大的考驗，因此提升CSEM三維反演算法的性能在長時間內(nèi)仍將是國內(nèi)外電磁法學(xué)者的研究重點之一.

反演的性能在很大程度上取決于正演計算.對當(dāng)前主流的正演方法包括有限差分、有限體積和有限單元法來說，大型稀疏線性方程組的求解是非常關(guān)鍵的一步.在早期計算機內(nèi)存的限制下，迭代法是唯一可行的求解方程的方法，這是由于迭代解法只需計算和存儲矩陣與向量的乘積，所需內(nèi)存較少，運算速度較快，例如以上作者均采用了Krylov子空間迭代法.然而迭代解法的收斂性與系數(shù)矩陣的譜性質(zhì)息息相關(guān)(Saad, 2003)，而模型的電導(dǎo)率變化范圍大和網(wǎng)格尺寸的非均勻性都會使得系數(shù)矩陣的條件數(shù)非常大，從而可能導(dǎo)致迭代法的收斂性不佳、計算時間長，甚至不收斂.隨著大型稀疏矩陣分解算法的持續(xù)優(yōu)化(Amestoy et al., 2001, 2006；Schenk and G?rtner, 2006)，以及計算機硬件技術(shù)的快速發(fā)展，使得在從工作站到并行集群上利用直接解法求解電磁法三維正演中的大型線性方程組成為可能.與迭代解法相比，對于同一線性方程組，直接解法需要消耗更多的計算內(nèi)存，計算時間一般也會更長.盡管如此，直接解法相對于迭代解法有三個優(yōu)點：一是求解精度更高(韓波等, 2015)；二是更加穩(wěn)定，即其求解時間和精度基本不受系數(shù)矩陣的條件數(shù)影響，對強電性差異和含多尺度構(gòu)造的模型同樣適用；三是直接解法將線性方程組的求解過程分解為系數(shù)矩陣的分解和前代-回代(forward-backward substitution)求解過程(Streich, 2009)，其中系數(shù)矩陣的分解雖然消耗大量的計算機內(nèi)存和計算時間，但對于單一頻率多個場源位置的問題，只需進行一次矩陣分解加上多次計算量極小的替代過程即可，而不像迭代解法對每個頻率每個場源都要單獨求解，因此直接解法特別適合具有多個場源的CSEM三維模擬(Oldenburg et al.，2013).正因為這些優(yōu)點，基于矩陣分解的直接解法最近幾年開始應(yīng)用于CSEM三維反演計算(Grayver et al.，2013；Schwarzbach and Haber, 2013).

正演中線性方程組的求解使用直接解法還是迭代解法，對反演中最優(yōu)化算法的選擇有直接影響.目前三維反演中廣泛應(yīng)用的NLCG法(Commer and Newman, 2008；翁愛華等，2012；劉云鶴和殷長春，2013)和QN法(Plessix and Mulder, 2008；Schwarzbach and Haber, 2013；翁愛華等，2015；趙寧等，2016)，盡管模型的每一迭代更新非常容易計算，但由于僅利用了目標(biāo)函數(shù)的一階導(dǎo)數(shù)(即梯度)信息，只具有線性收斂性(Wright and Nocedal, 1999)，模型往往需要更新很多次才能收斂.相較之下，高斯牛頓法(Grayver et al.，2013； Oldenburg et al，2013)由于同時利用了目標(biāo)函數(shù)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)(即Hessian矩陣)信息，表現(xiàn)出近似二次收斂性(Nocedal and Wright, 1999)，所需的模型更新次數(shù)往往遠少于NLCG和QN.當(dāng)正演使用直接解法時，反演中每更新一次模型必然需要進行新的矩陣分解，即模型更新的次數(shù)越多，需要的矩陣分解次數(shù)也越多；由于矩陣分解要消耗大量的內(nèi)存和時間，為減少反演迭代中矩陣分解的次數(shù)，相比NLCG和QN，高斯牛頓法是更好的選擇.并且，每次高斯牛頓迭代的法方程可以使用預(yù)條件共軛梯度法(PCG)求解，不僅可避免顯式求解和存儲靈敏度矩陣，減小計算量；而且還能夠重復(fù)利用正演計算階段所得到的矩陣分解結(jié)果，加快反演過程.

本文嘗試?yán)米钚聝?yōu)化的矩陣直接求解器，結(jié)合高斯牛頓最優(yōu)化算法，實現(xiàn)相對可靠、穩(wěn)定、高效且同時適用于陸地和海洋勘探場景的頻率域CSEM三維反演算法.首先簡要介紹了CSEM三維正演算法；接著詳細介紹了反演算法的細節(jié)：包括目標(biāo)函數(shù)的構(gòu)建、靈敏度矩陣計算、法方程的求解、正則化函數(shù)及正則化參數(shù)的選取.最后通過典型三維理論模型的反演測試來驗證所開發(fā)的算法的有效性.

2 三維正演計算

CSEM三維正演基于電場總場的Helmholtz方程為

(1)

其中，E為電場強度(V·m-1)，ω為角頻率(rad/s)，σ為介質(zhì)電導(dǎo)率(S·m-1)，μ為磁導(dǎo)率(H·m-1)，Js為外部激發(fā)場源的電流密度(A·m-2).使用擬態(tài)有限體積法(mimeticfinitevolume,MFV)對(1)式離散化.MFV方法的最大特點是能夠?qū)B續(xù)的微分算符進行精確的離散模擬，并確保離散化后的矢量場仍然滿足其對應(yīng)的連續(xù)形式的矢量性質(zhì)及物理特性(HymanandShashkov, 1999a,b)，有效避免了偽解的產(chǎn)生.對于正交規(guī)則網(wǎng)格及各向同性介質(zhì)來說，該離散化方法與常用的交錯網(wǎng)格有限差分法相一致(Smith, 1996).但MFV適用范圍更廣，其對于非正交網(wǎng)格同樣適用(HymanandShashkov, 1997).特別是對于具有高度電性差異及各向異性介質(zhì)，MFV能夠自然地得到對稱的離散化形式(Hymanetal., 1997；HaberandRuthotto, 2014)，這對于復(fù)雜三維介質(zhì)的模擬十分有效.關(guān)于MFV的具體離散化過程的推導(dǎo)可參見以上幾篇文獻.

采用交錯網(wǎng)格對(1)式進行離散化得到待求解的線性方程組為

(2)

其中系數(shù)矩陣A為大型稀疏正定對稱復(fù)矩陣.本文通過調(diào)用基于多波前分解算法的MUMPS (Amestoy et al., 2001；Amestoy et al., 2006)線性運算庫對A進行LDLH分解后，再求解方程(2)，即直接解法.對同一個模型的某一個頻率，A是固定的，場源位置的改變僅僅會改變右端項q，因此只需進行一次矩陣分解，使得多個場源位置的CSEM正演計算量與單一場源相當(dāng).此外，直接解法還有精度高、穩(wěn)定性好等優(yōu)點，具體可參考Streich(2009)、Oldenburg等(2013)和韓波等(2015).

3 三維反演理論

3.1 目標(biāo)函數(shù)的構(gòu)建

由于反演問題的不適定性，為獲得穩(wěn)定的最優(yōu)解，CSEM三維反演的目標(biāo)函數(shù)通常采用Tikhonov正則化方式來構(gòu)建(Tikhonov and Arsenin, 1977；Constable et al., 1987)，目標(biāo)函數(shù)可定義為

(3)

其中φd(m)為數(shù)據(jù)擬合差函數(shù)，φm(m)為模型正則化函數(shù)，β為Tikhonov正則化參數(shù).為保證反演中模型參數(shù)的物理意義，模型參數(shù)m設(shè)為m=lnσ.數(shù)據(jù)擬合差函數(shù)可表示為

(4)其中dobs為觀測的電磁場分量，Cd為數(shù)據(jù)協(xié)方差矩陣.P為觀測矩陣，其將正演計算得到的網(wǎng)格采樣點處的電磁場分量投影到測點處，觀測矩陣采用“體積加權(quán)”的線性插值方式得到(韓波，2015).模型正則化函數(shù)為

(5)

其中mref為反演計算中參考模型，包括模型的先驗信息.Wm為模型平滑矩陣.對于最小結(jié)構(gòu)模型泛函，Wm為單位矩陣；對于平滑結(jié)構(gòu)泛函，Wm為一階差分矩陣.本文采用平滑結(jié)構(gòu)泛函來構(gòu)建模型正則化函數(shù)(Lelièvre and Oldenburg, 2009)，公式為

(6)

其中Wx、Wy、Wz分別為x、y、z方向上的一階差分矩陣.系數(shù)αx、αy、αz控制著模型結(jié)構(gòu)在不同方向上的平滑程度.

對目標(biāo)泛函(3)在當(dāng)前模型變量處進行二階Taylor展開，并忽略高階項，得到高斯-牛頓法(Gauss-Newton, GN)的法方程(normal equation)為

(7)

(8)

(9)

其中J(m) 為CSEM的靈敏度矩陣.

3.2 法方程的求解

為獲得模型更新，在每次GN迭代中，都需求解(7)式的法方程.直接求解(7)式需要顯式地計算和存儲靈敏度矩陣以及海森矩陣.在梯度類反演算法中，靈敏度矩陣扮演著重要的角色.對于CSEM反演，利用隱式微分方法(Rodi and Mackie, 2001；Egbert and Kelbert, 2012)，靈敏度矩陣可表示為

(10)

其中P為(4)式中的觀測矩陣，A(m)為(2)式中的正演系數(shù)矩陣，u為正演得到的電磁場值.G(m, u)為

(11)

從(10)(11)式中可以看到，靈敏度矩陣的求解需要利用到所有發(fā)射場源(包括所有頻率及發(fā)射場源位置)的電磁響應(yīng)，顯式地求解涉及到大量正演計算.對于多頻率多場源的三維CSEM反演來說，需要消耗大量的計算內(nèi)存以及運算時間，現(xiàn)有計算資源往往無法滿足要求.為避免顯式的計算和存儲靈敏度矩陣，通常采用預(yù)條件共軛梯度法(preconditioned conjugate gradient，PCG，Barrett et al., 1994)來迭代求解(7)式(Haber et al., 2007；Oldenburg et al., 2013；Schwarzbach and Haber, 2013).在PCG算法中，只涉及到求解和存儲靈敏度矩陣與模型向量的乘積，以及靈敏度矩陣轉(zhuǎn)置與數(shù)據(jù)向量的乘積(林昌洪等，2012).圖1給出了求解靈敏度矩陣與模型向量乘積、以及靈敏度矩陣轉(zhuǎn)置與數(shù)據(jù)向量乘積的算法，相關(guān)計算公式的詳細推導(dǎo)，讀者可參考林昌洪等(2012)，本文不再贅述.具體的算法計算細節(jié)可見附錄A.從圖1中可以看出靈敏度矩陣(或轉(zhuǎn)置)與向量的乘積相當(dāng)于求解一次正演問題，因此每次PCG迭代需要求解兩次正演問題.由于采用直接解法進行正演計算，在同一次GN迭代中，矩陣分解的結(jié)果可以多次重復(fù)利用.因此，相比于迭代法，直接解法在反演計算中的優(yōu)勢會更加明顯(Oldenburg et al., 2013).

當(dāng)采用PCG來求解法方程(7)時，GN反演過程的計算消耗將主要集中在PCG的迭代上.為減少反演的計算時間，一個可行的策略是降低法方程求解的精度來減少每次GN迭代中所需的PCG迭代次數(shù)，這就是所謂的非精確高斯牛頓法(Inexact Gauss-Newton，IGN, Nocedal and Wright,1999).對于IGN來說，只要法方程的殘差(或目標(biāo)函數(shù)的梯度)隨著PCG迭代不斷降低，該算法的收斂性就能得到保證(Rieder, 2001,2005).由于反演的非線性特征，對于線性化迭代反演來說，這種非精確求解不僅減小了計算量，而且避免了每次GN迭代中對法方程的過度求解(oversolving，Eisenstat and Walker, 1996)，使得在反演搜索過程中能有機會跳出局部極小陷阱.另外，PCG迭代求解具有正則化效應(yīng)，并且PCG的迭代次數(shù)扮演著正則化參數(shù)的作用(van den Doel and Ascher, 2007, 2012)，使得PCG在求解法方程過程中具有很好的穩(wěn)定性.在三維反演中，一般PCG迭代次數(shù)取為十幾次到幾十次(Oldenburg et al., 2013).

對于PCG迭代中預(yù)條件子的選擇主要有兩種方式.一種方式是在每次GN迭代中，首先利用擬牛頓法(如L-BFGS)來近似得到目標(biāo)泛函Hessian矩陣的逆，將其作為PCG迭代的預(yù)條件因子(Oldenburg et al., 2013)；另一種方式是在每次GN迭代中，先求解線性方程

(12)

3.3 模型更新

在利用PCG求得法方程(7)的解后，需要進行線搜索來獲得每次GN迭代的模型更新量，公式為

(13)

步長α控制著模型更新大小，步長α的精確搜索需要求解關(guān)于α的單變量最優(yōu)化問題φ(mk+αδmk).對于三維反演來說，精確線搜索往往因計算消耗巨大而無法進行.因此一般采用非精確線搜索方式來確定更新步長α.最常見的策略是采用(14)式的弱Wolfe條件(NocedalandWright,1999)來進行非精確線性搜索，得到更新步長α.公式為

(14)

其中常數(shù)c一般取10-4.

3.4 正則化策略及正則化參數(shù)的選取

為了求解(3)式的最優(yōu)化問題，需要在反演迭代過程中選取合適的正則化參數(shù).由于正則化參數(shù)起著平衡數(shù)據(jù)擬合和模型結(jié)構(gòu)的作用，因此，正則化參數(shù)的選取必須兼顧反演的穩(wěn)定性以及先驗?zāi)Ｐ蛯Ψ囱萁Y(jié)果的影響.一方面正則化參數(shù)必須設(shè)置的充分大來保證反演的穩(wěn)定性；另一方面必須設(shè)置的足夠小來減小先驗?zāi)Ｐ蛯Ψ囱莸阉鲙淼钠?Siripunvaraporn, 2012；Grayver et al., 2013).對于Tikhonov正則化來說，可以通過數(shù)據(jù)擬合差和模型范數(shù)之間的差異準(zhǔn)則(discrepancy principle)來確定最優(yōu)的正則化參數(shù).在實際反演中一般通過采用所謂的松弛法(relaxation/cooling scheme)來選取正則化參數(shù)(Newman and Alumbaugh, 1997；Rodi and Mackie, 2001)：在反演的初期階段，選取較大的正則化參數(shù)來確保反演的穩(wěn)定性；隨著反演的進行，逐漸減小正則化參數(shù)來降低模型正則化項的權(quán)重，使反演集中在數(shù)據(jù)的擬合上，直到達到最佳的數(shù)據(jù)擬合差.本文利用松弛法來自動選取GN迭代過程中正則化參數(shù)，其更新方式為

(15)

其中γ為衰減因子，β0為初始正則化參數(shù)，N為當(dāng)前GN迭代數(shù).

算法1:計算靈敏度矩陣J(m)與模型向量v的乘積J(m)v1.計算:q=G(m,e)v2.求解正演問題:A(m)u=q3.計算:J(m)p=-Pu算法2:計算靈敏度矩陣轉(zhuǎn)置J(m)T與數(shù)據(jù)向量w的乘積J(m)Tw1.計算:q=PTw2.求解伴隨正演問題:A(m)Tu=q3.計算:J(m)Tw=-G(m,e)Tu

圖1 靈敏度矩陣及其轉(zhuǎn)置與向量乘積算法

Fig.1 Algorithms for matrix-vector products of Jocabian and its transpose with a vector

4 理論模型合成數(shù)據(jù)反演測試

為了驗證所開發(fā)的CSEM三維反演算法的有效性，本文設(shè)計了兩個三維模型來分別模擬陸地及海洋CSEM測量情況，并對這兩個三維模型的合成數(shù)據(jù)進行了三維反演測試.本文所有計算均在一個小型并行機上完成.該小型并行機群系統(tǒng)配置有8個計算節(jié)點，每個計算節(jié)點含有2個4核Intel○RXeon○RE5410型處理器，主頻為2.33 GHz，每個計算節(jié)點配置有32G內(nèi)存.操作系統(tǒng)為CentOS 5.5.

圖2 陸地CSEM三維模型平面(a) 水平切片圖z=800 m； (b) 垂直剖面圖y=2000 m.三角為發(fā)射場源位置，圓點為接收測站位置.Fig.2 Planar views of 3D model for land CSEM survey(a) Horizontal slice at z=800 m； (b) Cross-section at y=2000 m. Triangles and circles indicate receivers and transmitters, respectively.

4.1 陸地三維模型

首先考慮陸地CSEM測量情況，設(shè)計了如圖2所示的三維雙棱柱體地電模型.其中高阻和低阻棱柱體尺寸均為1000 m × 2000 m × 500 m，頂部埋深為500 m，電阻率分別為100 Ωm和1 Ωm.兩個棱柱體均埋藏在電阻率為10 Ωm的均勻半空間中.為得到高精度理論合成數(shù)據(jù)，模型核心區(qū)域被剖分為60 ×40 ×40的細密網(wǎng)格.為降低邊界條件對計算結(jié)果的影響，在模型邊界外各增加8個網(wǎng)格，網(wǎng)格增長因子為2.另外，計算模型區(qū)域包含10層空氣層，空氣層電阻率設(shè)為108Ωm.整個正演計算的模型區(qū)域共劃分成76×56×58的網(wǎng)格.觀測系統(tǒng)包括20個發(fā)射站，場源間距為800 m，80個接收站，測點間距為400 m.發(fā)射站和測點均布設(shè)在地表處，如圖2所示.每個發(fā)射站均布設(shè)沿x方向極化的長度為100 m的水平接地導(dǎo)線源，發(fā)射頻率為0.25 Hz和1 Hz.正演計算得到所有收發(fā)距大于600 m的電場Ex分量(包括實部和虛部)作為理論合成數(shù)據(jù)，然后對理論合成數(shù)據(jù)添加3%的隨機噪聲得到反演觀測數(shù)據(jù).

為測試反演算法的穩(wěn)定性，避免反演偏差，反演網(wǎng)格采用不同于合成數(shù)據(jù)所用的正演計算網(wǎng)格.反演核心區(qū)域劃分為40×30×40的網(wǎng)格，其中場源處網(wǎng)格進行了局部加密.同樣地，在核心區(qū)域外添加邊界網(wǎng)格和空氣層，整個反演計算區(qū)域劃分成52×42×56的網(wǎng)格.反演初始正則化參數(shù)設(shè)置為0.01，衰減因子γ為2.GN迭代中的PCG迭代次數(shù)設(shè)為20次.目標(biāo)擬合差設(shè)為1.反演初始模型及參考模型均采用電阻率為10 Ωm的均勻半空間.反演迭代中，空氣層電阻率保存不變.整個反演共耗時約12.5小時.

圖3給出了三維反演所得到的最后地電模型結(jié)果.從圖3a和圖3b可以看出，反演結(jié)果較好地刻畫出異常高低阻棱柱體的形態(tài).其中對于低阻異常體的恢復(fù)要明顯優(yōu)于高阻異常體，這主要是由于CSEM對于高低阻異常體的靈敏度及分辨率不同引起的.圖3c給出了反演結(jié)果的三維圖.對于低阻異常體來說，無論是異常體的電阻率還是異常體的位置分布，反演結(jié)果均與理論模型較為一致；對于高阻異常體，反演結(jié)果的最大電阻率數(shù)值約為50 Ωm，與真實模型有一定差距，另外，恢復(fù)的高阻體異常分布范圍要小于真實模型分布區(qū)域.

圖4a為反演過程中數(shù)據(jù)擬合差的變化.經(jīng)過9次GN迭代，反演數(shù)據(jù)擬合差從最初的12.78逐漸收斂到最終的1.01，顯示了GN反演算法優(yōu)良的收斂性及穩(wěn)定性.圖4b給出了反演合成數(shù)據(jù)與反演初始模型響應(yīng)數(shù)據(jù)及反演最終迭代模型響應(yīng)數(shù)據(jù)的振幅比率統(tǒng)計圖.從中可以看出，反演結(jié)果對于所有頻點和測點的數(shù)據(jù)總體擬合較為一致，最終的數(shù)據(jù)擬合誤差降為所添加的隨機誤差范圍之內(nèi)(3%).與之對應(yīng)，圖5顯示了頻率為1 Hz時，不同發(fā)射位置的所有測點的觀測數(shù)據(jù)(圖5a—c)、觀測數(shù)據(jù)與初始模型數(shù)據(jù)振幅比(圖5d—f)、以及觀測數(shù)據(jù)與反演最終模型數(shù)據(jù)振幅比(圖5g—i).從圖中可以看出，所有測點的數(shù)據(jù)擬合均勻一致，并未出現(xiàn)較明顯的數(shù)據(jù)“過擬合”或“欠擬合”現(xiàn)象，表明GN算法具有良好的一致性.

4.2 海洋三維模型

接著考慮海洋CSEM測量情況.為測試所開發(fā)的反演算法對于復(fù)雜模型的有效性，本文將地震勘探中常用于偏移成像算法測試的三維鹽穹體(salt dome)模型(Aminzadeh et al., 1997)進行修改來模擬海洋可控源勘探情況.整個鹽穹模型區(qū)域大小為6 km×6 km×2.8 km(含海水層)，如圖6所示.鹽穹體的電阻率設(shè)為100 Ωm.均勻海底沉積層電阻率為1 Ωm，海水層厚度為1 km，電阻率為0.33 Ωm.采用三維測網(wǎng)數(shù)據(jù)采集方式來模擬真實海洋可控源測量情況，如圖7所示.整個測量網(wǎng)格包括12條正交的測量線，每條測線長為6 km，測線間距為1 km.發(fā)射偶極子采用水平電偶極子(HED)，偶極方向與測線方向一致，位于海底上方50 m處.發(fā)射偶極子每隔300 m向海底發(fā)射頻率為0.25 Hz、0.75 Hz和1.25 Hz的電磁信號，共形成252個發(fā)射源位置.36個布設(shè)于海底的接收站在測區(qū)內(nèi)呈均勻分布，同時接收電場x分量和y分量，測點間距為1 km.為得到用于反演的合成數(shù)據(jù)，與上述陸地三維模型剖分方式類似，正演網(wǎng)格剖分為92×92×54.正演計算得到所有場源的電場Ex分量和Ey分量.考慮到電磁數(shù)據(jù)的冗余情況，本文只利用與發(fā)射源偶極方向一致的電場分量(即對于沿x方向極化的偶極源只利用Ex分量，而對于沿y方向極化的偶極源只利用Ey分量)作為反演數(shù)據(jù)，從而有效減少存儲需求和反演計算量.對理論合成數(shù)據(jù)添加5%的隨機噪聲得到反演的觀測數(shù)據(jù).為避免反演中對振幅較小的數(shù)據(jù)的過度擬合，數(shù)據(jù)誤差限設(shè)為5×10-15V·m-1(Constable and Weiss, 2006)，總共形成23100個海洋可控源電場數(shù)據(jù).反演網(wǎng)格采用不同于正演計算的網(wǎng)格，網(wǎng)格剖分為72×72× 46.反演初始正則化參數(shù)設(shè)置為0.01，為避免正則化因子下降過快，衰減因子γ為1.5.GN迭代中的PCG迭代次數(shù)設(shè)為15次.目標(biāo)擬合差設(shè)為1.反演迭代中，空氣層和海水層電導(dǎo)率均保持不變.

圖3 陸地CSEM三維反演結(jié)果(a) 水平切片圖z=800 m； (b) 垂直剖面圖y=2000 m； (c) 三維剖面圖(電阻率大于20 Ωm或電阻率小于5 Ωm).圖中三角為發(fā)射場源位置，圓點為接收測站位置，黑框表示真實模型位置.Fig.3 Inversion results of land CSEM survey(a) Horizontal slice at z=800 m； (b) Cross-section at y=2000 m； (c) 3D view. Cells with resistivity between 5 Ωm and 20 Ωm are cut-off. Triangles and circles indicate receivers and transmitters, respectively. Squares outline true locations of anomaly objects.

圖4 陸地可控源三維反演數(shù)據(jù)擬合統(tǒng)計(a) 數(shù)據(jù)擬合差RMS隨迭代變化圖，紅色虛線為目標(biāo)擬合差； (b) 反演合成數(shù)據(jù)與初始模型響應(yīng)數(shù)據(jù)(藍色)和反演最終模型響應(yīng)數(shù)據(jù)(紅色)振幅比率統(tǒng)計圖，紅色虛線為3%的數(shù)據(jù)誤差區(qū)間.Fig.4 Fitting statistics of 3D inversion(a) Data RMS misfit versus GN iteration for land CSEM survey. Red dashed line shows the desired data misfit at RMS=1；(b) Histogram of the amplitude ratios between observed data used for inversion and predicated data of initial iteration (blue) and of final iteration (red). The red dashed line indicates 3% data error interval.

圖5 陸地CSEM頻率為1 Hz時，(a—c)不同發(fā)射位置的所有測點的觀測數(shù)據(jù)振幅、(d—f)觀測數(shù)據(jù)與初始模型數(shù)據(jù)振幅比、(g—i)觀測數(shù)據(jù)與反演最終模型數(shù)據(jù)振幅比.紅色三角為發(fā)射場源位置，黑色圓點為接收測站位置.Fig.5 (a—c) Amplitudes of observed data. (d—f) Amplitude ratios between observed data and predicated data of starting model. (g—i) Amplitude ratios between observed data and predicated data of finial model for different transmitter locations with 1Hz transmission frequency for land CSEM survey. Triangles and circles denote receivers and transmitters, respectively.

圖8給出了三維反演過程中數(shù)據(jù)擬合差的變化情況.整個反演共進行了18次GN迭代，反演數(shù)據(jù)擬合差從初始的17.13逐漸收斂到最后的1.07，顯示了GN反演算法對于復(fù)雜模型同樣具有優(yōu)良的收斂性和穩(wěn)定性.每次GN迭代耗時約1小時40分鐘.

圖9給出了鹽穹體模型三維反演的不同剖面結(jié)果.可以看出，反演總體上很好地恢復(fù)了鹽穹體模型在空間上的分布.對比反演結(jié)果中的高阻區(qū)域與鹽穹體模型的真實區(qū)域(圖中黑色曲線所圈定區(qū)域)，可以發(fā)現(xiàn)根據(jù)反演結(jié)果能夠有效地推斷出鹽穹體的結(jié)構(gòu)特征，并估計出相應(yīng)的邊界位置.另外對于鹽穹電阻率值的恢復(fù)，反演在鹽穹核心區(qū)域處可達到80 Ωm以上 (如圖9h,i)，展示了海洋可控源電磁法對于高阻目標(biāo)體優(yōu)良的分辨能力.另一方面，由于鹽穹模型的不規(guī)則特征，在某些區(qū)域內(nèi)，反演結(jié)果并未能較好地刻畫出真實的形態(tài)特征(如圖9a,g)，這可能與本文所采用的平滑模型正則化泛函有關(guān)，若使用其他形式的模型正則化泛函有可能效果更好，這不在本文討論范圍以內(nèi).

圖6 海洋可控源三維鹽穹體模型剖面(a) XZ剖面(y=3200 m)； (b) YZ剖面(x=2000 m).Fig.6 Planar views of 3D salt dome reservoir for marine CSEM survey(a) XZ section at y=3200 m； (b) YZ section at x=2000 m.

圖7 海洋可控源三維觀測系統(tǒng)平面其中紅色圓點為布設(shè)于海底的接收站，黑線為發(fā)射源的航行線，距海底50 m，藍線為鹽穹模型在深度為1750 m時的平面投影.Fig.7 Sketch showing data acquisition geometry for marine CSEM surveySurvey consists of thirty-two sea-bottom detectors denoted by red circles and twelve sail lines shown in black lines. The blue curve indicates the projection of the salt dome reservoir at the depth of 1750 m.

圖8 海洋可控源鹽穹體模型反演數(shù)據(jù)擬合差RMS隨迭代次數(shù)變化(紅色虛線為目標(biāo)擬合差)Fig.8 Data RMS misfit versus GN iteration for marine CSEM survey (The red dashed line shows the desired data misfit at RMS=1)

5 結(jié)論

本文實現(xiàn)了一種頻率域CSEM三維反演算法.正演使用擬態(tài)有限體積法對控制方程離散化，并使用直接矩陣分解法求解離散所得的大型線性系統(tǒng)方程；反演中采用高斯牛頓最優(yōu)化算法，并使用預(yù)條件共軛梯度法(PCG)求解每次GN迭代所得到的法方程.以上這些技術(shù)的結(jié)合使用，使得本文的三維反演算法準(zhǔn)確、穩(wěn)定且高效，適用于陸地和海洋CSEM勘探場景，并且特別適合于多場源CSEM問題的模擬.理論模型的反演測試結(jié)果證明了本文算法的這些優(yōu)點.

圖9 海洋可控源三維鹽穹體模型不同剖面反演結(jié)果(a,b,c) XZ剖面； (d,e,f) YZ剖面； (g,h,i) XY剖面.黑色曲線為鹽穹的真實邊界位置.Fig.9 Inversion results shown at different cross sections and depth slices for marine CSEM surveyThe first row shows images at different XZ cross sections, the middle row for YZ sections, and the third row indicates horizontal slices at different depths. The black curve indicates true locations of the salt dome reservoir.

隨著電磁勘探精度要求的提高，對地電模型的假設(shè)將越來越接近真實情況，此外數(shù)據(jù)規(guī)模也將越來越大，例如最新的海洋CSEM發(fā)射與接收一體式的拖曳式測量會帶來數(shù)量極多的發(fā)射源和測點.這些因素都對CSEM的三維反演形成了巨大的挑戰(zhàn)，因此需要不斷地將最前沿的科學(xué)計算理論加以應(yīng)用以提升CSEM三維反演算法的性能.本文算是在這方面做了一個較小的嘗試，但還有大量的工作有待完成.

致謝 感謝加拿大英屬哥倫比亞大學(xué)(UBC)Eldad Haber教授的指導(dǎo)及與GIF組內(nèi)同學(xué)進行的有益討論.另外特別感謝MUMPS線性運算庫的開發(fā)者們的無私貢獻.向匿名審稿人對本文提出的寶貴修改意見表示衷心感謝.

附錄A Jv和JTw的計算

Helmholtz方程離散后得到的線性方程組可表示為

(A1)

其中A為稀疏正定對稱復(fù)矩陣，e為網(wǎng)格采樣點處電場分量，q為場源項.為獲得測點處的電場值，定義Nd×Ne階觀測矩陣P，假定其與模型參數(shù)無關(guān)，Nd、Ne分別為觀測數(shù)據(jù)個數(shù)和電場采樣點個數(shù).則觀測數(shù)據(jù)可表示為

(A2)

則數(shù)據(jù)的靈敏度矩陣J為

(A3)

Nm為模型參數(shù)個數(shù).

對(A1)式取全微分可得

(A4)

將(A4)式代入到(A3)式中可得

(A5)

為求得靈敏度矩陣J與模型向量v的乘積，定義向量u為

(A6)

令向量y滿足

Ay=u，

(A7)

結(jié)合(A5—A7)式，可得到

(A8)

為求得靈敏度矩陣轉(zhuǎn)置JT與模型向量w的乘積，令向量z滿足

ATz=PTw，

(A9)

結(jié)合(A5)和(A9)式可得

Amestoy P R, Duff I S, L′Excellent J Y, et al. 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling.SIAMJournalonMatrixAnalysisandApplications, 23(1): 15-41. Amestoy P R, Guermouche A, L′Excellent J Y, et al. 2006. Hybrid scheduling for the parallel solution of linear systems.ParallelComputing, 32(2): 136-156.

Aminzadeh F, Brac J, Kunz T. 1997. 3-D salt and overthrust models. SEG.

Barrett R, Berry M, Chan T F, et al. 1994. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. 2nd ed. Philadelphia: SIAM.

Commer M, Newman G A. 2008. New advances in three-dimensional controlled-source electromagnetic inversion.GeophysicalJournalInternational, 172(2): 513-535.

Constable S. 2010. Ten years of marine CSEM for hydrocarbon exploration.Geophysics, 75(5): 75A67-75A81.

Constable S, Weiss C J. 2006. Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling.Geophysics, 71(2): G43-G51.

Constable S C, Parker R L, Constable C G. 1987. Occam′s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data.Geophysics, 52(3): 289-300. Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems.GeophysicalJournalInternational, 189(1): 251-267, doi: 10.1111/j.1365-246X.2011.05347.x.Eisenstat S C, Walker H F. 1996. Choosing the forcing terms in an inexact Newton method.SIAMJournalonScientificComputing, 17(1): 16-32.Grayver A V, Streich R, Ritter O. 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver.GeophysicalJournalInternational, 193(3): 1432-1446. Grayver A V, Streich R, Ritter O. 2014. 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO2storage formation.Geophysics, 79(2): E101-E114.

Haber E, Oldenburg D W, Shekhtman R. 2007. Inversion of time domain three-dimensional electromagnetic data.GeophysicalJournalInternational, 171(2): 550-564.

Haber E, Ruthotto L. 2014. A multiscale finite volume method for Maxwell's equations at low frequencies.GeophysicalJournalInternational, 199(2): 1268-1277.

Han B. 2015 Three-dimensional parallel forward modeling and inversion of frequency-domain CSEM data[Ph.D. thesis](in Chinese). Wuhan:China University of Geosciences.

HAN Bo, HU Xiang-Yun, HUANG Yi-Fan et al .2015.3-D frequency-domain CSEM modeling using a parallel direct solver.ChineseJ.Geophys.,58(8): 2812-2826,doi: 10.6038/cjg20150816.

Hu X Y, Peng R H, Wu G J, et al. 2013. Mineral exploration using CSAMT data: Application to Longmen region metallogenic belt, Guangdong Province, China.Geophysics, 78(3): B111-B119. Hyman J, Shashkov M, Steinberg S. 1997. The numerical solution of diffusion problems in strongly heterogeneous non-isotropic materials.JournalofComputationalPhysics, 132(1): 130-148.

Hyman J M, Shashkov M. 1997. Adjoint operators for the natural discretizations of the divergence, gradient and curl on logically rectangular grids.AppliedNumericalMathematics, 25(4): 413-442.

Hyman J M, Shashkov M. 1999a. Mimetic discretizations for Maxwell′s equations.JournalofComputationalPhysics, 151(2): 881-909.

Hyman J M, Shashkov M. 1999b. The orthogonal decomposition theorems for mimetic finite difference methods.SIAMJournalonNumericalAnalysis, 36(3): 788-818.

Ledo J. 2005. 2-D versus 3-D magnetotelluric data interpretation.SurveysinGeophysics, 26(5): 511-543.

Lelièvre P G, Oldenburg D W. 2009. A comprehensive study of including structural orientation information in geophysical inversions.GeophysicalJournalInternational, 178(2): 623-637. Lin C H, Tan H D, Shu Q, et al. 2012. Three-dimensional conjugate gradient inversion of CSAMT data.ChineseJ.Geophys. (in Chinese), 55(11): 3829-3838, doi: 10.6038/j.issn.0001-5733.2012.11.030.

Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data.ChineseJ.Geophys. (in Chinese), 56(12): 4278-4287, doi: 10.6038/cjg20131230.

Newman G A, Alumbaugh D L. 1997. Three-dimensional massively parallel electromagnetic inversion—I. Theory.GeophysicalJournalInternational, 128(2): 345-354.

Nocedal J, Wright S. 1999. Numerical optimization. New York: Springer.

Oldenburg D W, Haber E, Shekhtman R. 2013. Three dimensional inversion of multisource time domain electromagnetic data.Geophysics, 78(1): E47-E57.

Pidlisecky A, Haber E, Knight R. 2007. RESINVM3D: A 3D resistivity inversion package.Geophysics, 72(2): H1-H10.

Plessix R é, Mulder W A. 2008. Resistivity imaging with controlled-source electromagnetic data: depth and data weighting.InverseProblems, 24(3): 034012. Rieder A. 2001. On convergence rates of inexact Newton regularizations.NumerischeMathematik, 88(2): 347-365.

Rieder A. 2005. Inexact newton regularization using conjugate gradients as inner iteration.SIAMjournalonnumericalanalysis, 43(2): 604-622. Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion.Geophysics, 66(1): 174-187. Saad Y. 2003. Iterative methods for sparse linear systems. Philadelphia: SIAM.

Schenk O, G?rtner K. 2006. On fast factorization pivoting methods for sparse symmetric indefinite systems.ElectronicTransactionsonNumericalAnalysis, 23(1): 158-179. Schwarzbach C, Haber E. 2013. Finite element based inversion for time-harmonic electromagnetic problems.GeophysicalJournalInternational, 193(2): 615-634, doi: 10.1093/gji/ggt006. Siripunvaraporn W. 2012. Three-dimensional magnetotelluric inversion: an introductory guide for developers and users.Surv.Geophys., 33(1): 5-27.

Siripunvaraporn W, Egbert G, Uyeshima M. 2005. Interpretation of two-dimensional magnetotelluric profile data with three-dimensional inversion: synthetic examples.GeophysicalJournalInternational, 160(3): 804-814. Smith J T. 1996. Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator.Geophysics, 61(5): 1319-1324.

Streich R. 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy.Geophysics, 74(5): F95-F105.

Tikhonov A N, Arsenin V I A K. 1977. Solutions of Ill-Posed Problems. Washington, D. C., USA: W. H. Winston and Sons. van den Doel K, Ascher U M. 2007. Dynamic level set regularization for large distributed parameter estimation problems.InverseProblems, 23(3): 1271-1288. van den Doel K, Ascher U M. 2012. Adaptive and stochastic algorithms for electrical impedance tomography and DC resistivity problems with piecewise constant solutions and many measurements.SIAMJournalonScientificComputing, 34(1): A185-A205.

Weng A H, Li D J, Li Y B, et al. 2015. Selection of parameter types in Controlled Source Electromagnetic Method.ChineseJ.Geophys. (in Chinese), 58(2): 697-708, doi: 10.6038/cjg20150230.

Weng A H, Liu Y H, Jia D Y, et al. 2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients.ChineseJ.Geophys. (in Chinese), 55(10): 3506-3515, doi: 10.6038/j.issn.0001-5733.2012.10.034.

Yang D K, Oldenburg D W. 2012. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit.Geophysics, 77(2): B23-B34.

Zhao N, Wang X B, Qin C, et al. 2016. 3D frequency-domain CSEM inversion.ChineseJ.Geophys. (in Chinese), 59(1): 330-341, doi: 10.6038/cjg20160128.

附中文參考文獻

韓波. 2015. 頻率域可控源電磁法并行化三維正反演算法研究[博士論文]. 武漢: 中國地質(zhì)大學(xué)(武漢).

韓波, 胡祥云, 黃一凡等. 2015. 基于并行化直接解法的頻率域可控源電磁三維正演. 地球物理學(xué)報, 58(8): 2812-2826, doi: 10.6038/cjg20150816.

林昌洪, 譚捍東, 舒晴等. 2012. 可控源音頻大地電磁三維共軛梯度反演研究. 地球物理學(xué)報, 55(11): 3829-3838, doi: 10.6038/j.issn.0001-5733.2012.11.030.

劉云鶴, 殷長春. 2013. 三維頻率域航空電磁反演研究. 地球物理學(xué)報, 56(12): 4278-4287, doi: 10.6038/cjg20131230.

翁愛華, 李大俊, 李亞彬等. 2015. 數(shù)據(jù)類型對三維地面可控源電磁勘探效果的影響. 地球物理學(xué)報, 58(2): 697-708, doi: 10.6038/cjg20150230.

翁愛華, 劉云鶴, 賈定宇等. 2012. 地面可控源頻率測深三維非線性共軛梯度反演. 地球物理學(xué)報, 55(10): 3506-3515, doi: 10.6038/j.issn.0001-5733.2012.10.034.

趙寧, 王緒本, 秦策等. 2016. 三維頻率域可控源電磁反演研究. 地球物理學(xué)報59(1): 330-341, doi: 10.6038/cjg20160128.

(本文編輯 張正峰)

3D inversion of frequency-domain CSEM data based on Gauss-Newton optimization

PENG Rong-Hua1,2, HU Xiang-Yun1*, HAN Bo3

1InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China2DepartmentofEarth,OceanandAtmosphericSciences,UniversityofBritishColumbia,VancouverV6T1Z4,Canada3CollegeofMarineGeosciences,OceanUniversityofChina,Qingdao266100,China

The Controlled-source electromagnetic (CSEM) method in the frequency domain has evolved into an established technique for mineral resources prospecting and hydrocarbon exploration. An increasing trend is to carry out three-dimensional (3D) CSEM surveys as EM explorations now are increasingly conducted in complex geological environments in order to improve the spatial resolution of subsurface conductivity structure.

Frequency-domain CSEM； 3D inversion； Gauss-Newton； Direct solver； Preconditioned conjugate gradient

10.6038/cjg20160929.

國家自然科學(xué)基金(41274077，41474055)、國家重點基礎(chǔ)研究發(fā)展計劃項目(2013CB733200)、國家留學(xué)基金委(201406410020)和湖北省自然科學(xué)基金(2015CFA019)聯(lián)合資助.

彭榮華，男，1988年生，博士研究生，研究方向為電磁法三維正演與反演模擬.E-mail:prhjiajie@gmail.com

*通訊作者 胡祥云，男，1966年生，教授，博士生導(dǎo)師，主要從事電磁法的理論與應(yīng)用研究.E-mail：xyhu@cug.edu.cn

10.6038/cjg20160929

P631

2015-12-06，2016-04-21收修定稿

彭榮華， 胡祥云， 韓波. 2016. 基于高斯牛頓法的頻率域可控源電磁三維反演研究. 地球物理學(xué)報，59(9):3470-3481,

Peng R H, Hu X Y, Han B. 2016. 3D inversion of frequency-domain CSEM data based on Gauss-Newton optimization.ChineseJ.Geophys. (in Chinese),59(9):3470-3481,doi:10.6038/cjg20160929.