本發(fā)明涉及一種二度體重力異常計(jì)算方法，特別是任意截面形狀、任意密度分布的二度體重力異?？焖佟⒏呔扔?jì)算方法。

背景技術(shù)：

重力勘探以地球重力場(chǎng)為根據(jù)、密度差異為物理基礎(chǔ)，是一種研究地球構(gòu)造及尋找礦產(chǎn)資源的地球物理勘探方法，該方法效率高、操作簡(jiǎn)單、成本低、勘探深度大、實(shí)施過程沒有過多條件限制。重力異常正演計(jì)算是指根據(jù)密度分布計(jì)算重力異常，反演是指根據(jù)觀測(cè)重力值計(jì)算密度分布。正演是反演的基礎(chǔ)，正演計(jì)算的效率直接影響反演計(jì)算的效率，而正演計(jì)算精度直接影響反演結(jié)果的質(zhì)量。一直以來(lái)，學(xué)者們都非常重視重力異常的正演計(jì)算。實(shí)際情況中，常見這樣一類線性地質(zhì)體：其走向方向的尺度要遠(yuǎn)比垂直其走向方向的尺度大，諸如此類地質(zhì)體，實(shí)際場(chǎng)源分布便可以用走向方向無(wú)限延伸的二度體代替。

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

對(duì)于任意截面形狀、任意密度分布的二度體重力異常計(jì)算，一般采用數(shù)值方法。文獻(xiàn)(張嶺,郝天珧.基于Delaunay剖分的二維非規(guī)則重力建模及重力計(jì)算.地球物理學(xué)報(bào),2006,49(3):877-884.)公開了一種Delaunay剖分方式，將截面分為若干三角形，將二度體分解為若干三角棱柱的組合，利用解析公式計(jì)算變密度三角棱柱重力異常，并將其累加，這種方法保證了計(jì)算精度，但計(jì)算效率低；文獻(xiàn)(朱自強(qiáng),曾思紅,魯光銀等.二度體的重力張量有限元正演模擬.物探與化探,2010,34(5):668-671.)公開了一種計(jì)算二度體重力梯度張量的有限元方法；文獻(xiàn)(Reeder K,Louie J,Kent G.Efficient 2D finite element gravity modeling using convolution.SEG Technical Program Expanded Abstracts,2014:1254-1258.)公開了一種有限元方法和卷積算法相結(jié)合的二度體重力異常計(jì)算方法。采用有限元方法計(jì)算二度體重力異常，能夠較精確刻畫復(fù)雜二度體，計(jì)算精度高但計(jì)算效率很低。

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

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

本發(fā)明針對(duì)目前現(xiàn)有二度體重力異常計(jì)算方法不能同時(shí)保證計(jì)算效率和計(jì)算精度，無(wú)法滿足大規(guī)模重力異常密度精細(xì)反演成像、人機(jī)交互建模和解釋的需求問題，本發(fā)明其目的在于提出一種二度體重力異常計(jì)算方法，其是一種任意截面形狀、任意密度分布的二度體重力異?？焖佟⒏呔扔?jì)算方法。

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

一種二度體重力異常計(jì)算方法，包括以下步驟：

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

確定觀測(cè)點(diǎn)，觀測(cè)點(diǎn)坐標(biāo)(x_m,z₀)，基于觀測(cè)點(diǎn)建立包含所有目標(biāo)區(qū)域的矩形模型，確定矩形模型在x，z方向的起始位置，使得包含起伏地形的目標(biāo)區(qū)域完全嵌入在該矩形模型中；

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

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

第二步矩形模型重力異常計(jì)算；

第一步中給出的矩形模型，其重力異常計(jì)算公式為

式中，(x_m,z₀)表示觀測(cè)點(diǎn)坐標(biāo)，z₀為常值；L表示矩形模型在z方向剖分小矩形的個(gè)數(shù)；M表示矩形模型在x方向剖分小矩形的個(gè)數(shù)；(ξ_p,ζ_r)表示編號(hào)為(p,r)的小矩形幾何中心坐標(biāo)即在矩形模型中x方向?yàn)榈趐個(gè)且在z方向?yàn)榈趓個(gè)的小矩形的幾何中心坐標(biāo)；ρ(ξ_p,ζ_r)表示編號(hào)為(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加權(quán)系數(shù)。

進(jìn)一步地，第二步中，矩形模型其重力異常計(jì)算公式的解算方法如下：

a，根據(jù)觀測(cè)點(diǎn)坐標(biāo)(x_m,z₀)和小矩形幾何中心坐標(biāo)(ξ_p,ζ_r)，計(jì)算加權(quán)系數(shù)h(x_m-ξ_p,z₀-ζ_r)，其計(jì)算公式為

式中，γ表示萬(wàn)有引力常數(shù)，arctan()表示反余切函數(shù)運(yùn)算符，ln()表示自然對(duì)數(shù)運(yùn)算符；其它符號(hào)含義如下

X₁＝ξ_p-0.5Δx-x_m，X₂＝ξ_p+0.5Δx-x_m，Z₁＝ζ_r-0.5Δz-z₀，Z₂＝ζ_r+0.5Δz-z₀

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

b，在矩形模型中在x方向上位于同一行的所有小矩形為一層。采用一維離散卷積快速計(jì)算方法來(lái)計(jì)算一層矩形模型重力異常，其計(jì)算公式為

式中，表示第r層(r＝1,2,…,L)矩形模型在測(cè)線z₀產(chǎn)生的重力異常；(x_m,z₀)表示觀測(cè)點(diǎn)坐標(biāo)；

c，將各層矩形模型重力異常進(jìn)行累加，得到整個(gè)矩形模型的重力異常，即

在第二步中的步驟b中，采用一維離散卷積快速計(jì)算方法來(lái)計(jì)算一層矩形模型重力異常，其步驟為：

(1)將加權(quán)系數(shù)h(x₁-ξ_p,z₀-ζ_r)排列成向量t，記為

t＝[t₀,t₁,t₂,…,t_M-1,0,t_M-1,t_M-2,…,t₂,t₁]^T (5)

式中，矩陣元素t_i與加權(quán)系數(shù)h(x₁-ξ_p,z₀-ζ_r)存在關(guān)系

t_i＝h(x₁-ξ_i+1,z₀-ζ_r) (6)

(2)將第r層密度值ρ(ξ_p,ζ_r)(p＝1,2,…,M)排列成向量ρ，向量元素ρ_i與密度值存在關(guān)系

ρ_i＝ρ(ξ_i,ζ_r) (7)

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

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

(3)計(jì)算

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

(4)計(jì)算

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

(5)計(jì)算

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

(6)提取矩陣g^ext的前M行元素，構(gòu)成向量g，即為一維離散卷積計(jì)算結(jié)果，向量g中元素就是所求

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

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

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

(3)大規(guī)模正演計(jì)算時(shí)，算法不但計(jì)算效率和計(jì)算精度高，并且所需計(jì)算機(jī)內(nèi)存小。

附圖說明

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

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

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

圖4為重力異常計(jì)算值和理論值對(duì)比圖；

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

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

ρ：表示密度。

具體實(shí)施方式

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

參照?qǐng)D1，為本實(shí)施例一種二度體重力異常計(jì)算方法的流程圖，包括以下步驟：

1、復(fù)雜二度體模型表示：

首先，確定觀測(cè)點(diǎn)，觀測(cè)點(diǎn)坐標(biāo)(x_m,z₀)，基于觀測(cè)點(diǎn)建立包含所有目標(biāo)區(qū)域的規(guī)則矩形模型，確定矩形模型在x，z方向的起始位置，使得目標(biāo)區(qū)域(包含起伏地形)完全嵌入在該矩形模型中；

其次，根據(jù)實(shí)際問題需求即根據(jù)預(yù)定的針對(duì)目標(biāo)區(qū)域確定的小矩形剖分個(gè)數(shù)，將矩形模型均勻劃分成多個(gè)規(guī)則小矩形(如圖2所示)，確定小矩形的幾何尺寸Δx，Δz；

最后，根據(jù)已知的目標(biāo)區(qū)域的密度分布，對(duì)每個(gè)小矩形密度進(jìn)行賦值，根據(jù)小矩形和二度體截面的幾何關(guān)系，當(dāng)小矩形中心位于二度體截面內(nèi)時(shí)，賦值給所在位置的二度體密度，每個(gè)小矩形密度為常值，不同小矩形密度取值不同，以此刻畫任意截面形狀、任意密度分布的二度體。位于空氣部分的小矩形，其密度值設(shè)為零，以此刻畫起伏地形。

2、矩形模型重力異常計(jì)算：

步驟一中給出的矩形模型，其重力異常計(jì)算公式為

式中，(x_m,z₀)表示觀測(cè)點(diǎn)坐標(biāo)，z₀為常值；L表示矩形模型在z方向剖分小矩形的個(gè)數(shù)；M表示矩形模型在x方向剖分小矩形的個(gè)數(shù)；(ξ_p,ζ_r)表示編號(hào)為(p,r)的小矩形幾何中心坐標(biāo)即在矩形模型中x方向?yàn)榈趐個(gè)且在z方向?yàn)榈趓個(gè)的小矩形的幾何中心坐標(biāo)；ρ(ξ_p,ζ_r)表示編號(hào)為(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加權(quán)系數(shù)。

實(shí)現(xiàn)式(1)的計(jì)算，分為三個(gè)環(huán)節(jié)：

首先，根據(jù)觀測(cè)點(diǎn)坐標(biāo)(x_m,z₀)和小矩形幾何中心坐標(biāo)(ξ_p,ζ_r)，計(jì)算加權(quán)系數(shù)h(x_m-ξ_p,z₀-ζ_r)，其計(jì)算公式為

式中，γ表示萬(wàn)有引力常數(shù)，arctan()表示反余切函數(shù)運(yùn)算符，ln()表示自然對(duì)數(shù)運(yùn)算符；其它符號(hào)含義如下

X₁＝ξ_p-0.5Δx-x_m，X₂＝ξ_p+0.5Δx-x_m，Z₁＝ζ_r-0.5Δz-z₀，Z₂＝ζ_r+0.5Δz-z₀

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

其次，在矩形模型中在x方向上位于同一行的所有小矩形為一層，采用一維離散卷積快速計(jì)算方法來(lái)計(jì)算一層矩形模型重力異常，其計(jì)算公式為

式中，表示第r層(r＝1,2,…,L)矩形模型在測(cè)線z₀產(chǎn)生的重力異常；(x_m,z₀)表示觀測(cè)點(diǎn)坐標(biāo)；

最后，將各層矩形模型重力異常進(jìn)行累加，得到整個(gè)矩形模型的重力異常，即

步驟二中所述的一維離散卷積快速計(jì)算方法，其步驟為：

(1)將加權(quán)系數(shù)h(x₁-ξ_p,z₀-ζ_r)排列成向量t，記為

t＝[t₀,t₁,t₂,…,t_M-1,0,t_M-1,t_M-2,…,t₂,t₁]^T (5)

式中，矩陣元素t_i與加權(quán)系數(shù)h(x₁-ξ_p,z₀-ζ_r)存在關(guān)系

t_i＝h(x₁-ξ_i+1,z₀-ζ_r) (6)

(2)將第r層密度值ρ(ξ_p,ζ_r)(p＝1,2,…,M)排列成向量ρ，向量元素ρ_i與密度值存在關(guān)系

ρ_i＝ρ(ξ_i,ζ_r) (7)

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

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

(3)計(jì)算

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

(4)計(jì)算

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

(5)計(jì)算

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

(6)提取矩陣g^ext的前M行元素，構(gòu)成向量g，即為一維離散卷積計(jì)算結(jié)果，向量g中元素就是所求

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

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

密度為常值的矩形區(qū)域內(nèi)嵌一個(gè)密度均勻的圓形，圓心與矩形中心重合。矩形范圍為：x方向從-1000m到1000m，z方向從0m到1000m(z軸向下為正)；圓形半徑為400m。矩形密度為0g/cm³，圓形的密度為1g/cm³。將矩形剖分成10000×10000個(gè)大小相同的小矩形，計(jì)算高度為-50m測(cè)線(圖3中虛線所示)上的重力異常，計(jì)算點(diǎn)個(gè)數(shù)為10000。

重力異常算法利用Fortran語(yǔ)言編程實(shí)現(xiàn)，運(yùn)行程序所用的個(gè)人臺(tái)式機(jī)配置為：CPU為i7-2620，主頻為2.7GHz，內(nèi)存為32GB，四核八線程。運(yùn)行所需時(shí)間約為6秒，由此可見算法效率很高。重力異常計(jì)算值和理論值如圖4所示，從形態(tài)上看，兩者是一致的。相對(duì)誤差由理論值減去計(jì)算值得到差值的絕對(duì)值除以理論值得到(圖5所示)，對(duì)相對(duì)誤差進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果由表1給出，可知算法精度很高。

表1重力異常理論值和計(jì)算值相對(duì)誤差統(tǒng)計(jì)

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

以上包含了本發(fā)明優(yōu)選實(shí)施例的說明，這是為了詳細(xì)說明本發(fā)明的技術(shù)特征，并不是想要將發(fā)明內(nèi)容限制在實(shí)施例所描述的具體形式中，依據(jù)本發(fā)明內(nèi)容主旨進(jìn)行的其他修改和變型也受本專利保護(hù)。本發(fā)明內(nèi)容的主旨是由權(quán)利要求書所界定，而非由實(shí)施例的具體描述所界定。