長方體全張量磁梯度的兩種正演公式對比

鐘煬,管彥武,石甲強(qiáng),肖鋒

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春130026

0 引言

磁法勘探探測技術(shù)的發(fā)展過程大致分為總磁場標(biāo)量測量、磁場三分量測量以及全張量磁梯度測量3個階段[1-2]。當(dāng)前航空磁測主要利用組合式磁梯度計進(jìn)行大工區(qū)的全張量磁梯度測量，全張量磁梯度數(shù)據(jù)較磁標(biāo)量數(shù)據(jù)和三分量矢量數(shù)據(jù)具有抗干擾能力更強(qiáng)、分辨率更高、利于分辨場源磁化方向及小異常體邊界特征等優(yōu)點(diǎn)[2-7]。長方體是正演計算中常用的三度體模型[8-9]，因此推導(dǎo)長方體全張量磁梯度理論表達(dá)式具有重要的理論意義。

針對長方體模型磁場正演問題，近年來許多專家學(xué)者開展了大量研究。郭志宏等人首次指出此前未曾被人發(fā)現(xiàn)或重視的長方體ΔT場表達(dá)式中場值無法計算的解析奇點(diǎn)問題，并從直立長方體模型引力位出發(fā)，在求解積分過程中多次運(yùn)用湊因式法及等價變換法消除導(dǎo)致奇點(diǎn)的多項式，推導(dǎo)出東南下坐標(biāo)系的長方體ΔT場及其梯度場無解析奇點(diǎn)表達(dá)式[9]，但其推導(dǎo)過程在消除奇點(diǎn)的同時也增大了計算量；管志寧在郭志宏等人的研究基礎(chǔ)上給出了長方體磁場三分量的計算公式[8]；為簡化長方體模型正演計算公式，駱遙等人引入歐拉方程對長方體磁場理論表達(dá)式進(jìn)行重新推導(dǎo)，得出了形式更加統(tǒng)一、簡潔的長方體ΔT場及磁場三分量在上半無源空間的無解析奇點(diǎn)理論表達(dá)式[10]；針對前人理論公式僅能計算長方體模型上半無源空間，而無法適應(yīng)起伏地形條件下的正演問題，匡星濤等人利用變量替換法，并單獨(dú)考察讓磁場表達(dá)式無意義的奇點(diǎn)處的磁場值，推導(dǎo)出適用于整個無源空間的長方體ΔT場無解析奇點(diǎn)表達(dá)式[11]；隨著航空全張量磁梯度測量儀器不斷取得重大突破，全張量磁梯度正反演理論也得到了發(fā)展。其中規(guī)則形體的全張量磁梯度正演算法受到廣泛的重視，徐熠通過求解直立長方體模型引力位二階偏導(dǎo)數(shù)，將其帶入泊松公式，求解出長方體全張量磁梯度表達(dá)式[12]；干博、呂文杰分別根據(jù)駱遙等人推導(dǎo)出的長方體磁場三分量表達(dá)式，給出長方體模型全張量磁梯度表達(dá)式[13-14]；修春曉在駱遙等人推導(dǎo)出的長方體磁場三分量表達(dá)式基礎(chǔ)上，給出基于體剖分模型的全張量磁梯度公式及計算效率較高的基于網(wǎng)格節(jié)點(diǎn)物性的長方體全張量磁梯度理論表達(dá)式[15]。

考慮到地球物理勘探常在場源外的上半空間進(jìn)行，故本文在管志寧給出的東南下坐標(biāo)系和駱遙等人給出的北東下坐標(biāo)系的磁場三分量表達(dá)式的基礎(chǔ)上，進(jìn)一步推導(dǎo)了全張量磁梯度的計算公式，并進(jìn)行了對比分析。沿用郭志宏及駱遙等人的長方體模型參數(shù)[10-16]，給出這兩種計算公式的全張量磁梯度正演結(jié)果，它們完全相同。在三維物性反演中，反演的地質(zhì)體往往比較復(fù)雜，通常需要剖分地下模型空間，利用正演公式進(jìn)行迭代計算[17-19]，進(jìn)而擬合觀測異常，達(dá)到反演的目的。例如將地下地質(zhì)體剖分為100×100×10的長方體組合模型，那么在正演公式中每增加一步計算，會增加10萬次的計算[20]。因而有必要選擇更為簡潔的計算公式，為模型正演計算節(jié)約時間，提高反演效率。

1 東南下坐標(biāo)系長方體全張量磁梯度計算公式

全張量磁梯度是磁場強(qiáng)度3個分量(Bx，By，Bz)在空間直角坐標(biāo)系X、Y、Z的3個坐標(biāo)軸方向的變化率。即：

由于磁場強(qiáng)度的旋度為零，所以全張量磁梯度矩陣為對稱矩陣，即式(1)中，Bxy=Byx，Bxz=Bzx，Byz=Bzy；在無源空間中，磁標(biāo)量位滿足拉普拉斯方程，即U=0，故有Bxx+Byy+Bzz=0，所以在全張量磁梯度的9個分量中，只有5個獨(dú)立分量[3]。在實際應(yīng)用中，為了便于表達(dá)全張量磁梯度或驗證磁標(biāo)量位滿足拉普拉斯方程，通常需要計算出上三角矩陣中的6個磁梯度分量。

首先，建立如圖1所示東南下空間直角坐標(biāo)系(X′,Y′,Z′)[8]，其中，X′正軸指向地理東方向，Y′正軸指向地理南方向，Z′軸鉛直向下。

圖1 東南下坐標(biāo)系長方體模型Fig.1 Cuboid model in east-south-down coordinate system

在該坐標(biāo)系中，直立長方體上半無源空間磁場三分量理論表達(dá)式為[16]:

(2)

(3)

(4)

對式(2)～(4)分別沿x′，y′，z′方向求偏導(dǎo)數(shù)，求得在東南下坐標(biāo)系中的長方體全張量磁梯度表達(dá)式，即公式(5)～(10)：

(5)

(6)

(7)

(8)

(9)

(10)

2 北東下坐標(biāo)系長方體全張量磁梯度計算公式

北東下坐標(biāo)系，即NED(North East Down)坐標(biāo)系，在航空航天、測繪和勘探等領(lǐng)域中經(jīng)常使用該坐標(biāo)系[24]。建立NED坐標(biāo)系如圖2所示。

圖2 北東下坐標(biāo)系長方體模型Fig.2 Cuboid model in north-east-down coordinate system

駱遙等人給出在該坐標(biāo)系下直立長方體上半無源空間磁場三分量理論表達(dá)式[10]：

(11)

(12)

(13)

駱遙等人根據(jù)長方體重力場及其梯度滿足構(gòu)造指數(shù)為1的歐拉方程，并將Okabe給出的長方體重力場垂向一階導(dǎo)數(shù)[25]帶入其中，得到不含分子分母同時為零的解析奇點(diǎn)的引力位二階導(dǎo)數(shù)，并將其帶入泊松方程中，得到長方體無解析奇點(diǎn)的磁場三分量理論表達(dá)式。但推導(dǎo)出的引力位二階導(dǎo)數(shù)中反正切項中仍存在分母為零的情況，再利用單變量階梯函數(shù)在對稱區(qū)間的三重積分恒為零的性質(zhì)，對分母為零的反正切函數(shù)進(jìn)行換元，求得形式更加整齊、簡潔的長方體無解析奇點(diǎn)的磁場三分量表達(dá)式(11)～(13)。以上三式將模型角點(diǎn)到觀測點(diǎn)的距離移至積分限中，利于簡化計算，從而在編程時減少冗余的計算步驟。

對以上三式分別沿x，y，z方向求偏導(dǎo)數(shù)，即保持積分上下限不變，對式中的ξ，η，ζ求偏導(dǎo)數(shù)。求得在北東下坐標(biāo)系中的長方體磁場全張量理論表達(dá)式，見式(14)～(19)。對比呂文杰及干博所推導(dǎo)的長方體全張量磁梯度公式[13-14]，他們沒有將模型角點(diǎn)到觀測點(diǎn)的距離移至積分限中，故所求長方體全張量磁梯度公式顯得冗長。徐熠及修春曉所推導(dǎo)的長方體全張量磁梯度公式[12,15]中，均存在未進(jìn)行因式分解至最簡結(jié)果的項。

(14)

(15)

(16)

(17)

(18)

(19)

為方便描述，簡稱公式(14)～(19)為算法二。

這兩種算法都是先計算出長方體引力位二階導(dǎo)數(shù)，再根據(jù)泊松方程計算磁場三分量，進(jìn)一步求方向?qū)?shù)，獲得全張量磁梯度表達(dá)式。但是在求解引力位二階導(dǎo)數(shù)的過程中，兩者采用了不同的方法。算法一通過等價變換及湊因式法避免了前人[21-22]在求解積分過程中引入的奇點(diǎn)問題；算法二則是引用了Okabe給出的長方體重力場公式[25]，并帶入到構(gòu)造指數(shù)為1的歐拉方程中求解直立長方體引力位二階導(dǎo)數(shù)。因而兩種算法得到了形式不盡相同的長方體磁場三分量表達(dá)式，進(jìn)而全張量磁梯度的表達(dá)式也簡繁不同。

算法二即為北東下坐標(biāo)系下長方體全張量磁梯度理論表達(dá)式，對比算法一。兩種算法需要的參數(shù)個數(shù)相同，但算法二的公式形式相對簡潔；算法二采用的坐標(biāo)系更符合常規(guī)，且磁化偏角與地磁偏角的定義一致，便于使用；算法二設(shè)定的積分限在簡化計算方面的優(yōu)勢也繼續(xù)得到體現(xiàn)。為進(jìn)一步對比兩種算法計算量的差異，以同一長方體模型，計算觀測面上同一點(diǎn)為例，統(tǒng)計了兩種算法中每個磁梯度分量的計算量(表1)。其中，算法二全張量磁梯度公式加減法計算總量不到算法一的四分之一；乘除法計算總量約為算法一的一半。

表1 兩種算法全張量磁梯度計算次數(shù)

Table 1 Calculation times of full tensor magnetic gradient with two algorithms

FTMG分量加減法計算次數(shù)乘除法計算次數(shù)算法一算法二算法一算法二Bxx1926013284Bxy2143014544Bxz1933012644Byy1906013074Byz1893012444Bzz116609474總計1 094270751370

3 長方體全張量磁梯度模型計算

為檢驗在兩種坐標(biāo)系下導(dǎo)出的長方體全張量磁梯度理論表達(dá)式的正確性，對算法一、算法二取完全相同的長方體模型進(jìn)行正演計算，模型參數(shù)沿用郭志宏、駱遙等人的參數(shù)[10、16]。圖3a～3f為利用算法一計算得到的全張量磁梯度等值線圖，為方便對比，已轉(zhuǎn)換到常用的北東下坐標(biāo)系中顯示。坐標(biāo)轉(zhuǎn)換方式為：

模型參數(shù)為：磁化方向I=50°，A=30°;磁化強(qiáng)度M=1 A/m;地磁場方向I0=60°，A0=10°;長方體中心埋深x0=10 000 m，y0=10 000 m，z0=3 000 m;長方體長a=8 000 m，寬b=4 000 m，高c=1 000 m;網(wǎng)格間距為100 m×100 m;計算高度z=0 m;等值線單位為nT/m。

圖4為模型在北東下坐標(biāo)系中利用算法二計算得到的全張量磁梯度等值線圖。對比兩圖，可見同一長方體模型，兩種算法的全張量磁梯度等值線圖形態(tài)完全一致。為進(jìn)一步對比兩坐標(biāo)系下全張量磁梯度差異，用圖4減去圖3所對應(yīng)的張量值，并繪制成等值線圖(圖5)。

圖5表明，兩種計算公式得到的全張量磁梯度差值的量級均在10-17～10-15nT/m之間，相比于當(dāng)前精度較高的航空全張量磁梯度測量的實際觀測精度為10-2nT/m[26]，這兩種算法的計算結(jié)果可視為無差別。此外，為對比計算效率，筆者建立了一個地下剖分為100×100×10的長方體組合模型，測區(qū)范圍20 000 m×20 000 m，點(diǎn)線距均為100 m，X和Y方向網(wǎng)格點(diǎn)數(shù)均為201個。所用計算機(jī)CPU為Intel(R)Core(TM)i7-8750H，內(nèi)存(RAM)8.0GB，Matlab版本2018a，算法一耗時1 896.5 s；而算法二耗時1 703.7 s。比較程序計算用時，算法二較算法一節(jié)省約10.2%的時間，可見兩者計算量相差較大。主要原因是兩種計算方法各張量表達(dá)式差異，算法二表達(dá)式相對簡潔。

圖3 由算法一計算的長方體全張量磁梯度等值線圖Fig.3 Contour map of full tensor magnetic gradient of cuboid with algorithm 1

圖4 由算法二計算的長方體全張量磁梯度等值線圖(模型參數(shù)與圖3相同)Fig.4 Contour map of full tensor magnetic gradient of cuboid with algorithm 2

圖5 圖4與圖3對應(yīng)張量之差等值線圖Fig.5 Contour map of the difference between Fig.4 and Fig.3

4 結(jié)論

(1)分別推導(dǎo)了東南下和北東下兩種坐標(biāo)系下長方體模型的全張量磁梯度理論表達(dá)式，經(jīng)驗證兩種算法的計算結(jié)果完全相同。

(2)從推導(dǎo)過程及公式形式上對比了兩種算法，它們采用了不同的推導(dǎo)過程，故公式形式不盡相同，但獲得了相同的計算結(jié)果；算法二中坐標(biāo)系和磁偏角的定義更符合常規(guī)，便于使用；算法二計算效率高于算法一，在模型試算中，算法二節(jié)省10.2%的計算時間。