频率域位场处理和转换实验
6第六章 磁异常的转换处理
+ a2 xi2
− ΔT (xi )⎤⎦
=0
∑ ∂δ
∂a1
=
m
2 ⎡⎣a0
i=−m
+
a1xi
+
a2 xi2
− ΔT (xi )⎤⎦xi
=
0
22
∑ ∂δ
∂a1
=
m
2 ⎡⎣a0
i=−m
+
a1xi
+
a2 xi2
− ΔT (xi )⎤⎦xi2
=
0
从而可由上述方程组解出 a0。若取 m=2，可采用五点圆滑时，当数据是等间距，点距△x=1，
24
∫ ΔT(x, z) = - z ∞ ΔT (ξ ,0) dξ
π −∞ (ξ − x)2 + z 2
(6-3-3)
式中 ΔT (ξ ,0) 为剖面上各点的实测值。
若坐标原点位于计算点下方实测剖面上，延拓高度为一倍点距，设为 h。即(6-3-3)式中 x=0，z=-h。式则(6-3-3)成为
∫ ΔT(0,-h) = 1 ∞ h ΔT (ξ ,0)dξ
实验二延拓——精选推荐
实验二频率域位场处理和转换实验******专业：勘查技术与工程学号：**********指导老师：王万银，纪新林目录一、基本原理 (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倍，且在边缘部位有几个异常点。
MAGS介绍(2.0)
（3）、任意形状三度体面元法正演
（4）、复杂几何形体正演，包括有限 长水平圆柱体、椭球体、走向与下延有 限倾斜板状体正演计算。
CUG
Institute of Geophysics and Geomatics
5、平面磁测资料转换处理
（1）、滑动平均法，插值切割场法， 趋势分析法，差值场法，匹配滤波， 3D频率域转换处理系统等方法是剖 面（二度）方法的推广,可实现平面 资料的各种转换处理。
CUG
Institute of Geophysics and Geomatics
• 70年代以来在我国，磁法勘探开始推广应用计 算机，但是工作人员计算都带着资料上北京计 算中心，80年代开始出现微机，但是DOC方式 下的黑屏和西文显示，曾给我们解释人员带来 诸多不便，90年代以来计算机科学的发展， Windows可视化技术的出现，使磁法勘探数据 处理上了一个新台阶：全中文提示可视化，实 时显示计算结果和绘制图件，以及以人机交互 方式修改参数和计算结果等等，新的软件处理 系统和工作方式可以使解释人员高效快速完成 内业整理与解释，甚至可以在野外条件下及时 获得处理解释的结果和绘制成果图件、编写报 告、制作汇报多媒体。
CUG
Institute of Geophysics and Geomatics
• 磁法勘探的仪器已经在80年代中期以电子式的 质子磁力仪取代了机械式的仪器，仪器的分辨 率由原来的1—5nT提高到0.1nT，现在广泛使 用的磁力仪是MP4、ENVI MAG、G856、G858 型仪器，工作效率也由一天几十、上百个测点 提高到几百、上千个测点。 • 室内资料的整理计算与图示再也不靠手工计算 与手工绘制，可以全部由计算机完成，一个小 组一周、一个月的工作现在一个人一天之内就 可以完成。
基于奇异谱分析的重磁位场分离方法
基于奇异谱分析的重磁位场分离方法朱丹;刘天佑;李宏伟【摘要】奇异谱分析是一种近年兴起的时间序列分析方法,它利用降秩原理实现信号分离.该方法将数据空间投影到不同特征的子空间中,并用奇异值来表征这些子空间的性质,最后通过截取奇异值实现数据的重构.重磁位场分离可以看成一种多信号叠加的分离问题.不同特征的重磁异常具有不同特征的奇异谱,这是奇异谱分析用于解决位场分离问题的应用基础.本文通过建立理论模型,分析重磁异常的奇异谱特征,得出适用于重磁位场分离的最优参数选择方法,并与传统方法进行比较.对比发现,无论是横向叠加模型、垂向叠加模型还是斜向叠加模型,奇异谱分析都具有很好的分离效果.最后,将奇异谱分析用于鄂东南某矿区的重力资料处理中,实现弱异常的识别和分离.【期刊名称】《地球物理学报》【年(卷),期】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。
低纬度磁异常的转换与处理
收稿日期２００８—０４—１ ５；修回日期２００８—０７—２５． 基金项目 国家基础研究９７３计划（２００７ＣＢ４１１７０１）．８６３计划（２００６ＡＡ０９Ａ１０１—０２０１—０２）联合资助． 作者简介赵百民，男，博士，１ ９７６年生，辽宁辽中人，主要从事航磁数据处理与解释工作．（Ｅ—ｍａｉｌ：ｚｈａｏｂｍ２００１＠１ ６３．ｃｏｍ） ＊通讯作者郝天珧，女，１９５７年生，北京市人，１９８２年毕业于长春地质学院，研究员，主要从事海、陆油气盆地综合地球物理研究
，取数的线距与
点距相差也不能大，以１：１或２：１或２．５：ｌ为 宜，当为４：１或更大的比例时，也会产生某些畸变 而影响化极效果．对于大比例尺的地面磁测数据，线 距与点距之比往往为４：１或５：１甚至更大比例， 此时应采取合适的手段插值，以缩小取数的线距与 点距之比．
低纬度磁异常的转换与处理
赵百民ｈ２， 郝天珧Ｈ， 徐 亚１
（１．中国科学院地质－ｑ地球物理研究所，北京１０００２９； ２．中国国土资源航空物探遥感中心，北京１０００８３）
摘 要 我国的南海大部分海域位于磁赤道带附近，属于低磁纬度区域，以水平磁化为主，磁性体产生的△丁异常特
征与中国大陆广大地区的△丁异常特征差别较大，且南北跨度大，达两千余公里，通常的化极技术在该区域很难取得
（２）运用畸变规律，用不同的化极倾角化极，比 较其化极异常图，既注意等值线形态的变化，也注意 正负异常强度的变化，选择合理的化极结果．
（３）磁异常数据区边界应尽量接近于零，尤其是 起始线，由于傅里叶变换的积分区域为正负无穷大 之间，而实际的数据区域总是有限的，为了减少因数 据区域有限性而引起的畸变（即吉布斯效应），应尽 量使数据边界值接近于零．但由于实际条件的限制， 边界值有时不接近于零，频率域位场转换程序一般 都采用了加权处理的方法，也就是使数据从边界向 外延伸，用余弦函数使延伸部分的数据逐步下降为 零；快速傅里叶变换要求数据的线数和点数都是２ 的整数幂，当实际的线数、点数不满足这一要求时， 则把线数、点数扩充为２的整数幂，从数据区域边界 向外扩充的部分作为向外延伸并逐步下降为零的部 分．
物探工作方法及技术要求—卓凤吉
第十五页，共87页。
工作装置 应用目的
1、中间梯度： 普查、多条线观测，效率高
特点：受低阻屏蔽较小； 穿透深度大； 地形影响小。
2、联合剖面： 寻找陡产状的地质体 确定目标体的顶端位置
分辨率高、效率低 3、对称四极测深： 研究地层垂向变化，确定界
4、数据处理
数据处理目的是提取岩性、构造等与成矿有关的信息，包括反演（包 括正演拟合）等方法。保存数据处理后的文件，其文件格式为网格化 文件（GRD）。 5、图件编制
基本图件有剖面平面图和等值线平面图。剖面平面图和等值线平 面图上颜色。剖面平面图正值区用红色，最好用渐变颜色。等值 线平面图正值区等值线之间用红色(浅红 — 深红），即深红 — 浅红对应大正 — 小正。处理图件种类根据数据处理结果来确定 ，处理图件按照等值线平面图的方法上颜色.图件格式一律为 MAPGIS的点、线、区文件。 6、报告编写 报告包括工作方法技术、方法试验、各项精度、图件编制、数据处 理、电性特征、异常编号、异常综合分析等内容。异常属性按照数 据库要求执行。
报告包括工作方法技术、方法试验、各项精度、图件编制、数据处 理、磁性特征、异常编号、异常综合分析等内容。异常属性按照数 据库要求执行。
第二页，共87页。
二、激发极化法测量
1、工作设计
按照工作任务和目的编写设计，包括方法技术、各项精度、原 始资料、图件编制、数据处理、报告编写等内容。方法技术部 分应包括方法试验内容，其中包括装置类型、供电极距、接收 极距、供电周期、延迟时间、积分时间、频率选择、观测次数 、观测精度等内容。
6、物探工作情况如实反映在项目工作年、季、月报。
航磁数据处理
航磁数据位场转换处理及效果∆测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测航磁T面上的综合反映。
由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。
为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。
第一节位场转换处理及效果航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。
频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。
位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的WINDOWS系统下地球物理数据处理解释软件（GeoProbe Mager）及航空物探彩色矢量成图系统（AgsMGis）。
一、原平面化极处理化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。
测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。
这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。
化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。
从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。
作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。
磁法勘探06磁异常的处理与换算资料
第一节 磁异常的处理与换算的目的意义
应当指出，磁异常处理和转换时，有两个问题必须要明确： 1.应当合理的选择处理和转换的方法。由于转换、处理方法 较多，具有各自的特点、作用、适用条件，不应盲目的对各 种方法都使用一遍。应当认真分析磁异常特征、测区内地质、 物性情况及所要解决的地质问题，根据各个方法的功能和适 用条件，合理的选择若干种处理方法; 2.磁异常的处理和转换只是一种数学加工处理，它能使资料 中某些信息更加突出和明显。但不能获得在观测数据中不包 含的信息。数学变换只能改变异常的信噪比，而不能提供新 信息；因此，在应用各个方法时必须要注意到实际资料的精
15
第二节 磁异常的处理
1.剖面网格化
16
17
第三节 磁异常的空间转换
延拓是把原观测面的磁异常通过一定的数学方法换算到高 于或低于原观测面上,分为向上延拓与向下延拓。向上延拓 是一种常用的处理方法，它的主要用途是削弱局部干扰异 常，反映深部异常。我们知道，磁场随距离的衰减速度与 具磁性的地质体体积有关。体积大，磁场衰减慢；体积小， 磁场衰减快。对于同样大小的地质体，磁场随距离衰减的 速度与地质体埋深有关。埋深大，磁场衰减慢；埋深小， 磁场衰减快。因此小而浅的地质体磁场比大而深的地质体 磁场随距离衰减要快得多。这样就可以通过向上延拓来压 制局部异常的干扰，反映出深部大的地质体。
是很重要的。随着磁测量精度的不断提高，实测异常中所包含 的可靠信息也不断增加。如何有效地提取和利用这些信息，就 成为磁异常解释理论研究的重要课题。早在20世纪40、50年代， 诸如导数异常的计算，磁场解析延拓，化磁极等处理方法已相 继问世。到60、70年代，由于电子计算机的广泛应用，使磁异 常的处理和转换容易实现，从而其理论和方法得到了迅速的发 展，并不断得到完善。由于在实践中磁异常的转换和处理对提 高磁方法解决问题的能力和改善地质效果起到了应有的作用， 因此它已成为当今磁异常推断解释中不可缺少的重要环节。
首钢秘鲁马尔科纳铁矿区磁异常的数据处理
意的化极结果。
２．４ 双 曲 正 （余 ）弦 函 数 法 ［５］
通过给化极因 子 的 分 母 加 上 一 个 校 正 函 数，使
该函数沿经向增强而且与转换后的磁化方向有关。
这一函数的加入 既 可 压 制 高 频 干 扰，又 可 改 善 分 母
接 近 ０ 时 的 不 稳 定 性 ，根 据 这 一 思 想 ，采 用 双 曲 函 数 构建校正函数。
２．３ 阻 尼 因 子 法 ［４］
同 样 ，为 改 善 低 纬 度 地 区 化 极 时 ，化 极 因 子 分 母
趋于０而产生过 度 放 大 的 情 况，将 化 极 因 子 的 分 母
加 一 个 很 小 的 值 ，则 化 极 因 子 变 为 一 个 有 效 值 ，这 个
附加值起到了阻尼的作用。阻尼因子形式如下：
只存有矿区 航 磁 图 １ 张。 首 钢 收 购 该 铁 矿 矿 权 后， 数是连续的，且积分上下限为正负无穷，而通常我们
并未投入大范围 的 铁 矿 普 查 工 作，但 随 着 近 年 来 铁 所获得的磁异常 数 据 总 是 离 散 的，且 测 区 也 总 是 有
矿石价格的不断 上 涨，在 矿 权 区 内 进 行 资 源 地 质 勘 限的范围，因此 转 换 中 容 易 产 生 假 频 而 引 起 高 频 振
将常规化极运 算 应 用 于 低 纬 度 地 区，得 到 的 结 果会有较大 的 误 差，异 常 走 向 也 会 发 生 变 化。 因 此 在低纬度地区进 行 化 极 运 算，必 须 进 行 特 定 的 处 理 以 压 制 干 扰 和 畸 变 ，但 数 据 处 理 时 赋 值 复 杂 ，有 时 的 计 算 量 很 大 ，很 难 取 得 满 意 的 效 果 。 ２．１ 低 纬 度 区 化 极 的 不 稳 定 因 素 ［２］
基于频率域磁异常转换的断裂构造划分判据
基于频率域磁异常转换的断裂构造划分判据宁堃;谢涛;邵昌盛;邓子清【摘要】实践表明,利用磁异常可以较为准确的对断裂、区域深大断裂、褶皱等构造进行推断划分.实践中,通过对西藏某热液型铅锌多金属矿区实测的磁异常数据进行频率域垂向二阶导数和水平梯度模的转换处理,突出一些浅部地质体或边界引起的弱磁异常响应,结合实际地质情况和?T、垂向二阶导数、水平梯度模的断裂确定标志,推断划分了多条断裂构造,取得了较理想的效果.【期刊名称】《四川地质学报》【年(卷),期】2018(038)002【总页数】5页(P332-336)【关键词】磁异常;频率域;垂向二阶导数;水平梯度模;断裂构造【作者】宁堃;谢涛;邵昌盛;邓子清【作者单位】四川省核工业地质局二八二大队,四川德阳 618000;四川省核工业地质局二八二大队,四川德阳 618000;四川省核工业地质局二八二大队,四川德阳618000;四川省核工业地质局二八二大队,四川德阳 618000【正文语种】中文【中图分类】P631地球物理勘探中，用磁法可以圈定断裂破碎带，是因为断裂的产生改变了岩石的磁性，或者改变了地层的产状，或者沿断裂带伴有后期或同期岩浆活动，或者沿断裂两侧具有不同的构造特点。
一般情况下，在DT等值图上可以根据磁异常特点直接对断裂进行划分。
然而，在地质情况复杂的地区，迭加异常相互干扰，弱信号与误差并存，使异常变得十分复杂而难以解释，这就需要根据实际情况对磁异常进行精细化处理和转换，从而进一步提高解释推断的精度。
以西藏某热液型铅锌多金属矿的高精度磁法数据为例，介绍如何综合利用频率域磁异常垂向二阶导数和水平梯度模技术对断裂构造进行推断划分。
频率域磁异常转换首先要将磁异常原始数据经过傅氏变换转换为频谱，然后乘以相应的频响因子，再经过反傅氏变换得到结果，这种处理方法对磁场高频成分有突出和放大作用，它侧重于浅层近地表地质的磁效应而压制深层区域背景场的影响，从而突出浅部地质体引起的局部异常和解释地质体的边界（断裂、岩体、矿体边界等）。
光学系统的光学传递函数OTF测定方法理论(实验)研究 - 终稿
本科毕业设计（论文）光学系统的光学传递函数OTF 测定方法理论（实验）研究学 院_ 物理与光电工程学院__专 业_____ 光信息科学与技术_（光电显示与识别技术方向）年级班别________2010级(2)班__学 号_________3110008945______学生姓名___________林清贤___指导教师___________雷 亮____2014 年 4 月 28 日摘要光学传递函数是定量描述成像性能的完备函数。
但是对于实际的光电成像器件（如CCD器件），通过解析法建立这一函数的表达式又是非常困难的，因此光学传递函数的实测技术就显得尤为重要。
光学传递函数是一个客观的、准确的、定量的像质评价指标，并且其能够直接方便的测量，因此已经广泛应用于光学设计、加工、检测和信息处理中。
本文主要介绍了光学传递函数的性质及其测量原理分析，并对固有频率目标法和狭缝扫描法进行了实验研究。
我们采用光学显微镜作为待测量光学传递函数的光学系统，通过改变显微镜的放大倍数，比较分析放大倍数对调制传递函数（MTF）测量的影响，并比较两种测量方法的优劣。
实数傅立叶变换是整个实验中需要透彻理解和运用的数学概念，在此基础上理解离散傅立叶级数与MTF定义的理论依据，并由此建立数学模型。
由本文建立的理论模型出发，结合实验所测得的数据，最后得到了基本可靠的实验结果。
本文最终给出两种测量法对应的matlab程序、数值测量结果、实验测得的可靠的MTF实验结果撰写毕业论文主要内容。
关键字: 光学传递函数，傅立叶变换，固有频率目标法，狭缝扫描法AbstractThe optical transfer function is quantitatively describe the imaging performance of the complete function.But for the actual photoelectric imaging devices (such as CCD device), through the analytic method to establish the function of expression is very difficult.Therefore the measurement technique of optical transfer function is particularly important.Optical transfer function is an objective, accurate and quantitative image quality evaluation index, and it can directly and convenient measurement, therefore has been widely applied optics design, processing, testing and information processing.This paper mainly introduces the properties of the optical transfer function and its measuring principle, and the inherent frequency target and slit scan method has carried on the experimental study.We use optical microscope as for measuring optical transfer function of optical system, through changing the magnification of the microscope, comparative analysis of magnification of modulation transfer function (MTF) measurement, the influence of the merits of the two measuring methods are compared.Real Fourier transform is the need to thoroughly understand and apply in the experiment of mathematical concepts, on the basis of the understanding of discrete Fourier series and the theoretical basis of the definition of MTF, and thus to establish mathematical model.Set up by this article on the theory model, combined with the data measured in laboratory, the fundamental and reliable experiment results are obtained.Finally, the paper proposes two kinds of measurement method of the corresponding matlab program, the results of numerical measurement and reliable experimental measured MTF experimental results of writing graduation thesis main content.Keywords: Optical transfer function, Fourier transform, Natural frequency method; Slit scan method目录第一章绪论 (1)1.1 光学传递函数简介 (1)1.2 光学传递函数的发展 (1)1.2.1 光学传递函数的发展历史 (1)1.2.2 光学传递函数的发展现状和趋势 (2)1.3 光学传递函数的测量意义 (3)1.4 本论文的主要内容 (4)第二章光学传递函数的基本理论 (5)2.1 光学成像系统的一般分析 (5)2.1.1 透镜的成像性质 (5)2.1.2 光学成像系统的普遍模型 (8)2.1.3 两种类型的物体照明方式 (9)2.1.4 阿贝成像理论 (9)2.2 光学传递函数的概念 (10)2.3 光学传递函数的计算 (12)2.3.1 以物像频谱为基础的计算 (12)2.3.2 以点扩散函数为基础的计算 (13)2.3.3 线扩散函数与一维调制传递函数 (13)2.4 离散傅里叶级数与MTF定义的理论依据 (14)第三章光学传递函数的测量原理分析 (17)3.1 光学传递函数的测量方法综述 (17)3.2 实验中的两种测量方法原理分析 (18)3.2.1 固有频率目标法 (18)3.2.2 狭缝扫描法 (20)3.3 光学传递函数测量系统软件 (21)3.4 CCD对光学传递函数测量的影响分析 (22)第四章光学传递函数测量实验及实验结果分析 (23)4.1 实验平台的搭建 (23)4.2 固有频率目标法实验 (23)4.3 狭缝扫描法实验 (25)4.4 两种测量实验结果分析 (31)第五章总结与展望 (32)参考文献 (33)致谢 (34)第一章绪论1.1 光学传递函数简介在应用光学领域中，有一个大家一直所瞩目的问题，那就是对光学系统成像质量的评价。
频率域偶层位方法在直升机磁测数据处理中的应用
笔者 提 出了一套 频率 域偶 层 位 曲面位 场数 据处 理和转 换 方法 ， 能够 处理 大数 据量 的航 磁数 据 ， 具有
收 稿 日期 ：０７— ９—１ ２０ ０ ０ 基 金项 目 ： “ 五 ” 技 攻关 课 题 （ ｏ ｌ Ａ ０ Ａ一５ 资 助 国家 十 科 ２ｏ Ｂ ６ ９ ）
磁场与地面磁测结果对应很好 ， 而且 曲面化极结果 与已知地质 体也有 很好 的对应 ， 实了频率 域偶 层位 曲面位场 证
数据处理和转换方法实用性强 ， 可用于起伏测量条件下航磁 数据的处理工作 中。 关键词 ： 航磁测量 ； 数据处理 ； 率域 ； 频 偶层位 ； 直升机
中图 分 类 号 ： ６ １ Ｐ ３ 文献标识码 ： Ａ 文 章 编 号 ：１０ ０ ０—８ １ （０ ７ ０ ０ ７ ９ ８ ２０ ）６— ５ ７—０ ４
和方法 ¨ 。如依据等效源原理发展 了单层位和偶
层位 法等 不 同的等 效 源 法 ， 比较 有 代 表 性 的作 者有 侯 重初 、 连喜 、 万银 等 。但 是 由于 遇到 了一 些诸 徐 王 如大数 据 量 的处 理 和 转换 、 型 的第 一 类 Ｆｅｈｌ 大 ｒｄ ｏ ｍ
面化极、 延拓 、 垂向一 阶导数、 延拓到地表 等处理成 果 图件 ， 已知资料 对应 较好 ， 与 为该 区地质 解释 提供
量来 说就 不能 忽 略 了 ＿ 。因 此 ， 须 利 用 三 维 曲面 ９ Ｊ 必
２ 方 法 技 术
２ １ 频 率域 偶层 位 法基本 原 理 ．
偶层 位法 曾是 空 间域 的一种 基 于等效 源 的 曲面 位 场处理 和转 换 方法 ， 由于在空 间域 实 现 ， 故将 其称
位 场数据 处理 和转 换方 法才 能保 证数 据处 理结 果 的
(整理)地球物理勘探习题集.
习题全集第一章习题（重力）1．绘出下图中各点的引力、离心力和重力的方向。
2. 假定地球是一个密度均匀的正球体，位于球心处单位质点所受的引力应是多大？有人说，按牛顿万有引力定律，该处的引力应为无穷大(因为2lim r GM r ），对不对？为什么？3．“引起重力变化的因素就是引起重力异常变化的因素．”这种说法对吗？为什么？4．将地球近似看成半径为6370km 的均匀球体，若极地处重力值为9．82/m s ，试估取地球的总质且J 为多少吨？5．重力等位面上重力值是否处处相等7为什么?如果处处相等，等位面的形状如何?如果重力有变化，等位面的形状又有何变化？6．从我国最南边的南沙群岛(约北纬5)到最北边的黑龙江省漠河(北纬约54)，正常重力值变化有多大(用赫尔默特公式计算）？7．请用赫尔默特公式计算：1)两极与赤道间的重力差是多大？2)若不考虑地球的自转，仅是由于地球形状引起的极地与赤道间重力差为多少？8．假定沿某一剖面上各点的正常重力值其大小和方向皆相同，试示意绘出当地下存在有剩余密度小于零的球形地质体时，沿剖面各点的重力分布图。
9．请考虑如图2所示的两种剖面情况，能否在地面上观测出有相对变化的重力异常来?这对重力勘探的应用条件提供了什么启示？10．请给出图3中各剖而内研究对象(划斜线的部分)的剩余密度值。
11．“同一质量的地质体在各处产生的重力异常胶应该都一样。
”你能指出这一说法的正误吗?12．若有一剩余质量为50万吨的球形矿体(可当作点质量看)，当其中心埋深为100 m 时，请计算：1)在地面产生的异常极大值是多少?2)异常值为极大值的1/3的点距极大值距离为多少米？3)若该矿体与岩围密度分别为3．03/g cm 和2．53/g cm ，该矿体的实际质量是多少吨？13．解释下列名词：地磁要紊、国际地磁参考场、地磁图、IGRF ．14．试述地磁场随空间、时间变化的基本特征？15．磁偏角在全球有几处为不定值？为什么？16．绘图表示一个通过重心绕水平轴自由转动的磁针其水平轴分别平行于磁子午面、垂直于磁子午面由地磁南极向北极移动时，磁针的静止状态。
空间滤波的实验研究
空间滤波的实验研究于雪冰;王伟【摘要】运用空间域和频率域方法讨论了阿贝-波特实验和空间滤波实验中光波分别经过单透镜2f系统和4f系统时透镜的傅里叶变换和物体的频谱,探究了物体成像过程中的分频和合频过程,分析了空间滤波实验中透镜孔径对实验的影响.应用M atlab仿真模拟了低通滤波、方向滤波、显色滤波等图像处理技术,验证了夫琅禾费衍射下的巴比涅原理.【期刊名称】《物理实验》【年(卷),期】2018(038)001【总页数】6页(P8-13)【关键词】阿贝成像原理;空间滤波;傅里叶变换【作者】于雪冰;王伟【作者单位】东北师范大学物理学院,吉林长春130024;东北师范大学物理学国家级实验教学示范中心(东北师范大学) ,吉林长春130024【正文语种】中文【中图分类】O438傅里叶变换光学是光学领域的一个分支，它的形成导致了光学信息处理技术的兴起. 光学信息处理以其容量大、速度快、并行性等显著优点，在二维图像信息处理和识别等方面有重要应用. 空间滤波是最基本的光学处理操作之一，其基本原理是根据具体需要制作适当的空间滤波器，并将其放在光路中输入图像的频谱平面处，通过对输入图像的频谱进行调制完成某种处理过程，如低通、高通、带通、边缘增强、相关识别等. 其理论基础是傅里叶变换，实验基础为阿贝成像原理，了解相关理论对掌握光学信息处理技术起着至关重要的作用.1 实验原理1.1 空间频率与频谱任意周期结构的屏函数均可以展开为傅里叶级数. 傅里叶系数的集合反映了原函数各种频率成分所占的分量，通常称其为傅里叶谱，简称频谱. 频谱可以是连续谱，也可以是离散谱. 周期函数的频谱是离散谱，非周期函数的频谱是连续谱. 实际栅函数为准周期函数，其频谱介于连续谱与离散谱之间，而更具有离散谱的特征，称为准离散谱.根据傅里叶分析可知，频谱面上的光场分布与物的结构密切相关，原点附近分布着物的低频信息，即傅里叶低频分量；离原点较远处，分布着物的较高的频率信息，即傅里叶高频分量[1].1.2 阿贝成像原理阿贝成像原理如图1所示，从频谱角度看，阿贝成像原理的基本思想是把相干光照明下的透镜成像过程分为两步：1)物是一系列的不同空间频率信息的集合，通过物的衍射光在透镜后焦面(频谱面)上形成空间频谱，所以衍射起“分解”频谱即“分频”的作用；2)代表不同空间频率的各光束在像平面上相干叠加而形成物体的像，因此干涉起“综合”频率即“合频”的作用[2].图1 阿贝成像原理图1.3 透镜的相位变化功能设物体的复振幅透射率为t(x0,y0)，物体与透镜间的距离为d0. 使用振幅为A的单色平面波垂直照射物体，U(x0,y0)为紧靠物体后平面上复振幅分布，U1(x,y)为紧靠透镜前平面上复振幅分布，则有U(x0,y0)=At(x0,y0),(1)F{U0(x,y)}=AF{t(x0,y0)}=AT(fx,fy),F{U1(x,y)}=F{U0(x0,y0)}L(fx,fy),(2)空间频率为忽略常量相位延迟，则有得(3)后焦面上的复振幅分布为(4)由上面的分析可见：后焦面上的复振幅分布正比于物体的傅里叶变换,变换式前的二次相位因子使物体的相位因子产生相位弯曲.当d0=f时，即当物体位于透镜前焦面时：(5)这时相位弯曲完全消失，后焦面上的光场分布,即置于前焦面物体的频谱,是物体准确的傅里叶变换[3].1.4 空间滤波对图像产生的复杂波前的傅里叶分析，意味着将其复杂的衍射场分解为一系列不同方向、不同振幅的平面衍射波，特定方向的平面衍射波，作为载波，携带着特定空间频率的光学信息，并将其集中于夫琅禾费衍射场的相应位置，实现了分频. 因为物信息的空间频谱展现在透镜的后焦面即傅氏面上，故若在频谱面上安置不同结构的光阑，以提取或摒弃某些频谱，从而改变了原物频谱，再合成于物的共轭像面上即为输出图像，这就完成了改造图像的信息处理. 频谱面上的光阑起选频作用，常称为空间滤波器[1].1.5 巴比涅原理2个互补衍射屏在衍射场中某点单独产生的复振幅之和等于光波自由传播时该点的复振幅，称为巴比涅原理. 巴比涅原理给出的3个场之间是复振幅关系，其中相位差因素也会起作用，故一衍射屏在某处的衍射强度是亮的，其互补屏在该处的衍射强度不一定是暗的.2 实验2.1 探究分频、合频过程调整好光路后，用准直的氦氖激光照明带网格的“光”字板. 先将白屏放在傅里叶透镜的后焦平面前且靠近傅里叶透镜，可看见所成像中央有一轮廓较为清晰的正立、缩小的“光”，并且“光”的周围有很多重影，发现网格已经分离[图2(a)]. 当将白屏慢慢向后焦平面移动时，中央的“光”已经被分解成更多的正立“光”，且轮廓越来越模糊，“光”也变得越来越小[图2(b)]，直至频谱面时网格已经完全分离成点阵，中央没有“光”字，只有几个光强较强的点[图2(c)]. 继续将白屏远离频谱面后方移动，白屏上的点阵逐渐消失，慢慢合成不很清晰的倒立“光”字，将屏后移，出现轮廓分明的倒立、逐渐放大了的“光”字，同时点阵也慢慢扩展复合成网格的像[图2(d)]，说明此时发生了频率的合成. 当把白屏再向后移动时，发现倒立的“光”越来越大，其上网格也越来越清晰[2][图2(e)].(a) (b) (c)(d) (e) 图2 探究“分频”及“合频”过程2.2 方向滤波2.2.1 2f成像系统物为正交光栅，滤波器为可旋转狭缝. 用准直的氦氖激光照明光栅，后焦面上出现一系列准离散的衍射谱斑. 在后焦面上安置可以旋转的狭缝作为滤波器，以选取不同谱斑，从而可以观测到相应不同的输出图像. 如果频谱面上放置的狭缝沿纵向，则输出图像只显示横条纹；如果狭缝处于水平方位，则输出图像只显示竖条纹；如果狭缝取向倾斜，则输出图像显示为较密而且与狭缝取向垂直的斜条纹. 即改变狭缝方向，观测到像的延展方向总是与谱斑铺展方向正交，表明横向的谱斑携带的是纵向信息. 因为斜向谱斑的角间隔比水平或垂直铺展的角间隔要大，对应的基频较高，所以呈现于像平面上的斜向条纹较密，如图3～4所示.(a)狭缝竖直 (b)狭缝水平 (c)狭缝倾斜图3 方向滤波图(透镜为大孔径)(a)狭缝竖直 (b)狭缝水平 (c)狭缝倾斜图4 方向滤波图(透镜为小孔径)2.2.2 4f成像系统光路图如图5所示. 物为正交光栅，滤波器为可旋转狭缝. 用准直氦氖激光照明光栅，后焦面上出现一系列准离散的衍射谱斑；在后焦面上安置可旋转狭缝作为滤波器，以选取不同谱斑，从而可观测到相应不同的输出图像. 如果频谱面上放置的狭缝沿纵向，则输出图像只显示横条纹；如果狭缝处于水平方位，则输出图像只显示竖条纹；如果狭缝取向倾斜，则输出图像较为密集且与狭缝取向垂直的斜条纹.4f系统中的“后焦面”有双重身份，对L而言是物场的频谱面；对L′而言是物平面，其频谱面即为系统的输出平面. 4f成像系统中，前后2个透镜共焦组合是必要条件，其保证了前后2次波前变换均为纯净的傅里叶变换. 透镜的前后2个焦面是1对傅里叶变换面，在4f系统中，像场是一系列不同方向平面波的干涉场，而前半部分的物场被分解为一系列不同方向的平面衍射波，即为阿贝成像原理中的“一分一合”在4f系统中的特别体现.图5 4f成像系统原理图对比分析2f成像系统和4f成像系统透镜孔径不同时输出图像的差异，发现透镜孔径较大时对应的输出图像包含更多细节，这是因为较大的孔径可以收集到高频信息引起的大角度的衍射光，这些衍射光到达像平面时相干叠加出较多的细节. 而孔径较小时，高频信息引起的大角度衍射无法进入镜头，频谱面上缺少了高频谱，像面上丢失了高频信息，如图6～8所示.(a)狭缝竖直 (b)狭缝水平 (c)狭缝倾斜图6 方向滤波图(两透镜均为小孔径) (a)狭缝竖直 (b)狭缝水平 (c)狭缝倾斜图7 方向滤波图(透镜孔径为一大一小)(a)透镜为小孔径 (b)透镜为大孔径图8 不同孔径透镜对应的输出图像2.3 显色滤波实验装置如图9所示，采用白光作照明光源，频谱面上同时展现图像的时间频谱与空间频谱. 在频谱面上特定位置设置小孔滤波器，提取特定波长的空间频率成分，如图10所示.图9 显色滤波实验装置图10 截取不同频率滤波对应不同输出图像调整滤波孔位置，输出图像的色彩发生改变，说明色彩是人为指定的而非天然色. 2.4 利用空间滤波技术进行图像处理实验光路为2f成像系统，物为正交光栅，用孔径大小可调的圆孔光阑作滤波器. 实验中改变光阑孔径大小，观测输出图像的变化情况. 当逐渐缩小孔径时，观测到图像的边缘逐渐变得柔和，这是因为图像的傅里叶变换频谱中的低频分量反映图像的背景，高频分量反映图像的细节、边缘及其他尖锐跳跃，孔径缩小，使透过的高频分量减少，故边缘变得柔和，如图11所示.(a) (b) (c) 图11 圆孔光阑孔径逐渐减小对应的滤波图像2.5 低通滤波2.5.1 物为带有周期性网格的“光”字在2f成像系统中，物为带有周期性网格的“光”字，采用孔径很小的圆孔做低通滤波器，观察后焦面上的频谱分布，可以看到排成十字形的点阵. 逐步减小圆孔孔径，观察输出图像变化情况. 孔径较大时，像中存在网格结构，逐渐减小孔径，最终观测到没有网格的“光”字. 因为与网格相比，“光”字的空间频率较低，集中在光轴附近很小范围内，而孔径较小的圆孔只通过低频分量，故可将周期性网格消除，如图12所示.(a)孔径较大 (b)孔径较小图12 不同光阑孔径的低通滤波图像将小圆孔移至频谱面上中央亮点以外的亮点上时，在输出平面上仍能看到无网格的“光”字，只是较暗淡. 这说明当物为“光”与网格的乘积时，其傅里叶谱是“光”的谱与网格的谱的卷积，因此每个亮点周围都是“光”的谱，再作傅里叶变换就还原成“光”字[4]，如图13所示.(a) (b)图13 滤波孔不在频谱中心时对应的输出图像2.5.2 物为“大”字与正交光栅组合在2f成像系统中，将正交光栅与不透明的“大”字重叠放在物面上,选取孔径很小的圆孔作低通滤波器，观测到周期性网格被消除[5],如图14所示.2.6 利用Matlab模拟傅里叶变换根据傅里叶变换的性质，2个函数卷积的傅里叶变换等于傅里叶变换的乘积. 在频谱面上插入空间滤波器相当于频谱分布函数乘以空间滤波器滤波函数的复振幅透过率函数. 空间滤波的光学处理器的模拟系统简图如图15所示，通过计算机模拟仿真可以完成空间滤波实验[6].(a) (b)图14 低通滤波后的输出图像图15 空间滤波光学处理器的模拟系统简图2.6.1 物二维光栅的频谱将二维光栅作为物,则可在傅里叶面上观测到如16图所示的频谱分布.图16 二维光栅的频谱图2.6.2 低通滤波的模拟结果在计算机模拟中，用Photoshop软件画出带有周期性网格的“光”字图片代替物体，并保存为bmp格式. 通过Matlab编程对这幅图进行傅里叶变换得到相应的频谱分布. 这一步骤相当于实验中透镜所起的傅里叶变换作用. 图17所示为原图像及其频谱图分布.(a)未放置滤波器(b)放置滤波器后图17 三维频谱图2.7 巴比涅原理的探究在2f成像系统中，在物平面分别放置方孔和去除方孔的屏，观察二者后焦面上夫琅禾费衍射图样的区别. 在平行光照明时，其自由光场聚焦于透镜的后焦点，即轴外自由光场为零. 由巴比涅原理知，在平行光照明下，2个互补屏在后焦面上产生的夫琅禾费衍射强度分布是完全相同的，看起来是完全相同的衍射图样，不同的仅仅是像点的光强. 图18所示实验现象与理论符合[1].(a)衍射屏为方孔 (b)衍射屏为方孔的互补屏图18 巴比涅原理实验验证图3 结束语空间滤波是目前应用较为广泛的光学信息处理技术，其理论依据为阿贝成像原理[7]. 本文通过实验验证了阿贝的二次成像原理，通过改变频谱结构改变了输出图像的性质，并用Matlab对相关过程进行了模拟. 其中在实现显色滤波过程中，应用自制的小孔滤波器改变频谱结构，实现了对输出图像色彩的改变，该方法操作简便且得到了明显的实验现象.【相关文献】[1] 钟锡华. 现代光学基础[M]. 北京：北京大学出版社，2003.[2] 彭小兰，王红成. 阿贝成像原理中“分频、合频”的实验演示[J]. 东莞理工学院学报，2011,18(3):38-41.[3] 冯璐. 空间滤波实验中光路和傅里叶变换透镜孔径对实验的影响[J]. 物理与工程，2010,20(4):26-28,35.[4] 杨述武，孙迎春，沈国土，等. 普通物理实验(三、光学部分)[M]. 5版. 北京：高等教育出版社，2016:106.[5] 何钰. 阿贝成像原理和空间滤波实验及计算机模拟实验[J]. 物理与工程,2006，16(2)：19-23.[6] 谢嘉宁，赵建林. 光学空间滤波过程的计算机仿真[J]. 光子学报,2002，31(7)：847-850.[7] 张朝晖，刘国超. 阿贝成像原理和空间滤波实验[J]. 物理实验，2017,37(9):23-29.[8] 朱昊，曹良才，何庆声. 空间滤波与体全息光存储实验[J]. 物理实验，2014,34(9)：4-8.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《重磁资料处理与解释》实验二频率域位场处理和转换实验学院：地测学院专业名称：勘查技术与工程学生姓名：学生学号：指导老师：提交日期：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。