周智文，何水原，孟小紅，王俊，方圓，鄭師警，廖開訓

1 中國地質(zhì)大學(北京)地球物理與信息技術學院，北京 100083 2 中國地質(zhì)調(diào)查局 廣州海洋地質(zhì)調(diào)查局，廣州 510760

0 引言

位場向下延拓技術能夠增強深部場源高頻信息，是位場數(shù)據(jù)處理解釋的關鍵手段.但由于邊界效應、延拓算子不收斂、易受噪聲影響等問題，成為眾多學者研究的焦點.Bateman(1946)將Taylor級數(shù)法的思想運用在位場延拓方面，但沒有給出明確的延拓公式,徐世浙等(2006)利用穩(wěn)定的向上延拓替代向下延拓，將延拓過程用迭代修正的方式進行.隨后劉東甲(2009)在波數(shù)域中使用此方法，給出了波數(shù)域中向下延拓的公式.基于此思想，王彥國等(2011)將向下延拓因子進行Taylor級數(shù)展開，推導出了位場向下延拓的Taylor級數(shù)迭代法迭代通式.Wang等(2018)將下界面實測數(shù)據(jù)與泰勒級數(shù)迭代法相結合，推導出約束向下延拓方法，但該方法受噪聲影響較大.馬國慶等(2016)提出向下延拓的水平導數(shù)迭代法，雖然有效的增強了計算結果的穩(wěn)定性，但容易在二度體周圍出現(xiàn)高頻振蕩現(xiàn)象.張沖等(2017)基于求解常微分方程的多步思想進行向下延拓，引入向上延拓至不同高度的垂向一階導數(shù)，提出空間域重力場向下延拓的3階Adams-Bashforth公式法，但結果有較大的幅值衰減現(xiàn)象.Dampney(1969)通過最小二乘法反演等效源，最早使用等效源進行重力場的空間延拓.

等效源法與常規(guī)改造延拓因子的延拓方法不同，而是模擬源進行正演過程.常用于位場數(shù)據(jù)多種處理中，比如：Needham(1970)利用球坐標系構建等效點層.Bhattacharyya 和Chan(1977)利用磁偶層等效源進行磁異常轉(zhuǎn)換.Cordell和Grauch(1982)結合等效源與傅里葉變換將重磁實測數(shù)據(jù)插值網(wǎng)格化.Xia等(1993)基于實現(xiàn)觀測數(shù)據(jù)最大平滑為準則，提出了根據(jù)點距數(shù)據(jù)光滑度的單層最優(yōu)等效層的概念.Barnes和Lumley(2011)構建等效源對重力梯度數(shù)據(jù)進行分量轉(zhuǎn)換.陳濤和張貴賓(2019)利用等效源技術實現(xiàn)重力異常構建重力梯度.王澤慶等(2021)提出了一種塊體與質(zhì)點的兩層等效源設置方案.由于其具有良好的抗噪性與較高的延拓精度，同樣也可以運用在位場曲面延拓處理方面，徐世浙等(2004)聯(lián)合邊界單元和單層等效源進行重構且實現(xiàn)了航測數(shù)據(jù)向下延拓.黃翼堅等(2009)采用磁偶極子等效源實現(xiàn)曲面隨機點磁異常位場數(shù)據(jù)延拓等處理.閆輝等(2010)結合等效源方法提出三分量磁場延拓的遞推算法.李曉杰和王真理(2018)研究了正則化等效源的波數(shù)響應，提出迭代補償算法.孟波和賈真人(2020)提出了等效場源設置的優(yōu)化策略，提高了延拓的穩(wěn)定性.但當延拓深度過大時，導致通過改造延拓算子放大高頻信號的延拓方式受噪聲影響嚴重，等效源延拓方法抗噪性雖好但深部場源的信息經(jīng)過衰減后在觀測界面難以探測，從而在擬合構建等效源時缺少深部高頻信息，這也會使得幅值衰減導致結果不準確.

本文針對上述問題，研究了在求解等效源參數(shù)過程中利用下界面實測數(shù)據(jù)并結合Tikhonov正則化的約束等效源向下延拓技術，使得衰減信號得以恢復的同時對噪聲不敏感.通過理論推導和模型試驗與目前常用穩(wěn)定向下延拓方法——波數(shù)域迭代法以及未加約束的等效源法進行比較，證明該方法的可行性.并且依托“深水油氣近海底重磁高精度探測關鍵技術”項目最新研發(fā)的一套具備水下自動定深及姿態(tài)控制功能的深水重磁勘探拖曳系統(tǒng)所測水下數(shù)據(jù)作為約束測線，將該方法運用到實際數(shù)據(jù)處理中，得到的延拓結果與水下已知數(shù)據(jù)擬合效果相比對照方法最佳，驗證了該方法的實用性.

1 等效源約束實現(xiàn)向下延拓

1.1 等效源構建

為了建立實測位場信息與等效源的關系，可以將地下空間剖分為一系列立方體，將問題轉(zhuǎn)變?yōu)榫€性方程的求解，其矩陣形式為：

dobs=km，

(1)

其中：dobs∈RM是實際觀測位場數(shù)據(jù),m∈RM為代求模型參數(shù)，模型正演算子k∈RM×N.假設觀測數(shù)據(jù)點與地下等效源剖分塊體個數(shù)分別為M、N,將等效源固定在深度為h處，則方程(1)可以表示為:

(2)

剖分塊體正演算子為:

Ki,j=Gσ|||[(xmj-xdi)ln(ymj-ydi+ρ)+(ymj+ydi)ln(xmj-xdi+ρ)

(3)

式中，G是萬有引力常量，σ為塊體密度，此處設為1 g·cm-3.xdj、ydj、zdj代表第i個觀測點的坐標.積分限(cximax,cximax)、(cyimax,cyimax)、(czimax,czimax)分別對應xmi、ymi、zmi，代表第j個剖分矩形塊體最小與最大角點坐標，ρ=[(xmj-xdi)2+(ymj-ydi)2+(zmj-zdi)2]1/2，模型目標方程為:

(4)

其中ζ(r)為深度加權算子,m(r)為代求的模型參數(shù)，m0為參考模型，αs、αx、αy以及αz為各方程權重系數(shù)，通常將大小設定為αs?1、αz=αx=αy=1.將模型目標方程變換為矩陣形式為：

(5)

同樣，可以給出模型正演數(shù)據(jù)與真實觀測數(shù)據(jù)誤差方程的矩陣形式為：

φd=(km-dobs)T(km-dobs).

(6)

1.2 等效源約束原理

當延拓距離過大，深部的場源信息在上界面嚴重衰減導致難以觀測到時，可以利用界面下空間的實測值作為約束點(圖1).但下界面測點不一定匹配測網(wǎng)，所以首先需要根據(jù)坐標插值到與上界面觀測網(wǎng)格相同的網(wǎng)格點位置上，假設某約束點坐標為(xck,yck,zck)觀測值為dck.

圖1 等效源向下延拓示意圖Fig.1 Sketch map of equivalent source downward continuation method

圖2 約束點插值示意圖Fig.2 Sketch map of constrained point interpolation

(7)

在擬合觀測面位場數(shù)據(jù)時，將插值完成后的約束點也代入方程組中進行等效源求解，相當于同時擬合兩個界面的數(shù)據(jù)來構建等效源，方法如公式(8)所示：

(8)

相應的，加入約束點后的正演算子(3)變?yōu)楣?9)：

Ki,j=Gσ|||[(xmj-xck)ln(ymj-yck+ρ)+(ymj+ydk)ln(xmj-xck+ρ)

(9)

其中xck、yck、zck代表第k個約束點的坐標.積分限(cximax,cximax)、(cyimax,cyimax)、(czimax,czimax)分別對應xmi、ymi、zmi，代表第i個剖分矩形塊體最小與最大角點坐標，ρ=[(xmj-xck)2+(ymj-yck)2+(zmj-zck)2]1/2.

求解公式(8)是典型的不適定問題，在此引入Tikhonov正則化的方法，其原理是在目標函數(shù)中加入正則化項λ，使得方程結果具有唯一性，則目標函數(shù)為：

(10)

由于公式(5)可知，Wm代表模型參數(shù)各個方向的梯度矩陣之和，則目標函數(shù)(10)為Tikhonov一階泛函，便可得到矩陣運算解表達式為：

(11)

求取合適的λ對最終結果尤為重要，本文采用廣義交叉驗證的方法(GCV)進行計算，公式為：

(12)

為了求取分母中矩陣的跡，采用隨機估計的方法,其原理為利用公式(13)，對矩陣A求跡：

trace(A)=uTAu，

(13)

其中u為-1與1各50%機率的隨機向量，所以可以得到公式(12)的近似表達式為：

(14)

給出λ的取值范圍(本文中采用ln5～ln25)代入公式(14)求得最小值所對應的λ值，作為后續(xù)計算的正則化參數(shù).于是可以得到第n次迭代結果為：

(15)

利用共軛梯度法(CG)迭代求解公式(15)，當Δm滿足預設條件或者達到迭代上限時迭代結束，mn為最終所求模型參數(shù).

2 模型試驗

利用模型試驗驗證等效源約束延拓的可行性，并且將其結果與波數(shù)域迭代法以及普通等效源向下延拓結果進行對比.模型參數(shù)如表1所示.

表1 模型物性參數(shù)表Table 1 Physical properties of synthetic model

模型由兩個空間大小相同，密度和埋深各異的長方體組成，測線距、測點距均為1 m，共計100×100個測點，模型空間分布如圖3所示.

圖3 模型三維示意圖Fig.3 The 3D vision of synthetic model

在z=30 m處觀測面理論重力異常值(圖4a)中加入標準差為5 μGal的隨機噪聲(圖4b),之后將其作為初始面向下延拓30 m，其結果理論值應該如圖4c所示.

通過分析徑向?qū)?shù)功率譜(圖5)來估計等效源埋深，本文模型實驗截取三段徑向?qū)?shù)功率進行分段擬合，紅線與綠線代表的擬合段分別代表對應深部與淺部場源埋深上界面，黑色擬合段為噪聲.之后利用公式(16)計算對應頻段深度(郭良輝，2012)：

(16)

求得深度分別為38.8867 m、4.3377 m和0.4585 m，本文選取38 m設置為等效源深度，將此深度的地下空間剖分為200×200個規(guī)則立方體，每個立方體大小為1 m×1 m×5 m.利用此等效源去擬合z=30 m處加噪理論重力異常(圖4b)，并向下延拓30倍點距(30 m)，結果如圖6b所示.同時，選取圖4c中的測線數(shù)據(jù)作為約束進行相同距離的延拓，結果如圖6c所示.此外，為了驗證該方法的有效性，將原始數(shù)據(jù)用目前穩(wěn)定的波數(shù)域迭代法下延相同深度進行比較，比較其結果(圖6a).

將三種方法的延拓結果與z=0 m觀測理論值做差，得到其偏差值(圖6d、e、f)，可以看出當場源埋深過大時，高頻的位場信息由于衰減后難以在觀測面上體現(xiàn)，再加上噪聲的影響，所以單純利用向下延拓技術得到的結果會與理論值有較大的偏差.

圖5 重力異常徑向?qū)?shù)功率譜及其分段擬合藍色：重力異常功率譜；紅色：頻段擬合1；綠色：頻段擬合2；黑色：頻段擬合3.Fig.5 The radially averaged logarithm power spectrum of the gravity anomaly and its fitting by using piece-wise linearizationThe blue line shows the gravity anomaly power spectrum. The red line shows the fitted straight line in the first radial frequency section, green line for the second section and black line for the third section.

波數(shù)域迭代法在大點距向下延拓時會受到噪聲的干擾，導致結果較差.等效源法延拓可以得到與理論值相似的異常幅值，但異常形態(tài)差別較大.約束等效源法延拓所得結果異常形態(tài)與異常幅值更加接近真實值，并且具有很好的抗噪性.從偏差圖6e、f能更加直觀的得到，等效源法延拓偏差達到了-40～30 μGal,而約束等效源延拓法將圖6e中大部分偏差消除，并且將偏差幅值縮小到了-5～15 μGal.

提取三種方法延拓結果中如(圖4c)所示的中心剖面數(shù)據(jù)進行對比，可以清晰的看出由于噪聲的影響波數(shù)域迭代法結果(黑線)與真實值有較大偏差，約束等效源法(紅線)能更好的擬合真實值，并且噪聲的影響不大(圖7).

3 實測數(shù)據(jù)向下延拓

利用廣州海洋地質(zhì)調(diào)查局所屬“海洋六號”多功能地質(zhì)地球物理調(diào)查船為作業(yè)母船，搭載其最新研發(fā)的一套具備水下自動定深及姿態(tài)控制功能的深水重磁勘探拖曳系統(tǒng)，采用此系統(tǒng)測量了如圖8范圍中的水面與水下300 m處相互對應的四條測線的重力異常.

采用四段直線擬合實測數(shù)據(jù)徑向平均對數(shù)功率譜(圖9)，利用分段斜率估計各段深度由深到淺分別為32822 m、494.739 m、106.650 m和26.5627 m.

本文采用中間兩段的平均深度550 m設置為等效層深度，將該深度的空間剖分成由400×400個1 m×0.5 m×5 m矩形構成的等效層，之后利用該等效層正演結果擬合水面實測重力異常來求取等效源參數(shù)，最后進行300 m向下延拓與水下實測數(shù)據(jù)進行對比.

分別利用波數(shù)域迭代法、等效源法以及約束等效源法向下延拓300 m與水下實測重力異常進行對比(圖10)，約束測線為水下測線抽稀之后的數(shù)據(jù).由圖10b可以看出，迭代法在大距離延拓時噪聲也不可避免的被放大從而影響其精確性，在異常形態(tài)上與實測數(shù)據(jù)有較大差別.等效源法延拓結果(圖10c)呈現(xiàn)出與實測數(shù)據(jù)相似的異常形態(tài)但其幅值相比實測值偏小，異常形態(tài)更多體現(xiàn)出趨勢性，相比實測值缺少一些細節(jié)信息.有約束等效源法延拓結果(圖10d)獲得了與實測值相似的異常形態(tài)與幅值，噪聲影響很小，相比未約束時的結果體現(xiàn)出更多的細節(jié)信息.

圖6 不同方法延拓結果及其偏差圖(a) 波數(shù)域迭代法延拓結果； (b) 等效源法延拓結果； (c) 約束等效源延拓結果； (d) 波數(shù)域迭代法延拓偏差； (e) 等效源法延拓偏差； (f) 約束等效源法延拓偏差.Fig.6 Continuation results of different methods and their deviation diagramsThe results of wave number domain iteration method (a) and the difference between it and the theoretical value (d). The results of equivalent source downward continuation method (b) and the difference between it and the theoretical value (e). The results of equivalent source constrained downward continuation method (c) and the difference between it and the theoretical value (f).

圖7 延拓結果剖面對比圖藍線：下延理論異常；黑線：波數(shù)域迭代法延拓結果；綠線：等效源法延拓結果；紅線：約束等效源法延拓結果.Fig.7 Comparison of different downward continuation methods result profileThe blue line shows the theoretical value of downward continuation. The black line is the result of wave number domain iteration method. The green line represents the result of equivalent source method and the red line is for the result of equivalent source constrain method.

在三種方法的延拓結果中找到與水下四條測線位置相對應的重力異常值并與實測值進行比較，如圖11所示，能更加直觀的看出，波數(shù)域迭代法(黑線)延拓結果異常趨勢與實測值相似但均小于實測值，此外，實測值幅值起伏均在5 mGal之內(nèi)，由于受噪聲影響嚴重延拓結果呈現(xiàn)出高低值震蕩的現(xiàn)象，幅值變化達到20 mGal.而未約束的等效源法(綠線)延拓結果受噪聲影響小，幅值變化范圍以及異常形態(tài)與趨勢都與實測值相似，不過在幅值上均小于真實值.而約束等效源法(紅線)延拓值相比上兩種方法更加接近真實值，幅值變化范圍以及異常形態(tài)與趨勢都能更好的擬合實測值，加入了水下測線作為約束而構建的等效源能在延拓結果中體現(xiàn)出更多衰減掉的高頻信息，從圖中可以看到多處實測值(藍線)有所體現(xiàn)但等效源法延拓結果(綠線)中未出現(xiàn)的信息，在加入約束方法后(紅線)得以體現(xiàn).

圖8 實測重力異常圖(a) 水面實測重力異常； (b) 水下300 m實測重力異常.Fig.8 Measured gravity anomaly map(a) The measured gravity anomaly at sea surface； (b) The measured gravity anomaly at depth of 300 m under the sea surface.

圖9 重力異常徑向?qū)?shù)功率譜及其分段擬合藍色：重力異常功率譜；黑色：頻段擬合1；綠色：頻段擬合2；紅色：頻段擬合3；黃色：頻段擬合3.Fig.9 The radially averaged logarithm power spectrum of the gravity anomaly and its fitting by using piece-wise linearizationThe blue line shows the gravity anomaly power spectrum. The black line shows the fitted straight line in the first radial frequency section, green line for the second section and red line for the third section, the yellow line fitted the last frequency section.

4 結論

本文為實現(xiàn)大距離穩(wěn)定向下延拓引入了約束等效源向下延拓方法，闡述了等效源構建中加入約束的公式.通過模型試驗驗證了該方法的可行性，并且與波數(shù)域迭代向下延拓法和未約束等效源向下延拓法進行比較，體現(xiàn)出了該方法的優(yōu)越性.最后將該技術運用于實際船測數(shù)據(jù)向下延拓，并且與水下實測重力異常進行比較，獲得了良好的匹配效果.

與直接構建延拓算子的方法相比，約束等效源下延法的本質(zhì)是利用源在目標深度的一個正演過程，不會像波數(shù)域迭代這一類的下延方法，需要構建一個濾波算子通過放大信號的方式來放大高頻信息，所以幾乎不會受到噪聲的影響.

圖10 不同方法實測數(shù)據(jù)向下延拓結果圖(a) 水下300 m實測重力異常； (b) 波數(shù)域迭代法延拓結果； (c) 等效源法延拓結果； (d) 約束等效源法延拓結果.Fig.10 Downward continuation results of measured data from different methods(a) The measured gravity anomaly at depth of 300 m under the sea surface; (b) The result of wave number domain iteration method; (c) The result of equivalent source method; (d) The result of equivalent source constrain method.

圖11 實測數(shù)據(jù)向下延拓結果剖面對比圖(a) 測線1延拓結果； (b) 測線2延拓結果； (c) 測線3延拓結果； (d) 測線4延拓結果.藍線：水下實測數(shù)據(jù)；黑線：波數(shù)域迭代法結果；綠線：等效源法延拓法；紅線：等效源約束延拓法.Fig.11 Comparison of different downward continuation methods result of measure data profile(a), (b), (c), (d) represent four measure lines downward continuation results respectively. The blue line shows the measure value. The black line is the result of wave number domain iteration method. The green line represents the result of equivalent source method and the red line is for the result of equivalent source constrain method.

當下延距離過大時，地下空間的密度異常體所產(chǎn)生的高頻信息，由于衰減距離過大導致在觀測面上難以被分辨，所以無論是利用濾波算子放大還是等效源構建正演的方式向下延拓，結果都會與真實值有一定偏差.本文所提出的方法，將地下空間的實測值作為約束點，在解算等效源模型參數(shù)的過程中進行約束，使得最終延拓結果更加接近真實值，在模型試驗與實際數(shù)據(jù)處理中都有很好的體現(xiàn).但地下空間剖分塊體的個數(shù)也會直接影響延拓精度，理論上剖分塊體越精細效果越好，但這樣會增加運算量，后續(xù)的研究會在如何盡可能細化的剖分更多的等效塊體與滿足計算速度之間找到平衡，使得該方法更加準確且高效.