韓 利 (西安石油大學地球科學與工程學院，陜西 西安 710065)

段瑞鋒 (陜西省地質(zhì)礦產(chǎn)勘查開發(fā)局物化探隊，陜西 西安 710043)

白 茹 (中冶陜壓重工設備有限公司，陜西 富平 711711)

位場各階垂向?qū)?shù)ISVD算法及其在Matlab中的實現(xiàn)

韓 利 (西安石油大學地球科學與工程學院，陜西 西安 710065)

段瑞鋒 (陜西省地質(zhì)礦產(chǎn)勘查開發(fā)局物化探隊，陜西 西安 710043)

白 茹 (中冶陜壓重工設備有限公司，陜西 富平 711711)

垂向?qū)?shù)在場分離及其確定場源體邊界等方面有重要作用。介紹了計算位場各階垂向?qū)?shù)的ISVD算法，并將其在Matlab中實現(xiàn)，給出了計算位場各階垂向?qū)?shù)ISVD算法的部分Matlab源代碼。研究表明，利用Matlab可以方便快速地實現(xiàn)ISVD算法并較為準確地獲得位場各階垂向?qū)?shù)。同時，與在頻率域和空間域計算相比，ISVD算法計算的位場各階垂向?qū)?shù)特別是高階垂向?qū)?shù)，具有較好的穩(wěn)定性。

位場；各階垂向?qū)?shù)；ISVD算法；Matlab

在位場數(shù)據(jù)的處理解釋中，垂向?qū)?shù)方法在場分離及確定場源體邊界等方面具有重要作用[1-2]。綜合二階垂向?qū)?shù)(integrated second vertical derivative，ISVD)算法是M.Fedi等在2001年提出并應用于實際資料處理解釋中的計算位場各階垂向?qū)?shù)的新方法[3]，能夠有效避免在頻率域計算各階垂向?qū)?shù)時傅里葉變換本身固有的吉布斯效應以及導數(shù)算子的高通濾波作用以及在空間域計算時以泰勒級數(shù)和拉普拉斯方程為基礎的計算公式產(chǎn)生的截斷誤差。Matlab是一種科學的數(shù)值計算編程語言，具有簡單易學、函數(shù)庫功能強大和界面友好等特點， 其主要功能包括數(shù)值分析、矩陣運算、信號處理和圖形顯示。利用上述功能編制Matlab程序，可以快速實現(xiàn)位場數(shù)據(jù)的計算處理[4-5]。下面，筆者對位場各階垂向?qū)?shù)ISVD算法及其在Matlab中的實現(xiàn)進行了研究。

1 方法原理

1.1在頻率域計算位場標量位

(1)

式中，u、v分別為X方向和Y方向的頻率。

1.2在空間域計算位場垂向二階導數(shù)

ISVD算法在空間域用位場的2個水平二階導數(shù)依據(jù)拉普拉斯方程計算位場垂向二階導數(shù)，因此計算垂向二階導數(shù)的關鍵在于計算位場的2個水平二階導數(shù)。 以計算位場x方向水平二階導數(shù)為例，設坐標原點在計算點上，并令：

(2)

應用最小二乘法對曲線進行擬合，當m=2時，用5個點的位場異常值計算的公式如下：

(3)

2 程序流程

頻率域中的位場數(shù)據(jù)處理流程包括讀入數(shù)據(jù)并擴邊、二維傅里葉變換、在頻率域計算、二維傅里葉反變換及計算位場垂向一階導數(shù)和數(shù)據(jù)縮邊并保存。以ISVD算法計算位場垂向一階導數(shù)為例，詳細說明其Matlab程序。

2.1讀入數(shù)據(jù)并擴邊

假設布格重力異常數(shù)據(jù)為BG.dat，即經(jīng)過插值計算的等間距網(wǎng)格數(shù)據(jù)。為了減小二維傅里葉變換的邊界效應，應對邊界進行處理。利用Matlab以矩陣為單位進行數(shù)據(jù)處理的特點，將異常值讀入到矩陣g(n,m)中，然后按照余弦公式：

(4)

將矩陣行列數(shù)擴邊至2的整數(shù)冪，擴邊后知陣temp1(n2,m2)滿足邊界值為零[6]。

2.2二維傅里葉變換

在Matlab中可以利用內(nèi)建函數(shù)fft2快速地對異常值矩陣temp4(n2,m2)實現(xiàn)二維傅里葉變換，然后利用內(nèi)建函數(shù)fftshift進行移頻，即將零頻移到矩陣中心位置，由此得到頻譜矩陣為glfs(n2,m2)。該段程序如下：

glfs=fftshift(fft2(temp4))

2.3在頻率域計算

在頻率域進行位場轉換，即異常值頻譜矩陣glfs(n2,m2)乘以垂向積分算子qv，以獲得轉換場的頻譜矩陣vfs(n2,m2)，其中ax、ay分別為x、y方向網(wǎng)格數(shù)據(jù)間隔，du、dv分別為x、y方向基頻。該段程序如下：

du=1/(m2*ax);dv=1/(n2*ay);

vfs=zeros(n2,m2);

for i=1:1:n2

for j=1:1:m2

if i～=(n2/2+1)||j～=(m2/2+1)

qv=(2*pi*(((j-(m2/2+1))*du)^2+((i-(n2/2+1))*dv)^2)^0.5)^(-1);

vfs(i,j)=g1fs(i,j)*qv;

vfs((n2/2+1),(m2/2+1))=vfs((n2/2+1),(m2/2+1))-vfs(i,j);

end

end

end

2.4二維傅里葉反變換及計算位場垂向一階導數(shù)

對轉換場頻譜矩陣vfs(n2,m2)進行移頻、二維傅里葉反變換和取實，即可得到轉換場矩陣vl(n2,m2)，其轉換場即為布格重力異常的位，求其2個水平二階導數(shù)矩陣vlxx(n2,m2)和矩陣vlyy(n2,m2)，然后依據(jù)拉普拉斯方程求其垂向二階導數(shù)矩陣vlzz(n2,m2)，即為布格重力異常的垂向一階導數(shù)。該段程序如下：

vl=real(ifft2(ifftshift(vfs)));vlxx=zeros(n2,m2);vlyy=zeros(n2,m2);vlzz=zeros(n2,m2);

for i=3∶1∶n2-3

for j=3∶1∶m2-2

vlxx(i,j)=2*(vl(i,j+2)+vl(i,j-2)-vl(i,j))/(7*ax^2)-(vl(i,j+1)+vl(i,j-1))/(7*ax^2);

vlyy(i,j)=2*(vl(i+2,j)+vl(i-2,j)-vl(i,j))/(7*ay^2)-(vl(i+1,j)+vl(i-1,j))/(7*ay^2);

vlzz(i,j)=-(vlxx(i,j)+vlyy(i,j));

end

end

2.5數(shù)據(jù)縮邊并保存

這是數(shù)據(jù)擴邊的逆過程，只取矩陣vlzz(n2,m2)中間與原始數(shù)據(jù)對應的部分，即可得到ISVD算法計算的布格重力異常垂向一階導數(shù)矩陣vzz(n,m)。然后將矩陣轉換為dat文件格式，并保存。該段程序如下：

vzz_temp=zeros(n,m);

vzz_temp=vlzz(cl1∶n2-cl2,rl1∶m2-rl2);

vzz_temp2=reshape(rot90(vzz_temp,-1),n*m,1);

vzz=cat(2,temp1(∶,1∶2),vzz_temp2);

save vzz.dat vzz -ascii;

程序代碼中n、m、ax、ay分別為已知參數(shù)，以上代碼在Matlab7.0中檢驗通過。

3 模型對比試驗分析

單個棱柱體模型長20km，寬5km，頂部深度4km，底部深度10km，密度差為0.3×10-3kg/cm3，為了更加貼近實際地質(zhì)情況，給原始重力異常值增加異常值1.25%的高斯噪音，即信噪比為80，采樣間距0.2km。單個地質(zhì)體模型如圖1所示。

圖1 單個地質(zhì)體模型

分別應用ISVD算法以及在頻率域和空間域計算模型重力異常的各階垂向?qū)?shù)，計算結果如圖2所示。從圖2可以看出：①對于重力異常垂向一階導數(shù)，3種方法計算的結果相差不大。②在頻率域計算的垂向二階導數(shù)由于垂向?qū)?shù)算子的高通濾波作用，將高頻噪音放大。在空間域采用羅森巴赫Ⅱ式計算的重力異常垂向二階導數(shù)與ISVD算法計算的結果相差較小，這是因為應用ISVD算法計算位場二階垂向?qū)?shù)時，也是在空間域計算，只是兩者采用的公式不同。③在頻率域計算重力異常三階導數(shù)時，高頻噪音被嚴重放大，有效信號完全淹沒在噪音當中。在空間域計算重力異常垂向三階導數(shù)，由于截斷誤差的傳遞疊加，噪音較為嚴重，效果也不理想，而ISVD算法計算的重力異常垂向三階導數(shù)，雖然也有噪音，但是其效果明顯好于前2種方法。

ISVD算法之所以比在頻率域和空間域計算位場各階垂向?qū)?shù)穩(wěn)定性要好，首先是因為位場在頻率域內(nèi)進行垂向積分，即乘以垂向積分算子，而垂向積分算子相當于圓滑濾波器，具有圓滑濾波作用，在一定程度上可以濾除部分高頻成分，從而提高計算垂向?qū)?shù)的穩(wěn)定性。另一方面是通過在空間域計算位場的2個水平二階導數(shù)，依據(jù)拉普拉斯方程計算垂向二階導數(shù)，有效地避免了在頻率域計算水平二階導數(shù)產(chǎn)生的震蕩，同時也避免了直接依據(jù)哈克公式或其它空間域垂向?qū)?shù)公式計算垂向二階導數(shù)產(chǎn)生的截斷誤差。

圖2 不同方法計算的各階垂向?qū)?shù)

[1]曾華霖.重力場與重力勘探[M].北京：地質(zhì)出版社，2009.

[2]蔡宗熹，陳維雄，姜蘭.位場數(shù)據(jù)求導精度的提高及其方法[J].物探化探計算技術，1991，13(1)：47-52.

[3]Fedi M，F(xiàn)lorio G.Detection of potential fields source boundaries by Enhanced Horizontal Derivative method[J].Geophysical Prospecting，2001，49：40-58.

[4]馮京，趙鐵虎，方中華.重力場上延計算及頻譜分析技術在Matlab中的應用[J].海洋地質(zhì)前沿，2011，21(3):63-68.

[5]肖鋒，孟令順，吳燕岡.在波數(shù)域計算一維重磁異常導數(shù)的Matlab語言算法[J].物探與化探，2008，32(3):316-320.

[6]王云專，王潤秋.信號分析與處理[M].北京：石油工業(yè)出版社，2008.

2013-01-23

韓利(1987-)，男，碩士生，現(xiàn)主要從事應用地球物理方面的研究工作。

P631.2

A

1673-1409(2013)10-0086-04

[編輯] 李啟棟