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

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

0 引言

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

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

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

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

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

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

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

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

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

(2)

(3)

(4)

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

(5)

(6)

(7)

(8)

(9)

(10)

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

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

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

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

(11)

(12)

(13)

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

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

(14)

(15)

(16)

(17)

(18)

(19)

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

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

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

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

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

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

3 長(zhǎng)方體全張量磁梯度模型計(jì)算

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

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

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

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

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

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

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

4 結(jié)論

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

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