有限元重力任意復(fù)雜地形校正方法研究

朱自強(qiáng)，邢澤峰，魯光銀

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

0 引言

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

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

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

1.1 邊值問題

采用有限元法計(jì)算重力異常時(shí)，通常將邊界選在遠(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)等價(jià)的變分問題[11]如式(2)所示。

δF(g)=0

(2)

其中：Ω為計(jì)算區(qū)域；Γ∞為無窮遠(yuǎn)邊界。

1.3 單元分析

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

(3)

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

(4)

如果三角形的某一條邊L落在Γ∞上，則

(5)

1.4 整體合成

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

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

(6)

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

Kg=P

(7)

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

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

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

2 正演算法精度驗(yàn)證

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

2.1 網(wǎng)格分剖

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

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

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

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

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

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

2.2 有限元正演算法精度驗(yàn)證

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

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

2.3 水平地形體模型驗(yàn)證

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

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

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

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

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

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

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

設(shè)計(jì)模型3，假設(shè)地下存在多個(gè)密度異常體，任意地形體沿水平面無限延伸，從負(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，此時(shí)地表所引起的重力響應(yīng)曲線隨著地形起伏產(chǎn)生相對應(yīng)的波動(dòng)，且與地形體起伏趨勢一致。地下存在密度異常體時(shí)，在異常體上方地表投影點(diǎn)附近重力響應(yīng)曲線出現(xiàn)相對應(yīng)向上或向下異常突出，與地形體產(chǎn)生的重力響應(yīng)曲線混合，使地表所觀測到的重力異常曲線變得極為復(fù)雜，增加后期解釋工作的難度。

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

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

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

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

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

設(shè)計(jì)模型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ù)雜地形影響下多個(gè)密度異常體有限元正演模擬結(jié)果

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

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

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

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

4 結(jié)論

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

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