本發(fā)明涉及一種二度體重力異常計算方法，特別是任意截面形狀、任意密度分布的二度體重力異常快速、高精度計算方法。

背景技術：

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

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

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

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

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

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

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

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

第一步復雜二度體模型表示：

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

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

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

第二步矩形模型重力異常計算；

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

式中，(x_m,z₀)表示觀測點坐標，z₀為常值；L表示矩形模型在z方向剖分小矩形的個數；M表示矩形模型在x方向剖分小矩形的個數；(ξ_p,ζ_r)表示編號為(p,r)的小矩形幾何中心坐標即在矩形模型中x方向為第p個且在z方向為第r個的小矩形的幾何中心坐標；ρ(ξ_p,ζ_r)表示編號為(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加權系數。

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

a，根據觀測點坐標(x_m,z₀)和小矩形幾何中心坐標(ξ_p,ζ_r)，計算加權系數h(x_m-ξ_p,z₀-ζ_r)，其計算公式為

式中，γ表示萬有引力常數，arctan()表示反余切函數運算符，ln()表示自然對數運算符；其它符號含義如下

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方向上位于同一行的所有小矩形為一層。采用一維離散卷積快速計算方法來計算一層矩形模型重力異常，其計算公式為

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

c，將各層矩形模型重力異常進行累加，得到整個矩形模型的重力異常，即

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

(1)將加權系數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與加權系數h(x₁-ξ_p,z₀-ζ_r)存在關系

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

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

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

將向量ρ補零擴展成向量ρ^ext

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

(3)計算

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

(4)計算

式中，“.*”表示對應元素相乘運算；

(5)計算

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

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

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

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

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

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

附圖說明

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

圖2為復雜二度體模型表示；

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

圖4為重力異常計算值和理論值對比圖；

圖5為計算值與理論值的相對誤差；

圖中符號說明如下：

ρ：表示密度。

具體實施方式

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

參照圖1，為本實施例一種二度體重力異常計算方法的流程圖，包括以下步驟：

1、復雜二度體模型表示：

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

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

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

2、矩形模型重力異常計算：

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

式中，(x_m,z₀)表示觀測點坐標，z₀為常值；L表示矩形模型在z方向剖分小矩形的個數；M表示矩形模型在x方向剖分小矩形的個數；(ξ_p,ζ_r)表示編號為(p,r)的小矩形幾何中心坐標即在矩形模型中x方向為第p個且在z方向為第r個的小矩形的幾何中心坐標；ρ(ξ_p,ζ_r)表示編號為(p,r)的小矩形的密度值；h(x_m-ξ_p,z₀-ζ_r)表示加權系數。

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

首先，根據觀測點坐標(x_m,z₀)和小矩形幾何中心坐標(ξ_p,ζ_r)，計算加權系數h(x_m-ξ_p,z₀-ζ_r)，其計算公式為

式中，γ表示萬有引力常數，arctan()表示反余切函數運算符，ln()表示自然對數運算符；其它符號含義如下

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方向上位于同一行的所有小矩形為一層，采用一維離散卷積快速計算方法來計算一層矩形模型重力異常，其計算公式為

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

最后，將各層矩形模型重力異常進行累加，得到整個矩形模型的重力異常，即

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

(1)將加權系數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與加權系數h(x₁-ξ_p,z₀-ζ_r)存在關系

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

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

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

將向量ρ補零擴展成向量ρ^ext

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

(3)計算

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

(4)計算

式中，“.*”表示對應元素相乘運算；

(5)計算

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

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

下面對本發(fā)明方法的效果進行檢驗。

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

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

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

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

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

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