基于分區并行的2.5維海洋可控源電磁法反演研究

張迎利， 李 賡

(河南理工大學 物理與電子信息學院，焦作 454000)

0 引言

海洋可控源電磁法(Marine Controlled Source Electromagnetic method, MCSEM)廣泛地應用于海底油氣資源和天然氣水合物儲層探測研究[1-2]。隨著能源需求的不斷增長和陸上油氣資源儲備量的日益減少，海洋油氣資源探測逐漸成為能源爭奪的焦點。近些年來，MCSEM被廣泛應用于海底油氣資源探測，并已經取得了明顯的效果[3-4]。海洋電磁正演策略分為二次場和總場算法[5-6]，為克服場源奇異性問題，二次場算法需要利用一維解析公式計算背景場，而總場算法需要計算源項[7]。一般來說，對于簡單模型而言，二次場算法的計算精度高，但對于復雜模型，由于二次場算法很難確定背景模型，因此總場算法的適應性更強[8-9]。

正演通過使用非結構化網格，可以很好地適應起伏地形和地下復雜電性結構[10]。但是在大規模電磁數值模擬中，通過增加適當的網格來適應地下電性結構是提高計算效率的關鍵。近些年，自適應有限元方法針對以上問題給出了很好的解決方案。其中二維電磁場數值模擬中基于非結構化網格的自適應有限元法已經取得了很好的應用效果[11-13]。

海洋可控源電磁反演研究中，Abubakar等[14]基于高斯-牛頓方法提出了一種快速、有效的2.5維反演算法;Chen等[15]采用非線性共軛梯度算法實現了MCSEM 2.5維反演;李剛等[16]采用有限差分法實現了2.5維頻率域MCSEM正演，并通過直接法求解方程組;趙東東等[17]基于有限元方法對2.5維MCSEM進行了正演模擬，引入網格加密收縮算法提升了數值的計算效率;陳光源等[18]利用非線性共軛梯度法實現了2.5維海洋可控源電磁反演,并探討了不同反演參數對計算速度的影響;Key等[12]實現了頻率、發射源并行的2.5維自適應有限元OCCAM反演算法，提高了反演計算速度，但仍需要計算、存儲整個反演域上的靈敏度矩陣。基于光滑模型約束的OCCAM反演方法，反演過程穩定，但結果不能清晰刻畫地下電性結構的邊界[19]。為了解決光滑模型約束的問題，Portniaguine等[20]提出了聚焦反演成像方法，其方法的優點在于能夠得到較為清晰地地電結構圖像，但反演過程不穩定，依賴初始電性模型的選擇。在MCSEM實際資料解釋過程中，反演計算量大、內存需求量高[21]。

筆者運用二維海洋可控源電磁資料處理，正演采用基于總場的自適應有限元方法，其中網格剖分采用非結構化三角單元。為提高計算速度，降低內存使用量，我們將反演區域進行分解計算。反演過程中利用地震資料和光滑模型進行約束反演，得到背景電阻率，然后將其作為初始模型進行模型空間約束反演，從而獲得工區內高阻區域分布。

1 正演模擬

1.1 控制方程

在二維海洋地電模型中，設走向方向為x、y方向軸水平向右，z方向軸垂直向下，假設時諧因子為e-iωt，則電場和磁場滿足的控制方程為[11]

(1)

式中:E為電場強度；ω為角頻率；μ0為真空中的磁導率；H為磁場強度；σ為介質電導率；J為外加源電流密度。

將式(1)中電磁場分量沿著構造走向(x方向)做傅里葉變換[22]

(2)

(3)

(4)

(5)

(6)

(7)

(8)

1.2 自適應有限元

對于簡單模型，依據經驗對網格進行剖分，可以得到精度較高的數值解；而對于復雜模型來說，需要把有限元網格剖分得足夠細，但網格過細會影響到計算效率。在自適應有限單元法中，通過后驗誤差估計方法獲得每個網格單元上的誤差分布，選擇性的加密誤差較大的部分單元，可以快速提高解的精度，同時避免全局網格加密所帶來的計算時間和內存使用量的急劇增加。

(9)

式(9)的離散形式表示為

B(un,v)=F(v)

(10)

其中，un為有限元數值解，u為精確解，令全局誤估計算子εn≈u-un。由于精確解u未知，因此εn不能直接求解，但是可以通過近似誤差函數計算。

B(εn,v)=B(u-un)=F(v)-B(un,v)

(11)

定義誤差函數[12]

(12)

圖1 自適應有限元正演算法流程圖Fig.1 Adaptive finite element forward algorithm chart

圖2 1D模型示意圖Fig.2 1D model

1.3 正確性驗證

為驗證正演算法正確性，設計如圖2所示海洋模型，使用本文正演算法計算其正演響應，將計算結果與解析解[24]進行對比。其中，海水層厚度為1 km，電阻率為0.3 Ω·m，在海底下方1 km處有一厚度為0.1 km的高阻儲層，電阻率為100 Ω·m；發射源置于海底上方0.05 km處，坐標為(0 km,0 km,0.95 km)，發射電流為1 A，發射頻率為0.3 Hz；30個海底接收點等間距分布于y=0.5 km～15 km的范圍內。將本文正演計算的電磁響應(第六次自適應加密)與解析解進行對比，結果如圖3(a)所示，可以看出數值解與解析解吻合較好。圖3(b)為相對誤差曲線，可以看到，在自適應加密初期相對誤差很大，隨著自適應次數增加，相對誤差逐漸下降到1%以下，表明本文的正演方法具有較高的計算精度。

2 基于背景模型約束的OCCAM反演

反演過程中，如果知道某區域地球物理場的背景信息，那么可以將其作為約束增加到反演中[25]，但是背景信息往往不易直接獲得。為了得到背景信息并且使其值更加逼近真實，利用已知的地震層位信息對反演區域進行剖分，將每個層位的電阻率作為一個參數參與反演，得到工區的電性參數背景(對應公式(13)，記為第一步)

圖3 數值解與解析解對比Fig.3 Comparison of numerica and analytical solutions(a)電場Ey振幅響應曲線;(b)不同迭代次數下相對誤差曲線

(13)

然后選擇它與模型m之差的二范數定義為模型空間的目標函數，結合數據空間的目標函數可形成新的目標函數，對新的目標函數繼續實施反演直至收斂(公式(14)，記為第二步)。

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

(14)

對于給定的初始模型mk，為了使總目標函數最小，需要求取目標函數的梯度，令mU=0，可得到新的迭代模型mk+1[18]：

對式(13)有

(15)

對式(14)有

mk+1=[λPTP+(WJk)TWJk]-1×

(16)

(17)

反演迭代過程中，通過計算模型的修改量以確定新的模型mk+1，從而計算、判斷模型的正演響應與觀測數據之間的擬合差是否達到設置要求，實現目標函數最小化求解。過程如下：首先利用層位信息建立反演初始模型，移除層位中模型粗糙度，并設置反演電阻率的左右邊界范圍為10-1Ω·m～103Ω·m，對其整體實施反演，直到RMS收斂為止。其次，用第一步反演結果作為初始模型，設置對應的P值，在目標區域范圍內采用矩形剖分，擴展區三角剖分，最大縱橫比設置為10。實施第二步反演，直到RMS不再下降為止。

3 分區并行方式

OCCAM反演方法需要計算和存儲整個反演域上的靈敏度矩陣，從而耗費大量計算時間和內存。根據航空footprint概念，某測點響應大小與其下方一定范圍內地電模型有關[26]。因此，為了提高反演計算效率，我們嘗試將footprint引入到海洋可控源電磁中。如圖4所示，計算了頻率為0.3 Hz的10個接收機為一組和發射源(O)分別位于-20 km、-5 km、10 km處的歸一化靈敏度示意圖。通過對MCSEM影響范圍進行分析可知，歸一化靈敏度值主要集中測點下方所在的有效區域范圍內，虛線范圍外歸一化靈敏度值相對很小。由于接收站附近“敏感區域”遠小于整個測區，因此在海洋電磁數據反演時，可先基于“敏感區域”大小對測區進行單元劃分，從而減小反演模型規模，在完成各測區單元的電磁數據反演后，組合得到反演區域內總的電性分布。

圖4 歸一化靈敏度示意圖Fig.4 Normalized sensitivity diagram(a)-20 km;(b)-5 km;(c)10 km

圖5 模型一Fig.5 Model I(a)模型示意圖;(b)傳統方法反演結果;(c)分區并行方式反演結果;(d)傳統方法、分區并行方式反演RMS對比

4 反演算例

4.1 模型一

為了測試分區并行方式的有效性，設計了如圖5(a)所示的二維地電模型。參數設置如下：空氣電阻率為109Ω·m；海水層厚度為1 km，電阻率0.3 Ω·m；海底巖層電阻率1 Ω·m，在海底下方0.6 km處有一個8 km×0.2 km的高阻薄層，電阻率大小為50 Ω·m；電偶源布置在海底上方0.05 km處，沿測線方向等間距分布100個(-25 km ≤y≤25 km)，發射頻率為0.3 Hz，發射電流1A；100個海底接收點等間距分布于y=-24 km～25 km的范圍內。正演計算得到電場Ey分量的幅值、相位數據，加入4%高斯噪聲，作為反演的擬合數據。反演初始模型由空氣-海水-巖層組成，其中空氣、海水層不參與反演，巖層初始電阻率為1 Ω·m。圖5(b)、5(c)展示了單節點28個處理器反演迭代20次結果，由圖可以發現二者反演結果很相似，異常體的位置與實際情況基本吻合，并且電阻率值和真實的電阻率值很接近。圖5(d)展示了單節點28個處理器反演迭代20次RMS(數據擬合誤差)曲線，由圖可以得出，二者初始RMS均為17，經過20次迭代最終RMS均收斂到1.0左右。為了體現分區并行方式的優越性，計算了單節點28個處理器反演迭代20次的耗時，傳統方法需要約29小時，分區并行方式反演需13.5小時左右，前者比后者的反演耗時增加了約兩倍。此外，分區并行比傳統方法反演的內存減少1/3左右。

圖6 反演并行結果Fig.6 Inversion of parallel results(a)不同處理器核數反演耗時;(b)分區并行方式加速比曲線

加速比是衡量并行處理性能的重要指標之一，因此，首先定義加速比為

SP=T1/TP

(18)

其中:SP為加速比;p為并行計算所用處理器的核數；T1為串行計算所花時間；Tp為p個處理器并行計算所花時間。

為了表示效率的提升情況，在圖5(a)所示模型基礎上，計算了不同處理器核數(1、2、6、12、18、24、28)反演迭代20次消耗CPU時間以及相應的加速比。并行計算環境為一個擁有28個Intel(R) Xeon (R) CPU E5-2620 v4 @ 2.10GHz的處理器構成的集群上，操作系統為SUSE Linux Enterprise Server 11 (x86_64)，內存128GB。圖6(a)為傳統方法、分區并行方式反演迭代20次消耗時間對比。由圖6可以看出，處理器核數相同時，傳統方法比分區并行方式反演耗時增加了兩倍左右。圖6b展示了加速比與并行計算所用處理器核數之間的關系，隨著計算規模的逐漸增加，計算效率有很大提升，同時取得了很好的加速比。

4.2 模型二

為了測試背景模型約束的OCCAM反演方法的有效性，設計了如圖7(a)所示的模型。其中，海底地形和地層信息來自實際數據資料，每層的電阻率均以各向同性的方式填充。空氣的電阻率為1×109Ω·m；海水電阻率為0.3 Ω·m；海底以下有一高阻薄層，寬度約6 km，深度為3.2 km～3.4 km，電阻率大小為60 Ω·m；電偶源布置在海底上方0.05 km處，發射頻率為0.1 Hz，發射電流為1A；15個接收點布設在海底，范圍為-13 km～6 km，間隔為1.35 km；正演計算電場Ey分量的幅值、相位數據，加入4%高斯噪聲，作為反演的擬合數據；反演初始模型由空氣-海水-巖層組成，空氣、海水層不參與反演，巖層初始電阻率為1 Ω·m半空間。圖7(b)展示了光滑反演(式(13))結果，可以看出，光滑反演結果中異常體邊界不夠清晰，并且異常體的縱向范圍有所拉伸。圖7(c)、圖7(d)為背景模型約束的OCCAM反演結果，由圖7(c)可以得出，整體電性結構與真實電性結構很接近，但靠近海底的兩層信息沒有明顯反映出來，是因為它們的(海底～2.5 km、2.5 km～2.9 km)真實電阻率值很接近。由圖7d可以看出，實施第二步反演后(式(14))，靠近海底的兩層信息略微呈現了出來，將其與光滑反演結果(圖7(b))對比可以發現，整體電性結構以及電阻率變化趨勢與真實情況基本一致，異常體電性邊界結構不僅更加清晰，而且位置更加準確、形態更加飽滿。圖7(e)、圖7(f)為RMS隨迭代次數的變化曲線，由圖7(e)可得，初始RMS約為13，經過14次迭代RMS降到1.9左右；由圖7(f)可以看出，經過第二步反演后RMS最終收斂到1.2左右。

圖7 模型二Fig.7 Model II(a)模型示意圖;(b)光滑反演結果;(c)第一步反演結果;(d)第二步反演結果;(e)、(f)RMS曲線

5 結論

筆者采用的反演策略與傳統反演方法進行對比，表明反演區域分解方法可以提升反演過程中內存的利用率和計算效率。為了使反演結果得到一個合理的解釋，利用模型約束信息首先對合成數據進行光滑反演，得到工區的電性參數背景，然后對模型施加背景約束繼續進行反演。實驗結果表明了方法的有效性、穩定性。

致謝

感謝河南理工大學高性能計算中心支持，感謝審稿專家提出的寶貴意見。