邰振華， 柴 琳， 黃德智， 商宇航

(黑龍江科技大學(xué) 礦業(yè)工程學(xué)院， 哈爾濱 150022)

0 引 言

在基礎(chǔ)地質(zhì)研究、礦產(chǎn)勘查、工程與環(huán)境調(diào)查等領(lǐng)域，長(zhǎng)方體是模擬磁性場(chǎng)源的重要模型，常用于基礎(chǔ)轉(zhuǎn)換、邊界檢測(cè)等數(shù)據(jù)處理方法的模型實(shí)驗(yàn)，是三維反演的核心構(gòu)件[1]。用于模型實(shí)驗(yàn)時(shí)，長(zhǎng)方體的正演值作為待處理數(shù)據(jù)，若存在解析奇點(diǎn)，可能導(dǎo)致實(shí)驗(yàn)效果不佳，如邊界檢測(cè)出現(xiàn)間斷、虛假特征等。三維形態(tài)、物性反演需要將地下空間剖分成長(zhǎng)方體單元，以其正演值與觀測(cè)值的擬合度調(diào)控迭代進(jìn)程[2]。然而，長(zhǎng)方體正演結(jié)果的解析奇點(diǎn)恰好位于角點(diǎn)，導(dǎo)致迭代運(yùn)算的收斂性差甚至無(wú)法收斂。

在常規(guī)磁異常的正演方面，郭志宏等[3]最早指出，長(zhǎng)方體的常規(guī)理論表達(dá)式在上半無(wú)場(chǎng)源空間存在無(wú)法計(jì)算的解析奇點(diǎn)，推導(dǎo)了新的磁異常與梯度正演公式，降低了奇點(diǎn)影響，但增加了計(jì)算量。駱遙等[4]利用歐拉方程重新推導(dǎo)了長(zhǎng)方體磁異常與梯度的無(wú)解析奇點(diǎn)理論表達(dá)式，簡(jiǎn)化了公式構(gòu)造，優(yōu)化了計(jì)算量。Luo等[5]推導(dǎo)了均勻磁化多面體的磁異常與梯度正演公式，可用于傾斜長(zhǎng)方體的正演，但必須預(yù)先給定法線方向、用于坐標(biāo)變換的9個(gè)旋轉(zhuǎn)角度，增加了設(shè)計(jì)階段的工作量。Kuang等[6]推導(dǎo)了適用于起伏地形的長(zhǎng)方體磁異常無(wú)解析奇點(diǎn)正演公式。

相比于常規(guī)磁異常，磁張量的分辨率更高，能提供更豐富、精細(xì)的地質(zhì)信息[7]。研發(fā)基于磁張量的數(shù)據(jù)處理新技術(shù)，突破地質(zhì)體位置、形狀、物性等反演的精度瓶頸，已成為磁法勘探的研究熱點(diǎn)[8]。與常規(guī)磁異常類似，磁張量的數(shù)據(jù)處理方法測(cè)試、三維反演依賴長(zhǎng)方體模型。因此，研究長(zhǎng)方體磁張量無(wú)解析奇點(diǎn)正演公式，能夠推動(dòng)磁張量三維反演技術(shù)的發(fā)展。

在前人的研究基礎(chǔ)上，筆者推導(dǎo)了直立長(zhǎng)方體的磁張量無(wú)解析奇點(diǎn)正演公式，給出了簡(jiǎn)潔的坐標(biāo)變換方法，建立了傾斜長(zhǎng)方體的磁張量無(wú)解析奇點(diǎn)正演方法。利用模型實(shí)驗(yàn)，驗(yàn)證了正演結(jié)果的準(zhǔn)確性，分析了地磁傾角、偏角對(duì)磁張量的影響。

1 有解析奇點(diǎn)的磁張量正演公式

在常規(guī)正演理論中，長(zhǎng)方體磁場(chǎng)的x、y、z方向分量依次為Hax、Hay、Za?；诘卮艌?chǎng)基本理論，分別對(duì)Hax、Hay、Za求取x、y、z方向的一階偏導(dǎo)數(shù)，可以得到長(zhǎng)方體的磁張量正演公式，即

(1)

(2)

(3)

(4)

(5)

(6)

ξi=(-1)ia-(x-x0) ，

ηj=(-1)jb-(y-y0)，

ζk=(-1)kc-(z-z0)，

式中：V——磁力位；

μ0——真空磁導(dǎo)率；

J——磁化強(qiáng)度；

ξi、ηj、ζk、Ri,j,k——距離關(guān)系。

L、M、N——總磁化強(qiáng)度的方向余弦；

(x0，y0，z0)——長(zhǎng)方體中心點(diǎn)坐標(biāo)；

(x，y，z)——計(jì)算點(diǎn)坐標(biāo)；

a、b、c——長(zhǎng)方體各邊長(zhǎng)的一半。

當(dāng)計(jì)算點(diǎn)與長(zhǎng)方體角點(diǎn)的垂直投影點(diǎn)重合時(shí)，ξ=0、η=0，式(1)、(2)、(4)存在分母為零的部分，計(jì)算結(jié)果趨于無(wú)窮大，形成了解析奇點(diǎn)。

2 無(wú)解析奇點(diǎn)的磁張量正演算法

2.1 直立長(zhǎng)方體的無(wú)解析奇點(diǎn)正演公式

駱遙等[4]重新推導(dǎo)了直立長(zhǎng)方體的Hax、Hay、Za正演公式。在此基礎(chǔ)上，依據(jù)地磁場(chǎng)基本理論，分別推導(dǎo)Hax在x、y、z方向的一階偏導(dǎo)數(shù)，可以導(dǎo)出直立長(zhǎng)方體磁張量Vxx、Vxy、Vxz的無(wú)解析奇點(diǎn)正演公式，結(jié)合磁張量的等價(jià)互換性質(zhì)，得出

(7)

(8)

(9)

同理，分別推導(dǎo)Hay在y和z方向一階偏導(dǎo)數(shù)，Za在z方向一階偏導(dǎo)數(shù)，得出

(10)

(11)

(12)

分析式(7)～(12)可知，不存在分母為零的情況，整個(gè)計(jì)算區(qū)域不存在解析奇點(diǎn)。同時(shí)，新推導(dǎo)的磁張量正演公式繼承了駱遙等研究成果的簡(jiǎn)潔屬性，計(jì)算量較小。

2.2 傾斜長(zhǎng)方體的無(wú)解析奇點(diǎn)正演算法

對(duì)于走向與x、y方向不平行或者傾角不為90°、0°的傾斜長(zhǎng)方體，磁張量正演可以借助坐標(biāo)變換的方法加以實(shí)現(xiàn)。若長(zhǎng)方體走向與x方向的夾角為α，則可以將原坐標(biāo)系xyz圍繞z軸逆時(shí)針旋轉(zhuǎn)α角，形成新的坐標(biāo)系xr1yr1zr1，如圖1所示。

圖1 圍繞z軸旋轉(zhuǎn)示意Fig. 1 Schematic of rotation around z-axis

xyz坐標(biāo)與xr1yr1zr1坐標(biāo)的變換關(guān)系為

(13)

將式(13)轉(zhuǎn)化為矩陣關(guān)系式

(14)

同理，若長(zhǎng)方體傾角不為90°、0 °，可以繼續(xù)圍繞x軸逆時(shí)針旋轉(zhuǎn)β或圍繞y軸逆時(shí)針旋轉(zhuǎn)γ，建立新的坐標(biāo)系?；诰仃嚨某朔ㄟ\(yùn)算，最終的矩陣關(guān)系式為

(15)

(16)

在坐標(biāo)系xryrzr下，傾斜長(zhǎng)方體轉(zhuǎn)化為直立長(zhǎng)方體，ξi=(-1)ia-xr；ηj=(-1)jb-yr；ζk=(-1)kc-zr，可以利用式(7)～(12)計(jì)算磁張量正演值。而后，將坐標(biāo)系xryrzr下的磁張量正演值轉(zhuǎn)化回坐標(biāo)系xyz，轉(zhuǎn)化關(guān)系式為

(17)

式中：Vxx～Vzz——傾斜長(zhǎng)方體在坐標(biāo)系xyz下的磁張量正演結(jié)果；

相比于Luo等[5]的坐標(biāo)轉(zhuǎn)換方法，文中方法只需給定與長(zhǎng)方體產(chǎn)狀直接相關(guān)的3個(gè)旋轉(zhuǎn)角度，可有效降低模型設(shè)計(jì)階段的工作量。

3 模型實(shí)驗(yàn)

3.1 直立長(zhǎng)方體模型

為了檢驗(yàn)直立長(zhǎng)方體磁張量無(wú)解析奇點(diǎn)正演公式的正確性，設(shè)計(jì)如表1所示的單一直立長(zhǎng)方體模型，其中，θ為地磁傾角，θD為地磁偏角，J為磁化強(qiáng)度。令z軸向下，觀測(cè)面是z=0的平面，利用式(7)～(12)計(jì)算各方向的磁張量無(wú)解析奇點(diǎn)正演結(jié)果(圖2)。圖2中，磁張量分布特征符合長(zhǎng)方體的異常變化情況，非對(duì)稱分布(圖2f的低值圈閉)主要源于地磁傾角與磁偏角的影響。利用式(1)、(2)、(4)計(jì)算出有解析奇點(diǎn)的Vxx、Vxy、Vyy(圖3)。對(duì)比圖2a～c與圖3a～c，除4個(gè)解析奇點(diǎn)外，式(1)、(2)、(4)與式(7)、(8)、(10)的計(jì)算結(jié)果完全一致。此外，以表1模型產(chǎn)生的磁異常為基礎(chǔ)，利用磁張量轉(zhuǎn)換方法計(jì)算磁張量[9]，與無(wú)奇點(diǎn)正演結(jié)果的差異見(jiàn)表2。考慮到磁張量轉(zhuǎn)換過(guò)程存在頻譜畸變與轉(zhuǎn)換誤差[10]，張量轉(zhuǎn)換與文中無(wú)奇點(diǎn)正演公式的計(jì)算結(jié)果基本相符。上述實(shí)驗(yàn)結(jié)果說(shuō)明，文中推導(dǎo)的直立長(zhǎng)方體磁張量無(wú)解析奇點(diǎn)正演公式是正確的。表2中均方差、均方差占比計(jì)算公式為

圖2 直立長(zhǎng)方體磁張量無(wú)解析奇點(diǎn)的正演結(jié)果Fig. 2 Forward results of upright cuboid magnetic tensor without analytical singularities

圖3 直立長(zhǎng)方體磁張量有解析奇點(diǎn)的正演結(jié)果Fig. 3 Forward results of upright cuboid magnetic tensors with analytical singularities

表2 斜磁化時(shí)無(wú)奇點(diǎn)正演結(jié)果與張量轉(zhuǎn)換結(jié)果的誤差Table 2 Errors between singularity-free forward results and tensor conversion results in inclined magnetization

(18)

(19)

式中：RMSE——均方差；

PRMSE——均方差占比；

FM——磁張量正演結(jié)果；

TC——磁張量轉(zhuǎn)換結(jié)果；

|FM|max——FM絕對(duì)值的最大值；

K、J——x、y方向計(jì)算點(diǎn)總數(shù)。

3.2 地磁傾角與磁偏角的影響

依據(jù)表1，將地磁傾角改為90°、磁偏角改為0°，利用式(7)～(12)計(jì)算垂直磁化長(zhǎng)方體的磁張量，與斜磁化正演結(jié)果(圖2)的差值見(jiàn)圖4。由圖4可知，對(duì)于同一地質(zhì)體，垂直磁化與斜磁化的磁張量差異較大，主要表現(xiàn)在張量的特征點(diǎn)(線)位置，如Vxx中的長(zhǎng)方體邊界位置、Vxy中的長(zhǎng)方體角點(diǎn)位置。地磁傾角與磁偏角能夠引起張量的幅值弱化、曲線形態(tài)變化，故化磁極同樣是磁張量數(shù)據(jù)處理的必備環(huán)節(jié)。

圖4 斜磁化與垂直磁化的直立長(zhǎng)方體磁張量差值Fig. 4 Differences between clino-magnetized and perpendicular-magnetized upright cuboid magnetic tensors

3.3 傾斜長(zhǎng)方體模型

3.3.1 垂直磁化傾斜長(zhǎng)方體模型

為了檢驗(yàn)文中推導(dǎo)的傾斜長(zhǎng)方體磁張量無(wú)解析奇點(diǎn)正演算法的正確性，設(shè)計(jì)如表3的單一傾斜長(zhǎng)方體模型，θd為長(zhǎng)方體傾角。令z軸向下，觀測(cè)面是z=0的平面，分別利用張量轉(zhuǎn)換方法與文中方法計(jì)算各方向的磁張量。分析圖5可知，垂直磁化下，文中方法與張量轉(zhuǎn)換方法計(jì)算的傾斜長(zhǎng)方體磁張量差值很小，主要集中在磁張量變化較大的位置，與表4的誤差統(tǒng)計(jì)相互印證?？紤]到張量轉(zhuǎn)換自身的誤差，推導(dǎo)的傾斜長(zhǎng)方體磁張量正演公式是正確的。

圖5 張量轉(zhuǎn)換與文中方法計(jì)算的傾斜長(zhǎng)方體磁張量差值Fig. 5 Differences of inclined cuboid magnetic tensors calculated by tensor transformation method and this method

表3 傾斜長(zhǎng)方體模型參數(shù)Table 3 Parameters of inclined cuboid model

表4 垂直磁化時(shí)文中方法與張量轉(zhuǎn)換計(jì)算的傾斜長(zhǎng)方體磁張量差異Table 4 Differences of perpendicular-magnetized inclined cuboid magnetic tensors calculated by this method and tensor transformation method

3.3.2 斜磁化傾斜長(zhǎng)方體模型

在不改變傾斜長(zhǎng)方體參數(shù)(表3)的情況下，地磁傾角改為75°、磁偏角改為15°(斜磁化條件)，利用文中方法計(jì)算的磁張量如圖6所示。與表3對(duì)比，磁張量正演結(jié)果的特征點(diǎn)(線)與長(zhǎng)方體的角點(diǎn)(邊界)對(duì)位，在張量Vxz、Vyz、Vzz中表現(xiàn)最明顯；受地磁傾角與磁偏角影響，磁張量的極值圈閉呈現(xiàn)非對(duì)稱的分布特征，軸線方向與地磁偏角對(duì)應(yīng)。如表5所示，與張量轉(zhuǎn)換方法的計(jì)算結(jié)果對(duì)比，PRMSE介于1.490%～2.990%。上述實(shí)驗(yàn)結(jié)果表明，文中方法計(jì)算的斜磁化傾斜長(zhǎng)方體磁張量在等值線形態(tài)、幅值兩方面均正確。表2與表5的對(duì)比結(jié)果說(shuō)明，張量轉(zhuǎn)換方法的計(jì)算精度同樣受地磁傾角與磁偏角的影響，化磁極與張量轉(zhuǎn)換的協(xié)同運(yùn)算才能得到較好的轉(zhuǎn)換結(jié)果。

圖6 文中方法計(jì)算的斜磁化傾斜長(zhǎng)方體磁張量Fig. 6 Magnetic tensors of clino-magnetized inclined cuboid calculated by this method

表5 斜磁化時(shí)文中方法與張量轉(zhuǎn)換計(jì)算的傾斜長(zhǎng)方體磁張量差異Table 5 Differences of clino-magnetized inclined cuboid magnetic tensors calculated by this method and tensor transformation method

4 結(jié)束語(yǔ)

在均勻磁化長(zhǎng)方體磁場(chǎng)理論表達(dá)式的基礎(chǔ)上，推導(dǎo)了直立長(zhǎng)方體的磁張量有解析奇點(diǎn)與無(wú)解析奇點(diǎn)正演公式，建立了傾斜長(zhǎng)方體的磁張量無(wú)解析奇點(diǎn)正演算法。實(shí)驗(yàn)結(jié)果表明，文中推導(dǎo)的磁張量無(wú)解析奇點(diǎn)正演算法是正確的；地磁傾角與磁偏角能夠影響磁張量的幅值與分布，是降低張量轉(zhuǎn)換精度的不利因素，故化磁極是磁張量數(shù)據(jù)精細(xì)化處理的必要環(huán)節(jié)。