下载文档原格式
实验二频率域位场处理和转换实验******专业：勘查技术与工程学号：**********指导老师：王万银，纪新林目录一、基本原理 (1)1.1空间域延拓原理 (1)1.2波数域延拓原理 (1)二、输入/输出数据格式设计 (2)1、输入数据 (2)2、输出数据 (2)三、总体设计 (2)1、程序流程图 (2)2、参数说明 (3)四、测试结果 (3)五、结论及建议 (7)附录 (8)一、基本原理1.1空间域延拓原理延拓分为：向上延拓和向下延拓。
向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。
对于二度体（令z 向下为正）：U(x,z)=-zπ U (ξ,0)(ξ−x)2+z 2d ξ +∞−∞（z<0） （1）其中U （ξ，0）为剖面上各点的实测值。
Z 为延拓的高度，即转换后的平面和观测平面的距离，由于z 坐标向下为正，因此z<0。
空间域的延拓实际是积分的计算。
向下延拓是由实测磁场向磁源方向延拓。
其计算公式和（1）相同。
1.2波数域延拓原理设场源位于z=H 平面以下（H>0），则磁场为在z=H 平面以上对x 、y 、z 连续的函数，若z=0观测平面的磁场T （x ，y ，0）为已知，则由外部的狄利克莱问题：T （x ，y ，z ）=−z2π T （ξ，η，0）[(x −ξ)2+(y −η)2+z 2]3/2+∞−∞+∞−∞d ξd η （2）对其进行变量x ，y 进行傅里叶变换成S T u ，v ，z = T （x ，y ，z ）e −2πi(ux +vy )dxdy ∞−∞∞−∞ （3） 因此：S T u ，v ，z = S T u ，v ，0 e 2π(u2+v 2)1/2z（4）在波数域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。
二、输入/输出数据格式设计1、输入数据（1）控制变量名：’cmd.txt’（2）观测面高度之差：height=3.3m或height=-3.3m（3）低高度网格化数据：从文件“A20_mag.grd”输入（4）高高度网格化数据：从文件“A53_mag.grd”输入（5）特征值：eigval=1.701411∗1038（6）输出文件名：out.grd2、输出数据通过修改参数，运行两次程序，向上延拓结果及向下延拓结果均输入“output.grd”先做向上延拓，将结果利用surfer画图后，再做向下延拓三、总体设计1、程序流程图图3.1 程序流程框图2、参数说明mpoint,line 点数，线数height 观测面高度差xmin,xmax,ymax,ymin x方向最小值、最大值，y方向最大值、最小值eigval 特征值m0,m1,m2,m3 x方向扩边点位n0,n1,n2,n3 y方向扩边点位Field 场值Field_real,Field_image 傅里叶变换中场值的实部，虚部Fconti_real Fconti_image 延拓因子实部、虚部四、测试结果1、原场值图：图4.1 低高度场值图图4.2 高高度场值2图4.3 向上延拓结果图4.4 向下延拓结果分析：将底高度向上延拓3.3m 后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了，将高高度向下延拓3.3m 后，异常部位没变，但是场值绝对值变大了近5倍，且在边缘部位有几个异常点。
基于奇异谱分析的重磁位场分离方法朱丹;刘天佑;李宏伟【摘要】奇异谱分析是一种近年兴起的时间序列分析方法,它利用降秩原理实现信号分离.该方法将数据空间投影到不同特征的子空间中,并用奇异值来表征这些子空间的性质,最后通过截取奇异值实现数据的重构.重磁位场分离可以看成一种多信号叠加的分离问题.不同特征的重磁异常具有不同特征的奇异谱,这是奇异谱分析用于解决位场分离问题的应用基础.本文通过建立理论模型,分析重磁异常的奇异谱特征,得出适用于重磁位场分离的最优参数选择方法,并与传统方法进行比较.对比发现,无论是横向叠加模型、垂向叠加模型还是斜向叠加模型,奇异谱分析都具有很好的分离效果.最后,将奇异谱分析用于鄂东南某矿区的重力资料处理中,实现弱异常的识别和分离.【期刊名称】《地球物理学报》【年(卷),期】2018(061)009【总页数】12页(P3800-3811)【关键词】奇异谱分析;重磁位场分离;降秩理论;最优参数;鄂东南地区【作者】朱丹;刘天佑;李宏伟【作者单位】中国地质大学(武汉)地球物理与空间信息学院,武汉430074;中国地质大学(武汉)地球物理与空间信息学院,武汉430074;中国地质大学(武汉)数学与物理学院,武汉430074【正文语种】中文【中图分类】P6310 引言奇异谱分析(Singular Spectrum Analysis, SSA)是一种近年兴起的时间序列分析方法,最早由Broomhead和King(1986)提出.自提出以来，SSA分析被广泛应用于多领域的信号处理中.SSA分析是信号去噪和预测的一种方法.该方法从Karhumen-Loeve分解理论的基础上发展而来(Vautard and Ghil，1989；Vautard et al.，1992).SSA分析是将原信号变换成Hankel矩阵，再对Hankel矩阵进行分解和重构，从而实现信号和噪声的分离.SSA分析能够实现信号和噪声分离的依据是它们的具有不同特征的奇异谱.Read等(1993)在SSA的基础上提出了用于多道信号处理的MSSA方法，Oropeza和Sacchi(2011)提出基于随机奇异值分解(Randomized Singular Value Decomposition，RSVD)的MSSA方法，Huang等(2016)提出提升去噪效果的阻尼MSSA方法.在地球物理领域，Sacchi(2009)、Oropeza和Sacchi(2011)、Kreimer和Sacchi(2012)、Chiu(2013)、Gan等(2015)、Huang 等(2016)利用SSA分析对一维和多维地震信号进行去噪和重建.重磁场是由具有密度与磁性差异的不同规模、不同深度、不同形状地质体共同引起，即为不同尺度、不同幅值异常的叠加.多种异常的混叠，给目标地质体的反演和解释带来困难.如何从混叠重磁场中分离出目标地质体引起的异常，是重磁勘探的研究方向之一.早期人们采用滑动平均和多项式拟合的方法实现不同尺度重磁异常的分离.在后来的重磁信号处理中，人们用频率域的概念描述重磁场.通过Fourier变换将重磁场由空间域变换到频率域，用振幅和相位来描述重磁场的特征.通常情况，浅部地质体产生的重磁场频率高，深部地质体产生的重磁场频率低.Spector和Grant(1970)推导了航磁异常的能谱公式，提出利用能谱分析粗略估计块状体埋深的方法，并在1975年提出匹配滤波法分离重磁场.文百红和程方道(1990)提出插值切割法分离磁异常.Mickus等(1991)首次将最小曲率法应用于美国某地区的布格重力异常分离中.随着20世纪70年代小波分析方法的提出，Fedu和Quarta(1998)提出一种基于离散小波变换的重磁位场分离方法.这些位场分离方法存在不同的问题，如窗口大小、拟合阶数、小波基的选择、区域场与局部场能谱斜率位置的确定等.本文应用奇异谱分析分离重磁位场，通过建立理论地球物理模型，分析重磁异常的奇异谱特征.通过理论模型计算分析，总结最优参数选取方法，提出适用于位场多信号分离的分解和重构方法，并对比SSA分析和传统场源分离方法的应用效果，最后将该方法应用到鄂东南某矿区赤铁矿重力异常的识别提取中.1 SSA和MSSA的基本原理实际工作中，剖面数据是一维的，平面数据经过网格化后可以表示成二维矩阵的形式.剖面和平面重磁数据分别采用SSA和MSSA进行分析.1.1 SSA设有剖面重磁数据x=[x1,x2,…,xN]，将该重磁数据以窗口M重新排列如下：(1)矩阵X称为时滞矩阵(Vautard and Ghil, 1989)，其自协方差矩阵表示如下：(2)Tx是Toeplitz矩阵.c是序列x的协方差：(3)计算Tx的特征值向量λ和对应的特征向量矩阵E，特征值向量表示如下：λ=[λ1,λ2,…,λM]，|λ1|≥|λ2|≥…≥|λM|,(4)特征向量矩阵表示如下：(5)其中E1,E2,…,EM是列向量，分别表示矩阵Tx的M个特征向量，并与λ1,λ2,…,λM对应.同时，矩阵Tx的奇异值向量σ可以表示如下：σ =[σ1,σ2,…,σM](6)通常，将奇异值向量σ称为x的奇异谱，其中明显大于0的奇异值称为有效奇异值.在实际情况中，信号奇异谱σ中的元素均大于0，其中较小的奇异值十分接近于0，而有效奇异值可以看成是非接近于0的奇异值.通过矩阵X和特征向量矩阵E求得权矩阵A(王解先等，2013)：A =[A1,A2,…,AM]=X′·E(7)权矩阵A反映了特征向量矩阵E在矩阵X中权重，那么通过权矩阵A和特征向量矩阵E可以重构出一维重磁数据x的M个不同尺度的细节细节可以用(8)式表示：(8)(8)式中分别为大小为(1×N)的向量，其中可由以下求得：(9)不同尺度的细节对应于奇异值σi，奇异值σi越大，与一维重磁数据x的近似程度越高，并且一维重磁数据那么一维重磁数据x的近似可以用前K个的和表示：(10)1.2 MSSA设有大小为Nr×Nc平面重磁数据如下：(11)通过F建立窗口大小Lr×Lc的轨迹矩阵Fi,j(Read, 1993; Oropeza and Sacchi, 2011; Huang et al., 2016)：(12)令Kr=Nr-Lr+1和Kc=Nc-Lc+1，那么1≤i≤Kr，1≤j≤Kc.构建Hankel矩阵Hi，Hi表示如下：(13)然后，结合Hankel矩阵Hi，构建块Hankel矩阵W，矩阵W表示如下：(14)矩阵W可以表示成如下形式：W=UΣV′,(15)其中，U和V是正交矩阵，Σ是对角矩阵.为了方便计算，令W是N×N的方矩阵，那么U和V是大小为N×N的正交矩阵，对角矩阵Σ可以表示成以下形式：(16)其中，|Σ1|≥|Σ2|≥…≥|ΣN|.那么矩阵W的奇异值向量σ表示如下：(17)令ΣK为对角矩阵Σ的截断位置，那么矩阵W的奇异值分解可以表示如下：(18)其中，ΣS=diag(Σ1,Σ2,…,ΣK)，ΣM=diag(ΣK+1,ΣK+2,…,ΣN)，US是U的第1到K列，UM是U的第K+1到N列， VS是V的第1到K列，VM是V的第K+1到N列.那么，块矩阵W的近似可以表示如下：(19)最后，通过对近似矩阵采取平均逆投影变换(Golyandina et al., 2007)就能够得到平面重磁数据F的近似1.3 实现步骤通常情况，信号的时滞矩阵是低秩的，而噪声会增加矩阵的秩.分解的过程可以看成时滞矩阵向低维子空间的投影，这些子空间的性质可以用奇异值来表征.重构的过程则是对奇异值分类，并选取特定奇异值所对应的子空间进行逆投影.分解和重构的过程实际上是矩阵降秩的过程.SSA和MSSA的实现步骤可以归纳为以下三点：(1) 原始数据的重新排列.利用一维窗口将一维重磁数据x重新排列成时滞矩阵X，X是Hankel矩阵.二维重磁数据F则采用二维窗口重构块Hankel矩阵W.(2) 矩阵的奇异值分解.对时滞矩阵X的自协方差矩阵或块Hankel矩阵W进行奇异值分解，得到奇异值和对应的特征向量.(3)重磁数据的重构.在合适位置将奇异值截断，并用截断的奇异值重构位场数据.2 重磁场的奇异谱特征下面我们分析不同尺度(不同埋深)、不同幅值重磁异常及噪声的奇异谱特征.(1) 幅值相同、尺度不同(埋深不同)的重磁异常的奇异谱特征建立四个埋深不同磁化强度不同的模型，由四个模型分别产生四组幅值相同、尺度不同的磁异常如图1所示，磁异常点数是800.分别对这四组磁异常做SSA分析，SSA分析窗口是400(采用较大窗口可以得到更多奇异值，其谱特征反映的更细致，窗口400是最大窗口，以下其他3种情况分析的点数窗口大小相同).从不同尺度磁异常的奇异谱中可以看出，随着磁异常尺度的增大，奇异值下降速度变快，有效奇异值的个数增多.这表明尺度大(埋深大)的重磁异常其奇异谱曲线陡峭，主要集中在奇异值较大的部分，尺度小(埋深小)的重磁异常其奇异谱曲线趋平缓，其特征与Fourier分析的功率谱相似.(2) 尺度相似(埋深相同)、幅值不同的重磁异常的奇异谱特征建立四个埋深相同磁化强度不同的模型，由四个模型分别产生四组尺度相似、幅值不同的磁异常如图2所示.从不同幅值磁异常的奇异谱中可以看出，随着磁异常幅值的增大，奇异值下降速度变快，但有效奇异值的个数几乎不变.说明相同埋深的地质体奇异谱具有相似的特征，时滞矩阵的秩几乎不变.(3) 尺度相同、幅值相同、水平位置不同的重磁异常的奇异谱特征对于同一模型，水平位置不同所产生的三组尺度相同、幅值相同的磁异常如图3所示.从图中可以看出，同一异常，分布在不同位置，其奇异谱相同.说明地质体的水平位置不影响奇异谱的分布.(4) 高斯白噪声的奇异谱特征若有三组不同强度的噪声如图4所示.从图中可以看出，三组不同强度噪声的奇异值在横坐标上都有分布.说明噪声的时滞矩阵是高秩的，噪声强度只决定奇异值的大小，而不会改变秩的大小.通过上面的理论模型计算分析可以知道，重磁场的奇异谱包含尺度(埋深)和幅值(物性)的信息，其中尺度决定奇异谱的下降速度和时滞矩阵的秩，幅值决定奇异谱的下降速度但不改变时滞矩阵的秩.对于重磁场，异常尺度越大，奇异谱下降速度越快，秩越低；异常幅值越高，奇异谱下降速度越快，但时滞矩阵的秩不变.噪声的奇异值在横坐标上都有分布，说明噪声的时滞矩阵是高秩的.异常的尺度(埋深)和幅值(物性)决定了奇异谱的特征，这是利用奇异谱进行场源分离的基础.3 参数讨论构建时滞矩阵的窗口与重构信号所选择的奇异值是SSA分析的两个参数，下面讨论在位场分离中如何选择这些参数.3.1 构建时滞矩阵的窗口选择时滞矩阵的自协方差矩阵Tx的大小由窗口决定，若窗口取M，矩阵Tx的奇异值个数也为M，这表示我们可以将总场分解成M个细节.因此，窗口的选择直接影响区域场和局部场的分离效果.在信号分析中，通常窗口取M/2(Golyandina et al., 2007; Chiu, 2013; Huang et al., 2016).通过理论模型计算分析发现，由于分离目标和信号特征的差异，这种窗口选择方式并不适用于解决重磁位场分离问题(图5c 给出了窗口与均方差的关系).图5是垂向叠加模型的奇异谱分析，剖面点数是240，点距50 m，化极磁异常零值点之间宽度约1200 m，最优窗口在25个点左右，当窗口选择10个点时，区域场和局部场没有完全分离，当窗口选择90个点时，区域场过拟合，局部场出现波动.图1 幅值相同、尺度不同模型的奇异谱特征Fig.1 The singular spectrum features of same amplitude and different scale models图2 幅值不同、尺度相同模型的奇异谱特征Fig.2 The singular spectrum features of different amplitude and same scale models图3 水平位置不同模型的奇异谱特征Fig.3 The singular spectrum features of different horizontal position models图4 噪声模型的奇异谱特征Fig.4 The singular spectrum features of noise models由图5可知，窗口选取过小会导致异常分离不完全，重构的区域场中仍包含过多局部场的信息，但是当窗口选取过大时，重构得到的区域场过拟合，导致分离得到的局部场产生波动.在重磁位场分离中，SSA分析的窗口取决于所要分离目标异常的尺度.当目标异常尺度较大时，应该选取较大的窗口进行SSA分析，当目标异常尺度较小时，则应选取较小的窗口进行SSA分析.通过理论模型计算发现，最优窗口的宽度近似等于待分离目标异常零值线的宽度.但是，窗口的改变对分离效果的影响是渐变的，因此最佳窗口可以不是某一确定的值，而是一个区间，窗口在这个区间内取值，得到的分离效果是相似的.图5 过小和过大窗口的理论模型奇异谱分析(a) 分离局部场； (b) 分离区域场； (c) 窗口对分离结果的影响.Fig.5 Singular spectrum analysis of synthetic model of undersize and oversize windows(a) The separated local fields; (b) The separated regional fields; (c) The effects of the windows on the results separated by SSA.3.2 重构阶数的选择很多学者利用SSA分析去噪时，采用均值截断法、二分法等方法重构信号(Oropeza and Sacchi, 2011; 王解先, 2013)，这些重构方法能否用于重磁位场分离？为了探究重磁位场分离的重构方法，建立理论叠加模型，该模型引起的磁异常如图6右上所示，其中总场是区域场、局部场和噪声的叠加.对理论模型的总场、区域场、局部场和噪声分别做奇异谱分析.奇异谱分析窗口是400(采用较大窗口可以得到更多奇异值，其谱特征反映的更细致)，截前60个奇异值如图6所示.为了阐述方便，这里令总场、理论区域场、理论局部场和噪声的奇异值分别为σi、、和，其中1≤i≤400.图6中，总场的奇异谱包含了理论区域场、局部场和噪声的奇异谱特性，其中总场奇异值σ1、σ2与理论区域场奇异值、相似，总场奇异值σ3,…,σ30与理论局部场奇异值,…,相似，总场奇异值σ31,…,σ400则与噪声的奇异值,…,相似.根据上述对应关系，将总场的奇异谱分解成[σ1,σ2]、[σ3,…,σ30]、[σ31,…,σ400]三个部分，这三个部分分别是区域场、局部场和噪声的奇异谱的主要成分.可以用它们去重构区域场、局部场和噪声(采用公式(9)).由上可知，总场的奇异谱包含其各成分的奇异谱特征，这为奇异值的分类提供参考.在实际问题中，我们可以根据奇异值的大小和下降趋势来进行划分(如图6中Regional field和Local field).一般来说，奇异值大且快速下降的部分反映的是区域场的特征，奇异值较小且缓慢下降的部分反映的是局部场的特征，奇异值近似为0且分布在整个横坐标的部分反映的是噪声的特征.4 奇异谱分析与传统方法的对比下面以不同的理论模型对比SSA分析和传统方法的分离效果，并从应用的角度讨论SSA分析在解决重磁位场分离问题中的可行性.图6 重构阶数的选择方法Fig.6 The setting of orders for reconstruction4.1 理论模型试验建立地球物理模型如图7b、图8b和图9b所示.图7b是横向叠加模型的磁异常场，共有800个数据点，点距是50 m，其中局部场由模型1(5000×10-3A·m-1)、模型2(2500×10-3A·m-1)和模型3(1500×10-3A·m-1)正演构成，区域场由模型4(2000×10-3A·m-1)正演构成.图8b是垂向叠加模型的磁异常场，总共有241个观测点，点距是50 m，其中局部场由模型5(3000×10-3A·m-1)正演构成，区域场由模型6(5000×10-3A·m-1)正演构成.图9b是斜向叠加模型的磁异常场，总共有241个观测点，点距是50 m，其中局部场由模型7(1000×10-3A·m-1)正演构成，区域场由模型8(4000×10-3A·m-1)正演构成.图7c和7d是窗口大小为400(以Model 4为分离目标)的横向叠加模型SSA分析结果.该模型的奇异谱(图7e)可以依据奇异值的大小和下降速度分成2个部分，第一个部分是σ1和σ2，剩余奇异值是第二个部分.这两个部分的奇异值由大到小，下降速度由快到慢.容易知道，奇异谱的第一个变化过程由区域场引起，那么采用奇异值σ1和σ2对应的细节重构信号就能够得到分离的区域场(图7c)，采用剩余奇异值对应的细节重构信号就能够得到分离的局部场(图7d).从图中可以看出，奇异谱分析在横向叠加模型上具有良好的分离效果.图7 横向叠加模型SSA分析Fig.7 Horizontal superposition model separation using SSA图8 垂向叠加模型SSA分析Fig.8 Vertical superposition model separation using SSA图9 斜向叠加模型SSA分析Fig.9 Diagonal superposition model separation using SSA图8c、8d是窗口大小为25(以Model 5为分离目标)的垂向模型SSA分析结果.该模型的奇异谱(图8e)能够分成两个部分，第一个部分是σ1，剩余奇异值是第二个部分，这两个部分分别对应地质体6和地质体5产生的异常.那么，可以将奇异值σ1对应的细节作为区域场(图8c)，将剩余奇异值对应的细节重构作为局部场(图8d).从分离结果可以看出，奇异谱分析能够实现垂向叠加模型的分离.图9c、9d是窗口大小为50(以Model 7为分离目标)的斜向模型SSA分析结果.该模型的奇异谱(图9e)能够分成两个部分，第一个部分是σ1、σ2和σ3，剩余奇异值是第二个部分，这两个部分分别对应地质体8和地质体7产生的异常.那么，可以将奇异值σ1、σ2和σ3对应的细节作为区域场(图9c)，将剩余奇异值对应的细节重构作为局部场(图8d).从分离结果可以看出，奇异谱分析能够实现斜向叠加模型的分离.4.2 传统场源分离方法的分离本小节主要利用小波分析法(Mallat and Hwang，1992；侯遵泽和杨文采，1997；杨文采等，2001)、匹配滤波法、插值切割法(文百红和程方道，1990)、多次迭代趋势分析法(李春芳，2011)、最小曲率法(纪晓琳等，2015；王万银等，2009)这五种传统重磁位场分离方法对4.1节中的模型进行分离.这些分离方法的参数都是在多次实验并与理论模型对比得到最好的分离效果下确定的.由于每种方法所针对的异常特点不同，每种模型采用效果最好的两种传统方法进行对比，分别见图10、图11和图12.图10 横向叠加模型传统方法实验Fig.10 Traditional methods separation of horizontal superposition model图11 垂向叠加模型传统方法实验Fig.11 Traditional methods separation of vertical superposition model图12 斜向叠加模型传统方法实验Fig.12 Traditional methods separation of diagonal superposition model4.3 方法对比重磁位场分离方法的优劣可以从两个方面判别，一是分离效果的好坏，二是参数选取的难易.重磁位场分离效果的好坏可以用理论场与分离场的均方差表示，均方差越小，表示分离场与理论场更接近.参数选取的难易程度又包括： 1 最优参数是否容易确定； 2 选取参数所带来的误差对分离结果的影响程度.通常，分离效果是判断方法好坏的先行标准，当方法的分离效果满足一定标准时，才有讨论参数选取难易程度的意义.表1 各方法分离理论模型的均方差Table 1 The MSE of synthetic models separated by different methods均方差(nT)横向叠加模型垂向叠加模型斜向叠加模型奇异谱分析法9.712.777.99小波分析法25.2610.309.69匹配滤波法25.9312.5913.69插值切割法21.455.6813.71迭代趋势分析法20.12--最小曲率法--7.82注：“-”表示该方法分离效果差，未列出.从表1可以看出，奇异谱分析能够用于三种不同模型的异常分离，并具有很好的分离效果.小波分析、匹配滤波和插值切割均能用于三种模型的异常分离，但均方误差大于奇异谱分析.迭代趋势分析对横向叠加模型具有较好的分离效果，但是无法分离垂向和斜向叠加模型.最小曲率法对局部异常尺度较小的斜向叠加模型具有较好的分离效果，但是无法分离局部异常尺度较大的横向和垂向叠加模型.在参数选取的方面，小波分析存在如何根据哪一阶细节重构区域场与局部场的问题；匹配滤波存在功率谱估计场源深度误差较大的问题；插值切割需要设置插值半径和迭代次数，这两个参数对结果都有较大影响；多次迭代趋势分析参数只需设置迭代次数；最小曲率法需要设置迭代步长和迭代次数，两个参数对结果有一定影响；SSA分析需要先确定窗口大小，然后根据奇异谱分布确定重构阶数.奇异谱分析具有适用性强，分离效果好的优点，并通过奇异谱的变化，容易选择重构阶数.但是，在二维数据的奇异谱分析中，其时滞矩阵大小是NrNc(Nr-Lr)(Nc-Lc)×NrNc(Nr-Lr)(Nc-Lc)，因此在Nr和Nc较大的情况下，需要较多的计算时间.5 奇异谱分析在鄂东南某矿区的应用鄂东南矿集区以接触交代型铁矿床为主，它们主要形成于燕山期中酸性侵入岩与碳酸盐岩接触带附近.区内出露白垩纪砂岩粉砂岩和燕山期闪长岩、石英斑岩(图13)，地层北倾.区内铁矿以赤铁矿为主，伴有磁铁矿，其中赤铁矿相对围岩的剩余密度为0.55 g·cm-3，磁铁矿相对围岩的剩余密度为1.47 g·cm-3，高密度的赤铁矿和磁铁矿能够引起局部高重异常.从布格重力异常图(图14)中可以看出，布格重力异常值在-28 mGal到-22 mGal之间，由北东高过渡到南西低.矿体产生的重力异常相对较弱，难以直接从布格重力异常图中分辨.图13 研究区地质图Fig.13 Geological map of study area图14 研究区布格重力异常图Fig.14 The Bouguer gravity anomaly map of study area图15 布格重力异常奇异谱分布图Fig.15 The singular spectrum of Bouguer gravity anomaly由于中浅部地质体引起的重力异常尺度较小，因此采用10×10(点线距20 m×20 m)的窗口进行多道奇异谱分析，奇异谱如图15所示.根据奇异谱的特征分类，用奇异值σ1重构区域场(图16a)，奇异值σ2和σ3重构局部场(图16b)，其余的奇异值则视为噪声引起.区内岩体和沉积岩的接触带是成矿的有利位置，分离的局部场在接触带附近存在多个局部高重异常.西缘的局部高重异常位于闪长岩和沉积岩的接触带上，形状近圆形，极大值是0.08 mGal.东缘的局部高重异常位于石英斑岩、闪长岩和沉积岩的接触带上，该接触带具有成矿的地质条件.东缘局部高重异常形状近条带状，规模比西缘局部高重异常大，极大值是0.16 mGal.由于区内岩体和沉积岩密度低，因此推断这两个局部高重异常由中浅部的高密度地质体引起，该地质体可能是赤铁矿、磁铁矿或蚀变岩体.2014年施工的ZKI01孔位于西缘局部高重异常的中心位置，钻井见矿117 m厚，其中59 m到159 m处为赤铁矿，159 m到177 m处为磁铁矿.结合地质和钻井资料，圈定接触带东缘局部重力大于0.03 mGal的区域作为成矿远景区(见图13).图16 鄂东南某工区多道奇异谱分析结果(a) 区域场； (b) 局部场.Fig.16 The results processed by MSSA of work area in southeastern Hubei(a) Regional field; (b) Local field.6 结论奇异谱分析是一种近年兴起的时间序列分析方法，它利用降秩原理实现信号分离，其过程包括变换、分解和重构.分解的过程可以看成时滞矩阵向低维子空间的投影，这些子空间的性质可以用奇异值来表征.重构的过程则是对奇异值分类，并选取特定奇异值所对应的子空间进行逆投影.结合理论模型计算分析，我们发现区域场具有低秩高奇异值的特征，局部场具有低秩中低奇异值的特征，噪声具有高秩低奇异值的特征，这是利用降秩原理进行场源分离和去噪的基础.通过理论模型试验我们发现，最优窗口的宽度近似等于待分离目标异常零值线的宽度，重构阶数则根据奇异谱的特征来选择.结合理论模型和实际数据，该方法能够实现重磁位场的分离，并具有适用性强、分离效果好的优点，但相比于传统方法，其计算效率稍低.致谢感谢《地球物理学报》编辑部对本文的支持，以及审稿专家的宝贵意见. ReferencesBroomhead D S, King G P. 1896. Extracting qualitative dynamics from experimental data. Physica D: Nonlinear Phenomena, 20(2-3): 217-236. Chiu S K. 2013. Coherent and random noise attenuation via multichannel singular spectrum analysis in the randomized domain. Geophysical Prospecting, 61(1): 1-9.Fedu M, Quarta T. 1998. Wavelet analysis for the regional-residual andlocal separation of potential field anomalies. Geophysical Prospecting,46(5): 507-525.Gan S W, Chen Y K, Zi S H, et al. 2015. Structure-oriented singular value decomposition for random noise attenuation of seismic data. Journal of Geophysics and Engineering, 12(2): 262-272.Gao J J, Sacchi M, Chen X H. 2013. A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions. Geophysics, 78(1): V21-V30.Golyandina N E, Usevich K D, Florinsky I V. 2007. Filtering of digital terrain models by two-dimensional singular spectrum analysis. International Journal of Ecology & Development, 8(7): 81-94.Hou Z Z, Yang W C. 1997. Wavelet transform and multi-scale analysis on gravity anomalies of China. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 40(1): 85-95.Huang W L, Wang R Q, Chen Y K, et al. 2016. Damped multichannel singular spectrum analysis for 3D random noise attenuation. Geophysics, 81(4): V261-V270.Ji X L, Wang W Y, Qiu Z Y. 2015. The research to the minimum curvature technique for potential field data separation. Chinese Journal of Geophysics (in Chinese), 58(3): 1042-1058, doi: 10.6038/cjg20150329. Kreimer N, Sacchi M D. 2012. A tensor higher-order singular value decomposition for prestack seismic data noise reduction and interpolation. Geophysics, 77(3): V113-V122.Li C F. 2011. The research on the separating methods for potential field in。
航磁数据位场转换处理及效果∆测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测航磁T面上的综合反映。
由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。
为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。
第一节位场转换处理及效果航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。
频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。
位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的WINDOWS系统下地球物理数据处理解释软件（GeoProbe Mager）及航空物探彩色矢量成图系统（AgsMGis）。
一、原平面化极处理化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。
测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。
这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。
化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。
从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。
作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。
基于频率域磁异常转换的断裂构造划分判据宁堃;谢涛;邵昌盛;邓子清【摘要】实践表明,利用磁异常可以较为准确的对断裂、区域深大断裂、褶皱等构造进行推断划分.实践中,通过对西藏某热液型铅锌多金属矿区实测的磁异常数据进行频率域垂向二阶导数和水平梯度模的转换处理,突出一些浅部地质体或边界引起的弱磁异常响应,结合实际地质情况和?T、垂向二阶导数、水平梯度模的断裂确定标志,推断划分了多条断裂构造,取得了较理想的效果.【期刊名称】《四川地质学报》【年(卷),期】2018(038)002【总页数】5页(P332-336)【关键词】磁异常;频率域;垂向二阶导数;水平梯度模;断裂构造【作者】宁堃;谢涛;邵昌盛;邓子清【作者单位】四川省核工业地质局二八二大队,四川德阳 618000;四川省核工业地质局二八二大队,四川德阳 618000;四川省核工业地质局二八二大队,四川德阳618000;四川省核工业地质局二八二大队,四川德阳 618000【正文语种】中文【中图分类】P631地球物理勘探中，用磁法可以圈定断裂破碎带，是因为断裂的产生改变了岩石的磁性，或者改变了地层的产状，或者沿断裂带伴有后期或同期岩浆活动，或者沿断裂两侧具有不同的构造特点。
一般情况下，在DT等值图上可以根据磁异常特点直接对断裂进行划分。
然而，在地质情况复杂的地区，迭加异常相互干扰，弱信号与误差并存，使异常变得十分复杂而难以解释，这就需要根据实际情况对磁异常进行精细化处理和转换，从而进一步提高解释推断的精度。
以西藏某热液型铅锌多金属矿的高精度磁法数据为例，介绍如何综合利用频率域磁异常垂向二阶导数和水平梯度模技术对断裂构造进行推断划分。
频率域磁异常转换首先要将磁异常原始数据经过傅氏变换转换为频谱，然后乘以相应的频响因子，再经过反傅氏变换得到结果，这种处理方法对磁场高频成分有突出和放大作用，它侧重于浅层近地表地质的磁效应而压制深层区域背景场的影响，从而突出浅部地质体引起的局部异常和解释地质体的边界（断裂、岩体、矿体边界等）。
《重磁资料处理与解释》实验二频率域位场处理和转换实验学院：地测学院专业名称：勘查技术与工程学生姓名：学生学号：指导老师：提交日期：2018年1月9日二0一八年一月目录1 基本原理 (2)1.1位场的方程 (2)1.2二维傅里叶变换及卷积性质 (2)（1）傅里叶变换 (2)（2）卷积性质 (2)1.3频率域位场延拓原理 (3)2 输入/输出数据格式设计 (3)2.1 输入数据格式设计 (3)2.2 输出数据格式设计 (3)2.3 参数文件数据格式设计 (3)3 总体设计 (4)3.1频率域位场处理与转换的一般步骤 (4)3.2软件总体设计结果流程图 (4)4 测试结果 (5)4.1 测试参数 (5)（1）向上延拓 (5)（2）向下延拓 (5)4.2 测试结果 (6)5 结论及建议 (7)附录：源程序代码 (8)1 基本原理1.1位场的方程由场论知识可知，位场方程分为两大类：有源的Possion 方程()02≠∇U ，以及无源的Laplace 方程()02=∇U 。
Laplace 方程的第一边值问题()1|f U S =通常为Dirichlet 问题，第二边值问题⎪⎭⎫⎝⎛=∂∂2|f n U s 通常称为Nueman 问题。
若P 点在S 平面内称为内部问题，反之称为外部问题。
由唯一性定理可知，Dirichlet 的内部和外部问题的解是唯一的，而Nueman 内部问题的解不是唯一的，有一常数差，但其外部问题解是唯一的。
外部问题的解的唯一性的原因：。
0;0=∂∂=∞→∞→r r nUU 无源区域位场可以表示为：ds n G W n W G p W ⎰⎥⎦⎤⎢⎣⎡∂∂-∂∂=π41)( （1-1）()()()()()[]()()z y x h W d d z y x W z z y W -=-+-+--=⎰⎰+∞∞-+∞∞-ξξηεηεξηεξηεπξ,,*,,,,2,,x 23222（1-2）1.2二维傅里叶变换及卷积性质（1）傅里叶变换[]⎰⎰+∞∞-+∞∞-+-==dxdy y x g y x g F v u G e vy ux i )(2),(),(),(π （1-3）[]⎰⎰+∞∞-+∞∞-+-==dudv v u G v u G y x g eF vy ux i )(21),(),(),(π （1-4）（2）卷积性质()()[]()()v u P v u G y x p y x F ,*,,*,g = （1-5）()()[]()()y x p y x v u P v u G F ,*,g ,*,1=- （1-6）1.3频率域位场延拓原理当已知实测平面的异常时，换算场源以外的异常称之为延拓，分为向上延拓和向下延拓。
半空间狄利克莱问题解析解：()[]()[]()[]()[]()z v u ey x W F z y x h F y x W F z y x W F -+-=-=ζπζζζ222,,,,,,,, （1-7）其中：ez v u )(222-+-ζπ称为延拓因子，ζ为计算面Z 坐标 ,Z 轴向下为正方向，()[]ζ,,y x W F 为计算面频率域位场，()[]z y x W F ,,为延拓面的频率域位场。
2 输入/输出数据格式设计2.1 输入数据格式设计观测面位场数据保存在filename_obser 文件中，为.grd 格式。
计算延拓误差时的精确场值文件保存在filename_obser2中，为.grd 格式。
2.2 输出数据格式设计实际计算得到的数据保存在filename_output1和文件filename_output2中，为.grd 格式。
2.3 参数文件数据格式设计将所要读取的参数保存在一个文件中，该文件名变量为cmdfile ，字符串变量，长度不超过80，全路径名。
在该文件中保存的参数如下：filename_obser:低高度观测面位场数据文件 filename_output1：向上延拓后位场数据文件 filename_output2：向下延拓后位场数据文件 filename_obser2:高高度观测面数据文件factor_m ：扩边因子distance1:向上延拓的高度（m/z 轴向下为正方向） distance2:向下延拓的高度（m/z 轴向下为正方向）3 总体设计3.1频率域位场处理与转换的一般步骤()[][][][][]即可。
去掉扩边部分，输出步第；进行反变换得到步第；计算步第；以及方向转换因子、延拓因子计算导数因子步第；对扩边后的数据进行步第进行扩边处理；对网格化的数据步第),,(:6),,(),,(:5),,(),,(F :4),,,,( )Y(),,,,D(:3,,W FFT :2),,(:1121z y x Q z y x Q F F z y x Q f Y D y x W F z y x Q M M t t u,v f z u,v,t t t u,v y x F y x W n -=⋅⋅⋅=''-⇒ζζζζ3.2软件总体设计结果流程图此次程序采用IPO 结构设计，首先通过读取cmd 文件，得到相关输入参数：观测面位场数据文件名变量、延拓后位场数据文件名变量、延拓后准确位场数据文件名变量、扩边因子、延拓的高度（m/z 轴向下为正方向）；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。
下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来计算延拓因子并且将扩完边的数据进行快速二维傅里叶变换后在频率域与延拓因子相乘；最后进行快速二维傅里叶反变换并且去除扩边部分后输出。
总体设计见表1。
4 测试结果4.1 测试参数（1）向上延拓原始场值数据保存在’’a20_mag.grd”中，向上延拓3.3m，延拓后理论数据保存在“a53_mag.grd”中,延拓后的数据保存在output1.grd中。
网格大小为27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.（2）向下延拓原始场值数据保存在’’a53_mag.grd”中，向下延拓3.3m，延拓后理论数据保存在“a20_mag.grd”中延拓后的数据保存在output2.grd中。
网格大小为27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有关参数保存在filename.cmd文件中，如下:A20_mag.grdA53_mag.grdoutput1.grdoutput2.grd3.3 -3.34.2 测试结果图1 低高度观测面(a20_mag)位场等值线图图2 高高度观测面(a53_mag)位场等值线图图3 高高度观测面(a53_mag)向下延拓位场等值线图图4 低高度观测面(a20_mag)向上延拓位场等值线图图5 运行结果图5 结论及建议由测试结果可知，向上延拓可以对浅部异常进行压制，且误差较小；向下延拓可突出浅部异常，且延拓结果误差较大。
通过本次频率域位场处理与转换设计实验，我收获颇丰，同时也感触很深！首先，刚开始我对空间域转频率域以及位场延拓的概念与处理是有不理解的，但是经过老师耐心细致的讲解和同学们之间的帮助，最终克服了这些困难。
其次，加深了我对于位场延拓的作用和具体延拓步骤的理解。
感谢老师孜孜不倦的教导与同学不厌其烦的讲解。
谢谢！附录：源程序代码!******************************************************************** **********! 程序功能：实现频率域位场延拓! cmd文件参数:! cmdfile:存放有关参数的文件名变量! filename_obser:低高度观测面位场数据文件! filename_output1：向上延拓后位场数据文件! filename_output2：向下延拓后位场数据文件! filename_obser2:高高度观测面数据文件! factor_m：扩边因子! distance1:向上延拓的高度（m/z轴向下为正方向）! distance2:向下延拓的高度（m/z轴向下为正方向）! .grd文件参数：! N_point,N_line：点数(x方向)、线数(y方向)! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! Ur:初始观测面场值! 扩边参数：! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! 延拓参数：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! Fr：延拓观测面的信号的实部! Fi：延拓观测面的信号的虚部! U：延拓后准确场值! W(m,n)：径向圆频率! Y(m,n)：延拓因子! 误差参数：! error:延拓后场值与真实场值的均方误差!******************************************************************** **********Program plyytimplicit nonereal eigvalparameter(eigval=3.701411*1e5)character*(80)cmdfile,filename_obser,filename_output1,filename_output2,filename _obser2real,allocatable::Ur(:,:),Ui(:,:),y(:,:),Fr(:,:),Fi(:,:),U(:,:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,n1,n2integer factor_mreal xmin,xmax,ymin,ymax,dx,dyreal distance1,distance2,errorcmdfile='filename.cmd'callread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,& filename_obser2,distance1,distance2)call read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n),U(1:m,1:n),W(1:m ,1:n))call input_grd(Ur,filename_obser,m1,m2,n1,n2,m,n)call input_grd(U,filename_obser2,m1,m2,n1,n2,m,n)!低高度向上延拓call expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)call FFT2(Ur,Ui,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,W)call calculate_Y(m,n,W,y,distance1)call calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)call FFT2(Fr,Fi,m,n,1)calloutput_grd(filename_output1,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax ,Ymin,Ymax)call calculate_error(distance1,Fr,U,m1,m2,n1,n2,m,n,error)!高高度向下延拓call expand_cos_2D(m1,m2,m,n1,n2,n,U,Ui)call FFT2(U,Ui,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,W)call calculate_Y(m,n,W,y,distance2)call calculate_FmulY(m,n,U,Ui,y,Fr,Fi)call FFT2(Fr,Fi,m,n,1)calloutput_grd(filename_output2,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax ,Ymin,Ymax)call calculate_error(distance2,Fr,Ur,m1,m2,n1,n2,m,n,error)deallocate(Ur,Ui,y,Fr,Fi,u,w)end program!******************************************************************** ********! 子程序：read_cmd! 功能：读取参数文件! 输入参数说明：! cmdfile：参数文件名! 输出参数说明：! filename_obser：存放低高度观测面的数据文件! filename_output1：存放向上延拓观测面的数据文件! filename_output2：存放向下延拓观测面的数据文件! filename_obser2：存放高高度观测面的数据文件! distance1:向上延拓的高度（m/z轴向下为正方向）! distance2:向下延拓的高度（m/z轴向下为正方向）! factor_m：扩边因子!******************************************************************** ********subroutineread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,& filename_obser2,distance1,distance2)implicit nonecharacter*(*)cmdfile,filename_obser,filename_output1,filename_output2,filename_ obser2integer factor_mreal distance1,distance2character*80 stropen(10,file=cmdfile,status='old')read(10,*)filename_obserread(10,*)filename_obser2read(10,*)filename_output1read(10,*)filename_output2read(10,*)factor_mread(10,*)distance1read(10,*)distance2close(10)end subroutine read_cmd!******************************************************************** *******! 子程序：read_grd! 功能：从原始观测.grd文件中读取相关参数! 输入参数说明：! filename_obser：输入文件名! 输出参数说明：! N_point,N_line：点数、线数! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值!******************************************************************** *******subroutine read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax) implicit nonecharacter*(*)filename_obserinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=filename_obser,status='old')Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)end subroutine read_grd!******************************************************************** ******! 子程序：calculate_mn! 功能：确定扩边数据点号位置! 输入参数说明：! factor_m: 扩边比例因子（>1.0）! a,b：点数、线数! 输出参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)!******************************************************************** ******subroutine calculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicit noneinteger a,b,m,n,m1,m2,n1,n2integer mtemp,mu,nuinteger factor_mmtemp=aDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENm=a*2ELSEmu=int(log(float(a))/0.693147+factor_m)m=2**muEND IFm1=1+(m-a)/2m2=m1+a-1pausemtemp=bDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENn=b*2ELSEnu=int(log(float(b))/0.693147+factor_m)n=2**nuEND IFn1=1+(n-b)/2n2=n1+b-1pauseend subroutine calculate_mn!******************************************************************** *****! 子程序：INPUT_GRD! 功能：读取grd文件中的数据! 输入参数说明：! filename_obser：输入文件名! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! 输出参数说明：! A：存放输出数据的二维数组名!******************************************************************** *****SUBROUTINE INPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)character*(*)filename_obserinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,kdo j=1,n,1do i=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=filename_obser,status='old')read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*) ((A(i,j),i=m1,m2),j=n1,n2)Close(20)END SUBROUTINE INPUT_GRD!******************************************************************** *****! 子程序：expand_cos_2D! 功能：二维扩边子程序并为信号虚部赋值! 输入参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! 输出参数说明：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部!******************************************************************** *****subroutine expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicit noneinteger m,n,m1,m2,n1,n2real Ur(1:m,1:n),Ui(1:m,1:n)real,allocatable::u(:),r(:)integer j,i,kallocate(u(1:m))do j=n1,n2,1do i=1,m,1u(i)=Ur(i,j)enddocall expand_cos_1d(1,m1,m2,m,u(1)) do i=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n))do i=1,m,1do j=1,n,1r(j)=Ur(i,j)enddocall expand_cos_1d(1,n1,n2,n,r(1)) do j=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)do i=1,mdo j=1,nUi(i,j)=0enddoenddoend subroutine expand_cos_2D!******************************************************************** ******! 子程序：expand_cos_1d! 功能：一维扩边子程序! 输入参数说明：! n0,n3：扩边后数据起点位置和终点位置! n1,n2：实际数据起点位置和终点位置! feild(i),(i=n1,n1+1,...,n2)：实际数据! 输出参数说明：! field(i),(i=n0,...,n1-1)：扩边后左边的数据! field(i),(i=n2+1,...,n3)：扩边后右边的数据!******************************************************************** ******Subroutine expand_cos_1d(n0,n1,n2,n3,Field)Real Field(n0:n3)pi=3.141592654Field (n0)=(Field (n1)+Field (n2))/2.0Field (n3)=Field (n0)i1=n0i2=n1DO i=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doi1=n2i2=n3DO i=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doEnd subroutine expand_cos_1d!******************************************************************** ****! 功能：FFT2! 功能：复数组2-D快速Fourier变换! 输入参数说明：! m0,m3：x方向的起点和终点! n0,n3：y方向的起点和终点! field：输入信号(需要赋值给Freal，实部)! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Freal：信号的实部! Fimage：信号的虚部(对于实信号而言，赋值为0)! 对应频率分布说明：! 数据Freal(m,n)和Fimage(m,n)对应的频率分布位置为:! m 方向：0,1,.......,m/2-1,m/2,-(m/2-1),......,-1! n 方向：0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!******************************************************************** ****SUBROUTINE FFT2(Freal,Fimage,m,n,nf)implicit noneINTEGER m,n,nfREAL Freal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE::Treal(:),Timage(:)integer nmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOPDO i=1,m,1IF (n.ne.1) THENdo j=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)end docall FFT(Treal,Timage,n,nf)do j=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)end doEND IFEND DODO j=1,n,1IF(m.ne.1) THENdo i=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)end docall FFT(Treal,Timage,m,nf)do i=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)end doEND IFEND DOdeallocate(Treal,Timage,STAT=ierr)END SUBROUTINE FFT2!****************************************************************** ! 子程序：FFT! 功能：复数组1-D快速Fourier变换! 输入参数说明：! Xreal(n): 输入数据实部! Ximage(n): 输入数据虚部! N: 点数（N 必须为2 的幂次方）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Xreal(n): 变换后频谱实部! Ximage(n): 变换后频谱虚部! 对应频率分布说明：! 数据Xreal(n)和Ximage(n)对应的频率分布位置为:! 0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!***************************************************************** SUBROUTINE FFT(Xreal,Ximage,n,nf)implicit noneINTEGER n,nfREAL Xreal(1:n),Ximage(1:n)integer nu,n2,nu1,k,k1,k1n2,l,i,ibitrreal f,p,arg,c,s,treal,timagenu=int(log(float(n))/0.693147+0.001)n2=n/2nu1=nu-1f=float((-1)**nf)k=0DO l=1,nu,1DO while (k.lt.n)do i=1,n2,1p=ibitr(k/2**nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k) thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageend ifEND DOIF(nf.ne.1) THENdo i=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)end doEND IFEND SUBROUTINE FFTFUNCTION IBITR(J,NU)implicit noneinteger ibitr,j,nuinteger j1,itt,i,j2j1=jitt=0do i=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doEND FUNCTION IBITR!******************************************************************** *******! 子程序：cal_dxdy! 功能：计算x,y方向的点距! 输入参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! N_point,N_line：点数(x方向)、线数(y方向)! 输出参数说明：! dx,dy：x,y方向的点距!******************************************************************** *******subroutine cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicit nonereal xmin,xmax,ymin,ymaxinteger N_POINT,N_LINEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEend subroutine cal_dxdy!****************************************************************** ! 子程序：WAVE2D! 功能：计算2D径向圆频率W! 输入参数说明：! dx: x 方向点距! dy: y 方向线距! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! 输出参数说明：! W(m,n): 径向圆频率!****************************************************************** SUBROUTINE WAVE2D(m,n,dx,dy,W)implicit noneINTEGER m,nREAL dx,dyREAL W(1:m,1:n),U(1:m),V(1:n)real pi,delx,delyinteger midm,midn,i,j,xx,yypi=3.1415926midm=m/2+1midn=n/2+1delx=pi*2.0/float(m)/dxdely=pi*2.0/float(n)/dydo j=1,n,1yy=jif(yy.gt.midn) yy=yy-nv(j)=dely*(yy-1)end dodo i=1,m,1xx=iif(xx.gt.midm) xx=xx-mu(i)=delx*(xx-1)end dodo j=1,n,1do i=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j))end doend doEND SUBROUTINE WAVE2D!******************************************************************* ! 子程序：calculate_Y! 功能：计算延拓因子Y! 输入参数说明：! distance:延拓的高度（m/z轴向下为正方向）! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! W：径向圆频率! 输出参数说明：! Y：延拓因子!******************************************************************* subroutine calculate_Y(m,n,W,y,distance)implicit noneinteger m,nreal pi,distancereal y(1:m,1:n),W(1:m,1:n)real i,jpi=3.1415926do i=1,m,1do j=1,n,1Y(I,J)=1*exp(-1*w(i,j)*distance)enddoenddoend subroutine calculate_Y!******************************************************************** *******! 子程序：calculate_FmulY! 功能：将延拓因子与原信号做乘积! 输入参数说明：! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! Y：延拓因子! Ur：初始信号的实部! Ui：初始信号的虚部! 输出参数说明：! Fr：延拓信号的实部! Fi：延拓信号的虚部!******************************************************************** *******subroutine calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)implicit noneinteger m,nreal Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n)real i,jdo i=1,mdo j=1,nFr(i,j)=Ur(i,j)*y(i,j)Fi(i,j)=Ui(i,j)*y(i,j)enddoenddoend subroutine calculate_FmulY!******************************************************************** **************************! 子程序：output_grd! 功能：按照grd格式输出延拓后的场值! 输入参数说明：! N_point,N_line：点数(x方向)、线数(y方向)! filename_output：输出文件的文件名! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! x_min,x_max,y_min,y_max：x/y的最小值、最大值! A：输出数组!******************************************************************** **************************subroutineoutput_grd(filename_output,N_point,N_line,m,n,m1,m2,n1,n2,eigval,A,Xmin,Xmax,Y min,Ymax)character*(*)filename_outputinteger N_point,N_lineinteger m,n,m1,m2,n1,n2real Xmin,Xmax,Ymin,Ymaxreal A(1:m,1:n)real amin,amax,i,jDO j=n1,n2,1do i=m1,m2,1if(ABS(A(i,j)).lt.eigval) thenamin=MIN(amin,A(I,j))amax=MAX(amax,A(I,j))endifEND DOOPEN(20,file=filename_output,status='unknown',form='formatted')write (20,'(A)')'DSAA'write (20,*)m2-m1+1,n2-n1+1write (20,*)Xmin,Xmaxwrite (20,*)Ymin,Ymaxwrite(20,*)amin,amaxDo j=n1,n2,1write(20,*) (A(i,j),i=m1,m2)End doClose(20)end subroutine output_grd!******************************************************************** **************************! 子程序：calculate_error! 功能：计算均方误差! 输入参数说明：! distance:延拓的方向! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! Fr,U：待求均方误差的两数组! 输出参数说明：! error:延拓后场值与真实场值的均方误差!******************************************************************** **************************subroutine calculate_error(distance,Fr,U,m1,m2,n1,n2,m,n,error)implicit noneinteger m1,m2,n1,n2,m,nreal Fr(1:m,1:n),U(1:m,1:n)real distancereal i,jerror=0.0do j=n1,n2do i=m1,m2error=error+(U(i,j)-Fr(i,j))**2enddoenddoerror=sqrt(error)if(distance>0.0)thenwrite(*,*)'请输入向上延拓的均方误差' elsewrite(*,*)'请输入向下延拓的均方误差' end ifwrite(*,*)errorend subroutine。
