matlab 波数域变换,在波数域计算一维重磁异常导数的Matlab语言算法

第 32卷第 3期 物 　探 　与 　化 　探 Vol. 32,No. 3　 　2008年 6月 GEOPHYSICAL & GEOCHEM ICAL EXPLORATION Jun. , 2008　在波数域计算一维重磁异常导数的 Matlab语言算法 肖锋 ,孟令顺 ,吴燕冈 (吉林大学 地球探测科学与技术学院 ,吉林 长春 　130026) 摘 要 :利用 Matlab内建的快速傅氏变换函数可以方便地在波数域计算重磁异常导数。介绍了基于 Matlab语言的波数域求导的算法 ,给出了程序源代码 ,讨论了一些有助于提高计算精度的编程技巧。通过模型试验和数据分析 , 发现在计算垂向导数时波数域求导算法的精度比傅氏级数的精度有明显改善 ;而水平导数的计算 , 2种方法的精度相当。在某区钾盐勘探中 ,用该方法处理高精度重力剖面数据 ,取得了较好效果。 关键词 :波数域 ;重磁异常 ;导数 ;Matlab语言 中图分类号 : P631　　　文献标识码 : A　　　文章编号 : 1000 - 8918 (2008) 03 - 0316 - 05 　　在重磁数据的处理中 ,导数的计算非常重要。它是压制区域场、圈定局部场、分离叠加异常的常用方法 [1 ] 。它还广泛地应用于其他方法中 ,如 :重力归一化总梯度法、界面反演的 Oldenburg2Parker算 法、 Tilt导数法等等。重磁异常导数的计算可以在空间域进行 ,也可以在波数域进行。空间域的导数计算方法由于编程复杂、精度较低等原因 ,现在已经很少使用了。波数域的导数计算方法在很多文献中都给出 [2 - 5 ] ,但是具体到如何编程实现还很少涉及。实际编程过程中 ,还要解决一些问题才能实现 ,比如单位换算等。这里着重介绍程序设计 ,目的是搭建一个桥梁将程序和该算法连接起来。 Matlab语言是一种在科研人员中已得到广泛使用的程序设计语言 ,用它来计算重磁异常导数非常方便。在这里给出波数域的算法和详尽的 Matlab 计算程序。文中的代码读者可直接利用。 1　算法的数学原理 重磁异常导数的计算在很多文献中被归纳为位场数据的转换 [1, 5 ]。对于计算一维 (剖面 )重磁异常导数而言 ,就是分别计算实测重磁异常对不同方向上的偏导数。建立直角坐标系 xOz, x轴方向代表剖面 (或测线 )方向 ,向右为正 , z轴正方向向下 ,代表垂直地表向下的方向。实测重磁异常用 T来表示 , 假设有 N 个测点 ,序号从 1, 2, ⋯, N (Matlab语言中表示数组元素的角标也是从 1开始的 )。第 i个测点表示为 Ti ( i =1, 2, ⋯, N ) ,测点间距为 Xstep ,单位是 m。基本原理是先将异常数据做傅氏变换 ,用它 的谱乘以对应的波谱算子 ,再经过傅氏反变换得到导数。实测重磁异常 T经过快速傅氏变换获得它的波谱gT,在波数域中求 n阶垂向导数主要依据为 nT 5zn = r ngT。 (1)其中 : r代表径向圆波数 ,对于一维的情况 , r等于 x 方向的波数 u。对于不同测点 Ti,对应不同的波数 r i。相应的波数域 n阶垂向导数算子是 r n。一维情 况下 ,波数域中对 x方向的求导公式如下 : nT 5xn = ( iu) ngT。 (2)对应的波数域 n阶水平导数算子是 ( iu) n。 2　程序流程 波数域中的重磁异常导数计算过程大致可以分为 5个部分 :扩边、傅氏变换、波谱计算、傅氏反变换、缩边 (图 1)。下面以重力剖面数据求一阶垂向导数和一阶水平导数为例 ,详细说明程序