matlab一维外延差值,在波数域计算一维重磁异常导数的Matlab语言算法_肖锋

第32卷第3期物　探　与　化　探V o l.32,N o.3 2008年6月G E O P H Y S I C A L＆G E O C H E M I C A LE X P L O R A T I O N J u n.,2008　在波数域计算一维重磁异常导数的M a t l a b语言算法

肖锋,孟令顺,吴燕冈

(吉林大学地球探测科学与技术学院,吉林长春　130026)

摘要:利用M a t l a b内建的快速傅氏变换函数可以方便地在波数域计算重磁异常导数。介绍了基于M a t l a b语言的波数域求导的算法,给出了程序源代码,讨论了一些有助于提高计算精度的编程技巧。通过模型试验和数据分析,发现在计算垂向导数时波数域求导算法的精度比傅氏级数的精度有明显改善;而水平导数的计算,2种方法的精度相当。在某区钾盐勘探中,用该方法处理高精度重力剖面数据,取得了较好效果。

关键词:波数域;重磁异常;导数;M a t l a b语言

中图分类号:P631 文献标识码:A 文章编号:1000-8918(2008)03-0316-05

在重磁数据的处理中,导数的计算非常重要。它是压制区域场、圈定局部场、分离叠加异常的常用方法[1]。它还广泛地应用于其他方法中,如:重力归一化总梯度法、界面反演的O l d e n b u r g-P a r k e r算法、T i l t导数法等等。重磁异常导数的计算可以在空间域进行,也可以在波数域进行。空间域的导数计算方法由于编程复杂、精度较低等原因,现在已经很少使用了。波数域的导数计算方法在很多文献中都给出[2-5],但是具体到如何编程实现还很少涉及。实际编程过程中,还要解决一些问题才能实现,比如单位换算等。这里着重介绍程序设计,目的是搭建一个桥梁将程序和该算法连接起来。

M a t l a b语言是一种在科研人员中已得到广泛使用的程序设计语言,用它来计算重磁异常导数非常方便。在这里给出波数域的算法和详尽的M a t l a b 计算程序。文中的代码读者可直接利用。

1　算法的数学原理

重磁异常导数的计算在很多文献中被归纳为位场数据的转换[1,5]。对于计算一维(剖面)重磁异常导数而言,就是分别计算实测重磁异常对不同方向上的偏导数。建立直角坐标系x O z,x轴方向代表剖面(或测线)方向,向右为正,z轴正方向向下,代表垂直地表向下的方向。实测重磁异常用T来表示,假设有N个测点,序号从1,2,…,N(M a t l a b语言中表示数组元素的角标也是从1开始的)。第i个测点表示为T i(i=1,2,…,N),测点间距为X s t e p,单位是m。基本原理是先将异常数据做傅氏变换,用它的谱乘以对应的波谱算子,再经过傅氏反变换得到导数。实测重磁异常T经过快速傅氏变换获得它的波谱T,在波数域中求n阶垂向导数主要依据为

n T

z n

=r n T。(1)其中:r代表径向圆波数,对于一维的情况,r等于x 方向的波数u。对于不同测点T i,对应不同的波数r i。相应的波数域n阶垂向导数算子是r n。一维情况下,波数域中对x方向的求导公式如下:

n T

x n

=(i u)n T。(2)对应的波数域n阶水平导数算子是(i u)n。

2　程序流程

波数域中的重磁异常导数计算过程大致可以分为5个部分:扩边、傅氏变换、波谱计算、傅氏反变换、缩边(图1)。下面以重力剖面数据求一阶垂向导数和一阶水平导数为例,详细说明程序的设计。磁法剖面数据的处理与之类似。

2.1　读入数据并扩边

假设原始重力观察数据已经按列存入文本文件,名为:r a w.t x t的记事本文件。其中第一列为沿剖面的坐标,即x坐标(单位:m),第二列为重力值(单位:g.u.),两列之间用空格分隔,而且原始数据经过插值计算,等间距分布。为了减弱F F T的边界效应,应该对边界进行处理。基本原理都是将边界数据归零,以减少间断点。程序采用的是将原始数据读入到数组A中,然后按照余弦公式将数据长度

收稿日期:2007-09-14

基金项目:国家自然科学基金(40774059)