張志厚 ，廖曉龍 ，姚 禹 ，范祥泰 ，路潤(rùn)琪

（西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756）

重磁位場(chǎng)向下延拓技術(shù)可以將海平面磁測(cè)數(shù)據(jù)大深度向下延拓至海底，將航空重磁數(shù)據(jù)延拓至地形線，不僅突出了有用信息，還可以彌補(bǔ)條件惡劣地區(qū)重磁場(chǎng)資料的不足，并且重磁位場(chǎng)的向下延拓運(yùn)算是許多重磁數(shù)據(jù)處理的核心步驟[1-5]. 此外，重磁位場(chǎng)的向下延拓技術(shù)在軍事領(lǐng)域具有非常重要的作用[6-9].

重磁位場(chǎng)向下延拓的主要方法從類型上大體可以分為泰勒級(jí)數(shù)法[10-13]、等效源法[14-20]、積分迭代法[1-4,21-23]、Tikhonov正則化法[24-29]、Landweber正則化法[26,28,30]、奇異值分解法[26]等.

泰勒級(jí)數(shù)法的主要缺點(diǎn)是更高階次的導(dǎo)數(shù)計(jì)算會(huì)導(dǎo)致數(shù)據(jù)的高頻成分被無(wú)限放大，曾小牛等[31]對(duì)求導(dǎo)的頻譜函數(shù)進(jìn)行了改進(jìn)，但階次越多誤差也被疊加得越多，最終導(dǎo)致下延誤差反而很大；等效源法的主要缺點(diǎn)是計(jì)算時(shí)間長(zhǎng)，且延拓效果受到“源”的類型、個(gè)數(shù)、深度及其組合的影響.

位場(chǎng)向下延拓的積分迭代法[21,32]計(jì)算效果好，諸多學(xué)者研究了其魯棒性等性能[33-36]，該方法在數(shù)學(xué)上也被稱為對(duì)稱核第一類Fredholm積分方程的逐次逼近解法[37]. Tikhonov正則化法相比Landweber迭代法具有最好的綜合下延性能[26,28]. 然而Tikhonov正則化法存在正則化參數(shù)的選取問題.

空間域的位場(chǎng)向下延拓實(shí)質(zhì)上是求解一個(gè)第一類Fredholm型積分方程，該方程最突出的特性是其“不適定”性. 向下延拓方程離散化后的系數(shù)矩陣為實(shí)對(duì)稱矩陣. 本文首先證明了該系數(shù)矩陣為對(duì)稱的雙重Toeplitz系統(tǒng)矩陣 （blak-Toeplitz-Toeplitzblock, BTTB）；其次假定在其正定的條件下，采用Barzilai-Borwein （BB）法[38-39]求解該線性方程組，并對(duì)迭代步長(zhǎng)進(jìn)行限制，保證了假定系數(shù)矩陣為正定情況下算法的收斂性，實(shí)現(xiàn)了位場(chǎng)的向下延拓計(jì)算，在迭代計(jì)算中使用BTTB與向量乘積的快速算法[40]對(duì)大數(shù)據(jù)進(jìn)行計(jì)算，一定程度上節(jié)約了計(jì)算成本，并對(duì)該方法進(jìn)行了理論模型無(wú)噪聲數(shù)據(jù)與實(shí)際資料的檢驗(yàn)，以及對(duì)比分析了該方法的收斂速度與位場(chǎng)向下延拓的積分迭代法的收斂速度.

1 空間域位場(chǎng)向下延拓

假設(shè)位場(chǎng)觀測(cè)面為笛卡兒直角坐標(biāo)系o-xyz的水平面 （z=0 ），該平面上的位場(chǎng)為f0(x,y)，計(jì)算z=h（h> 0）平面的位場(chǎng)fh(x,y)，且該平面為觀測(cè)面以下至場(chǎng)源以上空間（z軸向下）.

根據(jù)向上延拓公式，可得

表示成褶積的形式為

2 位場(chǎng)延拓系數(shù)矩陣的性質(zhì)

陳龍偉等[41]較為詳細(xì)地證明了位場(chǎng)向下延拓系數(shù)矩陣為對(duì)稱的分塊Toeplitz矩陣，本文提供一種簡(jiǎn)潔的證明方法. 已知f0(x,y) 通過(guò)式（1）求解fh(x,y)即為向下延拓，由于已知場(chǎng)與待求場(chǎng)皆為離散化數(shù)據(jù)，且f0(x,y) 與fh(x,y)離散化網(wǎng)格規(guī)格一致，因此將

式（1）離散化[42]為

式中：Δx、Δy分別為方向x、y的離散化網(wǎng)格間距；M、N分別為方向x、y的離散化網(wǎng)格點(diǎn)總數(shù)，m=1,2,3,···,M，n=1,2,3,···,N.

由于離散化網(wǎng)格間距固定，為了方便表示，將f0(mΔx,nΔy) 與fh(iΔx,jΔy) 分 別 寫 成f0(m,n)與fh(i,j) ，將式（4）中的f0和fh重排為列向量，分別記為

式（4）用矩陣表示為

式中：A為位場(chǎng)向下延拓的系數(shù)矩陣.

根據(jù)式（4），位場(chǎng)向下延拓的系數(shù)為

借鑒“格架函數(shù)”概念[43]，發(fā)現(xiàn)Q也具有平移等效性、互換對(duì)稱性，等價(jià)關(guān)系為

式（8）左右兩端分別代入式（7）即可得證. 根據(jù)式（8），4位數(shù)組Q(m,n,i,j)的前兩個(gè)參量都是常數(shù)1，故其可表示為二維數(shù)組，即將Q(m,n,i,j)簡(jiǎn)記為Q|m?i|+1,|n?j|+1.Q為四維數(shù)組，以行列式的形式排列即可得到：

可以看出：A為對(duì)稱矩陣，大小為MN×MN，內(nèi)部包含了N2個(gè)分塊的對(duì)稱矩陣，大小為M×M，即矩陣A為雙重Toeplitz系統(tǒng)矩陣.

3 BB法

BB法是一種特殊的最速下降法[44]，假設(shè)A正定，式（7）采用BB法求解，位場(chǎng)向下延拓BB法的迭代格式為

式中：k為迭代次數(shù)為第k次向下延拓的計(jì)算值；t(k)為第k次的迭代步長(zhǎng)；p(k)為第k次迭代待延拓場(chǎng)向上延拓值與已知場(chǎng)的差值向量

迭代停機(jī)準(zhǔn)則：當(dāng) (p(k))Tp(k)<τ，即Cτ 時(shí)可停止計(jì)算，其中 τ為給定的適當(dāng)小正數(shù)，C為常數(shù)，此時(shí)

若式（11）中的t(k)為常數(shù)時(shí)，上述迭代方法即為位場(chǎng)向下延拓的積分迭代法. 而向下延拓BB法相當(dāng)于變步長(zhǎng)的積分迭代法；王順杰等[36]采用不動(dòng)點(diǎn)原理證明了位場(chǎng)向下延拓積分迭代法的收斂性，但迭代步長(zhǎng)的取值范圍在0～2. 因此在求解式（4）時(shí)，將其迭代步長(zhǎng)約束在0～2范圍內(nèi)可確保向下延拓BB法的魯棒性. 考慮到迭代過(guò)程中必然會(huì)存在計(jì)算數(shù)據(jù)量大的問題，采用對(duì)稱Toeplitz矩陣與向量相乘的快速算法[40]對(duì)式（11）中的和進(jìn)行計(jì)算，有效地提高計(jì)算效率.

4 模型試驗(yàn)

4.1 二維模型

二維模型為無(wú)限延伸的水平圓柱體，正演計(jì)算兩個(gè)不同高度觀測(cè)線的理論重力異常，其中觀測(cè)線與模型體垂直，對(duì)計(jì)算結(jié)果和收斂速度進(jìn)行分析.

圓柱體模型的半徑為500 m，線密度 λ=1000kg/m3，模型的中心埋藏深度為3 km. 計(jì)算點(diǎn)距為100 m，計(jì)算了z=0 和z=2km的理論重力異常. 將z=0理論異常分別采用BB法與積分迭代法向下延拓至z=2km. 理論重力異常與不同方法向下延拓后的結(jié)果如圖1所示.

z=2km 理論重力異常的最大值與最小值分別為10.477 × 10?5、0.064 × 10?5m/s2；z=0理論重力異常的最大值與最小值分別為3.492 × 10?5、0.182 ×10?5m/s2. 結(jié)果表明：1） 迭代10次BB法計(jì)算結(jié)果的最大、最小值分別為8.645 × 10?5、?0.015 × 10?5m/s2，積分迭代法計(jì)算結(jié)果的最大、最小值分別為7.715 ×10?5、?0.164 × 10?6m/s2；2） 迭代50次BB法計(jì)算結(jié)果的最大、最小值分別為9.542 × 10?5、?0.173 × 10?6m/s2，積分迭代法計(jì)算結(jié)果的最大、最小值分別為9.132 ×10?5、?0.171 × 10?6m/s2. 由此可見，相同迭代次數(shù)下，BB法的計(jì)算結(jié)果與理論重力異常更為接近.

為了定量地評(píng)價(jià)本文所提方法的計(jì)算精度、穩(wěn)定性，采用不同迭代次數(shù)計(jì)算結(jié)果的均方誤差來(lái)評(píng)估其效果，均方誤差為

式中：fcon,l和ftheo,l分別為第l個(gè)計(jì)算點(diǎn)的延拓值和理論值；L為計(jì)算的點(diǎn)數(shù).

圖2為位場(chǎng)向下延拓的BB法與積分迭代法在不同迭代次數(shù)下計(jì)算結(jié)果的均方誤差. 可見兩種方法都具有較好的收斂性，但向下延拓BB法的收斂速度更快.

圖1 理論重力異常與不同方法向下延拓結(jié)果對(duì)比Fig. 1 Theoretical gravity anomaly compared with downward continuation results of different methods

圖2 重力模型計(jì)算結(jié)果的均方誤差Fig. 2 Mean square error calculated by gravity model

4.2 三維模型

三維模型為組合圓球型磁異常體，通常被用于位場(chǎng)向下延拓效果的檢驗(yàn)[45-47]. 首先通過(guò)正演計(jì)算兩個(gè)觀測(cè)面的磁異常，然后采用不同計(jì)算方法將距離場(chǎng)源較遠(yuǎn)處的平面磁異常向下延拓至另一平面，通過(guò)均方誤差來(lái)評(píng)估計(jì)算效果. 三維組合圓球形磁異常體參數(shù)見表1. 表中：(x0,y0,h0)為圓球型異常體中心坐標(biāo)；r0為圓球異常體半徑；M0異常體磁化強(qiáng)度；I0和A0分別為地磁場(chǎng)傾角和剖面磁偏角；I′和A′分別為異常體磁化強(qiáng)度的傾角和磁北夾角. 根據(jù)圓球型磁異常體的正演計(jì)算解析式[48]，分別計(jì)算兩個(gè)不同觀測(cè)面上的磁異常，如圖3所示，其中虛擬觀測(cè)區(qū)網(wǎng)格點(diǎn)個(gè)數(shù)為401× 401，Δx與 Δy都為50 m.

表1 三維組合圓球型磁異常體參數(shù)Tab. 1 Magnetic anomaly parameters of 3D combination sphere

圖3 不同高度面上的磁異常等值線（單位：nT）Fig. 3 Magnetic anomaly isolines at different planes (unit： nT)

將位于z= 0高度的理論磁異常分別采用BB法以及積分迭代法下延至z= 1 km的平面，計(jì)算結(jié)果見圖4所示. 其結(jié)果除了邊部的零等值線略有偏差，其余部分都與理論場(chǎng)（圖3（b））高度吻合.

圖4 理論磁異常向下延拓后的計(jì)算結(jié)果（單位：nT）Fig. 4 Calculation results of theoretical gravity anomaly after downward continuation (unit： nT)

兩種方法計(jì)算結(jié)果的均方誤差如圖5所示，可以看出：隨著迭代次數(shù)的增加，均方誤差都在急劇衰減，BB法相比積分迭代法的收斂速度更快；在相同迭代次數(shù)的情況下，位場(chǎng)向下延拓BB法的計(jì)算結(jié)果的精度高于積分迭代法；同一收斂精度下，當(dāng)均方誤差在0.4 nT以下時(shí)，BB法需要迭代次數(shù)15次，而積分迭代法需要迭代35次以上，可見向下延拓BB法的計(jì)算速度提高了2倍.

圖5 三維磁場(chǎng)模型計(jì)算結(jié)果的均方誤差Fig. 5 Mean square error calculated by3D magnetic field model

此外，隨著迭代次數(shù)增加，位場(chǎng)向下延拓BB法與積分迭代法的計(jì)算結(jié)果精度趨近，兩種方法的本質(zhì)都是迭代法，迭代法求解反問題時(shí)都會(huì)出現(xiàn)所謂的“半收斂”現(xiàn)象[44]，即在迭代的早期階段，估計(jì)解可以得到很好的改善，但其長(zhǎng)期行為會(huì)出現(xiàn)“飽和”效應(yīng)，甚至?xí)呦驘o(wú)序[49-50]，故而兩種方法的計(jì)算精度極有可能會(huì)越來(lái)越相近.

5 計(jì)算速度

在主頻2.30 GHz、內(nèi)存1024 MB的計(jì)算機(jī)上，1024 × 1024個(gè)數(shù)據(jù)向下延拓迭代50次，BB法計(jì)算時(shí)間約為130 s. 另外，一般情況下依據(jù)經(jīng)驗(yàn)取一定的迭代次數(shù)作為停機(jī)準(zhǔn)則，BB法相比積分迭代法只需要較少的迭代次數(shù)便可達(dá)到相同的計(jì)算精度.

6 實(shí)際資料檢驗(yàn)

通過(guò)某銅鐵礦區(qū)實(shí)測(cè)網(wǎng)格化磁異常資料來(lái)評(píng)估下延計(jì)算效果，如圖6所示. 圖6（a）為已知實(shí)測(cè)的網(wǎng)格化磁異常，點(diǎn)距為50 m. 將該異常向上延拓10倍點(diǎn)距（向上延拓方法采用快速傅里葉變換方法），其結(jié)果如圖6（b）所示. 隨后再將該高度上的異常分別采用BB法與積分迭代法向下延拓500 m至原始高度，其結(jié)果如圖6（c）、（d）所示（兩種方法迭代次數(shù)相同），計(jì)算結(jié)果都與圖6（a）基本一致.

圖6 實(shí)例驗(yàn)證（單位：nT）Fig. 6 Validation by measured data (unit： nT)

原始場(chǎng)的最大、最小值分別為1518.3、?669.2 nT.對(duì)應(yīng)的，向下延拓BB法計(jì)算結(jié)果的最大、最小值分別為1507.7、?640.7 nT，積分迭代法計(jì)算結(jié)果的最大、最小值分別為1497.3、?623.5 nT. 用式（11）計(jì)算BB法和積分迭代法下延結(jié)果的均方誤差分別為24.4 nT和29.7 nT. 相比原始場(chǎng)的最大、最小值，兩種方法的均方誤差都較小，但BB法的精度更高.

為了進(jìn)一步評(píng)價(jià)兩種方法的計(jì)算精度，采用平均相對(duì)誤差 ε進(jìn)行評(píng)估，如式（13）所示.

用式（13）計(jì)算本文所提BB法和積分迭代法的平均相對(duì)誤差分別為6.1%、7.7%，可見在實(shí)際應(yīng)用中本文所提方法的效果優(yōu)于積分迭代法.

7 結(jié) 論

1） 本文首先證明了位場(chǎng)向下延拓的系數(shù)矩陣為對(duì)稱的雙重Toeplitz系統(tǒng)矩陣，該性質(zhì)可為今后的位場(chǎng)數(shù)值模擬節(jié)省成本.

2） 實(shí)現(xiàn)了一種位場(chǎng)向下延拓的計(jì)算方法，即位場(chǎng)向下延拓的BB法，該方法主要優(yōu)點(diǎn)的計(jì)算精度高、速度快. 相比普遍采用的Tikhonov方法，文中方法不存在正則化因子的選擇問題.