基于非結(jié)構(gòu)三角網(wǎng)格的海洋CSEM和MT二維聯(lián)合反演研究

艾正敏，葉益信，湯文武，陳曉，杜家明

(1.東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西 南昌 330013; 2.東華理工大學(xué) 地球物理與測控技術(shù)學(xué)院，江西 南昌 330013)

0 引言

海洋可控源電磁法(marine controlled-source electromagnetic method, MCSEM)近些年被廣泛應(yīng)用于探測海底油氣與研究海洋深部構(gòu)造等[1-4]。海洋電磁法的數(shù)值模擬也隨之得到了很大的重視，目前大多以二維和三維電磁場的數(shù)值模擬為主，而傳統(tǒng)的電磁場正演模擬大多使用的是規(guī)則網(wǎng)格，一些學(xué)者使用非結(jié)構(gòu)三角網(wǎng)格進(jìn)行了有限元正演模擬[5-10]，這為非結(jié)構(gòu)三角網(wǎng)格在地球物理電磁場計算中的廣泛應(yīng)用奠定了基礎(chǔ)。電磁場的數(shù)值模擬通常會產(chǎn)生一些大型的稀疏矩陣，而大規(guī)模的反演問題通常是欠定的，這意味著要用有限的參數(shù)求取更多的未知數(shù)，這也是反演解的非唯一性問題的原因。為了減輕這一問題并提高反演的分辨率，通常利用不同地球物理數(shù)據(jù)集的互補(bǔ)性質(zhì)。

在電磁探測方面，CSEM對淺層異常體較為敏感，而MT對深部低阻異常體和大范圍電性結(jié)構(gòu)的探測更為有效，二者呈現(xiàn)出互補(bǔ)特性。隨著高精度探測的需求，使用CSEM和MT合成數(shù)據(jù)聯(lián)合反演的結(jié)果進(jìn)行海底探測成像，將會提高成像結(jié)果的分辨率。即使提高了探測精度，反演解的非唯一性問題仍舊存在。有研究者提出在地震和MT聯(lián)合反演中加入有效的模型約束，使反演的解更實(shí)際、更穩(wěn)定[11]，或者是在地震和MT聯(lián)合反演中加入有效的巖石物性約束，結(jié)合模擬退火模型將電阻率和深度耦合在一起，提高聯(lián)合反演反演出真實(shí)物性參數(shù)的可能性[12]。上述方法考慮的是地震數(shù)據(jù)和MT數(shù)據(jù)的聯(lián)合反演，而關(guān)于CSEM和MT數(shù)據(jù)的聯(lián)合反演更多的是使用數(shù)據(jù)加權(quán)的方式，通過改變CSEM或MT數(shù)據(jù)的權(quán)重來得到更優(yōu)的反演結(jié)果，提高反演解的分辨率[13]。隨著高精度探測和反演解釋的需要，不同方法的聯(lián)合反演將成為地球物理定量解釋的重要工具，是未來地球物理學(xué)的一個發(fā)展方向[14-15]。

本文使用合成的CSEM和MT數(shù)據(jù)進(jìn)行二維聯(lián)合反演的模型計算與分析，正演部分采用自適應(yīng)非結(jié)構(gòu)有限元正演模擬，非結(jié)構(gòu)三角網(wǎng)格能夠更好地模擬起伏地形和真實(shí)的地質(zhì)構(gòu)造，同時能夠提高有限元解的精度。反演采用快速Occam算法[16-17]，通過對聯(lián)合反演數(shù)據(jù)權(quán)重的公式進(jìn)行推導(dǎo)，構(gòu)建CSEM與MT數(shù)據(jù)加權(quán)系數(shù)，控制不同數(shù)據(jù)集的權(quán)重，提高反演對整體數(shù)據(jù)的擬合精度，從而得到擬合度更高的反演結(jié)果。雖然控制不同數(shù)據(jù)集的權(quán)重可以得到擬合度更高的模型，但這并不意味著擬合度更高的模型就一定和真實(shí)模型更接近，因此，本文使用不同的CSEM和MT數(shù)據(jù)的權(quán)重比值來進(jìn)行聯(lián)合反演，對不同權(quán)重反演得到的結(jié)果進(jìn)行分析與討論，希望找到最優(yōu)的權(quán)重因子。

1 方法理論

1.1 自適應(yīng)有限元方法

地球物理電磁場的二維模擬問題通常采用有限單元法(finite element method，F(xiàn)EM)[18]，而應(yīng)用自適應(yīng)有限元方法進(jìn)行電磁場數(shù)值模擬，能夠自動細(xì)化網(wǎng)格并在不顯著增加計算時間的條件下提供可靠的計算結(jié)果，是一種能夠自動調(diào)整算法以改進(jìn)求解過程的數(shù)值方法[5-6,19]。自適應(yīng)有限元法主要包括偏微分方程求解、誤差估計、網(wǎng)格標(biāo)記和網(wǎng)格優(yōu)化4個步驟，而其中研究重點(diǎn)是誤差估計和自適應(yīng)網(wǎng)格的優(yōu)化。本文采用梯度恢復(fù)型后驗(yàn)誤差估計[20]，即從初始網(wǎng)格出發(fā)，通過有限元解的后驗(yàn)誤差估計因子來控制后續(xù)模型網(wǎng)格的剖分細(xì)化。對于某一個有限元網(wǎng)格單元的后驗(yàn)誤差估計可表示為：

(1)

(2)

1.2 反演算法

本文海洋CSEM與MT數(shù)據(jù)聯(lián)合反演方法采用的是快速Occam算法[16-17]，該算法實(shí)質(zhì)為帶光滑約束的最小二乘法，其泛函U表達(dá)式為：

U=μ(P‖m-m*‖2+‖Rm‖2)+

‖W(d-F(m))‖2。

(3)

式中：m為N維的模型參數(shù)矢量，一般為電阻率值；P為正則項對應(yīng)的權(quán)重系數(shù)矩陣；R為粗糙度算子矩陣；μ為拉格朗日乘子，用以平衡模型粗糙度和數(shù)據(jù)擬合誤差，當(dāng)μ取較大值時，反演以搜索光滑模型為主，反之則以搜索最小擬合誤差為主；W為關(guān)聯(lián)的對角加權(quán)矩陣；d為觀測數(shù)據(jù)矢量；F(m)為模型m對應(yīng)的正演響應(yīng)。

關(guān)于給定初始模型mk的函數(shù)，通過以下迭代方法實(shí)現(xiàn)目標(biāo)函數(shù)的最小化：

(4)

(5)

通過上述模型迭代計算公式，將上一次的計算結(jié)果作為下一次計算的初始模型，μ在第一次迭代中給出初始值，之后每次迭代都按照線性搜索方式計算得出。如果這種迭代算法最后收斂，那就得到了目標(biāo)函數(shù)U最小化的解，并且與初始模型mk不相關(guān)，最后的結(jié)果將是在滿足數(shù)據(jù)擬合誤差前提下最平緩的模型。

1.3 聯(lián)合反演數(shù)據(jù)權(quán)重

MT和CSEM數(shù)據(jù)的聯(lián)合反演可以改進(jìn)海洋勘探的結(jié)構(gòu)成像[13,22]。通過聯(lián)合反演彌補(bǔ)單獨(dú)反演勘探深度或勘探精度的不足，進(jìn)而得到比單獨(dú)反演效果更好的結(jié)果。然而，聯(lián)合反演存在的一個問題是當(dāng)每個子集的數(shù)據(jù)量存在很大差異時，怎么樣才能相同地擬合每個數(shù)據(jù)子集。本文通過對不同子集擬合誤差的加權(quán)來控制聯(lián)合反演的權(quán)重，從而得到更好的反演效果。例如，將一個n維數(shù)據(jù)劃分為具有維度n1和n2的MT和CSEM子集，對于n維數(shù)據(jù)向量d，有：

d=[d1d2]T，

(6)

然后可以將目標(biāo)函數(shù)中的誤差函數(shù)擴(kuò)展為：

‖W(d-F(m))‖2=‖W1(d1-F1(m))‖2+

‖W2(d2-F2(m))‖2，

(7)

(8)

‖α1W1(d1-F1(m))‖2+‖α2W2(d2-F2(m))‖2

(9)

(10)

根據(jù)上面的公式，可以進(jìn)一步擴(kuò)展到N個子集的情況。具有N個子集的數(shù)據(jù)向量可以表示為：

d=[d1d2…dN]T。

(11)

這種情況下的誤差函數(shù)或者說是擬合差的均方根值(RMS)的平方可以表示為：

(12)

雖然CSEM與MT子集的數(shù)據(jù)量和數(shù)據(jù)類型是不一樣的，CSEM子集包含電磁場的振幅和相位，MT子集包含TE、TM模式的視電阻率和相位，但是這種方法對將CSEM和MT數(shù)據(jù)集放入一個數(shù)據(jù)集中進(jìn)行反演仍將是有用的,同時，該方法也可用于大地電磁測深數(shù)據(jù)和地震數(shù)據(jù)的聯(lián)合反演[23-24]。

1.4 聯(lián)合反演流程

本文采用的海洋可控源電磁和大地電磁數(shù)據(jù)聯(lián)合反演計算流程見圖1。

圖1 聯(lián)合反演流程示意Fig.1 Schematic diagram of the joint inversion process

2 模型計算

2.1 簡單模型反演計算

為驗(yàn)證反演算法的可靠性，本文設(shè)置了一個二維水平海底模型并進(jìn)行反演計算。模型如圖2所示，在海水層厚為500 m的海底設(shè)置了兩個異常體，淺部異常體A的電阻率為500 Ω·m，大小為1 000 m×500 m，深部異常體B的電阻率為0.5 Ω·m，大小為2 000 m×1 000 m。MT與CSEM采用相同的接收點(diǎn)，從左往右依次設(shè)置了20個測點(diǎn)，間距500 m，MT的頻率范圍為0.001～1 Hz，取對數(shù)等間隔分布選取了20個頻點(diǎn)。CSEM發(fā)射源均勻分布在-5 000～5 000 m范圍內(nèi)，共20個，位于海底上方50 m的海水中，發(fā)射源頻率分別為1 Hz、5 Hz。

圖2 二維水平海底模型Fig.2 Two-dimensional model with flat seafloor

針對上述模型分別進(jìn)行CSEM與MT正演，網(wǎng)格采用算法中的自適應(yīng)網(wǎng)格剖分。以CSEM正演網(wǎng)格為例，如圖3所示，模型正演計算的初次細(xì)化后的網(wǎng)格剖分為由1 445個節(jié)點(diǎn)形成的2 849個三角單元，第12次細(xì)化后的最終網(wǎng)格剖分為由43 127個節(jié)點(diǎn)形成的86 139個三角單元?？梢钥闯?，由于測點(diǎn)附近的網(wǎng)格剖分對測點(diǎn)接收到的電磁場分量影響較大，所以在自適應(yīng)網(wǎng)格優(yōu)化過程中，對淺地表和測點(diǎn)附近的網(wǎng)格進(jìn)行了加密，而邊界和深部的網(wǎng)格并沒有進(jìn)行加密。

a—初次優(yōu)化網(wǎng)格，1 445個節(jié)點(diǎn)，2 849個三角單元; b—第12次自適應(yīng)細(xì)化后的最終網(wǎng)格部分，43 127個節(jié)點(diǎn)，86 139個三角單元; 示意圖只展示了模型中間的異常區(qū)域a—the first optimized mesh containing 1 445 nodes and 2 849 triangular elements; b—the final mesh after the 12th adaptive refinement, containing 43 127 nodes and 86 139 triangular elements; the diagram only shows the abnormal area in the middle part of the model圖3 自適應(yīng)正演網(wǎng)格剖分示意(以CSEM網(wǎng)格為例)Fig.3 Schematic diagram of adaptive forward meshing (taking CSEM mesh as an example)

對模型進(jìn)行有限元正演計算，CSEM正演得到的是Ey的幅值和相位的響應(yīng)值，MT正演得到的是TE、TM模式的視電阻率和相位的響應(yīng)值，對得到的響應(yīng)值添加4%的隨機(jī)高斯噪聲，則生成了反演計算所需的觀測數(shù)據(jù)。反演的初始模型包含空氣、海水和海底地層，其中空氣和海水電阻率不參與反演，給定地層的初始電阻率為1.0 Ω·m。計算區(qū)域所劃分的單元越多，所需要的計算量越大，因此只對目標(biāo)反演區(qū)域采用較細(xì)的網(wǎng)格剖分，其余部分則采用粗網(wǎng)格以減少不必要的計算量。如圖4所示，選取y范圍為-4.5～4.5 km海底到5 km深的區(qū)域，調(diào)用Triangle程序[25]對該區(qū)域進(jìn)行精細(xì)網(wǎng)格剖分，共生成3 756個三角單元，其余部分共剖分為440個三角形單元。

圖4 初始模型的反演網(wǎng)格剖分Fig.4 Inversion mesh division of the initial model

設(shè)定目標(biāo)均方根誤差為1.0，輸入初始模型，調(diào)用程序?qū)δＰ蛿?shù)據(jù)進(jìn)行反演計算。在個人電腦上，MT數(shù)據(jù)單獨(dú)反演經(jīng)過20次迭代得到圖5a所示結(jié)果，CSEM數(shù)據(jù)單獨(dú)反演經(jīng)過15次迭代得到圖5b所示結(jié)果，CSEM與MT數(shù)據(jù)聯(lián)合反演經(jīng)過21次迭代最終得到圖5c所示結(jié)果。

a—MT反演; b—CSEM反演； c—CSEM與MT聯(lián)合反演a—MT data inversion result; b—CSEM data inversion result; c—CSEM+MT data joint inversion result圖5 簡單模型反演結(jié)果Fig.5 Inversion results of the simple two-dimensional model

對比3個反演的結(jié)果可以看出，MT對深部異常體的反演效果較好，對淺部異常體無法進(jìn)行識別，CSEM對淺部的反演效果非常好，對深部的反演則存在一些偏差，對淺部異常體邊界和埋深的反映與模型基本一致。整體來看，聯(lián)合反演對淺部和深部的反演效果都很不錯，但淺部的反演出的電阻率沒有CSEM反演出的視電阻率更接近模型的真電阻率，盡管如此，聯(lián)合反演在效果上依舊比MT或CSEM反演更具有優(yōu)勢。

圖6為聯(lián)合反演得到的最終模型MT響應(yīng)與MT觀測數(shù)據(jù)的擬斷面對比，圖7為發(fā)射源在y=0處CSEM觀測數(shù)據(jù)與聯(lián)合反演的CSEM模型響應(yīng)的振幅和相位擬合結(jié)果?？梢钥闯?，模擬數(shù)據(jù)與真實(shí)數(shù)據(jù)的視電阻率和相位信息高度吻合，驗(yàn)證了反演結(jié)果的可靠性和準(zhǔn)確度，說明本文采用的聯(lián)合反演方法能有效識別不同深度的異常體。

圖6 MT觀測數(shù)據(jù)(上)與聯(lián)合反演正演響應(yīng)(下)的擬斷面圖對比Fig.6 Comparison of pseudo-section diagrams of MT observation data (top) and MT forward responses of joint inversion model (bottom)

圖7 發(fā)射源在y=0處CSEM觀測數(shù)據(jù)與聯(lián)合反演模型響應(yīng)的振幅和相位擬合Fig.7 Amplitude and phase fitting diagram of the CSEM observation data and model response of the joint inversion result with the emission source at y=0

表1給出了上述3種方法反演的最終耗時、均方根擬合差RMS及粗糙度。從表中可以看出，MT的反演耗時為157.7 min，CSEM的反演耗時為368.5 min，聯(lián)合反演的耗時為1 247.3 min，聯(lián)合反演時數(shù)據(jù)量增加了一倍，導(dǎo)致反演時間大大增加，在以后的研究中應(yīng)該考慮對聯(lián)合反演數(shù)據(jù)集的選取研究和進(jìn)行算法的優(yōu)化。

表1 反演耗時、擬合差、粗糙度及迭代次數(shù)統(tǒng)計

2.2 簡單模型不同權(quán)重聯(lián)合反演計算

假設(shè)權(quán)重因子q=αMT/αCSEM=α1/α2，圖2模型聯(lián)合反演中使用的權(quán)重因子為1，使用權(quán)重因子為1的好處是可以在CSEM和MT數(shù)據(jù)量相差不大時更好地擬合2個數(shù)據(jù)子集，但并不是在所有的情況下都選擇1為權(quán)重因子的。在一些數(shù)據(jù)量差異較大時，就需要找到最優(yōu)的權(quán)重因子并用于聯(lián)合反演。因此，本文使用不同的權(quán)重因子對圖2中的模型再次進(jìn)行聯(lián)合反演。

圖8是使用不同權(quán)重因子的聯(lián)合反演結(jié)果示意，可以看出：q值越小，即CSEM數(shù)據(jù)權(quán)重系數(shù)越大時，聯(lián)合反演對淺部高阻異常體探測效果更好，而當(dāng)q值越大，即MT數(shù)據(jù)權(quán)重系數(shù)更大時，聯(lián)合反演對淺部異常體的探測能力變?nèi)趿?，對深部異常體的探測效果并沒有很大改變。所以，可認(rèn)為權(quán)重因子q對聯(lián)合反演中的CSEM數(shù)據(jù)影響更大。

a—q=2；b—q=0.5；c—q=10；d—q=0.1圖8 不同權(quán)重因子q值時的聯(lián)合反演結(jié)果Fig.8 Joint inversion results with different q value

圖9為不同q值時聯(lián)合反演最終所得模型的RMS值，可以看出，權(quán)重因子q=0.1時的擬合差要明顯低于q=10的擬合差，而且聯(lián)合反演對MT數(shù)據(jù)的擬合差并沒有隨著q值的變大而變小，反而是隨著q值的變大而變大，這也是因?yàn)镸T權(quán)重太大導(dǎo)致的，在聯(lián)合反演時需要充分?jǐn)M合2個數(shù)據(jù)集，而不是過于偏向某個數(shù)據(jù)集。

圖9 不同q值時聯(lián)合反演最終RMS值Fig.9 The final RMS value of joint inversion with different q value

2.3 復(fù)雜模型反演計算

為進(jìn)一步驗(yàn)證反演算法的準(zhǔn)確性和可靠性，設(shè)計了1個海底二維帶地形的復(fù)雜模型。如圖10所示，海水電阻率為0.3 Ω·m，海底下方設(shè)置了4個異常體，Ⅰ號異常體是一個中間層的侵入體，電阻率為20 Ω·m；Ⅱ號與Ⅳ號異常體在2個不同的地層中，電阻率分別為500 Ω·m和0.2 Ω·m；Ⅲ號異常體位于地層淺部，電阻率為500 Ω·m。

圖10 二維復(fù)雜海洋模型Fig.10 Two-dimensional complex ocean model

在海底設(shè)置了25個觀測裝置(圖10中白色三角符號所示)，測點(diǎn)在1～25 km范圍沿海底地形均勻的分布。MT和CSEM的接收點(diǎn)位置相同，CSEM發(fā)射源位于海底上方50 m的海水中，范圍在0～25 km內(nèi)均勻分布共25個。MT觀測頻率范圍為10-4～0.1 Hz，取對數(shù)等間隔分布，共20個，CSEM發(fā)射頻率分別為0.1、1、10 Hz。對該模型進(jìn)行有限元正演計算，得到其視電阻率、振幅和相位響應(yīng)值，對響應(yīng)值添加4%高斯隨機(jī)噪聲，得到反演所需的觀測數(shù)據(jù)。反演初始模型包括空氣、海水和海底地層，其中空氣和海水電阻率為固定值不參與反演，給定地層的初始電阻率為1.0 Ω·m。

對初始模型進(jìn)行剖分，如圖11所示，根據(jù)接收點(diǎn)和發(fā)射點(diǎn)的位置分布對異常體或構(gòu)造可能存在的區(qū)域進(jìn)行精細(xì)網(wǎng)格剖分，其余部分則采用不加約束的粗網(wǎng)格剖分。精細(xì)網(wǎng)格的剖分可以保證計算的精度，采用粗網(wǎng)格則可以減少不必要的計算量。調(diào)用Triangle程序[25]對反演區(qū)域剖分，最終剖分得到 6 089 個三角單元。

圖11 初始模型的反演網(wǎng)格剖分Fig.11 Inversion mesh division of the initial model

設(shè)定目標(biāo)均方根誤差為1.0，輸入初始模型。MT觀測數(shù)據(jù)2 000個，CSEM觀測數(shù)據(jù)1 842個，總的觀測數(shù)據(jù)為3 842個，其中數(shù)據(jù)量比值NCSEM/NMT=0.921,對于MT與CSEM數(shù)據(jù)量比值接近1的情況，采用權(quán)重因子q=1，以便聯(lián)合反演能夠同等擬合2個觀測數(shù)據(jù)集，提高反演的分辨率。反演結(jié)果見圖12。

a—MT反演; b—CSEM反演； c—CSEM與MT聯(lián)合反演； d—誤差變化a—MT data inversion result; b—CSEM data inversion result; c—CSEM+MT data joint inversion result; d—RMS diagram圖12 反演結(jié)果Fig.12 Inversion result

圖12a是MT數(shù)據(jù)經(jīng)迭代17次的反演結(jié)果，圖12b是CSEM數(shù)據(jù)迭代24次反演結(jié)果，圖12c是CSEM+MT數(shù)據(jù)經(jīng)迭代18次的聯(lián)合反演結(jié)果。整體來看，MT反演能夠恢復(fù)較大的低阻異常體電阻率，但不能恢復(fù)較小和較大的高阻異常體，雖然能夠識別深部基巖，但基巖起伏面不夠明顯。CSEM反演結(jié)果對淺部高阻異常體和深部基巖電阻率的恢復(fù)比MT好一點(diǎn)，同時更好區(qū)分淺部的兩個不同地層，但對低阻異常體的電阻率和幾何形態(tài)恢復(fù)得較差。相對前兩種反演方法，聯(lián)合反演對異常體的位置反演得十分準(zhǔn)確，也沒有產(chǎn)生過多虛假異常，對淺部地層分界面能夠進(jìn)行較好的區(qū)分，同時基巖起伏界面也比較明確，與真實(shí)模型十分接近，這也說明該算法可用于復(fù)雜模型的聯(lián)合反演，進(jìn)一步證明了該算法的有效性。圖12d為反演過程中RMS隨迭代次數(shù)的變化曲線，如圖所示，反演經(jīng)過多次迭代，最終RMS都下降到1左右。

3 結(jié)論

本文探討和分析了基于非結(jié)構(gòu)三角網(wǎng)格的海洋CSEM和MT數(shù)據(jù)的聯(lián)合反演。從目標(biāo)函數(shù)出發(fā)，采用數(shù)據(jù)加權(quán)的方式以達(dá)到對CSEM和MT數(shù)據(jù)子集的迭代擬合；然后對模型進(jìn)行了正反演試算，發(fā)現(xiàn)聯(lián)合反演的效果比單一反演效果更好。聯(lián)合反演結(jié)合了CSEM與MT方法的各自優(yōu)勢，使用不同的權(quán)重因子q進(jìn)行聯(lián)合反演，能夠有效識別不同深度的異常體與構(gòu)造，且對地層電阻率恢復(fù)度較高，與真實(shí)模型更加接近，同時，反演生成模型的正演響應(yīng)與觀測數(shù)據(jù)非常吻合，反演迭代收斂穩(wěn)定，結(jié)果真實(shí)可靠，證明了本文算法的有效性。

本文在數(shù)據(jù)加權(quán)上只是單純地對不同數(shù)據(jù)集進(jìn)行加權(quán)，并沒有得到加權(quán)系數(shù)與探測深度和探測精度之間的關(guān)系，而且聯(lián)合反演的耗時較大，因此需要進(jìn)行聯(lián)合反演數(shù)據(jù)集的選取和算法優(yōu)化研究。目前尚缺少實(shí)測數(shù)據(jù)的聯(lián)合反演實(shí)例。對于后續(xù)實(shí)測數(shù)據(jù)的聯(lián)合反演過程，建議首先需要了解不同方法數(shù)據(jù)的探測深度和分辨率，通過幾次單獨(dú)的反演嘗試了解，然后使用適當(dāng)?shù)臄?shù)據(jù)權(quán)重對合成的CSEM和MT數(shù)據(jù)進(jìn)行聯(lián)合反演，數(shù)據(jù)權(quán)重可以根據(jù)數(shù)據(jù)量的比值適當(dāng)選取。