張迎利， 李 賡

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

0 引言

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

正演通過使用非結(jié)構(gòu)化網(wǎng)格，可以很好地適應(yīng)起伏地形和地下復(fù)雜電性結(jié)構(gòu)[10]。但是在大規(guī)模電磁數(shù)值模擬中，通過增加適當(dāng)?shù)木W(wǎng)格來適應(yīng)地下電性結(jié)構(gòu)是提高計(jì)算效率的關(guān)鍵。近些年，自適應(yīng)有限元方法針對以上問題給出了很好的解決方案。其中二維電磁場數(shù)值模擬中基于非結(jié)構(gòu)化網(wǎng)格的自適應(yīng)有限元法已經(jīng)取得了很好的應(yīng)用效果[11-13]。

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

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

1 正演模擬

1.1 控制方程

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

(1)

式中:E為電場強(qiáng)度；ω為角頻率；μ0為真空中的磁導(dǎo)率；H為磁場強(qiáng)度；σ為介質(zhì)電導(dǎo)率；J為外加源電流密度。

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

(2)

(3)

(4)

(5)

(6)

(7)

(8)

1.2 自適應(yīng)有限元

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

(9)

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

B(un,v)=F(v)

(10)

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

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

(11)

定義誤差函數(shù)[12]

(12)

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

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

1.3 正確性驗(yàn)證

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

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

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

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

(13)

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

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

(14)

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

對式(13)有

(15)

對式(14)有

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

(16)

(17)

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

3 分區(qū)并行方式

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

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

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

4 反演算例

4.1 模型一

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

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

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

SP=T1/TP

(18)

其中:SP為加速比;p為并行計(jì)算所用處理器的核數(shù)；T1為串行計(jì)算所花時(shí)間；Tp為p個(gè)處理器并行計(jì)算所花時(shí)間。

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

4.2 模型二

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

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

5 結(jié)論

筆者采用的反演策略與傳統(tǒng)反演方法進(jìn)行對比，表明反演區(qū)域分解方法可以提升反演過程中內(nèi)存的利用率和計(jì)算效率。為了使反演結(jié)果得到一個(gè)合理的解釋，利用模型約束信息首先對合成數(shù)據(jù)進(jìn)行光滑反演，得到工區(qū)的電性參數(shù)背景，然后對模型施加背景約束繼續(xù)進(jìn)行反演。實(shí)驗(yàn)結(jié)果表明了方法的有效性、穩(wěn)定性。

致謝

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