本發(fā)明涉及一種二度體磁場數(shù)值計(jì)算方法，特別是任意截面形狀、任意磁化率分布的二度體磁場快速、高精度數(shù)值計(jì)算方法。

背景技術(shù)：

磁法勘探以地球磁場為根據(jù)、磁化率差異為物理基礎(chǔ)，是一種研究地球構(gòu)造及尋找固體礦產(chǎn)資源的地球物理勘探方法，該方法效率高、操作簡單、成本低、勘探深度大、實(shí)施過程沒有過多條件限制。實(shí)際勘探工作中，常見這樣一類線性地質(zhì)體：其走向方向的尺度要遠(yuǎn)比垂直其走向方向的尺度大，諸如此類地質(zhì)體，實(shí)際場源分布便可以用走向方向無限延伸的二度體代替。

伴隨傳感器和測量技術(shù)的發(fā)展，尤其是航空磁測的出現(xiàn)，使得磁場測量在精度、空間分辨率和范圍上都有了顯著提高，為磁化率反演提供了大規(guī)模高精度、高分辨率數(shù)據(jù)。同時(shí)，隨著計(jì)算機(jī)軟硬件水平的不斷提高，人機(jī)交互建模、解釋也日益得到人們的重視，在地球物理勘探中發(fā)揮著越來越重要作用。人機(jī)交互建模能夠使人們通過直觀的方式對(duì)地質(zhì)體進(jìn)行建模，更容易融合地下結(jié)構(gòu)的先驗(yàn)信息。反演方法與人機(jī)交互建模、解釋方法相輔相成，將極大提高人們對(duì)地球內(nèi)部結(jié)構(gòu)的認(rèn)識(shí)。

磁場正演計(jì)算是指根據(jù)磁化率分布計(jì)算磁場，反演是指根據(jù)觀測磁場數(shù)據(jù)計(jì)算磁化率分布。正演是反演的基礎(chǔ)，正演計(jì)算效率直接影響反演計(jì)算效率，而正演計(jì)算精度直接影響反演結(jié)果的質(zhì)量。一直以來，學(xué)者們都非常重視磁場的正演計(jì)算。對(duì)于任意截面形狀、任意磁化率分布的二度體磁場正演計(jì)算，一般采用數(shù)值方法?；舅悸肥菍?fù)雜二度體剖分成許多規(guī)則二度體，然后利用規(guī)則二度體磁場的累加來逼近復(fù)雜二度體磁場。文獻(xiàn)(劉增群.任意多邊形截面二度體磁場正演程序.物探化探計(jì)算技術(shù),1991,13(4):357-360.)公開了一種磁化率為常值情況下任意多邊形截面二度體磁場計(jì)算方法，將二度體分解為若干臺(tái)階的組合，利用解析公式計(jì)算臺(tái)階磁場，并將其累加，這種方法保證了計(jì)算精度，但計(jì)算效率低，且只能計(jì)算磁化率為常值情況；文獻(xiàn)(張嶺,郝天珧.基于delaunay剖分的二維非規(guī)則重力建模及重力計(jì)算.地球物理學(xué)報(bào),2006,49(3):877-884.)公開了一種delaunay剖分方式，將截面分為若干三角形，將二度體分解為若干三角棱柱的組合，利用解析公式計(jì)算變密度三角棱柱重力異常，并將其累加。這種方法適用于任意截面形狀、任意磁化率分布情況下二度體磁場正演計(jì)算，該方法保證了計(jì)算精度，但計(jì)算效率低。

剖分方式和計(jì)算方法共同決定了磁場正演計(jì)算的效率和精度。計(jì)算效率和計(jì)算精度是一對(duì)矛盾體，目前已有的方法存在的最大問題是不能同時(shí)保證計(jì)算效率和計(jì)算精度，無法滿足大規(guī)模磁場數(shù)據(jù)磁化率精細(xì)反演成像、人機(jī)交互建模和解釋的需求。因此，尋找一種計(jì)算效率高、同時(shí)能保證計(jì)算精度的磁場計(jì)算方法具有重要的現(xiàn)實(shí)意義。

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

針對(duì)目前現(xiàn)有二度體磁場數(shù)值計(jì)算方法不能同時(shí)保證計(jì)算效率和計(jì)算精度，無法滿足大規(guī)模磁場數(shù)據(jù)磁化率精細(xì)反演成像、人機(jī)交互建模和解釋的需求問題，本發(fā)明提供一種二度體磁場數(shù)值計(jì)算方法。

本發(fā)明的技術(shù)方案是：

一種二度體磁場數(shù)值計(jì)算方法，包括以下步驟：

第一步復(fù)雜二度體模型表示：

根據(jù)二度體展布范圍，建立包含所有目標(biāo)區(qū)域的矩形模型，確定矩形模型在x，z方向的起始位置，使得包含起伏地形的目標(biāo)區(qū)域完全嵌入在該矩形模型中；

然后將該矩形模型均勻劃分成若干個(gè)規(guī)則小矩形，確定小矩形的幾何尺寸δx，δz；

最后根據(jù)二度體的磁化率分布，對(duì)每個(gè)小矩形磁化率進(jìn)行賦值，每個(gè)小矩形磁化率為常值，不同矩形磁化率取值不同，以此刻畫任意截面形狀、任意磁化率分布的二度體；將位于空氣部分的小矩形的磁化率值設(shè)為零，以此刻畫起伏地形。

步驟二、矩形模型磁場計(jì)算：

步驟一中給出的矩形模型，其磁場計(jì)算公式為

式中，(xm,z0)表示觀測點(diǎn)坐標(biāo)，z0為常值；l表示矩形模型其z方向的小矩形剖分個(gè)數(shù)；m表示矩形模型其x方向的小矩形剖分個(gè)數(shù)；(ξp,ζr)表示編號(hào)為(p,r)的小矩形幾何中心坐標(biāo)；mx(ξp,ζr)和mz(ξp,ζr)分別表示編號(hào)為(p,r)的小矩形的磁化強(qiáng)度的x分量和z分量；hx(xm-ξp,z0-ζr)和hz(xm-ξp,z0-ζr)分別表示對(duì)應(yīng)磁化強(qiáng)度的x分量和z分量的加權(quán)系數(shù)。

步驟二中，矩形模型，其磁場計(jì)算公式的解算方法如下：

s2.1，根據(jù)地球主磁場模型igrf，計(jì)算剖分小矩形中心點(diǎn)位置(ξp,ζr)的地球主磁場x分量tx(ξp,ζr)，z分量tz(ξp,ζr)。根據(jù)位置(ξp,ζr)處二度體的磁化率κ(ξp,ζr)，計(jì)算磁化強(qiáng)度的x分量mx(ξp,ζr)和z分量mz(ξp,ζr)

mx(ξp,ζr)＝κ(ξp,ζr)·tx(ξp,ζr)(2)

mz(ξp,ζr)＝κ(ξp,ζr)·tz(ξp,ζr)(3)

s2.2，計(jì)算加權(quán)系數(shù)hx(xm-ξp,z0-ζr)和hz(xm-ξp,z0-ζr)，其計(jì)算公式分別為

式(4)和式(5)中，μ0表示真空磁導(dǎo)率，arctan()表示反余切函數(shù)運(yùn)算符，ln()表示自然對(duì)數(shù)運(yùn)算符；其它符號(hào)含義如下

x1＝ξp-0.5δx-xm，x2＝ξp+0.5δx-xm，z1＝ζr-0.5δz-z0，z2＝ζr+0.5δz-z0

δx，δz表示小矩形幾何尺寸。

與現(xiàn)有計(jì)算二度體磁場方法相比，采用式(4)和式(5)給出的方法計(jì)算加權(quán)系數(shù)是本申請(qǐng)的一個(gè)創(chuàng)新點(diǎn)，其能夠提高復(fù)雜截面形狀、任意磁化率分布二度體磁場計(jì)算精度，原因是采用式(4)和式(5)，保證了式(1)是計(jì)算矩形組合模型磁場的精確公式。

s2.3，計(jì)算相對(duì)z方向而言某一層矩形組合模型的磁場，其計(jì)算公式為

式中，表示第r層(r＝1,2,…,l)矩形組合模型在高度面z0產(chǎn)生的磁場；(xm,z0)表示離散觀測點(diǎn)坐標(biāo)；

s2.4，將各層矩形組合模型磁場進(jìn)行累加，得到整個(gè)組合模型的磁場，即

在步驟s2.3中，采用一維離散卷積快速算法，計(jì)算式(6)中的兩個(gè)一維離散卷積即和其計(jì)算步驟為：

(1)將加權(quán)系數(shù)h(x1-ξp,z0-ζr)排列成向量t，記為

t＝[t0,t1,t2,…,tm-1,0,tm-1,tm-2,…,t2,t1]^t(8)

式中，矩陣元素ti與加權(quán)系數(shù)h(x1-ξp,z0-ζr)存在關(guān)系

ti＝h(x1-ξi+1,z0-ζr)(9)

(2)將第r層磁化強(qiáng)度值m(ξp,ζr)(p＝1,2,…,m)排列成向量m，向量元素mi與磁化強(qiáng)度值存在關(guān)系

mi＝m(ξi,ζr)

(10)

將向量m補(bǔ)零擴(kuò)展成向量m^ext

式中，0m×1表示m×1零向量；

(3)計(jì)算

式中，fft()表示一維快速傅里葉變換；

(4)計(jì)算

式中，“.*”表示對(duì)應(yīng)元素相乘運(yùn)算；

(5)計(jì)算

式中，ifft()表示一維快速傅里葉反變換；

(6)提取矩陣f^ext的前m行元素，構(gòu)成向量f，即為一維離散卷積計(jì)算結(jié)果。具體來講，當(dāng)向量t中h分別取hx，hz，而向量m分別對(duì)應(yīng)取mx，mz時(shí)，即可得到式中兩個(gè)一維離散卷積和的累加結(jié)果。

與現(xiàn)有二度體磁場計(jì)算方法相比，采用上述一維離散卷積快速算法是本申請(qǐng)的創(chuàng)新點(diǎn)之一，它保證了在海量剖分矩形個(gè)數(shù)情況下磁場計(jì)算的高效性。

本發(fā)明是一個(gè)有機(jī)整體，即在特定的模型表示方式條件下，建立矩形二度體磁場疊加模型，根據(jù)一種特殊的加權(quán)系數(shù)計(jì)算公式，采用一維離散卷積快速算法，實(shí)現(xiàn)了磁場計(jì)算在效率和精度上的統(tǒng)一。

與現(xiàn)有磁場計(jì)算技術(shù)相比，本發(fā)明具有以下優(yōu)點(diǎn)：

(1)模型表示方法簡單、靈活，很容易刻畫復(fù)雜截面形狀、任意磁化率分布復(fù)雜二度體以及起伏地形；

(2)能夠?qū)崿F(xiàn)任意截面形狀、任意磁化率分布情況下復(fù)雜二度體磁場的快速、高精度計(jì)算，可以滿足大規(guī)模磁場數(shù)據(jù)磁化率精細(xì)反演、人機(jī)交互建模和解釋的需求；

(3)剖分矩形的個(gè)數(shù)巨大時(shí)，算法不但計(jì)算效率和計(jì)算精度高，并且所需計(jì)算機(jī)內(nèi)存??；

(4)磁化強(qiáng)度可以分為感應(yīng)磁化強(qiáng)度和剩余磁化強(qiáng)度。本發(fā)明既可以計(jì)算感應(yīng)磁化強(qiáng)度磁場，也可以計(jì)算剩余磁化強(qiáng)度磁場，還可以計(jì)算兩者同時(shí)存在時(shí)的磁場，具體體現(xiàn)在對(duì)mx(ξp,ζr)和mz(ξp,ζr)賦值環(huán)節(jié)。本發(fā)明以感應(yīng)磁化強(qiáng)度磁場計(jì)算為例，示出了計(jì)算感應(yīng)磁化強(qiáng)度磁場的具體計(jì)算過程。

附圖說明

圖1為本發(fā)明的計(jì)算流程圖；

圖2為復(fù)雜二度體模型表示；

圖3為截面為圓形二度體模型圖；

圖4為磁場計(jì)算值和理論值對(duì)比圖；

圖5為計(jì)算值與理論值的相對(duì)誤差；

圖中符號(hào)說明如下：

κ：表示磁化率；

本發(fā)明目的的實(shí)現(xiàn)、功能特點(diǎn)及優(yōu)點(diǎn)將結(jié)合實(shí)施例，參照附圖做進(jìn)一步說明。

具體實(shí)施方式

為使本發(fā)明的目的、技術(shù)方案和優(yōu)點(diǎn)更加清楚，下面將結(jié)合附圖對(duì)本發(fā)明實(shí)施方式作進(jìn)一步地詳細(xì)描述。

參照?qǐng)D1，為一種二度體磁場數(shù)值計(jì)算方法的流程圖，包括以下步驟：

步驟一、復(fù)雜二度體模型表示：

根據(jù)二度體展布范圍，建立包含所有目標(biāo)區(qū)域的矩形模型，確定矩形模型在x，z方向的起始位置，使得包含起伏地形的目標(biāo)區(qū)域完全嵌入在該矩形模型中；

然后根據(jù)實(shí)際問題需求，將該矩形模型均勻劃分成若干個(gè)規(guī)則小矩形(如圖2所示)，確定小矩形的幾何尺寸δx，δz；

最后根據(jù)二度體的磁化率分布，對(duì)每個(gè)小矩形磁化率進(jìn)行賦值，每個(gè)小矩形磁化率為常值，不同矩形磁化率取值不同，以此刻畫任意截面形狀、任意磁化率分布的二度體；將位于空氣部分的小矩形的磁化率值設(shè)為零，以此刻畫起伏地形。

步驟二、矩形模型磁場計(jì)算：步驟一中給出的矩形模型，其磁場計(jì)算公式為

式中，(xm,z0)表示觀測點(diǎn)坐標(biāo)，z0為常值；l表示矩形模型其z方向的小矩形剖分個(gè)數(shù)；m表示矩形模型其x方向的小矩形剖分個(gè)數(shù)；(ξp,ζr)表示編號(hào)為(p,r)的小矩形幾何中心坐標(biāo)；mx(ξp,ζr)和mz(ξp,ζr)分別表示編號(hào)為(p,r)的小矩形的磁化強(qiáng)度的x分量和z分量；hx(xm-ξp,z0-ζr)和hz(xm-ξp,z0-ζr)分別表示對(duì)應(yīng)磁化強(qiáng)度的x分量和z分量的加權(quán)系數(shù)。

步驟二中，矩形模型，其磁場計(jì)算公式的解算方法如下：

s2.1，根據(jù)地球主磁場模型igrf，計(jì)算剖分小矩形中心點(diǎn)位置(ξp,ζr)的地球主磁場x分量tx(ξp,ζr)，z分量tz(ξp,ζr)。根據(jù)位置(ξp,ζr)處的磁化率κ(ξp,ζr)，計(jì)算磁化強(qiáng)度的x分量和z分量

mx(ξp,ζr)＝κ(ξp,ζr)·tx(ξp,ζr)(2)

mz(ξp,ζr)＝κ(ξp,ζr)·tz(ξp,ζr)(3)

s2.2，計(jì)算加權(quán)系數(shù)hx(xm-ξp,z0-ζr)和hz(xm-ξp,z0-ζr)，其計(jì)算公式分別為

式(4)和式(5)中，μ0表示真空磁導(dǎo)率，arctan()表示反余切函數(shù)運(yùn)算符，ln()表示自然對(duì)數(shù)運(yùn)算符；其它符號(hào)含義如下

x1＝ξp-0.5δx-xm，x2＝ξp+0.5δx-xm，z1＝ζr-0.5δz-z0，z2＝ζr+0.5δz-z0

δx，δz表示小矩形幾何尺寸。

s2.3，采用一維離散卷積快速算法，計(jì)算式(6)中的兩個(gè)一維離散卷積，得到一層(相對(duì)z方向而言)矩形組合模型的磁場，其計(jì)算公式為

式中，表示第r層(r＝1,2,…,l)矩形組合模型在高度面z0產(chǎn)生的磁場；(xm,z0)表示離散觀測點(diǎn)坐標(biāo)；

s2.4，將各層矩形組合模型磁場進(jìn)行累加，得到整個(gè)組合模型的磁場，即

在步驟s2.3中，采用一維離散卷積快速算法，計(jì)算式(6)中的兩個(gè)一維離散卷積和其計(jì)算方法為：

(1)將加權(quán)系數(shù)h(x1-ξp,z0-ζr)排列成向量t，記為

t＝[t0,t1,t2,…,tm-1,0,tm-1,tm-2,…,t2,t1]^t(8)

式中，矩陣元素ti與加權(quán)系數(shù)h(x1-ξp,z0-ζr)存在關(guān)系

ti＝h(x1-ξi+1,z0-ζr)(9)

(2)將第r層磁化強(qiáng)度值m(ξp,ζr)(p＝1,2,…,m)排列成向量m，向量元素mi與磁化強(qiáng)度值存在關(guān)系

mi＝m(ξi,ζr)

(10)

將向量m補(bǔ)零擴(kuò)展成向量m^ext

式中，0m×1表示m×1零向量；

(3)計(jì)算

式中，fft()表示一維快速傅里葉變換；

(4)計(jì)算

式中，“.*”表示對(duì)應(yīng)元素相乘運(yùn)算；

(5)計(jì)算

式中，ifft()表示一維快速傅里葉反變換；

(6)提取矩陣f^ext的前m行元素，構(gòu)成向量f，即為一維離散卷積計(jì)算結(jié)果。具體來講，當(dāng)向量t中h分別取hx，hz，而向量m分別對(duì)應(yīng)取mx，mz時(shí)，即可得到式中兩個(gè)一維離散卷積和的累加結(jié)果。

下面對(duì)本發(fā)明方法的效果進(jìn)行檢驗(yàn)。

為了說明本發(fā)明所提出的方法用于計(jì)算任意截面形狀、任意磁化率分布情況下復(fù)雜二度體重力異常時(shí)的效率和精度，設(shè)計(jì)了如下二度體組合模型(圖3所示)：

目標(biāo)區(qū)域有一個(gè)截面為圓形的二度體，目標(biāo)區(qū)域范圍為：x方向從-1000m到1000m，z方向從0m到1000m(z軸向下為正)；圓心與矩形中心重合，圓形半徑為400m；磁化率為0.01；目標(biāo)區(qū)域地球主磁場為50000nt，磁傾角為45度。將目標(biāo)區(qū)域剖分成10000×10000個(gè)大小相同的小矩形，計(jì)算高度為-50m測線(圖3中虛線所示)上的磁場，計(jì)算點(diǎn)個(gè)數(shù)為10000。

磁場算法程序利用fortran語言實(shí)現(xiàn)，運(yùn)行程序所用的個(gè)人筆記本電腦配置為：cpu為i7-4810，主頻為2.8ghz，內(nèi)存為32gb。運(yùn)行所需時(shí)間約為15秒，由此可見算法效率很高。磁場計(jì)算值和理論值如圖4所示，從形態(tài)上看，兩者是一致的。理論值減去計(jì)算值的差值取絕對(duì)值后除以理論值得到相對(duì)誤差(圖5所示)，對(duì)相對(duì)誤差進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果由表1給出，可知算法精度很高。

表1磁場理論值和計(jì)算值相對(duì)誤差統(tǒng)計(jì)

以上所述僅是本發(fā)明的優(yōu)選實(shí)施方式，本發(fā)明的保護(hù)范圍并不局限于上述實(shí)施例，凡屬于本發(fā)明思路下的技術(shù)方案均屬于本發(fā)明的保護(hù)范圍。應(yīng)當(dāng)指出，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明原理前提下的若干改進(jìn)和潤飾，這些改進(jìn)和潤飾也應(yīng)視為本發(fā)明的保護(hù)范圍。