【doc】重力位谱分析及重力异常导数换算新方法——余弦变换

【doc】重力位谱分析及重力异常导数换算新方法——余弦变换 重力位谱分析及重力异常导数换算新方法 ——余弦变换 第49卷第1期 2006年1月 地球物理 CHINESEJOURNALOFGEOPHYSICS Vo1.49.No.1 Jan.,2O06 张凤旭,盂令顺,张风琴等.重力位谱分析及重力异常导数换算新方法——余弦变换.地球物理,2OO6,49(1):244—248 ZhangFX,MengLS,ZhangFQ,eta1.Anewmethodforspectralanalysisofthepotentialfields ndconversionofderivativeofgravity. anomalie8:cosinetrai~onl'1.ChineseI,.Geop伽.(inChinese),2006,49(1):244—248 重力位谱分析及重力异常导数换算新方法 —— 余弦变换 张凤旭,孟令顺,张凤琴,刘财,吴艳冈,杜晓娟 1吉林大学地球探测科学与技术学院,长春130026 2吉林大学物理学院,长春130023 摘要为了提高重力异常导数换算的精度,真实有效地反映地质体的异常特征,提出用余弦变换计算异常导数 的新方法.给出并证明了两个定理,利用它们推导出重力位余弦谱一般表达式以及重力异常各阶导数计算公式,建 立了位场余弦谱分析理论.模型实验中发现,用Fourier变换计算的一阶导数与理论导数偏差很大,而余弦变换计 算的导数与理论异常导数拟合效果非常好,除边界几个数据因重力异常的有限截断产生的吉布斯效应残留使误差 较大外数据的计算精度均很高,误差为一O.o9%一5%. 关键词重力位余弦变换谱,异常导数,余弦变换,精度 文章编号0001—5733{:~oo6}ol一0244—05中图分类号P631收稿日期2004—11 —23,2OO5—10—12收修定稿 AnewmethodforspectralanalysisofthepotentialfieIdandconversionof ' derivativeofgravity-anomalies:cosinetransform ZHANG.Feng?Xu,MENGLing.Shun,ZHANGFeng.Qin 一,LIUCai,WuYan.Gang,DUXiao.Juall lGeo-ExplorationScienceandTechnologyInstitute,JilinUniversity,Changchun130026,China 2CollegeofPhysics,JilinUniversity,Changchun130023,China AbstractInordertoimprovetheaccuracyofderivativeconversionofgravityanomaliesandreflectanomaly characteristicofgeologicbodieseffectively,weproposeanewmethodofcalculatinganomalyderivativeusing cosinetransform.Twotheoremsareputforwardandprovedandthecommonexpressionofthecosinetransform spectrumofthegravitypotentialfieldandtheformulatocalculatederivativesofgravityanomaliesarededuced, S0thetheoryofcosine?transform?spectrumofthepotential— fieldisestablished,Inmodelexperiments.wefind thatthedeviationofthefirstderivativecalculatedbyFouriertransformisverylargetocomparingwiththe theoreticalderivative,butthefittingeffectofthederivativeofanomaliescalculatedbycosinetransfomisvcry g00d.Thecalculationaccuracyofdataareallveryhighexceptthaterrorsofseveraldata0ntheboundaryare largebecauseoftheresidualGibbuseffectinducedbyfinitetruncationofgravityanomalies.Errorsare — O.O9%,5%. KeywordsCosine—transformspectrumofgravitypotentialfield,Derivativeofanomaly,Cosinetransform, Accuracy 基金项目大庆探区外围中新生代断陷盆地群演化与油气远景项目(XQ.2007.07.05) 及吉林大学"九八五工程"项目资助. 作者i畸介张风旭,男,1969年生,讲师,1995年毕业于长春地质学院应用地球物理专业,现在吉林大学攻读博士学位,主要从事固体地球物理 学研究.E-mail:三hal】出@jlu.edu.cn 《}}}}}}f 1期张风旭等:重力位谱分析及重力异常导数换算新方法——余弦变换 1引言 在物探数据处理中,重磁异常导数换算具有重 要意义.这是因为重磁异常的导数在不同形状地质 体上有不同的特征,有助于对异常的划分和解释;它 能够突出浅而小的地质体异常特征而压制区域性深 部地质因素的影响,同时在一定程度上可以划分不 同深度和大小异常源产生的叠加异常,且适当地选 取导数的阶次,可提高异常的分辨率?]. 异常导数的换算大多在空问域和波数域中进 行.空间域主要理论依据为函数的微分定理,用最 小二乘法对实测值进行函数拟合来实现,还可利用 拉格朗日插值法?j,样条函数插值法H计算垂向导 数,但这些方法的计算精度都受到一定的限制.波 数域大多采用Fourier变换],利用谱分析原理来 实现,其理论较为完善,但由于函数非周期性因子以 及有限截断的影响,使其应用范围受到很大的限制, 而且计算的精度也不高. 1974年Ahmedeta1.提出了离散余弦变换 (DCT),使DCT在语音,图像编码以及数据压缩等方 面得到了广泛的应用'.但近30年来,余弦变换 在地球物理学方面的开发和应用还很不足,尤其是 利用余弦变换研究位场谱和计算重磁异常导数的问 题,在物探数据处理理论中还是一个盲区. 余弦变换与Fourier变换一样都属于正交变换. 但余弦变换不但避免了复数运算,而且它有着与K. L变换(Karhunen.LoSvetransform)相似的性能,能够 去除原信号的相关性,从而保留原信号的最大能 量n,同时随着DCT快速算法的研究和实现n, 使其计算的速度远远超过了Fourier变换.由此可见 余弦变换有其独特的优越性,必然能够提高数据处 理的速度和精度. 2位场余弦变换谱理论的建立 2.1两个定理的证明 用余弦变换计算异常的导数,目前在物探数据 处理理论中依然是一个盲区,因此本文提出并证明 下面两个有用定理,以便由此推导出重力位余弦谱 一 般表达式,并引出重力异常导数计算的基本公式. 定理1如果函数f(t),g(t)满足狄利希里条 件,且绝对可积,而且其中f(t)或g(t)至少有一个 函数为偶函数,则它们的褶积的余弦变换式为 c[f(t)*g(t)]fc(e.O)g.(e.O),(1) 式中C表示余弦变换,09=2厂为角频率. 证明:令F(e.O)=C[f(t)*g(t)],由于函数 f(t),g(t)满足狄利希里条件,且绝对可积,因此其 余弦变换存在,根据余弦变换的定义有 ?? F(e.O)=IIf(r)g(t—r)drcos(mt)dt.:: 设t1=t—r,则t=t1+r,dt:dt1,将它们代入上式得 F():-f-f,(r)g(t)d~eos,o(t+r)dt = -r-r厂(r)g(t,)dr(c.scutc.scr — sino~t1sinwr)dtI = Ig(f1)eosmtldtlIf(r)coso~rdr :: 一 Ig(1)sinwt1dtlIf(r)sino~rdr.:: 由于f(t)或g(t)至少有一个函数为偶函数,考虑积 分函数的奇偶性,上式的第二项为零,所以 F(e.O)C[,(t)g(t)]=fc(叫)g(e.O), 表达式1得证. 定理2如果函数. 厂(t)处处可微分任意多次, 且各阶导数的Fourier变换式存在,则f(t)各阶导数 的余弦变换式为 ,, f(一1)丁e.OC[f(t)]=1,3,5,?? c[(){,【 (一1)(uC[f(t)]n=2,4,6,… (2) 其中Cn/2[Jr(,)]= . fJr(,)c.s(号一,)d,表示相移 /2的余弦变换式.为区别于正弦变换而给出如上 余弦形式. 表达式2的推导过程比较简单,不做详细论述. 2.2重力位余弦谱一般表达式推导 设直角坐标系oxyz,轴垂直向下,地下有任意 密度P为 形体地质体,其空间范围定义为,剩余 ,,的函数,那么该地质体在计算点P(,Y,z) 的重力位表达式 (,Y,z)=GIP(,7,) ×F.(3)×iF' 如果把密度函数lD(,,)的定义域扩大到整 地球物理(ChineseJ.Geophys.)49卷 个三维空间,定义以外空间的』D=0,那么(3)式右式(2)导出,即 边就变成了一个褶积积分的形式,g(u,,)=(2~/_)g(u,,z),(11) 令,y'南,则有 V(,Y,z)=』D(,Y,)*R(,Y,z).(4) 显然,R(一,一Y,一z)=R(,Y,),R为偶函 数,因此(4)式满足表达式1的条件,由(1)式得 (H,,)=lDc(H,口,)Rc(H,口,),(5) 其中,",,分别为,Y,z方向的波数. 因为密度函数lD(,,)的定义域已经扩大到 整个三维空间,由微分定理可得 :++ .aVo =一 47t印(,Y,),(6) 对该式取余弦变换,由(5)式及表达式2的二阶水平 导数余弦变换知 C(V)=一47c2("+口+)(H,口,) =一 47t印c(H,口,),(7) 于是 (,,):.(8) (8)式即为位场的余弦谱计算公式,但在通常的 物探资料处理中,一般是在平面上观测并对资料做 谱分析,因此要找出为常量下的位场余弦变换谱 表达式.为此,对波数W做余弦逆变换,由式(1)得 (U,,z)=(H,,)』D(",,), 而 R(u,,)=-f-i—c.s(2wz)d G2I; .' 所以 z):f衙"一.(9) 上式就是重力位余弦谱一般表达式. (9)式对做/1,阶微分,则可推导出重力位垂向 /1,阶导数 : (2厢)(,).(10)a 2.3重力异常导数的余弦谱表达式 在重力资料数据处理过程中,接触的实测数据 大多是重力异常数据,其各阶导数均由重力异常导 出,下面给出三度体重力异常各阶导数余弦谱公式, 其中垂向导数由公式(10)给出,水平方向导数由公 g(,,)= g(,,)= 其中 (一)(27c)C1[g(,Y,z)](一)T(27c)[(,,z)] =1,3,5,…, (一1){(27tu)c[g(,Y,)] n=2,4,6,…, (12) (一1)(2xv)C[g(,y,z)] =1,3,5,…, (一1)号(27c)C[g(戈,y,)] n=2,4,6,…, (13) [g(]=-f-f,ycos【号 一 (u+uy)】dxdy. 令(11),(12)式中的=0就可以得到二度体导 数余弦谱公式. 3重力异常导数的实现 3.1离散余弦变换 Ahmedeta1.首先给出了一维离散余弦变换 的定义.给定一维数据序列{(17,):n=0,1,…, N一1},其离散余弦变换和逆变换定义为 Xc?:))c0s,)(JI})=?c(JI})?(n)c0s,a) ): n=o , (14 k (n)??=O c(JI})s,(14b) 式中,=0,1,2,…,N一1, c(JI}){?20,(14c)LIk0 仿照一维离散余弦变换,对于二维数据序列 {(i,):i=0,1,…,N一1;_『=0,1,…,M一1),其 二维离散余弦正,反变换分别定义为 Xc(m,n)=c(m,n)【??(f,)NM? ×COSCOS 】,×———一一I, (15a) 【c() }};f}l}}}}}}}}}}}}《ll}}}}}《}}l}l}》lll}_i;《" 1期张凤旭等:重力位谱分析及重力异常导数换算新方法——余弦变换247 ×c.sc.s ],(15b)×..————一...———芝—一J, C中m,i:0,1,…,N一1;n,=0,1,…,M一1; ccm,n={:?2,n麓..c-5c 3.2异常导数的数值实现 对于二度体重力异常垂向导数的计算,首先采 用一维离散余弦变换公式(14a)计算出重力异常余 弦谱,然后代入到公式(11)的二维公式中得到垂向 导数余弦谱,最后由余弦逆变换公式(14b)计算出垂 向导数;三度体异常垂向导数的计算步骤同上,其 正,反变换过程采用二维余弦变换公式(15a,b)计算. 对于水平方向导数谱,由于其各阶导数的余弦 变换均发生了相位的偏移,从表达式2中可以推出, 水平n阶导数,相位偏移为nn/2,因此不能笼统地 采用(14a)和(15a)进行计算.理论和实验分析证明, 当水平导数的阶数为偶数时,离散余弦正变换部分 采用(14a)和(15a)式进行;当水平导数的阶数为奇 数时,其离散余弦谱变换式部分采用下式进行计算, ?:舢)莹)sin,(1l})=?c(后)?(n),n=O (16) Xc(m,n)=c(m,n) LNM × 【sin ,) ×sin】 其中C(),C(m,n)取值同(14c),(15c). 把以上计算结果分别代入到(12),(13)式,得到 相应的二度,三度体重力异常水平导数余弦谱,然后 利用公式(14b),(15b)对其实施余弦逆变换便可计 算出n阶水平方向导数. 应该说明的是,尽管(16),(17)式采用了正弦表 达式,但它们完全不同于离散正弦变换n,在计算 中应给予区别. 4模型计算与精度分析 前文从理论上推导并证明了利用余弦变换计算 重力异常各阶导数的可行性.下面以无限长水平圆 柱体为例,利用Fourier变换,余弦变换计算的水平 圆柱体重力异常垂向,水平一阶导数(g,g)同其理 论一阶导数进行对比分析,研究用余弦变换计算的 重力异常一阶导数的精度及其实用价值. 采用余弦变换计算重力异常导数同Fourier变 换一样也会产生吉布斯效应和皱波现象.吉布斯效 应是由于函数的突然截断而产生的,利用镶边函数 可以弥补该效应产生的误差;皱波现象是由于异常 导数谱变换过程中的褶积产生的,选择适当的离散 数据取样间隔便可以消除这种现象,如果取样间隔 过大,应采用适当的数据插值方法来缩小数据间隔. 下面用余弦变换和Fourier变换计算的异常导数,均 采取了同样的镶边函数和采样间隔,以便在同等的 条件下分析二者的计算精度. 计算图1的模型为无限长水平圆柱体,其模型参 数为:半径R:1.0kin;圆柱体中轴线坐标(.,d),. =1Okm,d=2.0km;剩余密度p=O.5g/cm3.图1的 理论垂向和水平一阶导数分别采用如下公式计算: 'o)=,(18) 枷)=.(19) 图1给出了两种方法计算的与理论垂向,水平 一 阶导数对比分析图,图中可以清晰地看到,利用余 弦变换计算的异常导数(图1中c)与理论异常导数 n拟合效果非常好,而利用Fourier变换计算的异常 导数b尽管与理论导数曲线走向相似,但其偏离程 度非常大. 由于Fourier变换计算的异常导数偏差很大,其 误差精度分析的意义不大,因此表1只给出采用余 弦变换计算的无限长水平圆柱体的g和g的误差 表1利用余弦变换计算的无限长水平圆柱体g.g,精度 Table1Accuracyofg,gofinfinite cylindercalculatedbycosinetraII1n 248地球物理(ChineseJ.Geophys.)49卷 e 舍0 - 20 O246810121416l820 O2468101214161820 图1不同方法计算的无限长水平圆柱体 一 阶导数对比分析图 (a)垂向一阶导数;(b)水平一阶导数. 0理论值;bFourier变换;c余弦变换. Fig.1Contrastanalysisofverticalfirst-orderderivatives ofnilinfinitecylindercalculatedbydifferentmethods 分析.从表中可以看出,除边界几个数据因重力异 常的有限截断产生的吉布斯效应残留使误差较大 外,数据计算精度均很高,误差为一0.09%,5%. 5结论 本文从理论上证明并推导了用余弦变换计算重 力异常导数的公式,模型实验证实了方法的可靠性 和实用性.通过精度对比分析,发现用余弦变换计 算的重力异常导数具有很高的精度,它与理论导数 的拟合程度是Fourier变换所不能比拟的,因此余弦 变换完全可以替代Fourier变换计算异常的导数.该 方法的提出,对进一步提高物探数据处理的精度有 着积极的意义. 参考文献(References) [1]罗孝宽,郭绍雍.应用地球物理教程.北京:地质出版社, 1991.139—154 LuoXK.GuoSY+AppliedGeophysicsCourse(inChinese). Beijing:GeologicalPublishingHouse,1991.139,154 [2]王谦身,安玉林,张赤军等.重力学.北京:地震出版社,2003. 196—290 WangQS,AnYL,ZhangCJ,eta1.Gravitology(inChinese). Beijing:SeismologicalPress,2003.196—290 【3]刘保华,张维冈,盂恩.重力异常垂向一阶导数的一种简便 算法.青岛海洋大学,1995,25(2):233—238 LiuBH,ZhangWG,MengEn.Asimplemethodforcalculatingthe firstverticalderivative.fournalofOceanUniversityofQingdao(in Chineae),1995.25(2):233—238 [4]汪炳柱.用样条函数法求重力异常二阶垂向导数和向上延拓 计算.石油地球物理勘探,1996,31(3):415—422 WangBZ.Computingtheverticalsecondderivativeandupward continuationofgravityanomalybysplinefunctionmethod.Oil GeophysicalProspecting(inChineae).1996,31(3):415—422 【5]吴宣志,刘光海,薛光奇等.富里叶变换和位场谱分析方法及 其应用.北京;测绘出版社,1987 WuXZ,nuGH,XueGQ,eta1.AnalysisofMethodand ApplicationonFourierTransformandPotentialFieldSpectrum(in Chinese).Beijing:SurveyingandMappingPress,1987 [6JKevinLMichus,JuanHomeroHinojosa.Thecompletegravity gradienttensorderivedfromtheverticalcomponentofgravity:A Fouriertransformtechnique.Journa/ofAppliedGeophysics,2001, 46:159一l74 l7]LourencoJS,MorrisonHF,Vectormagneticanomaliesderivedfrom measurementsofsinglecomponentofthefield.Geophysics,1973, 38:395—368 [8JAhmedNT,NatarajanT,RaoKR.Discretecosinetransform. IEEETrans.Comput.,1974.23(1):90—93 [9]RaoKFt,YipP.DiscreteCosineTransfrom:Algorithms, AdvantagesandApplications.NewYork:AcademicPress,1990 [10JDinsteini,RoseK,HeimanA.Variableblock—sizetransformimage coder.IEEETans.Comm.,1990,38(11):2073—2078 [11]胡广书.数字信号处理.北京:清华大学出版社,2003.363— 375 HuGS.DigitalSigr~Processing(inChinese).Beijing:TsinghHa UniversityPress.2003.363—375 [12]Cvetkovicz,PopovieMV.Newfast~orithmsforthecomputation disctetecosineandsinetransform.IEEE,卯Process., 1992,40(8):2oB3—2086 [13JVladimlrBritafidk.Aunifieddiscretecosineanddiseretesine transformcomputation.SignalProcessing,1995,43:333—339 I14JJainAK.AfastKarhunen—Loevetransformforaclassofstochastic processes.1EEETrans.Commun..1976,7.4:1023—1029 (本文编辑汪海英) l;;l}?}ll}}l2ll}};}f}}}}}}il《}l}}}}ll}}}}}}}}}}}l,《 ?如加如加mo加 一_l