航磁数据位场转换处理及效果
航磁∆T 测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测面上的综合反映。由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。
第一节 位场转换处理及效果
航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。
位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的
WINDOWS 系统下地球物理数据处理解释软件（GeoProbe Mager）及航空物探彩色矢量成图系统（AgsMGis ）。
一、原平面化极处理
化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。
作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。从项目组野外物性测量结果看，区内多数岩石以感磁为主，剩磁方向与感磁方向接近，符合化极的前提条件。
全区采用" 频率域偶层位变倾角磁方向转换方法" 实现磁场全变倾角化极。在观测面上建立笛卡尔直角坐标系，使x 轴志向磁北，z 轴垂直向下。假设观测场∆T 是
一分布在观测面下方z=h平面上的偶层磁荷面引起的。它在观测点P(x,y,z)处产生的磁位U 与磁场T 分别为
U =-1⎛1⎫⎰⎰M ∆ ⎪d ξd η (1) 4π⎝r ⎭
∆T =- (2) u t ∇U 0
(3)
设
⎛∂∂∂⎫∆= , , ⎪ ⎝∂x ∂y ∂z ⎭
t =t 0+∆t 0 (4)
l =l 0+∆l 0 (5)
式中M 表示偶层磁荷面的磁化强度矢量； t 0、l 0表示均为常矢量研究区t 、l 的平
均值； ∆t 0、∆l 0表示他们的变化值。
将(4)式代入(2)式并进行傅氏变换得：
F [∆T ]=-u g 0)t 0{( [F ]U +[ F ∆ t ∇]} U 0 (6)
移项得
F [U ]=-
式中
⎫ 1⎧1∇U ]⎬⎨F [∆T ]+F [∆t 0 g t 0⎩u 0⎭ (7) g =(i 2πu , i 2πv , 2πf ), f =对(1)式两端作傅氏变换得：
F [U ]=-1⎡g ⋅F [M ]⋅4π⎣-2πf (h -)z e /⎤⎦f (8) 同样的，将(5)式代入(3)式进行傅氏变换后代入(8)式得：
F [U ]=-1{(g ⋅l )⋅F [M ]+g ⋅F [M ⋅∆l ]}⋅e 4πf 00-2πf (h -z )
移项得
F [M ]=-14πf ⋅e 2πf (h -z )F [U ]+g ⋅F [M ⋅∆l 0] g ⋅l 0{}(9)
（7）式、(9)式即为频率域变倾角磁方向转换的两个基本公式。已知观测场∆T ，可应用（7）式计算F [U ]，再将F [U ] F [U ]代入(9)式中计算F [M ]。再把F [M ]代入地磁极处的F [T ]极，即可实现变倾角化极处理。F [T ]极与F [M ]的关系如下。
在地磁北极有：t 0=l 0=(0,0,1), ∆t 0=l 0=0，由(7) 式、(9)式可得：
F [M ]极=π
-2πf (h -u 0f e )z (10) [F ]M
经傅立叶反变换后可得化磁极长Z ⊥
πf (h -)z Z ⊥=F -1{πu (11) f -2e M [F ]}0
假设要换算的场为T p ，其磁场方向单位矢为t ，磁化方向单位矢为l ，则只要把
他们代入(7) ，令∆t 0、∆l 0偏差为零，即得F [U ]的初值F 0[U ]：
F [U ]=-1
1 F ⎡T p ⎤⎣⎦u 0g ⋅t 0 对F 0[U ]反变换求得初值U 0，把U 0代入(7)式得F 1[U ]，如此反复迭代，直到求得的U 值之差小于给定的标准为止。求得F [U ]后，类似地取F 0[U ]为
F [M ]=--4πf F [U ]⋅e 2πf (h -z ) 0g ⋅l 0
反变换求得M 0，把M 0代入(9)式，求得F 1[M ]，如此反复可最后求得F [M ]，并代入(11)式求得化磁极磁场。
由于本区处于xxx°xxx ′~xxx°xxx ′，属于中纬度地区，斜磁化能够产生一定影响，对原磁场数据进行化极处理后，在垂直磁化的条件下，磁异常的形态以及磁异常与磁性体的关系都比较简单，便于进行地质解释。对比航磁∆T 等值线平面图和航磁∆T 化极等值线平面图，航磁化极处理作用非常明显（图1）：局部异常整体向北偏移，表明通过化极处理，使异常回归到磁性地质体上方；减小或消除了由于斜磁化而引起的多数局部异常正负异常伴生现象，为进一步圈定岩体边界创造条件；使异常带及梯度带更加明显，有利于揭示出不同地质体的分布与形态，对圈定各种不同类型的断裂、确定磁性体的性质及边界具有重要的意义。
图1 航磁∆T 化极处理效果对比图
a-航磁∆T 等值线平面图；b-航磁∆T 化极等值线平面图
二、化极垂向导数处理
航磁局部异常通常是叠加在区域背景场上的次级异常，在原始航磁或化极航磁等基础图件中表现并不明显，需要通过一定的数学处理手段来突出其特征。垂向导数处理是解决这个问题的一种有效手段，它反映了磁场在垂直方向上的梯度变化，在增强由浅部磁性体引起的局部异常、压制长波区域场有很强的功能，可以突出在总场图上不明显的细节，并能分解横向叠加异常，理论上导数的次数越高，这种分辨能力就越强。
磁异常垂向导数换算公式如下：
如果令S zx (x , y , z )、S zy (x , y , z )、S zz (x , y , z )及S zzz (x , y , z )、S zyy (x , y , z )、S zzz (x , y , z )分别为Z a (x , y , z )对x 、y 、z 的一阶导数及二阶导数的频谱，则有微分定力易于得到：
(S zx (u , v , z )=2πiuS z (u , v ,0)e 2πu 2+v 2)1/ z 2
(S zy (u , v , z )=2πiuS z (u , v ,0)e
S 2πu 2+v 2)1/2 z 2πu 2+v 2(u , v , z )=2πu 2+v 2
zz ()1/2S z (u , v ,0)e ()1/2 z
同理，可以写出：
⎧2⎪(2πiu )⎧S zxx ⎫⎪⎪⎪2⎪S u , v , z =)⎨(2πiv )⎨zyy ⎬(⎪⎪⎪⎩S zzz ⎭⎪⎡2πu 2+v 2
⎪⎢⎣⎩⎫⎪ ⎪2π(u 2+v 2)1/2 z ⎪S u , v ,0e )⎬z (⎪1/22⎤⎪⎥⎪⎦⎭
n ()由此可知，求磁场的n 阶垂向导数的频谱，应乘上的导数因子为⎡2π(u 2+v 2)1/2⎤；⎢⎣⎥⎦
而求磁场沿x 方向或y 方向的n 阶水平导数的频谱，应乘上导数因子为(2πiu )n 或
n (2πiv )。
如果求磁场的m 阶垂向导数、n 阶沿x 方向水平导数、l 阶沿y 方向的导数的
(n +l +m )Z a (x , y , z )的频谱) ，应乘上的导数因子为 频谱(即求∂
∂x n y l ∂z m
(2πiu )n (2πiv )l ⎡2πu 2+v 2⎢⎣()1/2⎤ (12) ⎥⎦m
航磁∆T 垂向一阶导数已经广泛地应用于磁异常的解释，它能区分相邻磁性体异常，减少其相互叠加的影响，并把叠加在背景场中的局部异常分离出来，是压制区域场，圈定局部异常，分离叠加异常的常用方法。在实际磁场转换处理中，由于垂向一次导数相当于高通滤波器，在突出高频异常的同时，也突出了测量、磁场调平等干扰误差。对本区化极场的数据进行压制干扰垂向一阶导数处理，处理后的图件与原磁场图相比（图2a 、b ），突出了浅部磁性体信息，而压制了深层区域场的影响。该处理也消除或减弱了局部异常之间的叠加和干扰现象。因此，航磁∆T 化极垂向一阶导数处理在提取强背景场中的弱缓异常，圈定局部异常、火山
构造、划分构造边界等方面具有重要作用。
根据厚板状磁性体异常公式，垂向二阶导数的零值线为磁性体边界位置。因此，航磁∆T 化极垂向二阶导数处理的主要目的是利用航磁异常垂向二次变换率来圈定磁性体的范围和边界。本区航磁∆T 垂向二阶导数处理是在化极处理的基础上，对化极后的网格数据采用频率域位场转换方法求取磁异常沿垂直方向上的二次变换率，并编制了航磁∆T 垂向二阶导数等值线平面图（图2c ）。在理论上，经垂向二阶导数处理后，区域场得到了进一步的压制，很大程度上消除了深部磁性体的影响，使得磁性体的范围和边界更加明显，仅供参考使用。
图2 航磁∆T 化极垂向一阶导数处理效果对比图
a-航磁∆T 等值线平面图；b-航磁∆T 化极垂向一阶导数等值线平面图；
c-航磁∆T 化极垂向二阶导数等值线平面图；
三、化极0°方向水平一阶导数处理
化极0°方向水平一阶导数处理的目的是突出异常在东西向的线性特征，分辨东西方向上构造线的展布，以准确的划定浅层构造、断裂构造，以便推断区内的构造格架。
磁异常水平导数换算公式如下：
如果令S zx (x , y , z )、S zy (x , y , z )、S zz (x , y , z )及S zxx (x , y , z )、S zyy (x , y , z )、S zzz (x , y , z )分别为Z z (x , y , z )对x 、y 、z 的一阶导数及二阶导数的频谱，则有微分定力易于得到：
(S zx (u , v , z )=2πiuS z (u , v ,0)e 2πu 2+v 2)1/ z 2
(S zy (u , v , z )=2πiuS z (u , v ,0)e
S zz 2πu 2+v 2)1/2 z 2πu 2+v 2(u , v , z )=2π(u 2+v 2)1/2(S z (u , v ,0)e )1/2 z
同理，可以写出：
⎧2⎪(2πiu )⎧S zxx ⎫⎪⎪⎪2⎪⎨S zyy ⎬(u , v , z )=⎨(2πiv )⎪⎪⎪⎩S zzz ⎭⎪⎡2πu 2+v 2⎢⎪⎣⎩⎫⎪ ⎪2π(u 2+v 2)1/2 z ⎪⎬S z (u , v ,0)e ⎪1/22⎤⎪⎥⎭⎪⎦
n ()由此可知，求磁场的n 阶垂向导数的频谱，应乘上的导数因子为⎡2π(u 2+v 2)1/2⎤；⎢⎥⎣⎦n 而求磁场沿x 方向或y 方向的n 阶水平导数的频谱，应乘上导数因子为(2πiu )或
n (2πiv )。
设l 是实测平面上任一方向，它与x 轴的夹角为α，则有：
∂S T (x , y , z )∂S T (x , y , z ) ∂S T (
x , y , z )
∂l =cos a ∂x +sin a ∂y
两边作傅氏变换并应用微分定理，得知
S n (u , v , z )=i (π2u c o a +s πv 2a s S T n (u v z ) , , )i (13) 利用(13)式即可实现磁场的频率域方向导数计算，当a =0，代入(13)式即可求得0°方向水平导数。
Sn(u,v,z)=2πiu Sr(u,v,z) Sn (u , v , z )=2πiu Sr (u , v , z ) (14)
航磁∆T 化极0°方向水平导数处理结果显示（图3），局部域近东西向的线性异常特征及弧形异常特征都非常明显，为该区划分浅层构造、近东西向断裂构造等提供依据。
图3 航磁∆T 化极0°方向水平向导数处理效果对比图
a-航磁∆T 化极等值线平面图；c-航磁∆T 化极0°方向水平导数等值线平面图
四、向上延拓处理
磁场向上延拓处理就是将原观测面上的磁场值向上换算到另一个高度面上。随着上延高度的增加，磁性体引起的异常幅度按指数规律衰减。衰减最快的为浅部局部磁性体引起的高频异常成分，而具有一定延伸的大规模磁性体引起的低频异常成分衰减较慢。可见，向上延拓处理起到压制浅部小规模磁性体异常而突出深大地质体异常的作用。
设场源位于z=H平面一下（H>0），则磁场在z=H平面以上是对x 、y 、z 的连续函数。若z=0观测平面上的磁场T (x , y ,0)为已知，可以得到向上延拓公式为
T (x , y , z )=T (ξ, η,0) -z ∞∞d ξd η2π⎰-∞⎰-∞⎡x -ξ2+y -η2+z 2⎤3/2)()⎣(⎦
1-z ⋅2πx 2+y 2+z 2(14) 由褶积积分公式可知，上式为T (x , y ,0)与()3/2关于变量(x , y )二维褶积。
空间域的褶积与频率域的乘积相对应。下面分别求T (x , y ,0)及-z
2πx 2+y 2+z 2的傅立叶变换，设T (x , y , z )对于变量(x , y )的傅立叶变换为S T (u , v , z ), 有
S T (u , v , z )=∞∞T (x , y , z )e -2πi (ux +vy )dxdy (15) ⎰-∞⎰-∞
则
S T (u , v , 0)=⎰∞
-∞-∞⎰∞-2πi (u +T e x (x , y ), 0)v y d y (16) d x
利用上式可以由已知的T (x , y ,0)求出其频谱S T (u , v ,0)。进一步求
傅立叶变换，应用Erdelyi （1954）给出的积分变换表可以得到：
-z 2πx 2+y 2+z 2的()3/2⎰⎰∞∞-z 2πx 2+y 2+z 2-∞-∞-2πi (ux +vy )dxdy =e 2πu 2+v 21/2 z () (17)
当z
上式对于z≤0成立。
T (x , y ,0)是S T (u , v , z )的反傅立叶变换，即
S T (u , v , z )=S T (u , v ,0)e 2πu 2+v 21/2 z () (18)
T (x , y , z )=⎰∞⎰∞S T (x , y ,0)e 2π(u +v )
-∞-∞2221/ z e 2πi (ux +vy )dudv (19)
(19)式即为向上延拓的频谱表达式。
通过磁场向上换算，相当于加大了观测面与场源的距离，可以使局部小规模异常随换算高度的增加而减小，而深部规模较大的磁性体所产生的异常更加凸出。为了了解深源磁性体的特征及航磁异常随高度衰减变化特征，判断磁性体的埋深及延伸情况，在化极基础上进行了0.5km 、1.0km 、3.0km 、5.0km 四种不同高度的向上延拓处理。通过不同高度的向上延拓，消除了高频磁异常的干扰，使得磁场面貌逐渐单调，达到了突出低频区域异常的目的，对了解深源磁性体的特征和基底构造具有一定的地质意义。
图4 航磁∆T 化极上延处理效果对比图
a-航磁∆T 化极上延0.5km 等值线平面图；b-航磁∆T 化极上延1.0km 等值线平面图；
c-航磁∆T 化极上延3.0km 等值线平面图；b-航磁∆T 化极上延5.0km 等值线平面图；
对比测区上延高度磁场图可以看出（图4）：化极上延0.5km 后，高频干扰异常被压制，有意义的局部异常基本保留；化极上延1.0km 后，规模较小的局部异常衰减得很快，中等规模的异常明显突出，区域磁场面貌反映得更加清晰；化极上延3km 后，由于测区覆盖较浅，而引起局部异常的磁性地质体延伸有限，高频异常几乎消失；化极上延5km 后，有效的压制了浅部磁性体引起的异常，突出了深源低频磁异常。因此，化极上延0.5km 或1.0km 磁场图，对研究本区区域构造、划分隐伏岩体非常有效；上延3.0km 后仍然存在的磁异常则反映出了规模较大、延伸较深或埋深较大的磁性地质体；上延5.0km 后反映的深源低频磁异常，对于确定磁性基岩、深大断裂及区域构造格架有着重要意义。
五、化极匹配滤波求取局部及区域场
区域场和局部场的分离问题是航磁数据处理的主要内容之一，对实际资料的解释有重要意义。利用匹配滤波算法可对航磁数据进行区域场和局部场的分离，进一步达到突出浅部异常或突出深部异常的目的。
我们可以假设局部异常是由许多下延较小的磁形体引起的，由场的等效原理可知，这类磁性体可以用随机分布的偶极源组成的等效层代替。设偶极源等效层的深度为d 2，源的偶极磁距为mdp (ξ, η, d 2)，其傅氏正变换为mdp (u , v , d 2)；并设观测面高度为z ，则以求出偶极等效层所产生的场的振幅谱为（为简洁式中略去u 0/4π，并不影响结果）
A dp (
u , v , z )=4πdp (u , v , d 2)e -2π(
d 2-z dp (u , v , d 2)e -2π(
d 2-z (20) 显然，振幅A dp (u , v , d
2)
在假定区域场是由许多下延伸很大的磁性体所引起，它们可以用由随机分布的电极等效层所代替。设点极等效层深度为d 1，电极源的磁荷密度为m p (ξ, η, d 1)，其傅氏正变换为m p (u , v , d 1)，则由点极等效层引起场的振幅谱为
A p (u , v , z )=2πM p (u , v , d 1)e -2π(
d 1- (21)
当不考虑干扰场时，设实测场为区域场与局部场之和，它的振幅谱为
A (u , v , z )=2πM p (u , v , d 1)e -2π(
d 1-z 4πM dp (u , v , d 2)e -2π(
d 2-z (22)
为了从实测场中提取区域场，可以设计一个只对区域场响应的滤波器，即与区域场匹配的滤波器。令滤波器的频率响应函数为）Φ(u , v , z ), 为使它仅对区域场响应，即要：
A (u , v , )z ⋅Φ(u , , v =)z p z (23) (A , u ), v
将(21)、 (22)式代入 (23)式可得：
Φ(u , v , z )
=p ⎢⎣ (24) ⎥⎦
实际上此时Φ(u , v , z )与z 无关，应为Φ(u , v )。求出频率相应Φ(u , v )后，由(23)式求出A p (u , v , z ), 在利用无相移的傅氏反变换（即相位谱不变），就获得了区域场。
同样，如果从实测场中提取局部场，也可以设计个只对局部异常匹配的滤波器，用同样的方法也可以得到它的频率响应函数为
Φ(
u , v )=1
⎡M ⎢1-2⎢⎣2-d 1) (25) ⎤⎥⎥⎦
利用匹配滤波分离场时, 其关键是是否具有良好的径向对数能谱曲线, 因为滤波因子是通过径向对数能谱曲线获得的。径向平均对数功率谱经过数据拟合得到的直线斜率与点距间的关系式:
h =-k ⋅∆x 4π (25)
式中, ∆x 表示点距, k 表示直线斜率, h 为区域场与局部场分离的最佳深度。
航磁数据处理中，匹配滤波法是分离局部与区域异常场的重要手段之一，因为匹配滤波器是已知相关滤波器，它要求二者具有明显差异的波数成分，该工区的航磁异常特征满足这个条件（图5），浅部大多为带状、块状或零散分布的火山岩、侵入岩体，提取区域场时，它是一个低通滤波器，这与数学解析向上延拓不同之处在于它有一个较为复杂的类似于汉宁滤波器的窗函数。在提取高频成分时，它不会放大导致高频成分的振荡效应，因为高通时的滤波器渐近线为1。与向上延拓相比，该方法简单易行，而且还能够获取分析局部场与区域场的相关窗口。通过匹配滤波算法可以分离浅部异常信息，通过不同的匹配因子，可以逐层剥离异常，为地下异常区域的分层解析提供了技术手段。
图5 航磁∆T 化极匹配滤波处理效果对比图
a-航磁∆T 化极匹配滤波浅源异常提取图；b-航磁∆T 化极匹配滤波深源异常提取图；
六、剩余异常剖面平面图
对于那些频率高、幅值低的航磁异常，在航磁∆T 等值线平面图上受网格化取数圆滑滤波影响以及成图精度的制约，往往显示不清或被漏掉；而在航磁∆T 平面剖面图上往往叠加在较大异常的背景场上，不容易识别出来。为了从剖面上将一些局部异常和较微弱的短波异常从区域背景中分离提取出来，采用空间域非线性滤波法对测区剖面数据进行非线性滤波处理计算，提取航磁剩余异常。
非线性滤波方法是根据剖面异常曲线的曲率变化，按一定的异常宽度窗口拟合计算出剖面的区域背景异常值，然后将其从原始剖面异常中减去，所获得的剖面
异常即为相应的剩余异常。该方法可以有效的提取指定波长的剩余磁异常，具体计算公式如下和迭代步骤如下（图6）：
步骤1 用剖面实测异常作为第一次圆滑的原始数据，用(26)、(27)、(28)式估计给定滤波窗口宽度WD 两倍内异常曲线弯曲方向的特征值。
(26) S k 1=T F K m 1-0.5 ⎡⎣T F (K m )⎤⎦+T F (K )
()
S k 2=T F (K )-0.5 ⎡T F K m 1⎤+T F K p 1 (27)
⎣⎦
()()
(28) S k 3=T F K p 1-0.5 ⎡⎣T F (K )⎤⎦+T F K p 2()()
步骤2 如果S k1、S k2、S k3同时满足下面条件：
图6 非线性滤波方法原理示意图
S k 3*S k 2
表明K m1、K 、K p1三点之间，异常曲线弯曲方向相同，拥有同一个区域异常。
由(29)、(30)、(31)式计算出K m1、K 、K p1三点区域异常值，然后向前移动窗口，回到步骤1。
T F K m 1=T F K m 1-S k 1.0.5 (29)
()()
T F (K )=T F (K )-S k 2.0.5 (30)
T F K p 1=T F K p 1-S k 3.0.5 (31)
()()
步骤3 如果S k1、S k2、S k3不同时买足步骤2的条件，则向前移动窗口，步长等于滤波窗口宽度，回到步骤1，指导整条剖面计算完成为止。
步骤4 以点距为单位，逐步减少窗口宽度，用新得到剖面异常代替原始异常，重复步骤1、2、3，共重复W D 次，使滤波窗口内每一个异常点都得到处理，提高滤波效果。
步骤5 为保证计算精度，重复上述步骤1、2、3、4过程2~5次，获得剖面区域异常。
步骤6 用原始剖面异常值减去由上面过程获得的区域异常值，得到剖面剩余异
常值。
为了获得较好的滤波效果，根据本区地质解释的需要，并分析对比不同滤波异常窗口处理结果信息，经过使用对比，本区最后选用滤波异常半宽度≥1000m ，沿剖面逐线完成对航磁弱缓异常的提取。从航磁弱缓异常处理效果看，消除了背景场的影响，分离了叠加异常，局部异常的形态更加完整清晰，并对弱小异常起到增强作用。如图7中，经过弱缓异常处理后，弱小磁异常在数条测线上反映明显，异常边界也更加清晰。
图7航磁、剩余异常综合图
a-航磁∆T 剖面平面图；b-航磁∆T 剩余异常剖面平面图
第二节 2.5D 人机交互解释
目前对起伏地形条件下航磁数据进行定量正反演计算主要采用最优化算法，在原起伏测量面上直接进行解释。本次研究以最优化理论为基础，使用二维半模型（2.5D ）来模拟二度和三度地质体。
为了提供模型正演计算的速度，很多专家对基于二度半多边形截面楞柱体重磁正反演公式做了推导简化工作。该区采用了有直接积分推导出的公式（姚长利，1998）。如图8所示，一个多边形界面2.5D 棱柱体，设其密度为α，磁化强度为M ，则在任意点P （x,y,z ）引起的磁场三分量为：
图8 2.5D 棱柱体模型及坐标系统
H ax (P )=-∑sin ϕi (M x I li +M y I 2i +M z I 3i ) (32)
i =1
N
(33) H ay (P )=-∑⎡⎣(M x sin ϕi -M z cos ϕi )I 2i -M y (sin ϕi I 1i -cos ϕi I 3i )⎤⎦
i =1
N
Z a (P )=∑sin ϕi (M x I 1i +M y I 2i +M z I 3i ) (34)
i =1
N
其中， I ji =P ji (Y 2)-P ji (Y 1), j =1,2,3。
P 1i (y )=cos ϕi ln
⎛R i +y u y u y ⎫
-sin ϕi arctan i +1-arctan i ⎪ R i +1+y w i R i +1w i R i ⎭⎝P 2i (y )=ln
u i +R i
u i +1+R i +1
P 3i (y )=sin ϕi ln
⎛R i +y u y u y
+cos ϕi arctan i +1-arctan i
R i +1+y w i R i +1w i R i ⎝
⎫
⎪⎭
总场异常：
其中：
u i =x i cos ϕi +Z i sin ϕi , u i +1=x i +1cos ϕi +Z i +1sin ϕi R i =x i 2+y 2+z i 2
∆T (P P c o s 0I c o s 0D +)=H a x ()
a x
H ()P c o i +n D (a )Z 0s I s 0
s i n I P 0
(35)
()
1/2
, R i +1=x i 2+1+y 2+z i 2+1
()
1/2
ϕi =arctan
z i +1-z i
, w i =-x i sin ϕi +z i cos ϕi
x i +1-x i
M x =M cos I cos D ，M y =M cos I sin D , M x =M sin I
以上各式中：i 为棱柱体角点序号，N 为棱柱体的边数，I 0、D 0为地磁场倾角、偏角，I 、D 为磁化强度方向的倾角、偏角。
依据上述公式，中国国土资源航空物探遥感中心自主研发的高精度重磁剖面解释系统（GMVPS ），软件中加入高程参量以适合各种起伏地形条件，结合最新计算机技术、自动反演技术及多种异常转换技术，实现了复杂地形条件下磁异常的定量解释（图9）。
该技术具有下列明显的技术优势：用一系列有限长的组合多边形棱柱体逼近三度体进行正反演，较以往用多边形二度体拟合剖面磁异常更接近实用；用起伏面作为异常起算面，较以往水平面更适合复杂地形地区；集成了模型可视化编辑、剖面数据处理、磁化强度和模型实时正反演等功能。这里的“可视化”是指磁性地质体模型在计算机屏幕上始终以图形或图象实体出现，并可直接对其进行修改、反演等操作，其形态、物性变化过程也是图形化的；“实时”是指模型的变化与其引起的异常曲线的响应变化几乎是同步的。
图9高精度重磁剖面解释系统（GMVPS ）界面
航磁数据位场转换处理及效果
航磁∆T 测量数据是不同深度、不同形态、规模的磁性地质体磁场信息在观测面上的综合反映。由于场的叠加效应，使得某些具有一定地质意义的异常变得复杂，在原始图件上很难识别，给地质解释工作带来了难度。为了提高对航磁异常的分辨能力，突出更多有用信息，根据测区航磁异常特征和地质解释需要，对原始测量数据进行了原平面化极、上延、垂向一阶导数以及剩余异常提取等几种位场转换处理。
第一节 位场转换处理及效果
航磁平面网格数据位场转换处理采用表达式简单、运算速度快捷的频率域算法，进行化极、导数换算、解析延拓等处理。频率域转换的过程是：首先对异常资料进行傅立叶正变换，以得到异常资料的频谱；而后把异常的频谱和与转换相应的频率相应函数点积，得到处理后异常的频谱；最后对处理后异常的频谱进行傅立叶反变换，从而得到处理后的异常。
位场转换处理使用的软件是中国国土资源航空物探遥感中心自主开发的
WINDOWS 系统下地球物理数据处理解释软件（GeoProbe Mager）及航空物探彩色矢量成图系统（AgsMGis ）。
一、原平面化极处理
化极，即化磁极，就是把斜磁化异常转变为垂直磁化异常，相当于在磁北极观测异常。测区处于中纬度地区，由于倾斜磁化的影响，造成磁异常中心不是正好对应在地质体的正上方，而是相对于地质体的中心向南部产生一定的偏移。这对于确定磁性地质体的空间位置、形态、分布范围以及对磁异常的定性定量解释均带来一定的困难。化极可用于消除由于非垂直磁化引起的异常不对称性，在剩磁很小或感磁远大于剩磁且两者方向一致的情况下，将实测的斜磁化异常转化为垂直磁化异常，这样可以较为准确的确定异常的场源位置，提高异常解释的定位精度。从而使异常形态简化，并与磁性体位置保持一致，有利于圈定磁性体边界和走向。
作化极处理时要注意剩磁的影响，化极处理一般都假定磁化方向与地磁场方向一致，对于那些剩磁远远大于感磁且剩磁方向与地磁场方向不一致的磁性体就不符合这一假设条件，特别是测区中的火山岩分布区，由于剩磁较大会出现磁场畸变现象，使用时应注意甄别。从项目组野外物性测量结果看，区内多数岩石以感磁为主，剩磁方向与感磁方向接近，符合化极的前提条件。
全区采用" 频率域偶层位变倾角磁方向转换方法" 实现磁场全变倾角化极。在观测面上建立笛卡尔直角坐标系，使x 轴志向磁北，z 轴垂直向下。假设观测场∆T 是
一分布在观测面下方z=h平面上的偶层磁荷面引起的。它在观测点P(x,y,z)处产生的磁位U 与磁场T 分别为
U =-1⎛1⎫⎰⎰M ∆ ⎪d ξd η (1) 4π⎝r ⎭
∆T =- (2) u t ∇U 0
(3)
设
⎛∂∂∂⎫∆= , , ⎪ ⎝∂x ∂y ∂z ⎭
t =t 0+∆t 0 (4)
l =l 0+∆l 0 (5)
式中M 表示偶层磁荷面的磁化强度矢量； t 0、l 0表示均为常矢量研究区t 、l 的平
均值； ∆t 0、∆l 0表示他们的变化值。
将(4)式代入(2)式并进行傅氏变换得：
F [∆T ]=-u g 0)t 0{( [F ]U +[ F ∆ t ∇]} U 0 (6)
移项得
F [U ]=-
式中
⎫ 1⎧1∇U ]⎬⎨F [∆T ]+F [∆t 0 g t 0⎩u 0⎭ (7) g =(i 2πu , i 2πv , 2πf ), f =对(1)式两端作傅氏变换得：
F [U ]=-1⎡g ⋅F [M ]⋅4π⎣-2πf (h -)z e /⎤⎦f (8) 同样的，将(5)式代入(3)式进行傅氏变换后代入(8)式得：
F [U ]=-1{(g ⋅l )⋅F [M ]+g ⋅F [M ⋅∆l ]}⋅e 4πf 00-2πf (h -z )
移项得
F [M ]=-14πf ⋅e 2πf (h -z )F [U ]+g ⋅F [M ⋅∆l 0] g ⋅l 0{}(9)
（7）式、(9)式即为频率域变倾角磁方向转换的两个基本公式。已知观测场∆T ，可应用（7）式计算F [U ]，再将F [U ] F [U ]代入(9)式中计算F [M ]。再把F [M ]代入地磁极处的F [T ]极，即可实现变倾角化极处理。F [T ]极与F [M ]的关系如下。
在地磁北极有：t 0=l 0=(0,0,1), ∆t 0=l 0=0，由(7) 式、(9)式可得：
F [M ]极=π
-2πf (h -u 0f e )z (10) [F ]M
经傅立叶反变换后可得化磁极长Z ⊥
πf (h -)z Z ⊥=F -1{πu (11) f -2e M [F ]}0
假设要换算的场为T p ，其磁场方向单位矢为t ，磁化方向单位矢为l ，则只要把
他们代入(7) ，令∆t 0、∆l 0偏差为零，即得F [U ]的初值F 0[U ]：
F [U ]=-1
1 F ⎡T p ⎤⎣⎦u 0g ⋅t 0 对F 0[U ]反变换求得初值U 0，把U 0代入(7)式得F 1[U ]，如此反复迭代，直到求得的U 值之差小于给定的标准为止。求得F [U ]后，类似地取F 0[U ]为
F [M ]=--4πf F [U ]⋅e 2πf (h -z ) 0g ⋅l 0
反变换求得M 0，把M 0代入(9)式，求得F 1[M ]，如此反复可最后求得F [M ]，并代入(11)式求得化磁极磁场。
由于本区处于xxx°xxx ′~xxx°xxx ′，属于中纬度地区，斜磁化能够产生一定影响，对原磁场数据进行化极处理后，在垂直磁化的条件下，磁异常的形态以及磁异常与磁性体的关系都比较简单，便于进行地质解释。对比航磁∆T 等值线平面图和航磁∆T 化极等值线平面图，航磁化极处理作用非常明显（图1）：局部异常整体向北偏移，表明通过化极处理，使异常回归到磁性地质体上方；减小或消除了由于斜磁化而引起的多数局部异常正负异常伴生现象，为进一步圈定岩体边界创造条件；使异常带及梯度带更加明显，有利于揭示出不同地质体的分布与形态，对圈定各种不同类型的断裂、确定磁性体的性质及边界具有重要的意义。
图1 航磁∆T 化极处理效果对比图
a-航磁∆T 等值线平面图；b-航磁∆T 化极等值线平面图
二、化极垂向导数处理
航磁局部异常通常是叠加在区域背景场上的次级异常，在原始航磁或化极航磁等基础图件中表现并不明显，需要通过一定的数学处理手段来突出其特征。垂向导数处理是解决这个问题的一种有效手段，它反映了磁场在垂直方向上的梯度变化，在增强由浅部磁性体引起的局部异常、压制长波区域场有很强的功能，可以突出在总场图上不明显的细节，并能分解横向叠加异常，理论上导数的次数越高，这种分辨能力就越强。
磁异常垂向导数换算公式如下：
如果令S zx (x , y , z )、S zy (x , y , z )、S zz (x , y , z )及S zzz (x , y , z )、S zyy (x , y , z )、S zzz (x , y , z )分别为Z a (x , y , z )对x 、y 、z 的一阶导数及二阶导数的频谱，则有微分定力易于得到：
(S zx (u , v , z )=2πiuS z (u , v ,0)e 2πu 2+v 2)1/ z 2
(S zy (u , v , z )=2πiuS z (u , v ,0)e
S 2πu 2+v 2)1/2 z 2πu 2+v 2(u , v , z )=2πu 2+v 2
zz ()1/2S z (u , v ,0)e ()1/2 z
同理，可以写出：
⎧2⎪(2πiu )⎧S zxx ⎫⎪⎪⎪2⎪S u , v , z =)⎨(2πiv )⎨zyy ⎬(⎪⎪⎪⎩S zzz ⎭⎪⎡2πu 2+v 2
⎪⎢⎣⎩⎫⎪ ⎪2π(u 2+v 2)1/2 z ⎪S u , v ,0e )⎬z (⎪1/22⎤⎪⎥⎪⎦⎭
n ()由此可知，求磁场的n 阶垂向导数的频谱，应乘上的导数因子为⎡2π(u 2+v 2)1/2⎤；⎢⎣⎥⎦
而求磁场沿x 方向或y 方向的n 阶水平导数的频谱，应乘上导数因子为(2πiu )n 或
n (2πiv )。
如果求磁场的m 阶垂向导数、n 阶沿x 方向水平导数、l 阶沿y 方向的导数的
(n +l +m )Z a (x , y , z )的频谱) ，应乘上的导数因子为 频谱(即求∂
∂x n y l ∂z m
(2πiu )n (2πiv )l ⎡2πu 2+v 2⎢⎣()1/2⎤ (12) ⎥⎦m
航磁∆T 垂向一阶导数已经广泛地应用于磁异常的解释，它能区分相邻磁性体异常，减少其相互叠加的影响，并把叠加在背景场中的局部异常分离出来，是压制区域场，圈定局部异常，分离叠加异常的常用方法。在实际磁场转换处理中，由于垂向一次导数相当于高通滤波器，在突出高频异常的同时，也突出了测量、磁场调平等干扰误差。对本区化极场的数据进行压制干扰垂向一阶导数处理，处理后的图件与原磁场图相比（图2a 、b ），突出了浅部磁性体信息，而压制了深层区域场的影响。该处理也消除或减弱了局部异常之间的叠加和干扰现象。因此，航磁∆T 化极垂向一阶导数处理在提取强背景场中的弱缓异常，圈定局部异常、火山
构造、划分构造边界等方面具有重要作用。
根据厚板状磁性体异常公式，垂向二阶导数的零值线为磁性体边界位置。因此，航磁∆T 化极垂向二阶导数处理的主要目的是利用航磁异常垂向二次变换率来圈定磁性体的范围和边界。本区航磁∆T 垂向二阶导数处理是在化极处理的基础上，对化极后的网格数据采用频率域位场转换方法求取磁异常沿垂直方向上的二次变换率，并编制了航磁∆T 垂向二阶导数等值线平面图（图2c ）。在理论上，经垂向二阶导数处理后，区域场得到了进一步的压制，很大程度上消除了深部磁性体的影响，使得磁性体的范围和边界更加明显，仅供参考使用。
图2 航磁∆T 化极垂向一阶导数处理效果对比图
a-航磁∆T 等值线平面图；b-航磁∆T 化极垂向一阶导数等值线平面图；
c-航磁∆T 化极垂向二阶导数等值线平面图；
三、化极0°方向水平一阶导数处理
化极0°方向水平一阶导数处理的目的是突出异常在东西向的线性特征，分辨东西方向上构造线的展布，以准确的划定浅层构造、断裂构造，以便推断区内的构造格架。
磁异常水平导数换算公式如下：
如果令S zx (x , y , z )、S zy (x , y , z )、S zz (x , y , z )及S zxx (x , y , z )、S zyy (x , y , z )、S zzz (x , y , z )分别为Z z (x , y , z )对x 、y 、z 的一阶导数及二阶导数的频谱，则有微分定力易于得到：
(S zx (u , v , z )=2πiuS z (u , v ,0)e 2πu 2+v 2)1/ z 2
(S zy (u , v , z )=2πiuS z (u , v ,0)e
S zz 2πu 2+v 2)1/2 z 2πu 2+v 2(u , v , z )=2π(u 2+v 2)1/2(S z (u , v ,0)e )1/2 z
同理，可以写出：
⎧2⎪(2πiu )⎧S zxx ⎫⎪⎪⎪2⎪⎨S zyy ⎬(u , v , z )=⎨(2πiv )⎪⎪⎪⎩S zzz ⎭⎪⎡2πu 2+v 2⎢⎪⎣⎩⎫⎪ ⎪2π(u 2+v 2)1/2 z ⎪⎬S z (u , v ,0)e ⎪1/22⎤⎪⎥⎭⎪⎦
n ()由此可知，求磁场的n 阶垂向导数的频谱，应乘上的导数因子为⎡2π(u 2+v 2)1/2⎤；⎢⎥⎣⎦n 而求磁场沿x 方向或y 方向的n 阶水平导数的频谱，应乘上导数因子为(2πiu )或
n (2πiv )。
设l 是实测平面上任一方向，它与x 轴的夹角为α，则有：
∂S T (x , y , z )∂S T (x , y , z ) ∂S T (
x , y , z )
∂l =cos a ∂x +sin a ∂y
两边作傅氏变换并应用微分定理，得知
S n (u , v , z )=i (π2u c o a +s πv 2a s S T n (u v z ) , , )i (13) 利用(13)式即可实现磁场的频率域方向导数计算，当a =0，代入(13)式即可求得0°方向水平导数。
Sn(u,v,z)=2πiu Sr(u,v,z) Sn (u , v , z )=2πiu Sr (u , v , z ) (14)
航磁∆T 化极0°方向水平导数处理结果显示（图3），局部域近东西向的线性异常特征及弧形异常特征都非常明显，为该区划分浅层构造、近东西向断裂构造等提供依据。
图3 航磁∆T 化极0°方向水平向导数处理效果对比图
a-航磁∆T 化极等值线平面图；c-航磁∆T 化极0°方向水平导数等值线平面图
四、向上延拓处理
磁场向上延拓处理就是将原观测面上的磁场值向上换算到另一个高度面上。随着上延高度的增加，磁性体引起的异常幅度按指数规律衰减。衰减最快的为浅部局部磁性体引起的高频异常成分，而具有一定延伸的大规模磁性体引起的低频异常成分衰减较慢。可见，向上延拓处理起到压制浅部小规模磁性体异常而突出深大地质体异常的作用。
设场源位于z=H平面一下（H>0），则磁场在z=H平面以上是对x 、y 、z 的连续函数。若z=0观测平面上的磁场T (x , y ,0)为已知，可以得到向上延拓公式为
T (x , y , z )=T (ξ, η,0) -z ∞∞d ξd η2π⎰-∞⎰-∞⎡x -ξ2+y -η2+z 2⎤3/2)()⎣(⎦
1-z ⋅2πx 2+y 2+z 2(14) 由褶积积分公式可知，上式为T (x , y ,0)与()3/2关于变量(x , y )二维褶积。
空间域的褶积与频率域的乘积相对应。下面分别求T (x , y ,0)及-z
2πx 2+y 2+z 2的傅立叶变换，设T (x , y , z )对于变量(x , y )的傅立叶变换为S T (u , v , z ), 有
S T (u , v , z )=∞∞T (x , y , z )e -2πi (ux +vy )dxdy (15) ⎰-∞⎰-∞
则
S T (u , v , 0)=⎰∞
-∞-∞⎰∞-2πi (u +T e x (x , y ), 0)v y d y (16) d x
利用上式可以由已知的T (x , y ,0)求出其频谱S T (u , v ,0)。进一步求
傅立叶变换，应用Erdelyi （1954）给出的积分变换表可以得到：
-z 2πx 2+y 2+z 2的()3/2⎰⎰∞∞-z 2πx 2+y 2+z 2-∞-∞-2πi (ux +vy )dxdy =e 2πu 2+v 21/2 z () (17)
当z
上式对于z≤0成立。
T (x , y ,0)是S T (u , v , z )的反傅立叶变换，即
S T (u , v , z )=S T (u , v ,0)e 2πu 2+v 21/2 z () (18)
T (x , y , z )=⎰∞⎰∞S T (x , y ,0)e 2π(u +v )
-∞-∞2221/ z e 2πi (ux +vy )dudv (19)
(19)式即为向上延拓的频谱表达式。
通过磁场向上换算，相当于加大了观测面与场源的距离，可以使局部小规模异常随换算高度的增加而减小，而深部规模较大的磁性体所产生的异常更加凸出。为了了解深源磁性体的特征及航磁异常随高度衰减变化特征，判断磁性体的埋深及延伸情况，在化极基础上进行了0.5km 、1.0km 、3.0km 、5.0km 四种不同高度的向上延拓处理。通过不同高度的向上延拓，消除了高频磁异常的干扰，使得磁场面貌逐渐单调，达到了突出低频区域异常的目的，对了解深源磁性体的特征和基底构造具有一定的地质意义。
图4 航磁∆T 化极上延处理效果对比图
a-航磁∆T 化极上延0.5km 等值线平面图；b-航磁∆T 化极上延1.0km 等值线平面图；
c-航磁∆T 化极上延3.0km 等值线平面图；b-航磁∆T 化极上延5.0km 等值线平面图；
对比测区上延高度磁场图可以看出（图4）：化极上延0.5km 后，高频干扰异常被压制，有意义的局部异常基本保留；化极上延1.0km 后，规模较小的局部异常衰减得很快，中等规模的异常明显突出，区域磁场面貌反映得更加清晰；化极上延3km 后，由于测区覆盖较浅，而引起局部异常的磁性地质体延伸有限，高频异常几乎消失；化极上延5km 后，有效的压制了浅部磁性体引起的异常，突出了深源低频磁异常。因此，化极上延0.5km 或1.0km 磁场图，对研究本区区域构造、划分隐伏岩体非常有效；上延3.0km 后仍然存在的磁异常则反映出了规模较大、延伸较深或埋深较大的磁性地质体；上延5.0km 后反映的深源低频磁异常，对于确定磁性基岩、深大断裂及区域构造格架有着重要意义。
五、化极匹配滤波求取局部及区域场
区域场和局部场的分离问题是航磁数据处理的主要内容之一，对实际资料的解释有重要意义。利用匹配滤波算法可对航磁数据进行区域场和局部场的分离，进一步达到突出浅部异常或突出深部异常的目的。
我们可以假设局部异常是由许多下延较小的磁形体引起的，由场的等效原理可知，这类磁性体可以用随机分布的偶极源组成的等效层代替。设偶极源等效层的深度为d 2，源的偶极磁距为mdp (ξ, η, d 2)，其傅氏正变换为mdp (u , v , d 2)；并设观测面高度为z ，则以求出偶极等效层所产生的场的振幅谱为（为简洁式中略去u 0/4π，并不影响结果）
A dp (
u , v , z )=4πdp (u , v , d 2)e -2π(
d 2-z dp (u , v , d 2)e -2π(
d 2-z (20) 显然，振幅A dp (u , v , d
2)
在假定区域场是由许多下延伸很大的磁性体所引起，它们可以用由随机分布的电极等效层所代替。设点极等效层深度为d 1，电极源的磁荷密度为m p (ξ, η, d 1)，其傅氏正变换为m p (u , v , d 1)，则由点极等效层引起场的振幅谱为
A p (u , v , z )=2πM p (u , v , d 1)e -2π(
d 1- (21)
当不考虑干扰场时，设实测场为区域场与局部场之和，它的振幅谱为
A (u , v , z )=2πM p (u , v , d 1)e -2π(
d 1-z 4πM dp (u , v , d 2)e -2π(
d 2-z (22)
为了从实测场中提取区域场，可以设计一个只对区域场响应的滤波器，即与区域场匹配的滤波器。令滤波器的频率响应函数为）Φ(u , v , z ), 为使它仅对区域场响应，即要：
A (u , v , )z ⋅Φ(u , , v =)z p z (23) (A , u ), v
将(21)、 (22)式代入 (23)式可得：
Φ(u , v , z )
=p ⎢⎣ (24) ⎥⎦
实际上此时Φ(u , v , z )与z 无关，应为Φ(u , v )。求出频率相应Φ(u , v )后，由(23)式求出A p (u , v , z ), 在利用无相移的傅氏反变换（即相位谱不变），就获得了区域场。
同样，如果从实测场中提取局部场，也可以设计个只对局部异常匹配的滤波器，用同样的方法也可以得到它的频率响应函数为
Φ(
u , v )=1
⎡M ⎢1-2⎢⎣2-d 1) (25) ⎤⎥⎥⎦
利用匹配滤波分离场时, 其关键是是否具有良好的径向对数能谱曲线, 因为滤波因子是通过径向对数能谱曲线获得的。径向平均对数功率谱经过数据拟合得到的直线斜率与点距间的关系式:
h =-k ⋅∆x 4π (25)
式中, ∆x 表示点距, k 表示直线斜率, h 为区域场与局部场分离的最佳深度。
航磁数据处理中，匹配滤波法是分离局部与区域异常场的重要手段之一，因为匹配滤波器是已知相关滤波器，它要求二者具有明显差异的波数成分，该工区的航磁异常特征满足这个条件（图5），浅部大多为带状、块状或零散分布的火山岩、侵入岩体，提取区域场时，它是一个低通滤波器，这与数学解析向上延拓不同之处在于它有一个较为复杂的类似于汉宁滤波器的窗函数。在提取高频成分时，它不会放大导致高频成分的振荡效应，因为高通时的滤波器渐近线为1。与向上延拓相比，该方法简单易行，而且还能够获取分析局部场与区域场的相关窗口。通过匹配滤波算法可以分离浅部异常信息，通过不同的匹配因子，可以逐层剥离异常，为地下异常区域的分层解析提供了技术手段。
图5 航磁∆T 化极匹配滤波处理效果对比图
a-航磁∆T 化极匹配滤波浅源异常提取图；b-航磁∆T 化极匹配滤波深源异常提取图；
六、剩余异常剖面平面图
对于那些频率高、幅值低的航磁异常，在航磁∆T 等值线平面图上受网格化取数圆滑滤波影响以及成图精度的制约，往往显示不清或被漏掉；而在航磁∆T 平面剖面图上往往叠加在较大异常的背景场上，不容易识别出来。为了从剖面上将一些局部异常和较微弱的短波异常从区域背景中分离提取出来，采用空间域非线性滤波法对测区剖面数据进行非线性滤波处理计算，提取航磁剩余异常。
非线性滤波方法是根据剖面异常曲线的曲率变化，按一定的异常宽度窗口拟合计算出剖面的区域背景异常值，然后将其从原始剖面异常中减去，所获得的剖面
异常即为相应的剩余异常。该方法可以有效的提取指定波长的剩余磁异常，具体计算公式如下和迭代步骤如下（图6）：
步骤1 用剖面实测异常作为第一次圆滑的原始数据，用(26)、(27)、(28)式估计给定滤波窗口宽度WD 两倍内异常曲线弯曲方向的特征值。
(26) S k 1=T F K m 1-0.5 ⎡⎣T F (K m )⎤⎦+T F (K )
()
S k 2=T F (K )-0.5 ⎡T F K m 1⎤+T F K p 1 (27)
⎣⎦
()()
(28) S k 3=T F K p 1-0.5 ⎡⎣T F (K )⎤⎦+T F K p 2()()
步骤2 如果S k1、S k2、S k3同时满足下面条件：
图6 非线性滤波方法原理示意图
S k 3*S k 2
表明K m1、K 、K p1三点之间，异常曲线弯曲方向相同，拥有同一个区域异常。
由(29)、(30)、(31)式计算出K m1、K 、K p1三点区域异常值，然后向前移动窗口，回到步骤1。
T F K m 1=T F K m 1-S k 1.0.5 (29)
()()
T F (K )=T F (K )-S k 2.0.5 (30)
T F K p 1=T F K p 1-S k 3.0.5 (31)
()()
步骤3 如果S k1、S k2、S k3不同时买足步骤2的条件，则向前移动窗口，步长等于滤波窗口宽度，回到步骤1，指导整条剖面计算完成为止。
步骤4 以点距为单位，逐步减少窗口宽度，用新得到剖面异常代替原始异常，重复步骤1、2、3，共重复W D 次，使滤波窗口内每一个异常点都得到处理，提高滤波效果。
步骤5 为保证计算精度，重复上述步骤1、2、3、4过程2~5次，获得剖面区域异常。
步骤6 用原始剖面异常值减去由上面过程获得的区域异常值，得到剖面剩余异
常值。
为了获得较好的滤波效果，根据本区地质解释的需要，并分析对比不同滤波异常窗口处理结果信息，经过使用对比，本区最后选用滤波异常半宽度≥1000m ，沿剖面逐线完成对航磁弱缓异常的提取。从航磁弱缓异常处理效果看，消除了背景场的影响，分离了叠加异常，局部异常的形态更加完整清晰，并对弱小异常起到增强作用。如图7中，经过弱缓异常处理后，弱小磁异常在数条测线上反映明显，异常边界也更加清晰。
图7航磁、剩余异常综合图
a-航磁∆T 剖面平面图；b-航磁∆T 剩余异常剖面平面图
第二节 2.5D 人机交互解释
目前对起伏地形条件下航磁数据进行定量正反演计算主要采用最优化算法，在原起伏测量面上直接进行解释。本次研究以最优化理论为基础，使用二维半模型（2.5D ）来模拟二度和三度地质体。
为了提供模型正演计算的速度，很多专家对基于二度半多边形截面楞柱体重磁正反演公式做了推导简化工作。该区采用了有直接积分推导出的公式（姚长利，1998）。如图8所示，一个多边形界面2.5D 棱柱体，设其密度为α，磁化强度为M ，则在任意点P （x,y,z ）引起的磁场三分量为：
图8 2.5D 棱柱体模型及坐标系统
H ax (P )=-∑sin ϕi (M x I li +M y I 2i +M z I 3i ) (32)
i =1
N
(33) H ay (P )=-∑⎡⎣(M x sin ϕi -M z cos ϕi )I 2i -M y (sin ϕi I 1i -cos ϕi I 3i )⎤⎦
i =1
N
Z a (P )=∑sin ϕi (M x I 1i +M y I 2i +M z I 3i ) (34)
i =1
N
其中， I ji =P ji (Y 2)-P ji (Y 1), j =1,2,3。
P 1i (y )=cos ϕi ln
⎛R i +y u y u y ⎫
-sin ϕi arctan i +1-arctan i ⎪ R i +1+y w i R i +1w i R i ⎭⎝P 2i (y )=ln
u i +R i
u i +1+R i +1
P 3i (y )=sin ϕi ln
⎛R i +y u y u y
+cos ϕi arctan i +1-arctan i
R i +1+y w i R i +1w i R i ⎝
⎫
⎪⎭
总场异常：
其中：
u i =x i cos ϕi +Z i sin ϕi , u i +1=x i +1cos ϕi +Z i +1sin ϕi R i =x i 2+y 2+z i 2
∆T (P P c o s 0I c o s 0D +)=H a x ()
a x
H ()P c o i +n D (a )Z 0s I s 0
s i n I P 0
(35)
()
1/2
, R i +1=x i 2+1+y 2+z i 2+1
()
1/2
ϕi =arctan
z i +1-z i
, w i =-x i sin ϕi +z i cos ϕi
x i +1-x i
M x =M cos I cos D ，M y =M cos I sin D , M x =M sin I
以上各式中：i 为棱柱体角点序号，N 为棱柱体的边数，I 0、D 0为地磁场倾角、偏角，I 、D 为磁化强度方向的倾角、偏角。
依据上述公式，中国国土资源航空物探遥感中心自主研发的高精度重磁剖面解释系统（GMVPS ），软件中加入高程参量以适合各种起伏地形条件，结合最新计算机技术、自动反演技术及多种异常转换技术，实现了复杂地形条件下磁异常的定量解释（图9）。
该技术具有下列明显的技术优势：用一系列有限长的组合多边形棱柱体逼近三度体进行正反演，较以往用多边形二度体拟合剖面磁异常更接近实用；用起伏面作为异常起算面，较以往水平面更适合复杂地形地区；集成了模型可视化编辑、剖面数据处理、磁化强度和模型实时正反演等功能。这里的“可视化”是指磁性地质体模型在计算机屏幕上始终以图形或图象实体出现，并可直接对其进行修改、反演等操作，其形态、物性变化过程也是图形化的；“实时”是指模型的变化与其引起的异常曲线的响应变化几乎是同步的。
图9高精度重磁剖面解释系统（GMVPS ）界面