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

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

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

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

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

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

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

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

垂向?qū)?shù)在場分離及其確定場源體邊界等方面有重要作用。介紹了計(jì)算位場各階垂向?qū)?shù)的ISVD算法，并將其在Matlab中實(shí)現(xiàn)，給出了計(jì)算位場各階垂向?qū)?shù)ISVD算法的部分Matlab源代碼。研究表明，利用Matlab可以方便快速地實(shí)現(xiàn)ISVD算法并較為準(zhǔn)確地獲得位場各階垂向?qū)?shù)。同時(shí)，與在頻率域和空間域計(jì)算相比，ISVD算法計(jì)算的位場各階垂向?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年提出并應(yīng)用于實(shí)際資料處理解釋中的計(jì)算位場各階垂向?qū)?shù)的新方法[3]，能夠有效避免在頻率域計(jì)算各階垂向?qū)?shù)時(shí)傅里葉變換本身固有的吉布斯效應(yīng)以及導(dǎo)數(shù)算子的高通濾波作用以及在空間域計(jì)算時(shí)以泰勒級(jí)數(shù)和拉普拉斯方程為基礎(chǔ)的計(jì)算公式產(chǎn)生的截?cái)嗾`差。Matlab是一種科學(xué)的數(shù)值計(jì)算編程語言，具有簡單易學(xué)、函數(shù)庫功能強(qiáng)大和界面友好等特點(diǎn)， 其主要功能包括數(shù)值分析、矩陣運(yùn)算、信號(hào)處理和圖形顯示。利用上述功能編制Matlab程序，可以快速實(shí)現(xiàn)位場數(shù)據(jù)的計(jì)算處理[4-5]。下面，筆者對位場各階垂向?qū)?shù)ISVD算法及其在Matlab中的實(shí)現(xiàn)進(jìn)行了研究。

1 方法原理

1.1在頻率域計(jì)算位場標(biāo)量位

(1)

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

1.2在空間域計(jì)算位場垂向二階導(dǎo)數(shù)

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

(2)

應(yīng)用最小二乘法對曲線進(jìn)行擬合，當(dāng)m=2時(shí)，用5個(gè)點(diǎn)的位場異常值計(jì)算的公式如下：

(3)

2 程序流程

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

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

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

(4)

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

2.2二維傅里葉變換

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

glfs=fftshift(fft2(temp4))

2.3在頻率域計(jì)算

在頻率域進(jìn)行位場轉(zhuǎn)換，即異常值頻譜矩陣glfs(n2,m2)乘以垂向積分算子qv，以獲得轉(zhuǎn)換場的頻譜矩陣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二維傅里葉反變換及計(jì)算位場垂向一階導(dǎo)數(shù)

對轉(zhuǎn)換場頻譜矩陣vfs(n2,m2)進(jìn)行移頻、二維傅里葉反變換和取實(shí)，即可得到轉(zhuǎn)換場矩陣vl(n2,m2)，其轉(zhuǎn)換場即為布格重力異常的位，求其2個(gè)水平二階導(dǎo)數(shù)矩陣vlxx(n2,m2)和矩陣vlyy(n2,m2)，然后依據(jù)拉普拉斯方程求其垂向二階導(dǎo)數(shù)矩陣vlzz(n2,m2)，即為布格重力異常的垂向一階導(dǎo)數(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ù)擴(kuò)邊的逆過程，只取矩陣vlzz(n2,m2)中間與原始數(shù)據(jù)對應(yīng)的部分，即可得到ISVD算法計(jì)算的布格重力異常垂向一階導(dǎo)數(shù)矩陣vzz(n,m)。然后將矩陣轉(zhuǎn)換為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中檢驗(yàn)通過。

3 模型對比試驗(yàn)分析

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

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

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

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

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

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

[2]蔡宗熹，陳維雄，姜蘭.位場數(shù)據(jù)求導(dǎo)精度的提高及其方法[J].物探化探計(jì)算技術(shù)，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]馮京，趙鐵虎，方中華.重力場上延計(jì)算及頻譜分析技術(shù)在Matlab中的應(yīng)用[J].海洋地質(zhì)前沿，2011，21(3):63-68.

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

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

2013-01-23

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

P631.2

A

1673-1409(2013)10-0086-04

[編輯] 李啟棟