本發(fā)明涉及地球物理勘探領(lǐng)域中地下介質(zhì)的電導(dǎo)率任意各向異性問題，提出了一種基于二次場矢量位、標(biāo)量位的海洋可控源電磁三維非結(jié)構(gòu)化有限元數(shù)值模擬方法，適用于幾何復(fù)雜、物性分布的航空電磁、井中電磁、陸面電磁等地球物理勘探方法數(shù)值模擬研究。

背景技術(shù)：

地球內(nèi)部的結(jié)構(gòu)和屬性是地球物理學(xué)研究的核心內(nèi)容。21世紀(jì)是地球科學(xué)進入對介質(zhì)，構(gòu)造和深層動力學(xué)過程的各向異性的深入研究的時代。隨著現(xiàn)代觀測技術(shù)的不斷進步以及認(rèn)識水平的不斷提高，各向異性問題逐漸引起了國內(nèi)外學(xué)者廣泛的重視，成為地球物理學(xué)研究的熱點。

海洋可控源電磁法是利用人工源電磁場研究地球電性結(jié)構(gòu)的一種地球物理勘探方法。近十年來，該方法發(fā)展迅速，已經(jīng)成為實現(xiàn)尋找油氣田、研究海洋深部構(gòu)造等勘探目的一種行之有效的勘探方法?，F(xiàn)有的海洋可控源電磁數(shù)值模擬方法主要有積分方程法(iem),有限差分法(fdm)，有限單元法(fem)等。

目前為止，對于電導(dǎo)率各向異性介質(zhì)的海洋可控源電磁數(shù)值模擬，大都是僅僅考慮三主軸電導(dǎo)率各向異性，這固然簡化問題，然而，由地質(zhì)運動所造成的地層的隆起、翻轉(zhuǎn)等引起的電導(dǎo)率各向異性往往并不局限于三主軸，甚至可能更為復(fù)雜，并且在實際探測中，復(fù)雜的海底地形與電導(dǎo)率各向異性的影響對實際數(shù)據(jù)的解釋可能會產(chǎn)生一定的偏差。

技術(shù)實現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題在于提供一種各向異性介質(zhì)的海洋可控源電磁法有限元正演方法，實現(xiàn)高效、快速的海洋電磁各向異性介質(zhì)的數(shù)值模擬。

本發(fā)明是這樣實現(xiàn)的，

一種各向異性介質(zhì)的海洋可控源電磁法有限元正演方法，該方法包括：

1)對研究區(qū)域進行非結(jié)構(gòu)化網(wǎng)格剖分，得到網(wǎng)格單元編號，節(jié)點編號和坐標(biāo)參數(shù)；

2)讀取設(shè)計的地電模型參數(shù)，包括背景層參數(shù)、頻率參數(shù)、參考電導(dǎo)率、三個歐拉旋轉(zhuǎn)角度、網(wǎng)格單元編號以及節(jié)點坐標(biāo)；

3)計算背景層相應(yīng)的一次場電磁各分量；

4)對所有網(wǎng)格單元進行循環(huán)計算單元系數(shù)矩陣，接著總體合成所有網(wǎng)格單元的系數(shù)矩陣，接著根據(jù)步驟3)背景層計算出來的一次電磁場計算線性方程組的右端項；

5)加載本質(zhì)邊界條件，求解線性方程組，得到各個節(jié)點的二次場矢量位及標(biāo)量位；

6)對二次場矢量位、標(biāo)量位進行求導(dǎo)，得到所有節(jié)點的電磁場各分量。

進一步地，將研究區(qū)域剖分成有限多個四面體單元e。

進一步地，該方法還包括設(shè)定一個參考電導(dǎo)率，該參考電導(dǎo)率三個非零對角線元素的電導(dǎo)率，接著設(shè)定三個歐拉旋轉(zhuǎn)角度，經(jīng)過三次歐拉旋轉(zhuǎn)，得到任意方向的電導(dǎo)率張量模型。

進一步地，接著從麥克斯韋方程組出發(fā)，引入基于coulomb規(guī)范下的磁矢量位a、標(biāo)量位ψ來表示電場、磁場，在二次場中，總磁矢量位和標(biāo)量分解成二次場與背景場之和，得到關(guān)于電磁場二次位表達(dá)式：

其中，為異常電導(dǎo)率，ap為一次場矢量位，as為二次場矢量位，ψp為一次場標(biāo)量位，ψs為二次場標(biāo)量位。

進一步地，將電磁場二次位表達(dá)式按照xyz軸依次展開，利用伽遼金方法對展開式進行加權(quán)，結(jié)合矢量恒等式以及散度定理，得到關(guān)于二次位的體積分方程組。

進一步地，步驟4)中線性方程組：

ku＝b

k為系數(shù)矩陣，u為求解域中各個節(jié)點待求的場值，b為右端項。

進一步地，步驟5)中求解線性方程組：ku＝b，加入狄利克雷邊界條件：(as,ψs)γ＝0。

進一步地，采用不完全lu分解的預(yù)條件因子的idr(s)迭代算法對線性方程組進行求解。

進一步地，步驟6)采用指數(shù)加權(quán)移動平均最小二乘法求解二次場矢量位、標(biāo)量位對x，y，z的偏導(dǎo)數(shù)。

本發(fā)明與現(xiàn)有技術(shù)相比，有益效果在于：

本發(fā)明提出了一種基于二次場矢量位、標(biāo)量位的非結(jié)構(gòu)化有限元數(shù)值模擬算法來模擬電導(dǎo)率任意各向異性的海洋可控源電磁問題。該方法既可以避免源點的奇異性的影響，又能滿足節(jié)點有限元對于場的連續(xù)性要求，同時采用非結(jié)構(gòu)化網(wǎng)格，有利于構(gòu)建幾何復(fù)雜地電模型，更好地模擬真實的地質(zhì)情況，因此，可以實現(xiàn)高效、快速的海洋電磁各向異性介質(zhì)的數(shù)值模擬。

本發(fā)明主要針對可控源電磁探測地下介質(zhì)中廣泛存在的電導(dǎo)率任意方向各向異性的問題，首先假設(shè)一個電導(dǎo)率張量及三個歐拉旋轉(zhuǎn)角度，通過進行三次歐拉旋轉(zhuǎn)，構(gòu)造了任意方向各向異性的電導(dǎo)率張量模型。并且首次推導(dǎo)了基于二次場磁矢量位、標(biāo)量位的電導(dǎo)率任意各向異性條件下的可控源電磁非結(jié)構(gòu)化有限元方程，可以有效避免源點奇異性的影響，提高數(shù)值解的精度，又能夠滿足節(jié)點有限元對于場的連續(xù)性的要求。為了更好模擬復(fù)雜地質(zhì)結(jié)構(gòu)，本發(fā)明采用非結(jié)構(gòu)化網(wǎng)格對研究區(qū)域進行離散剖分，可構(gòu)建復(fù)雜的地質(zhì)結(jié)構(gòu)及物性分布，同時避免了對全區(qū)域進行同等程度的加密，減少了計算量。針對常規(guī)迭代算法求解線性方程組收斂慢的問題，提出了采用不完全lu分解預(yù)處理方法與krylov子空間迭代算法中的idr(s)法相結(jié)合對線性方程組求解，有效地提高了計算效率。采用指數(shù)加權(quán)的滑動平均最小二乘法對二次場矢量位、標(biāo)量位進行求導(dǎo)，相對于常規(guī)的求導(dǎo)方法，進一步提高了導(dǎo)數(shù)求解的精度。因此，本發(fā)明可以實現(xiàn)高精度、快速的電導(dǎo)率呈任意各向異性條件下的海洋可控源電磁法數(shù)值模擬。

附圖說明

圖1是本發(fā)明實施例提供的非結(jié)構(gòu)化網(wǎng)格剖分單元的節(jié)點編號示意圖；

圖2是本發(fā)明實施例提供的任意各向異性介質(zhì)電導(dǎo)率張量的構(gòu)建相應(yīng)的三次歐拉旋轉(zhuǎn)示意圖；其中，圖2(a)沿x軸，圖2(b)沿y軸，圖2(c)沿z軸；

圖3是本發(fā)明實施例提供的任意各向異性介質(zhì)海洋可控源電磁法非結(jié)構(gòu)化有限元數(shù)值模擬流程圖；

圖4是本發(fā)明實施例建立的水平層狀電導(dǎo)率各向異性介質(zhì)模型；

圖5是本發(fā)明實施例各向異性層狀模型擬解析解與有限元數(shù)值解的對比；

圖6是本發(fā)明實施例建立的任意各向異性介質(zhì)中含有高阻異常體模型；

圖7是本發(fā)明實施例建立的任意各向異性介質(zhì)中含有高阻異常體模型非結(jié)構(gòu)化網(wǎng)格剖分示意圖；

圖8是本發(fā)明實施例建立的任意各向異性介質(zhì)中含有高阻異常體模型電磁場各分量幅值第一平面圖，其中圖8(a)、圖8(b)、圖8(c)為參考電導(dǎo)率不旋轉(zhuǎn)時電場振幅平面圖；8(d)、8(e)、8(f)、為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)30。時的電場振幅平面等值線圖圖；8(g)、8(h)、8(l)為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)45。時的電場振幅平面等值線圖；8(m)、8(n)、8(o)為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)60。時的電場振幅平面等值線圖；

圖9是本發(fā)明實施例建立的任意各向異性介質(zhì)中含有高阻異常體模型電磁場各分量幅值第二平面圖，圖9(a)、圖9(b)、圖9(c)分別為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)30。時的電場振幅平面等值線圖；圖9(d)、圖9(e)、圖9(f)為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)45。時的電場振幅平面等值線圖；圖9(g)、圖9(h)、圖9(l)為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)60。時的電場振幅平面等值線圖。

具體實施方式

為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點更加清楚明白，以下結(jié)合實施例，對本發(fā)明進行進一步詳細(xì)說明。應(yīng)當(dāng)理解，此處所描述的具體實施例僅僅用以解釋本發(fā)明，并不用于限定本發(fā)明。

參見圖1，本發(fā)明實施例提供的非結(jié)構(gòu)化網(wǎng)格剖分的四面體單元節(jié)點(1、2、3、4)編號示意圖。

采用現(xiàn)成的網(wǎng)格剖分軟件，應(yīng)用任意四面體網(wǎng)格單元對研究區(qū)域進行離散剖分，在電性分界處及源點區(qū)域進行加密，離激發(fā)源較遠(yuǎn)的地方剖分稀疏些，這樣便可避免了對全區(qū)域進行同等程度的加密，節(jié)省了計算內(nèi)存，減少了計算量。

參見圖3，任意電導(dǎo)率各向異性介質(zhì)可控源電磁有限元模擬方法流程圖，包括如下的步驟：

1)對研究區(qū)域進行非結(jié)構(gòu)化網(wǎng)格剖分，得到網(wǎng)格單元編號，節(jié)點編號和坐標(biāo)等參數(shù)；

2)讀取設(shè)計的地電模型參數(shù)，包括背景層參數(shù)、頻率參數(shù)、參考電導(dǎo)率、三個歐拉旋轉(zhuǎn)角度、網(wǎng)格單元編號、節(jié)點坐標(biāo)等；

3)計算背景層相應(yīng)的一次場電磁各分量；

4)對所有網(wǎng)格單元進行循環(huán)計算系數(shù)矩陣，接著總體合成所有單元的單元系數(shù)矩陣，同時，根據(jù)背景層計算出來的一次電磁場計算線性方程組的右端項。

5)加載本質(zhì)邊界條件，求解線性方程組，得到各個節(jié)點的二次場矢量位及標(biāo)量位；

6)對二次場矢量位、標(biāo)量位進行求導(dǎo)，得到所有節(jié)點的電磁場各分量。

電導(dǎo)率任意各向異性介質(zhì)的可控源電磁邊值問題地球物理勘探中，往往是采用低頻電磁信號，忽略位移電流,采用正弦諧變時間因子e^-(wt，maxwell′s方程可表示如下形式:

其中：w為角頻率，μ0為真空中的磁導(dǎo)率。.s為場源電流分布。為任意各向異性電導(dǎo)率張量，e表示的電場，h表示的磁場，如下所示：

為了計算電導(dǎo)率張量通常事先設(shè)定一個參考電導(dǎo)率σ4(其中三個對角線的元素表示三個主軸的電導(dǎo)率)，如下所示：

接著通過三次歐拉旋轉(zhuǎn)，參見圖2(a)、圖2(b)、圖2(c)，便可得到任意各向異性介質(zhì)的電導(dǎo)率張量：

其中旋轉(zhuǎn)矩陣為：

α為坐標(biāo)軸沿著x軸旋轉(zhuǎn)角度；β為坐標(biāo)軸沿著y軸旋轉(zhuǎn)角度；γ為坐標(biāo)軸沿著z軸旋轉(zhuǎn)角度。，此過程對應(yīng)步驟3。

引入基于coulomb規(guī)范下的磁矢量位a、標(biāo)量位ψ來表示電場、磁場，如下所示：

將方程(6)、(7)代入方程(1)、(2)得：

在二次場算發(fā)中，總磁矢量位和標(biāo)量為可分解成二次場與背景場之和，如下所示：

a＝ap+as(11)

ψ＝ψp+ψs(12)

其中，為背景電導(dǎo)率，為異常電導(dǎo)率，ap為一次場矢量位，as為二次場矢量位，ψp為一次場標(biāo)量位，ψs為二次場標(biāo)量位。

將方程(10)-(12)代入方程(8)(9)可得關(guān)于電磁場二次位表達(dá)式：

通過觀察方程(13)、(14)可知，方程右端項是關(guān)于一次位，而不含電流項。因此可避免源點的奇異性的影響。而一次位可通過層狀介質(zhì)或者均勻半空間為背景層求得?？紤]如下關(guān)系：

ep為一次電場。

方程(13)(14)可改寫成如下形式：

這樣就避免求解一次場標(biāo)量位的梯度，提高一次場的精度。

有限單元分析:

將方程(13)(14)按照xyz軸依次展開，可得如下方程：

利用伽遼金方法對方程(15)至(18)進行加權(quán)，結(jié)合矢量恒等式以及散度定理，最后可得關(guān)于二次位的體積分方程組：

其中，n為線性插值函數(shù)。

本發(fā)明中，將研究區(qū)域剖分成有限多個四面體單元e，參見圖1，在各個單元中，電導(dǎo)率是固定不變的，單元中的二次矢量位和標(biāo)量位呈線性變化：

線性插值函數(shù)定義如下所示：

其中v^e為單元體積，分別為插值函數(shù)的系數(shù)。

將(23)(24)方程代入(19)-(22)最后可得離散線性方程組：

ku＝b(25)

其中:

b^e＝(bx^e,by^e,b2^e,bψ^e)^t

u^e＝(asx^e,asy^e,as2^e,ψs^e)^t

為了求解方程(25)，還需要加入相應(yīng)的邊界條件。本發(fā)明采用狄利克雷邊界條件：(as,ψs)γ＝0(26)。

步驟5中線性方程組的求解：krlov子空間迭代算法因收斂速度快，求解精度高，而且穩(wěn)定性好等優(yōu)點而被廣泛用于求解大型稀疏線性方程組，特別是當(dāng)其與預(yù)處理技術(shù)結(jié)合能夠有效加快求解線性方程組的速度，為了進一步提高求解的效率，本發(fā)明采用不完全lu分解的預(yù)條件因子的idr(s)迭代算法對線性方程組進行求解。

步驟6電磁場各分量的求?。?/p>

電磁場二次場各分量可通過以下方程進行求解：

指數(shù)加權(quán)移動平均最小二乘法具體實現(xiàn)過程：

首先設(shè)單元中的二次位的線性描述：

fi＝t1xi+t2yi+t3zi+d(29)

這里:xi、yi、zi為所需要求解的二次位的網(wǎng)格節(jié)點x、y、z坐標(biāo)。

顯而易見，系數(shù)t1、t2、t3即為待求解的二次位的導(dǎo)數(shù)，即為電磁場二次場各分量，在此直接給出高斯加權(quán)的滑動平均最小二乘法所形成的方程，如下所示：

t＝(x^twx)^-1x^twf(30)，

這里

x＝(x,y,z,m),x＝(x1,x2,x3…xn)^t

y＝(y1,y2,y3…yn)^t,z＝(z1,z2,z3…zn)^t

m＝(1,1,1…1)^t

其中，x為表示網(wǎng)格節(jié)點三維空間坐標(biāo)數(shù)組，w為權(quán)函數(shù)，r表示網(wǎng)格節(jié)點到坐標(biāo)原點的距離。h為r最大值。

權(quán)函數(shù)

h＝max(r)；

這里β為權(quán)函數(shù)系數(shù)，一般取1。f為待求二次位導(dǎo)數(shù)點附近n個節(jié)點的二次位值，x、y、z為待求二次位的節(jié)點坐標(biāo)。

為了驗證本發(fā)明中方法的正確性。設(shè)計一個三層水平層狀沉積模型，如圖4所示。采用海水-沉積層兩層介質(zhì)為背景層。設(shè)上半空間中的空氣電導(dǎo)率為10^(-12)s·m-1。海水電導(dǎo)率為3.2s·m-1，深度為0.3km。距離海水底面深度1km有一水平方向無限延伸的各向異性高阻體，其厚度為100m，參考電導(dǎo)率為σc＝diag(0.01,0.01,0.025)s·m-¹，向y軸旋轉(zhuǎn)30。沉積層的電導(dǎo)率為1s·m-1。采用沿著x方向的水平電偶極子為激發(fā)源，其坐標(biāo)位于(0,0,970)。發(fā)射頻率為0.25hz。偶極矩為1。沿著z＝970處進行觀測。剖分區(qū)域大小(-4km≤x,y≤4km,-1km≤z≤5km)，網(wǎng)格單元總數(shù)1467492。在amdathlon(tm)x463quad-coreprocessor2.6ghz，內(nèi)存12g計算平臺上執(zhí)行。

圖5為在計算區(qū)域內(nèi)的有限元數(shù)值解與擬解析解的電場x分量幅值對比圖，擬解析解的來源為losethlo,ursimb.electromagneticfieldsinplanarlylayeredanisotropicmedia[j]。從圖5中可以看出，有限元數(shù)值解與擬解析解的數(shù)據(jù)兩者高度吻合，充分驗證了本發(fā)明關(guān)于電導(dǎo)率任意各向異性問題的電磁法數(shù)值模擬算法的正確性，可為其他電磁法在各向異性介質(zhì)方面的研究提供一種新的方法。

圖6為建立的任意各向異性介質(zhì)中含有高阻異常體模型，上半空間空氣電導(dǎo)率設(shè)為10^(-12)，海水深度為1km，電導(dǎo)率為3.3s·m-1。沉積層中有一高阻體，頂面距離海水底面1km，大小為2kmx2kmx0.1km，中心坐標(biāo)為(2000,0,2050)，電導(dǎo)率為0.01s·m-1。沉積層的參考電導(dǎo)率設(shè)為σc＝diag(2,1,1)。激發(fā)源為沿著x方向的水平電偶極子，發(fā)射頻率為0.25hz，偶極距為100，源點坐標(biāo)(0,0,950)。接收機位于海水深度z＝950處。非結(jié)構(gòu)化網(wǎng)格剖分如圖7所示，剖分區(qū)域大小(-4km≤x,y≤4km,-1km≤z≤4km)。令沉積層參考電導(dǎo)率分別沿著y軸、z軸旋轉(zhuǎn)0。、30。、45。、60。后研究電場三分量幅值的異常響應(yīng)特征。

圖8，圖9分別為沉積層參考電導(dǎo)率沿著y軸跟z軸旋轉(zhuǎn)不同角度后電場振幅平面分布圖。圖8中圖8(a)、8(b)、8(c)為參考電導(dǎo)率不旋轉(zhuǎn)時電場振幅平面圖；8(d)、8(e)、8(f)、為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)30。時的電場振幅平面等值線圖圖；8(g)、8(h)、8(l)為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)45。時的電場振幅平面等值線圖；8(m)、8(n)、8(o)為參考電導(dǎo)率沿著y軸旋轉(zhuǎn)60。時的電場振幅平面等值線圖。

圖9是本發(fā)明實施例建立的任意各向異性介質(zhì)中含有高阻異常體模型電磁場各分量幅值第二平面圖，9(a)、9(b)、9(c)分別為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)30。時的電場振幅平面等值線圖；9(d)、9(e)、9(f)為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)45。時的電場振幅平面等值線圖；9(g)、9(h)、9(l)為參考電導(dǎo)率沿著z軸旋轉(zhuǎn)60。時的電場振幅平面等值線圖。

從圖8可知，當(dāng)參考電導(dǎo)率張量沿著y軸旋轉(zhuǎn)的角度30。、45。時，電場ex分量的幅值在x、y方向都有所延伸，特別在y方向的延展更為明顯；而電場ey分量正好相反。相比于ex、ey分量，ez分量幅值變化幅度比較大，由一開始的對稱到明顯的不對稱，且右側(cè)的響應(yīng)值比左側(cè)的大。當(dāng)參考電導(dǎo)率張量沿著y軸旋轉(zhuǎn)60。時，電場ex分量幅值分布有所收縮；電場ey分量變化不大；而電場ez幅值分布明顯向外擴展，且左側(cè)的幅值變大，但右側(cè)幅值依然比左側(cè)大。通過觀察圖9，可知，參考電導(dǎo)率沿著z軸旋轉(zhuǎn)30。、45。時，電場ex分量分布發(fā)生傾斜向外擴張.ey分量的幅值隨著旋轉(zhuǎn)的角度增大而有所增大，且在x、y方向分布也有所延伸；ez分量幅值分布發(fā)生明顯的傾斜變化。當(dāng)參考電導(dǎo)率沿著z軸旋轉(zhuǎn)60。時，ex分量分布進一步發(fā)生擴展；ey分量的幅值分布有所收縮；ez分量幅值分布變成了斜對稱。因此，可知不同方向的電導(dǎo)率各向異性所引起的電場異常響應(yīng)與三主軸電導(dǎo)率各向異性的響應(yīng)特征有明顯的區(qū)別。

以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)所作的任何修改、等同替換和改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。