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

鐘煬,管彥武,石甲強,肖鋒

吉林大學地球探測科學與技術學院，長春130026

0 引言

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

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

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

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

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

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

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

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

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

(2)

(3)

(4)

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

(5)

(6)

(7)

(8)

(9)

(10)

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

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

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

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

(11)

(12)

(13)

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

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

(14)

(15)

(16)

(17)

(18)

(19)

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

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

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

表1 兩種算法全張量磁梯度計算次數

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

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

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

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

模型參數為：磁化方向I=50°，A=30°;磁化強度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;網格間距為100 m×100 m;計算高度z=0 m;等值線單位為nT/m。

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

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

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

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

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

4 結論

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

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