朱自強，邢澤峰，魯光銀

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083)

0 引言

我國為多山地區(qū)，地形起伏較大，而不同地形條件對重力勘察工作有很大的影響。此外，地形校正的精度對后續(xù)重力資料的進(jìn)一步處理解釋影響顯著[1-2]。目前，國內(nèi)、外學(xué)者針對重力地形校正這一課題已進(jìn)行了很多相關(guān)的研究工作，其中廣泛應(yīng)用于生產(chǎn)實踐的地形校正方法包括圓域地形校正法、方形域地形校正法以及表面積分法地形校正[3-5]，其中表面積分法由于其對地形近似程度高，往往可以取得令人滿意的校正精度。Bychkov等[6]利用線性解析近似的方法來進(jìn)行地形校正，獲得較好的地形校正結(jié)果；張偉等[7]提出基于三角網(wǎng)扣合的山區(qū)重力測量高精度近區(qū)地改算法，其計算精度可達(dá)“微伽級”；Cella[8]基于MATLAB開發(fā)地形校正處理程序，可較精確地實現(xiàn)地形校正，且可用于陸地海面等不同勘探環(huán)境。目前，我國的勘探環(huán)境越來越復(fù)雜，重力資料處理中，在地下地質(zhì)體密度變化較大且地下異常體引起的重力異常信號較微弱的情況下，如果校正項出現(xiàn)問題，將直接影響后期數(shù)據(jù)處理與解釋工作[9-10]。因此，開展任意復(fù)雜重力地形校正方法研究具有重要意義。

有限單元法因能夠準(zhǔn)確模擬復(fù)雜情況下地球物理場的分布[11]，被廣泛應(yīng)用于地球物理數(shù)值模擬研究工作。徐世浙等[12]通過有限單元法來求解重力場，其結(jié)果具有較高的精度和計算速率。筆者嘗試采用有限單元法來研究任意復(fù)雜地形條件下有限元重力地形校正問題。采用有限單元法進(jìn)行正演模擬時，不可避免的會遇到由于計算區(qū)域外地形體缺失而導(dǎo)致的重力異常問題。為此，這里提出一種重力補償方法，以解決由于計算區(qū)域外地形體缺失而導(dǎo)致的重力虛假異常。基于改進(jìn)的有限元地形體正演模擬算法，對任意地形體進(jìn)行重力正演模擬，從而獲得任意地形影響下在地表所引起的重力異常響應(yīng)曲線，在此基礎(chǔ)上開展地形校正，以獲取任意復(fù)雜地形條件下有限元重力地形校正結(jié)果。

1 改進(jìn)的有限元正演計算方法

1.1 邊值問題

采用有限元法計算重力異常時，通常將邊界選在遠(yuǎn)離目標(biāo)體的地方，如果區(qū)域的邊界取的足夠遠(yuǎn)，場源對邊界的影響可以近似看作質(zhì)量集中于截面質(zhì)心處質(zhì)線的影響[12]。求解區(qū)域內(nèi)的重力g就等同于求解下列邊值問題[11]：

(1)

其中:G是萬有引力常數(shù)；σ是密度差；θ為r與水平線的夾角；α是邊界外法線n與矢徑r的夾角。

1.2 變分問題

與邊值問題式(1)等價的變分問題[11]如式(2)所示。

δF(g)=0

(2)

其中：Ω為計算區(qū)域；?！逓闊o窮遠(yuǎn)邊界。

1.3 單元分析

對待求解區(qū)域Ω進(jìn)行三角單元剖分，在三角元中進(jìn)行線性插值。通過網(wǎng)格剖分，把積分區(qū)域剖分為許多三角形單元，從而可以把式(2)中的積分分解為各三角單元e上的積分。計算單元e的積分得

(3)

假設(shè)三角元e內(nèi)的密度ρ均勻，則有

(4)

如果三角形的某一條邊L落在?！奚?，則

(5)

1.4 整體合成

在每個三角元內(nèi)，將式(3)、式(4)和式(5)相加，再擴展由全體節(jié)點組成的矩陣或陣列，進(jìn)而全部單元相加得式(6)。

圖1 改進(jìn)的重力有限元正演算法流程

(6)

其中，g是全部節(jié)點組成的列向量。對式(6)求變分，并令其等于0，得線性方程組:

Kg=P

(7)

解線性方程組，即可得到各節(jié)點的重力異常g。

1.5 改進(jìn)的有限元地形體算法

基于上述有限元法對二度地形體模型進(jìn)行正演模擬，由于計算范圍的限制，只能考慮到預(yù)先設(shè)定計算范圍內(nèi)地形體在地表引起的重力響應(yīng)，導(dǎo)致計算區(qū)域外重力響應(yīng)數(shù)據(jù)缺失，從而使地形體重力響應(yīng)曲線產(chǎn)生畸變。針對上述問題，提出一種改進(jìn)的有限元地形體正演模擬算法，對缺失的重力響應(yīng)做出相對應(yīng)的重力補償，以消除重力響應(yīng)曲線畸變影響。改進(jìn)的算法流程見圖1。

2 正演算法精度驗證

對地下二度密度異常球體進(jìn)行有限元正演模擬，分析有限元正演模擬結(jié)果與理論重力異常間的相對誤差，以驗證有限元正演算法的精度。此外，利用改進(jìn)的有限元地形體正演算法對水平地形體模型進(jìn)行正演模擬，并將計算結(jié)果與理論曲線進(jìn)行對比，分析驗證該改進(jìn)的有限元地形體正演算法有效性。

2.1 網(wǎng)格分剖

本文中的網(wǎng)格剖分利用了Per-Olof Persson所提供的基于MATLAB所編寫的三角網(wǎng)格剖分開源代碼DISTMESH，可通過調(diào)整網(wǎng)格剖分函數(shù)各項參數(shù)來達(dá)到控制每個三角網(wǎng)格剖分范圍，網(wǎng)格密度以及三角網(wǎng)格形狀等網(wǎng)格參數(shù)的目的。因此，可基于MATLAB利用DISTMESH開源代碼進(jìn)行相關(guān)模型構(gòu)建，設(shè)置任意幾何形狀的地下異常體模型，并可實現(xiàn)每個網(wǎng)格參數(shù)獨立賦值。

圖2 模型1網(wǎng)格剖分圖

圖3 二度球體有限元正演模擬結(jié)果與理論曲線間擬合情況

圖4 正演模擬結(jié)果與理論曲線間相對誤差曲線

圖5 模型2網(wǎng)格剖分圖

圖6 水平地形體不同正演模擬方法效果對比

2.2 有限元正演算法精度驗證

設(shè)計模型1(圖2)，假設(shè)地下密度異常體為二度球體，對該二度球體進(jìn)行正演模擬，將正演模擬結(jié)果與理論曲線對比進(jìn)行誤差分析。模型1參數(shù)設(shè)置為，二度球體中心坐標(biāo)(0，-50)，半徑為25 m。假設(shè)該異常體密度均勻為2 g/cm3，背景物質(zhì)密度設(shè)置為0 g/cm3，模型網(wǎng)格剖分如圖2所示。

對模型1正演計算結(jié)果與其理論值進(jìn)行對比分析，從圖3中可看出，正演模擬結(jié)果和理論曲線形態(tài)相似度很高，擬合情況較好，其相對誤差最大值小于3×10-4(圖4)，推測其誤差主要為三角形模擬近似二度球體所引起的?？傮w來看，該有限元正演模擬程序穩(wěn)定性好，精度較高。

2.3 水平地形體模型驗證

設(shè)計模型2(圖5)，假設(shè)水平地形體沿水平面無限延伸，地形體垂向厚度為100 m，密度均勻且密度為1 g/cm3。模型2網(wǎng)格剖分如圖5所示。

對模型2正演計算結(jié)果與其理論值進(jìn)行對比分析，從圖6中可看出，未經(jīng)過重力補償?shù)挠邢拊匦误w重力響應(yīng)曲線與理論重力響應(yīng)曲線間存在明顯的誤差，且在地形體兩端誤差最大；而經(jīng)過重力補償后的有限元地形體正演模擬結(jié)果與理論重力響應(yīng)曲線擬合情況較好，誤差較小，證明了該改進(jìn)的重力有限元地形體正演模擬算法的有效性。

3 改進(jìn)的有限元任意復(fù)雜地形校正

利用上述改進(jìn)的有限元地形體正演模擬算法，計算任意地形影響下地下密度異常體在地表所引起的重力異常響應(yīng)，并進(jìn)行有限元地形校正。

3.1 有限元任意簡單地形校正

圖7 模型3網(wǎng)格剖分圖

圖8 任意簡單地形影響下多個密度異常體有限元正演模擬結(jié)果

設(shè)計模型3，假設(shè)地下存在多個密度異常體，任意地形體沿水平面無限延伸，從負(fù)無窮到-200，及200到正無窮范圍內(nèi)，地形體垂向厚度相同，設(shè)置為100 m。在-200至200范圍內(nèi)設(shè)置為任意起伏地形，該地形體密度均勻且密度為2 g/cm3，模型3參數(shù)設(shè)置：二度球體中心坐標(biāo)(-10，-30)，半徑為20 m；二度長方體左下角坐標(biāo)(-115，-30)，右上角坐標(biāo)為(-85,0)。二度長方體和二度球體密度均勻，與地形體相比其剩余密度分別設(shè)置為mg/cm3、ng/cm3(m,n為剩余密度的大小)。模型3網(wǎng)格剖分如圖7所示：

利用改進(jìn)的有限元地形體正演算法對上述任意簡單地形體模型進(jìn)行正演模擬，分析其在地表所引起的重力響應(yīng)曲線分布規(guī)律。圖8中，設(shè)置m=0，n=0，任意簡單地形體密度均勻且密度為2 g/cm3，此時地表所引起的重力響應(yīng)曲線隨著地形起伏產(chǎn)生相對應(yīng)的波動，且與地形體起伏趨勢一致。地下存在密度異常體時，在異常體上方地表投影點附近重力響應(yīng)曲線出現(xiàn)相對應(yīng)向上或向下異常突出，與地形體產(chǎn)生的重力響應(yīng)曲線混合，使地表所觀測到的重力異常曲線變得極為復(fù)雜，增加后期解釋工作的難度。

根據(jù)任意簡單地形影響下地下多個密度異常體正演模擬結(jié)果，進(jìn)一步進(jìn)行有限元地形校正，得到圖9中重力響應(yīng)曲線。基于有限元重力地形校正結(jié)果，再進(jìn)行中間層校正，得到圖10中地形及中間層校正后的重力響應(yīng)曲線。從圖10可明顯看出，不同剩余密度重力響應(yīng)曲線的極值點分布在x=-100及x=-10附近，與所設(shè)置模型參數(shù)相符，重力響應(yīng)曲線可更加直觀可靠地反映地下密度異常體的分布情況。

圖9 有限元地形校正后重力異常響應(yīng)曲線

圖10 有限元中間層校正后重力異常響應(yīng)曲線

圖11 模型4網(wǎng)格剖分圖

3.2 有限元任意復(fù)雜地形校正

設(shè)計模型4，任意復(fù)雜地形體沿水平面無限延伸，從負(fù)無窮到-200，以及200到正無窮范圍內(nèi)，地形體垂向厚度相同，設(shè)置為100 m。在-200至200范圍內(nèi)設(shè)置為任意復(fù)雜起伏地形。地形體呈層狀分布。各層密度不同；單一層內(nèi)物質(zhì)分布均勻，密度相同。模型4參數(shù)設(shè)置：地形體上層密度為2 g/cm3，下層密度設(shè)置為3 g/cm3；二度球體中心坐標(biāo)(-10，-30)，半徑為20 m；二度長方體左下角坐標(biāo)(-115，-30)，右上角坐標(biāo)為(-85,0)。二度長方體和二度球體密度均勻，與上層地形體相比其剩余密度分別設(shè)置為mg/cm3、ng/cm3(m,n為剩余密度的大小)。模型4網(wǎng)格剖分如圖11所示。

圖12 任意復(fù)雜地形影響下多個密度異常體有限元正演模擬結(jié)果

圖13 有限元任意復(fù)雜地形校正后重力異常響應(yīng)曲線

圖14 有限元任意復(fù)雜地形中間層校正后重力異常響應(yīng)曲線

利用改進(jìn)的有限元地形體正演算法對上述任意復(fù)雜地形體模型進(jìn)行正演模擬，得到地下多個異常體在地表所引起的重力響應(yīng)曲線(圖12)，當(dāng)m=n=0時，地下無密度異常體。當(dāng)?shù)叵麓嬖诙鄠€密度異常體時，復(fù)雜地形體與多個地下密度異常體相互影響，很大程度上增大了地表所觀測到重力響應(yīng)曲線復(fù)雜性。

根據(jù)任意復(fù)雜地形影響下地下多個密度異常體正演結(jié)果，進(jìn)行有限元地形校正。從圖14可看出，任意復(fù)雜地形影響下經(jīng)有限元地形校正及中間層校正后重力響應(yīng)曲線可更加直觀地反映地下密度異常體分布狀態(tài)，可明顯區(qū)分兩個異常體的位置分布，驗證了該有限元地形校正算法的有效性。

4 結(jié)論

利用筆者提出改進(jìn)的重力有限元任意地形模擬算法對任意地形影響下地下多個密度異常體進(jìn)行正演模擬，獲得地表處的重力響應(yīng)曲線，經(jīng)校正后的重力響應(yīng)曲線更加直觀地反映出地下密度異常體分布狀態(tài)，驗證了該有限元任意復(fù)雜地形校正算法的有效性。

基于有限元法可建立任意形狀、任意密度復(fù)雜地下模型。在保證計算精度的情況下，可快速計算出重力任意復(fù)雜地形校正結(jié)果，有限元重力任意復(fù)雜地形校正方法尤其適合于地下地質(zhì)體密度變化較大的復(fù)雜地形。