黃宗理, 郭佳

1 中國地質(zhì)科學(xué)院地球深部探測中心, 北京　100035 2 地下信息探測技術(shù)與儀器教育部重點實驗室(中國地質(zhì)大學(xué)，北京)， 北京　100083

?

二維位場在有源空間的解析延拓

黃宗理1, 郭佳2

1 中國地質(zhì)科學(xué)院地球深部探測中心, 北京100035 2 地下信息探測技術(shù)與儀器教育部重點實驗室(中國地質(zhì)大學(xué)，北京)， 北京100083

摘要位場的向下延拓對重、磁等位場資料的處理與解釋具有重要意義，但現(xiàn)有的延拓方法在位場向下接近場源時會產(chǎn)生劇烈震蕩，計算過程很不穩(wěn)定，而且理論上都不能越過場源. 本文在復(fù)數(shù)域里定義了一個具有位場和幾何雙重特性的新函數(shù)PFGR，提出了一套新的延拓方法，在理論上實現(xiàn)了位場在包含場源的下半空間里的解析延拓，經(jīng)過模型試算驗證，表明該方法可以延拓越過場源，準確確定場源奇點，計算過程穩(wěn)定，對測量中的隨機噪聲也具有較好的抗干擾能力.

關(guān)鍵詞位場； 解析延拓； 有源空間

1引言

位場延拓在重力、磁測等位場資料的處理和解釋中具有廣泛用途.因為向下延拓平面更接近于場源，場在場源附近的特性對于反演有重要意義，所以對位場向下延拓方法的研究一直是重磁數(shù)據(jù)處理方法研究中的熱點之一.方法研究沿空間域和頻率域兩個方向進行，但數(shù)學(xué)原理相同，都是由外部狄里希萊問題求出向上延拓的嚴格解，再導(dǎo)出向下延拓的解，并根據(jù)計算方法的不同，發(fā)展了多種多樣向下延拓的具體方法(Huestis et al., 1979, 1998；Fedi et al., 2002；梁錦文，1989；王興濤等，2004；徐世浙，2006；陳生昌和肖鵬飛，2007；劉東甲等，2009；閆輝等，2010).但現(xiàn)有的延拓方法在向下延拓接近場源時，都會產(chǎn)生劇烈震蕩，計算都十分不穩(wěn)定.楊文采(1997)曾從數(shù)學(xué)上證明 “位場向下延拓的最大深度是接近最淺的奇點，再深時解不存在”.因此目前所有的向下延拓方法，理論上都不能越過場源.

位場向下延拓不能越過場源的根本原因，是因為位場在場源奇點上趨于無窮大，位場函數(shù)在有源空間里不是解析函數(shù).針對這一問題，本文在復(fù)數(shù)域里定義了一個新函數(shù)，該函數(shù)在場源奇點上是有限值，并在包含場源的空間里解析，解決了位場在有源空間不解析給向下延拓造成的困擾.提出了一套新的延拓方法，可以越過場源，不會產(chǎn)生震蕩，在理論上實現(xiàn)了位場在有源空間的解析延拓.

2方法原理

二維水平圓柱體的復(fù)磁場強度F的復(fù)數(shù)表達式(吳功建等，1980；張秋生，1985)：

(1)

F的實部和虛部分別是磁場的垂直分量Z和水平分量H.

二維水平圓柱體的復(fù)磁位V的復(fù)數(shù)表達式：

(2)

(1)、(2)式中M為水平圓柱體的復(fù)磁矩；

ζ0=x0+iy0為水平圓柱體柱心復(fù)坐標；

z=x+iy為空間任意點的復(fù)坐標(見圖1).

圖1　復(fù)平面上水平圓柱體磁場的復(fù)坐標Fig.1　Complex coordinates of horizontal cylinder magnetic field on the complex plane

定義一個新函數(shù)R：

(3)

以(1)、(2)式代入R：

(4)

因為復(fù)磁位V和復(fù)場強F在無場源區(qū)是解析函數(shù)，解析函數(shù)的商也是解析函數(shù)(梁昆淼，1965；張秋生，1985)，所以R在無場源區(qū)域是解析的.在場源處z=ζ0，R=z-ζ0=0是有限值.很容易證明R滿足柯西-黎曼方程.

R=z-ζ0=(x-x0)+i(y-y0)=u+iv,

(5)

設(shè)水平圓柱的中心坐標 ζ0=20+20i，已知其在半徑50的圓環(huán)邊界上的R函數(shù)值：

(6)

按照以上推導(dǎo)，根據(jù)柯西公式由邊界上的值計算圓內(nèi)任意z點的R值：

(7)

實際計算得到R函數(shù)值，繪出R的等值線如圖2.

圖2　在半徑50 m的圓內(nèi)延拓出的|R|函數(shù)等值線Fig.2　Continuation contour of the R-Function in a circle with a radius of 50 m

(1) 圖2結(jié)果表示，可以根據(jù)邊界值計算出圓內(nèi)所有點上R的函數(shù)值，即R函數(shù)可以在有源空間內(nèi)解析延拓.R在圓柱體中心的ζ0點有極小值，R等值線為包圍場源中心ζ0的同心圓.隨著遠離場源，R逐漸增大.

(2) 由(4)式知，R等于空間觀測點z與場源奇點ζ0之間的距離，在復(fù)平面上，R的長度可用兩點復(fù)坐標之差來計算.R等值線圖形是描述空間觀測點與場源奇點長度關(guān)系的幾何圖形，R具有幾何屬性.

(3) 因R是由復(fù)磁位V和復(fù)磁場強度F相除計算得出的，故R可以反映V和F的位場屬性，比如R在空間任意點的函數(shù)值，應(yīng)反映所有場源產(chǎn)生的場在該點的疊加效應(yīng).

(4)R的幾何屬性和反映的位場屬性在以下的運算中都要使用，很難說清這個具有雙重屬性函數(shù)的物理意義.借用愛因斯坦在廣義相對論中把引力和空間彎曲結(jié)合在一起的思想，我們把R叫做PFGR函數(shù)(Potential field geometric representation).

3延拓計算方法

實際工作中只能測量得到地面或空中的磁場強度，不能測出地下閉曲線上的磁場值.要使用柯西公式沿閉曲線積分進行解析延拓，需將直角坐標系中的地面直線映射到閉合圓周上才能進行.

3.1進行三次保角映射

第一次映射通過公式(8)將Z平面映射為W1平面：

(8)

第二次映射通過公式(9)將W1平面映射為W2平面：

(9)

第三次映射通過公式(10)將W2平面映射為W3平面：

(10)

三次映射坐標對應(yīng)關(guān)系見圖3.

圖3　三次映射中Z， W1， W2 ，W3平面的對應(yīng)關(guān)系Fig.3　The corresponding relation of Z，W1 ，W2 and W3 in three mappings

三次映射可合并為一個統(tǒng)一的正映射公式：

(11)

至此已把Z平面映射為W平面，Z平面上的有限長線段AC映射為單位圓周，根據(jù)解析函數(shù)論中的邊界對應(yīng)原理，Z平面上的A→B→C→D→A包圍的三角形區(qū)域映射為W1，W2，W3上相應(yīng)點包圍的區(qū)域，在W(W3)平面上映射為單位圓的內(nèi)域.

按照(11)式計算出Z平面AC線段上測點Z(x,0),與W平面上單位圓周上對應(yīng)點eiφ.x與φ的對應(yīng)關(guān)系見圖4.

圖4　映射前后x與φ對應(yīng)關(guān)系Fig.4　The corresponding relation of x and φ before and after mapping

由圖4見映射前直線上均勻分布的測點x，映射為單位圓周上總體均勻分布的φ點，這對使用柯西公式進行數(shù)值積分是有利的.

由正映射公式(11)可推出由W到Z的反映射公式：

(12)

3.2構(gòu)建PFGR函數(shù)

根據(jù)地面上的實測場值構(gòu)建單位圓周上的PFGR函數(shù).已知地面上的V和F，則地面上的R函數(shù)為

(13)

所以

(14)

經(jīng)三次保角映射后，Z平面上的Z(x,0)點映射為W平面單位圓周上的Wz點，Wz=eiφ.Z平面上的ζ0點映射為W平面上的Wζ0點(見圖5)：

(15)

圖5　映射前后坐標點對應(yīng)圖Fig.5　The corresponding relation of coordinate points before and after mapping

則Z平面上的距離R映射為W平面上的距離WR：

(16)

WR是單位圓邊界上的PFGR函數(shù)值.

3.3解析延拓

已知單位圓邊界上的PFGR函數(shù)值，則W平面單位圓內(nèi)任意點m的PFGR函數(shù)值可按照柯西公式計算得出

(17)

以上計算的流程框圖如下(圖6)：

圖6　延拓計算流程Fig.6　Flowchart of analytical continuation

3.4模型驗證

設(shè)定水平圓柱的圓心坐標x0=y0=-100, ζ0=-100-100i,2M=1000cgs,代入(1)、(2)式，其理論復(fù)磁位V和復(fù)場強F為

(18)

(19)

根據(jù)地面上的V和F，采取上述方法向下半空間解析延拓，延拓結(jié)果見圖7.

圖7　水平圓柱體的解析延拓結(jié)果Fig.7　Result of analytical continuation for horizontal cylinder body

由圖7見，延拓計算很穩(wěn)定，延拓出的場源奇點等于設(shè)定的理論值，延拓計算精度很高. PFGR形態(tài)為包含場源圓心的同心圓，與磁位V的形態(tài)類似.

4任意形狀二度體的解析延拓

4.1延拓方法原理

物體的磁化是其內(nèi)部磁偶極子定向排列的結(jié)果，磁性體的外部磁場則是其內(nèi)部磁偶極子場的疊加，二維時則是偶極線磁場的疊加.因為水平圓柱體的外部磁場等效于偶極線的磁場，故任意形狀磁性體的磁場就可用無限多個水平圓柱體磁場的疊加來等效表示.

(20)

式中，Δζi=z-ζi，ζi為場源內(nèi)奇點(水平圓柱體柱心)的坐標.

(21)

考慮到均勻磁化時各個水平圓柱體的磁矩Mi相等：

(22)

(23)

則

·Δζ1Δζ2…Δζn.

(24)

(24)式中R是一系列(z-ζi)的組合，只要任意一個Δζi=z-ζi=0則R=0，即任意形狀的磁場源，當(dāng)觀測點坐標等于場源內(nèi)任意一個奇點的坐標時，PFGR函數(shù)都等于0.而當(dāng)所有的Δζi≠0時，|R|≠0.PFGR等值線表現(xiàn)為在所有的場源奇點上，R取極小值，不在場源奇點上，R值增大，任意形狀磁性體的PFGR等值線性質(zhì)與單個水平圓柱體場源時相同.

分析(24)式可見，式中所有自變量Δζi=z-ζi表示的都是測點至場源的距離，故R和單一圓柱體一樣，表現(xiàn)的也是距離參數(shù)，只是很難求出如x-ζi那樣的解析表達式.我們采用數(shù)值方法來計算.

4.2任意形狀二度體解析延拓的計算方法

在每個測點x下方放一個等效圓柱體，規(guī)定該點的磁場只是由這一個圓柱體產(chǎn)生的, 則該圓柱體的柱心位置就可由(14)式計算:

(25)

N個測點下有N個圓柱體，可以算出N個ζi，求出各x測點對應(yīng)的ζi，就可以按照單一圓柱體源同樣的方法進行解析延拓了.

這些圓柱體在測線下方很近的位置排列成一個密集的鏈條，只要測點足夠密，這個鏈條產(chǎn)生的磁場就與磁偶層等效.垂直磁化磁偶層面上的磁場是垂直方向的(實測磁場可先化極)，故規(guī)定x點的磁場只受其下方一個圓柱體的作用，在一定精度下是合理的.

求出的ζi是否正確需要驗證. 理論上x軸各點的磁場應(yīng)等于所有等效圓柱體場的疊加，設(shè)其磁位為V1，場強為F1，每個圓柱體的磁矩等于常數(shù)K，則:

(26)

(27)

將計算出的等效源磁場V1和F1和實測磁場V和F進行比較，如果二者在允許誤差范圍內(nèi)，則認為求出的ζi是正確的，可以繼續(xù)進行延拓；如果不相等，延拓的誤差可能增大.

4.3驗證實例1——傾斜板狀體的延拓

設(shè)定傾斜板狀體的上頂坐標為(50, -100)，下底坐標為(100, -150)，復(fù)磁矩2M=1000 cgs，先計算出其理論復(fù)磁位和復(fù)場強，再采用上述延拓方法進行延拓，結(jié)果如圖8.

圖8　板狀體的解析延拓結(jié)果Fig.8　Result of analytical continuation for tabular body

在圖8中，圖a的實線表示傾斜板狀體的理論磁場V和F的模，標記×表示由等效源計算出的磁場V1和F1的模，由圖可見二者吻合得很好.說明計算出等效源的ζ(i)是正確的.用這些ζ(i)進行解析延拓，延拓順利地越過場源，PFGR等值線確定的極小值范圍與傾斜板狀體重合，等值線圓滑，沒有出現(xiàn)震蕩，延拓精度較高.

4.4驗證實例2——三角棱柱體的延拓

設(shè)定三角棱柱體頂點坐標分別為(80, -100)、(60, -120)和(100, -120)，復(fù)磁矩2M=1000 cgs，計算出理論磁場V和F，延拓結(jié)果如圖9.

圖9　三角棱柱體的解析延拓結(jié)果Fig.9　Result of analytical continuation for triangular prism body

在圖9中，圖a的實線表示三角棱柱體的理論磁場V和F的模，標記×表示由等效源計算出的磁場V1和F1的模，由圖可見，等效源計算出的磁場與三角棱柱體理論磁場吻合很好，向下延拓計算穩(wěn)定，沒有出現(xiàn)震蕩，PFGR等值線確定的奇點位置與三角棱柱體重合，延拓精度較高.

5延拓影響因素分析

5.1保角映射時如何選擇測線長度的大小

由圖3知，測線長度AC=2a的大小決定著向下延拓三角形的深度和邊長.因為PFGR函數(shù)R=z-ζ，如果a過大，R中的z相應(yīng)也很大.前已證明R只在包含ζ的有限大空間里解析，z過大，ζ提供的有用信息相對減小，延拓誤差會增大；如果a過小，實測剖面長度太短，磁場參數(shù)V和F形態(tài)不完整，不能給出場源全部信息，誤差也會增大.試算結(jié)果表明，以測點距離為單位1，a選擇500較好.

5.2測量誤差的影響

實際工作中磁場參數(shù)V和F都存在測量誤差，我們對一個組合長方體模型的V和F添加10%的隨機噪聲，延拓結(jié)果見圖10.延拓結(jié)果表明，延拓方法的抗干擾能力較好，計算穩(wěn)定，沒有出現(xiàn)震蕩，延拓出的奇點位置也很準確.

圖10　存在10%隨機誤差時組合板狀體磁場的解析延拓Fig.10　Result of analytical continuation with 10% random error

5.3如何得到復(fù)磁位和復(fù)磁場強度

本文推導(dǎo)計算中用的是復(fù)磁位和復(fù)場強，驗證實例中使用的是垂直磁化的復(fù)磁場.實際工作中一般只測量磁場強度ΔT，或重力場強的垂直分量Δg，這些參數(shù)都可以在頻率域中通過快速傅里葉變換(FFT)互相轉(zhuǎn)換，簡單易行，本文不贅述.

6結(jié)論

(1) 針對現(xiàn)有方法向下延拓出現(xiàn)困難的根本原因是位場在有源空間里不解析，本文定義了一個新函數(shù)PFGR，其值等于復(fù)磁位和復(fù)磁場強度的商.理論推導(dǎo)和模型試驗表明，PFGR函數(shù)在包含場源的有限大空間里解析，在場源奇點處為有限值，可以用柯西公式對PFGR函數(shù)進行解析延拓.

(2) 利用PFGR函數(shù)的幾何特性和反映的位場特性，提出了一套新的向下延拓的數(shù)值計算方法，經(jīng)過理論模型試算表明，對于任意形狀的二度體，可以實現(xiàn)在包含場源的下半空間的解析延拓，PFGR等值線圈定了場源奇點.計算過程穩(wěn)定，沒有產(chǎn)生震蕩，計算速度很快.

(3) 模型試算表明，延拓方法具有較高的精度，對測量中出現(xiàn)的隨機誤差也有較好的抗干擾能力.

致謝衷心感謝董樹文研究員、姚長利教授、孟小紅教授和張永謙博士在本文研究工作中所給予的支持.

References

Chen S C, Xiao P F. 2007. Wavenumber domain generalized inverse algorithm for potential field downward continuation.ChineseJ.Geophys. (in Chinese), 50(6):1816-1822.

Fedi M, Florio G. 2002. A Stable downward continuation by using the isvd method.GeophysicalJournalInternational, 151(1):146-156.Huestis S P, Parker R L. 1979. Upward and downward continuation as inverse problems.GeophysicalJournalInternational, 57(1):171-188.Liang J W. 1989. Downward continuation of regularization methods for potential fields.ChineseJ.Geophys. (in Chinese), 32(5):600-608.Liang K M. 1965. Methods of Mathematical Physics (in Chinese). Beijing: Higher Education Press.

Liu D J, Hong T Q, Jia Z H, et al. 2009. Wave number domain iteration method for downward continuation of potential fields and its convergence.ChineseJ.Geophys. (in Chinese), 52:(6):1599-1605,doi:10.3969/j.issn.0001-5733.2009.06.022.

Huestis S P. 1998. The continuation inverse problem revisited.GeophysicalJournalInternational, 133(3):705-712.

Wang X T. 2004. A comparison of different downward continuation methods for airborne gravity data.ChineseJ.Geophys. (in Chinese), 47(6):1017-1022.Wu G J. 1980. Applied Geophysics—Magnetic Course (in Chinese). Beijing: Geological Publishing House.

Xu S Z. 2006. The integral-iteration method for continuation of potential fields.ChineseJ.Geophys. (in Chinese), 49(4):1176-1182.

Xue Q F. 1978. Field Theories (in Chinese). Beijing: Geological Publishing House.

Yan H, Xiao C H, Zhang Z Y, et al. 2010. Iterative method for continuation of three-component magnetic field.Chinese

Journal of Computational Physics. (inChinese), 27(05):705-710.YangWC. 1997.Inversetheoriesandmethodsofgeophysicalproblems(inChinese).Beijing:GeologicalPublishingHouse.

ZhangQS. 1985.FieldTheory(inChinese).Beijing:GeologicalPublishingHouse.

附中文參考文獻

陳生昌, 肖鵬飛. 2007. 位場向下延拓的波數(shù)域廣義逆算法. 地球物理學(xué)報, 50(06):1816-1822.

梁錦文. 1989. 位場向下延拓的正則化方法. 地球物理學(xué)報, 32(5):600-608.

梁昆淼．1965. 數(shù)學(xué)物理方法．北京：高等教育出版社.

劉東甲，洪天求，賈志海等. 2009. 位場向下延拓的波數(shù)域迭代法及其收斂性. 地球物理學(xué)報, 52:(6):1599-1605,doi:10.3969/j.issn.0001-5733.2009.06.022.

王興濤. 2004. 航空重力測量數(shù)據(jù)向下延拓方法比較. 地球物理學(xué)報, 47(6):1017-1022.

吳功建．1980. 應(yīng)用地球物理學(xué)-磁法教程．北京：地質(zhì)出版社.

徐世浙．2006. 位場延拓的積分-迭代法．地球物理學(xué)報, 49(4):1176-1182.

薛琴訪．1978. 場論．北京：地質(zhì)出版社.

閆輝,肖昌漢，張朝陽等. 2010. 三分量磁場延拓的遞推算法. 計算物理, 27(05):705-710.

楊文采．1997. 地球物理反演的理論與方法．北京：地質(zhì)出版社.

張秋生．1985. 場論．北京：地質(zhì)出版社.

(本文編輯汪海英)

基金項目國家深部探測技術(shù)與實驗研究專項(SinoProbe-08)資助.

作者簡介黃宗理，男，1947年生，高級工程師，現(xiàn)從事地球物理位場方法及應(yīng)用領(lǐng)域的研究工作. E-mail：Zlhuang5653@hotmail.com

doi:10.6038/cjg20160227 中圖分類號P631

收稿日期2014-09-16，2015-12-22收修定稿

Analytical continuation of the two-dimensional potential field in an active space

HUANG Zong-Li1, GUO Jia2

1SinoProbeCenter，ChineseAcademyofGeologicalSciences,Beijing100035,China2KeyLaboratoryofGeo-detection(ChinaUniversityofGeosciences,Beijing),MinistryofEducation,Beijing100083,China

AbstractDownward continuation of the potential field has a great significance in practices. However, the existing downward continuation methods cannot continue a potential field over the source and their numerical calculation is unstable. We define a new function PFGR which has both the potential field and geometry properties. We also present a novel continuation method which can analytically continue the potential field into the lower half-space including the source. This method can steadily continue the potential field over the source. And the test on synthetic data confirmed the reliability of this new method.

KeywordsPotential field； Analytical continuation； Active space

黃宗理, 郭佳. 2016. 二維位場在有源空間的解析延拓.地球物理學(xué)報，59(2):704-710,doi:10.6038/cjg20160227.

Huang Z L, Guo J. 2016. Analytical continuation of the two-dimensional potential field in an active space.ChineseJ.Geophys. (in Chinese),59(2):704-710,doi:10.6038/cjg20160227.