2017, Vol. 60
异常密度是地下地质体密度与围岩密度差，也称剩余密度；异常磁化强度可类似地定义，但它不能称为剩余磁化强度(因其与岩、矿石形成时的地磁场有关).深度d上方的地质体异常密度或异常磁化强度产生的位场异常(重力异常或磁异常)又称深度d上方的剩余异常，位场观测异常减去这种剩余异常称之深度d的区域异常，或称深度d的背景异常.通常分离出深度d的区域异常，将其从位场观测异常中减去，得到深度d上方的剩余异常.剩余异常也称局部异常.
要把任一深度d上方的剩余异常从叠加异常中分离出来，这是十分困难的，也是目前未能很好解决的难题(徐世浙等，2009).若能解决该问题，就能得到地下不同水平层的位场异常，进一步反演出各层的异常物性参数，即地下三维异常密度或异常磁化强度，这对位场资料的反演及解释有重要意义.
经几十年的研究积累，已产生多种位场异常分离方法.例如：圆周平均法(Griffin, 1949)、多项式拟合法(Agcos, 1951; Oldham and Sutherland, 1955；Rao et al., 1975；李庆春等，1994)、非线性滤波法(Naudy and Dreyer, 1968；Keating and Pinet, 2011)、匹配滤波法(Spector and Grant, 1970)、切割法(刘东甲，1984；程方道等，1987；文百红和程方道，1990)、延拓法(Jacobsen, 1987; Pawlowski, 1995; 许德树和曾华霖, 2000；Zeng et al., 2008; Meng et al., 2009)、优化滤波法(郭良辉等, 2012)、小波变换法(侯遵泽和杨文采，1997；Fedi and Quarta, 1998；杨文采等, 2001, 刘天佑等，2007；刘彩云等，2015).
以上方法，效果各异.根据模型和实际资料的试验，切割法是一种比较好的分离技术(徐世浙等，2009).在切割半径较小(如小于等于10倍点距)，且切割次数较多时，切割法分离效果较好，但切割半径较大时，异常畸变、假异常及高频干扰严重，切割法效果差(刘东甲等，20151)).比如，位场资料点距为50 m，要分离出深度2000 m上方的剩余异常，一般需要40倍点距的切割半径，切割法失效.
1) 刘东甲.2015.视密度和视磁化强度三维反演技术研究及其在深部找矿中的应用.安徽国土资源科技项目研究报告.
因此，要对位场资料进行定量反演和解释，需要创新位场异常分离方法.本文推导出一个新的八点圆周平均公式，由八点圆周平均公式的传递函数及其滤波特性，按递减半径的方式构造出一个用于分离区域异常的传递函数，该传递函数显著地改善了前者的滤波特性.由该传递函数给出空间域分离区域异常的递减半径线性迭代法，通过对该法进行改进，给出空间域分离区域异常的递减半径非线性迭代法，并用模型和实例进行试验.
2 一个新的八点圆周平均公式及递减半径迭代传递函数用圆周平均法对方格网位场(重、磁)数据进行分离，简单的四点圆周平均公式使分离的局部异常产生十字形畸变，因此，本节导出用来计算位场区域异常的一个新的八点圆周平均公式及其传递函数，研究八点圆周平均公式随圆周半径变化时的滤波特性；针对圆周平均法的不足，由该传递函数构造出递减半径迭代传递函数.
2.1 八点圆周平均公式设网格线距Δx等于点距Δy，x轴为纵轴，y轴为横轴。在以网格节点(x, y)为圆心、r为半径(暂且等于Δx)的圆周上取均布的八个点，其中仅四点为网格节点，另四个点的每点值由其邻近三节点插值，便导出由位场异常u(x, y)计算位场区域异常的一个新的八点圆周平均公式：
|
(1) |
式中系数
|
(2) |
公式(1) 被推广用于r等于网格距整倍数的情形，即圆周半径.
|
(3) |
则半径为r时的区域异常和局部异常分别为
|
(4) |
|
(5) |
利用δ函数，式(1)还可表示为褶积：
|
(6) |
其中
|
(7) |
用圆周平均公式计算位场区域异常，类似于数字图象处理领域(Gonzalez and Woods, 2008)的圆滑滤波(smoothing filter)，可把这个新的八点圆周平均公式看作以计算点(x, y)为中心、尺度为2r×2r邻域(或模板)的9个节点位场异常的加权平均，权系数可用矩阵表示为
|
(8) |
权系数之和等于1，即
|
(9) |
对式(1) 或式(6) 作Fourier变换，得
|
(10) |
其中U(kx, ky)为位场异常u(x, y)的波谱；U(kx, ky, r)为u(x, y, r)的波谱；八点圆周平均传递函数
|
(11) |
它是波数(kx, ky)的实偶函数，其参变量为圆周半径r=rj(式(3)).
对采样间隔为Δx=Δy的网格，折叠波数(又称Nyquist波数)为
|
(12) |
由微分学，H(kx, ky, rj)在波数平面的正方形区域{(kx, ky), |kx|≤kN，|ky|≤kN}存在极值点和鞍点，它们的坐标为
|
(13) |
l和m均为偶数时，有极大值H(kxl*, kym*, rj)= 4(b1+b2)=1；l和m均为奇数时，有极小值H(kxl*, kym*, rj)=-4(b1-b2)=－(2－
图 1(a—c)分别是Δx=Δy=100 m时由式(11) 计算绘制的3个半径(r1，r2，r3)对应的八点圆周平均传递函数，由图 1(a—c)及多个更大半径的计算结果可见H(kx, ky, rj)的滤波特性：(1) 因-(2-
|
图 1 传递函数 八点圆周平均(a), (b)和(c)；等半径迭代(d), (e)和(f)；递减半径迭代(g), (h)和(i). Fig. 1 Transfer function (a) H(kx, ky, r1); (b) H(kx, ky, r2); (c) H(kx, ky, r3); (d) Her (kx, ky, r1); (e) Her(kx, ky, r2); (f) Her (kx, ky, r3); (g) Hdr (kx, ky, r1); (h) Hdr (kx, ky, r2); (i) Hdr (kx, ky, r3). |
按等半径(equal radius)迭代方式构造r=rj时的一个传递函数，称之等半径迭代传递函数：
|
(14) |
其迭代次数即指数ner一般取6. 图 1(d—f)分别是3个半径(r1，r2，r3)对应的等半径迭代传递函数，其与图 1(a—c)相比，主瓣和旁瓣变窄，负的旁瓣得到压制，正的旁瓣依旧存在，即同八点圆周平均公式类似，计算出的区域异常中存在规则分布的假异常及高频干扰.
由于不同圆周半径对应的传递函数主瓣极值点位置相同而旁瓣极值点位置绝大多数不同，可按递减半径(decreasing radius)迭代方式尝试构造r=rj时的一个新的传递函数，称之递减半径迭代传递函数：
|
(15) |
它最后一项Hn1(kx, ky, r1)表示迭代到最小半径时等半径(r=r1)迭代n1次，目的在于利用H(kx, ky, r1)无正旁瓣(见图 1a)的优点，进一步改善Hdr(kx, ky, rj)的滤波特性.本文一般取n1=6，否则，将另给其值. 图 1g所示的Hdr(kx, ky, r1)由图 1a所示的传递函数自乘n1次，即Hn1(kx, ky, r1)；图 1h所示的Hdr(kx, ky, r2)由图 1b和图 1g所示的两个传递函数相乘，即H(kx, ky, r2)Hn1(kx, ky, r1)；图 1i所示的Hdr(kx, ky, r3)由图 1c、图 1b和图 1g所示的三个传递函数相乘，即H(kx, ky, r3)H(kx, ky, r2)Hn1(kx, ky, r1).
一系列的计算结果表明，Hdr(kx, ky, rj)显著地改善了H(kx, ky, rj)的滤波特性：(1)Hdr(kx, ky, rj)在零波数点等于1；(2)Hdr(kx, ky, rj)只有主瓣，几乎无旁瓣，类似一个低通滤波器；(3)rj越大，Hdr(kx, ky, rj)的主瓣宽度越小，即低通滤波器的截止波数越小；(4)Hdr(kx, ky, rj)的截止波数约等于π/(4rj)，与向上延拓传递函数exp(-rjk)的截止波数1/rj相近，这两种传递函数形状有差异(如图 2，ky=0, 
|
图 2 kx剖面传递函数 (a)递减半径迭代；(b)向上延拓. Fig. 2 Transfer functions at kx cross section (a) Decreasing radius linear iteration; (b) Upward continuation. |
式(15)Hdr(kx, ky, rj)的计算类似波数域迭代法(刘东甲等，2009；2012).由构造的Hdr(kx, ky, rj)，利用正、反Fourier变换，可按下式计算半径r=rj时位场的区域异常：
|
(16) |
如果直接利用构造的Hdr(kx, ky, rj)，由式(16) 分离位场的区域异常，其效果会显著优于式(1) 和(4) 的计算结果，但这属于波数域的级联滤波(cascade filtering)，仍然避免不了线性滤波的不足(见下文说明(4) 和(5))，需要把它变到空间域，以便进一步构造出更好的空间域迭代公式.
3.1 空间域递减半径线性迭代法取区域异常迭代初始值为位场异常值
|
(17) |
区域异常迭代公式
|
(18) |
其中系数Cn=1，Bn(x, y)为第n次迭代区域异常，Sn(x, y, r)为第n次迭代剩余异常，其算式为
|
(19) |
其中Bn(x, y, r)由式(1) 计算；r=rj－n=(j－n)Δx.
当迭代到第n=j－1次，即r=r1=Δx时，以r=r1连续迭代n1次，共迭代j+n1-1次，得到对应r=rj=jΔx时的区域异常
|
(20) |
和局部异常
|
(21) |
几点说明：
(1) 因为系数Cn=1，所以，式(18) 和式(19) 实质就是迭代公式
|
(22) |
(2) 从式(17) 到式(20) 的迭代与式(16) 等价.证明：用式(22) 进行第1次迭代后，在波数域中就有
|
(23) |
对
(3) 若由式(17) 到式(20) 的线性迭代按递增半径的方式进行，即以r=r1迭代n1次，再以r=r2，…，r=rj迭代，则这种递增半径迭代与递减半径迭代等价.这是由于成立
|
(4) 把式(23) 代入式(21) 的Fourier变换式中，得半径r=rj时的局部异常波谱：
|
(24) |
其中，局部异常传递函数
|
(25) |
因为递减半径线性迭代类似低通滤波，其传递函数Hdr(kx, ky, rj)在零波数点为1，所以，局部异常计算就类似高通滤波，其传递函数HL(kx, ky, rj)在零波数点为0，故局部异常波谱
|
(26) |
基于长方体磁异常波谱在零波数点为0，而长方体重力异常波谱在零波数点为长方体剩余质量的2πγ倍(γ是万有引力常数)，可认为上式对有限延伸磁性体的磁异常成立，对重力异常不成立，这会导致分离出的重力局部异常外侧出现反号假异常.
(5) 线性迭代公式(18) 每次迭代都减去迭代剩余异常值，其系数Cn恒定为1，未根据位场异常的变化来减小这个系数，这也会导致出现假异常且使分离出的异常畸变.
(6) 把简单的迭代公式(22) 表示为式(18) 和式(19) 这种较复杂的形式，有两个作用：其一是把区域异常迭代表现为多次减去迭代剩余异常的过程；其二是只要改变式(18) 中的Cn，就可以把这种线性迭代推广到下面的非线性迭代，以改善上述线性迭代法.
3.2 空间域递减半径非线性迭代法取区域异常迭代初始值为位场异常值
|
(27) |
区域异常迭代公式
|
(28) |
其中，Bn(x, y)是第n次迭代区域异常；Sn(x, y, r)是第n次迭代剩余异常，其算式为
|
(29) |
构造的非线性修正系数为
|
(30) |
式中，指数q取值范围0<q<2；Dn(x, y, r)是以计算点(x, y)为中心的四个中心差分绝对值的均值，按下式计算：
|
(31) |
式(28) 到式(31) 中r=rj－n=(j－n)Δx.当迭代到第n=j－1次，即r=r1=Δx时，以r=r1连续迭代n1次，共迭代j+n1-1次，得到对应r=rj=jΔx时的区域异常
|
(32) |
和局部异常
|
(33) |
从算法看，本节“递减半径非线性迭代”与上节“递减半径线性迭代”的唯一差别是，通过构造非线性修正系数Cn(x, y, r)，对迭代剩余异常Sn(x, y, r)进行修正，得修正后的迭代剩余异常Cn(x, y, r)·Sn(x, y, r)，使区域异常的计算变成多次从迭代区域异常中减去修正后的迭代剩余异常的过程.
非线性修正系数Cn(x, y, r)由Dn(x, y, r)和迭代剩余异常绝对值| Sn(x, y, r) |构造而成，且满足关系
|
(34) |
用Dn(x, y, r)反映迭代区域异常的变化，迭代区域异常不变化处Cn(x, y, r)=1，迭代区域异常变化大处，Cn(x, y, r)变小，目的在于压制假异常，减少异常畸变.
由于是非线性迭代，递减半径迭代与递增半径迭代不等价，且递减半径迭代优于递增半径迭代.
递减半径非线性迭代主要影响参数是圆周半径，较大半径对应较大深度上方各异常源产生的局部异常，因为该法简便快捷，可依据多个半径的局部异常和区域异常等值线图选取某一半径对应的异常分离结果.
顺便指出：若式(1) 中取b1=1/4，b2=0，则式(1) 为简单的四点圆周平均公式；取b1=b2=1/9，在式(1) 右边加一项b1u(x, y)，则式(1) 便成为9点滑动平均(moving average)公式，它们分别属于如下一组滤波器(Gonzalez and Woods, 2008；王一丁等. 2015) 中的h 1和h 5.
|
(35) |
把这组滤波器中任一个用于本文递减半径非线性迭代法，都能取得明显效果，比较各滤波器的递减半径迭代传递函数Hdr(kx, ky, rj)(由式(15) 计算)可知，h 4、h 5、h 6、h7、h 8优于h 2、h 3，而h 2、h 3优于h 1.
4 应用 4.1 应用于理论模型设有四个长方体场源，其中浅部三个小长方体作为局部场源，深部一个大长方体作为区域场源(图 3).小长方体源南北向长各1000、1200、1000 m，东西向宽各1000、1000、1000 m，厚各300、500、300 m，顶面埋深各300、500、300 m；大长方体源南北向长10000 m，东西向宽8000 m，厚6000 m，顶面埋深9000 m；剩余密度均为1.0 g·cm-3；磁化强度均为1.0 A·m-1，磁倾角均为45°，磁偏角均为0°.
|
图 3 长方体重、磁异常源模型 Fig. 3 Cuboid model for gravity and magnetic anomaly sources |
北为x坐标正向，东为y坐标正向.矩形网格参数为：256条线，每线256个点，线距Δx=100 m、点距Δy=100 m.
图 4和图 5分别为四场源模型(方框为其地表投影)的叠加重力异常Δg和磁异常ΔT；图 6a和图 7a分别为理论区域重力异常和理论局部重力异常；图 8a和图 9a分别为理论区域磁异常和理论局部磁异常.
|
图 4 模型重力异常Δg Fig. 4 Model gravity anomaly Δg |
|
图 5 模型磁异常ΔT Fig. 5 Model magnetic anomaly ΔT |
|
图 6 模型区域重力异常 (a)理论; (b)等半径非线性迭代; (c)递减半径非线性迭代. Fig. 6 Model regional gravity anomaly (a) Theory; (b) Equal radius nonlinear iteration; (c) Decreasing radius nonlinear iteration. |
|
图 7 模型局部重力异常 (a)理论; (b)等半径非线性迭代; (c)递减半径非线性迭代. Fig. 7 Model local gravity anomaly (a) Theory; (b) Equal radius nonlinear iteration; (c) Decreasing radius nonlinear iteration. |
|
图 8 模型区域磁异常 (a)理论; (b)等半径非线性迭代; (c)递减半径非线性迭代. Fig. 8 Model regional magnetic anomaly (a) Theory; (b) Equal radius nonlinear iteration; (c) Decreasing radius nonlinear iteration. |
|
图 9 模型局部磁异常 (a)理论; (b)等半径非线性迭代; (c)递减半径非线性迭代. Fig. 9 Model local magnetic anomaly (a) Theory; (b) Equal radius nonlinear iteration; (c) Decreasing radius nonlinear iteration. |
对四场源的重力异常(图 4)，分别用等半径非线性迭代法(即式(28) 按相等半径迭代ner次)和递减半径非线性迭代法计算了区域重力异常(图 6b和图 6c)和局部重力异常(图 7b和图 7c); 对四场源的磁异常(图 5)，分别用等半径非线性迭代法和递减半径非线性迭代法计算了区域磁异常(图 8b和图 8c)和局部磁异常(图 9b和图 9c).所取参数为：r=15Δx，ner=6，n1=6，对重力异常q=1.5，对磁异常q=0.3.
对比可见：两种方法分离异常(区域异常和局部异常)均能反映理论异常(区域异常和局部异常)，但等半径非线性迭代法的重力异常分离结果中存在异常畸变，其磁异常分离结果中除异常畸变外存在规则分布的假异常(图 8b和图 9b)；递减半径非线性迭代法明显优于等半径非线性迭代法.
图 10由图 5模型理论ΔT异常加方差为7.67 nT(ΔT最大值的5%)的高斯噪声，等半径非线性迭代分离出的区域磁异常(图 11a)中，存在规则分布的假异常及大量的噪声，递减半径非线性迭代分离出的区域磁异常(图 11b)与理论区域磁异常(图 8a)基本一致.
|
图 10 模型理论ΔT加高斯噪声 Fig. 10 Model theoretical ΔT with Gaussian noise |
|
图 11 模型区域磁异常 (a)等半径非线性迭代; (b)递减半径非线性迭代. Fig. 11 Model regional magnetic anomaly (a) Equal radius nonlinear iteration; (b) Decreasing radius nonlinear iteration. |
图 12是图 6—9南北主剖面(y=0) 上递减半径非线性迭代法、等半径非线性迭代法与理论值的对比.重力异常分离结果中，两种方法与理论值差值的均方根分别为0.0213 mGal和0.0192 mGal，中间局部异常较两侧局部异常更接近理论异常(图 12c)；磁异常分离结果中，两种方法与理论值差值的均方根分别为0.0843 nT和0.1396 nT，三个局部异常都接近理论异常(图 12d)；等半径非线性迭代法产生规则分布的假异常，剖面上表现为高频振荡(图 12b较图 12a更显著).
|
图 12
y=0剖面分离异常与理论异常的比较
(a)区域重力异常; (b)区域磁异常; (c)局部重力异常; (d)局部磁异常. t：理论曲线；dr：递减半径非线性迭代法；er：等半径非线性迭代法. Fig. 12 Comparisons between the theoretical anomalies and the separated anomalies by the two methods along profile y=0 (a) Regional gravity anomaly; (b) Regional magnetic anomaly; (c) Local gravity anomaly; (d) Local magnetic anomaly. t-Theory; dr-Decreasing radius nonlinear iteration; er-Equal radius nonlinear iteration. |
应该指出，图 3模型中区域场源埋深较大.随埋深减小，亦能分离出三个浅源产生的局部异常，但分离的局部异常与理论局部异常差异增大.
4.2 应用于泥河铁矿区泥河铁矿区位于安徽省庐江县泥河镇以南霍家院子一带.处于庐枞火山岩盆地西北部边缘，罗河—黄屯成矿带上，罗河铁矿东北侧3.5 km处.地表有薄层第四系覆盖，下伏地层由老到新为下白垩统砖桥组、双庙组、杨湾组.砖桥组和双庙组为火山碎屑岩及火山熔岩类岩性，分布全区，倾向北西，倾角平缓，总厚度600~800 m，下伏闪长玢岩岩体.杨湾组主要岩性为红色砂砾岩、泥岩等，习惯称之为“红层”，呈角度不整合覆盖于火山岩之上，分布于矿区西北部，厚度0~112 m，向西北方向增厚(吴明安等，2011).火山岩地层整体未形成明显的褶皱，主体为向北西倾的单斜构造(赵文广等，2011).
矿体埋深650~1000 m区间，产于闪长玢岩体与砖桥组接触带附近，矿体呈透镜状或似层状，其中铁矿主矿体产于闪长玢岩岩侵穹窿的顶部.铁硫矿体总体平面形态呈北东向不规则多边形(汪青松，2012)，本文把其叠放到图 13—18各图.其西南矿段以铁矿为主，局部也有硫铁矿体，主要矿石矿物为磁铁矿，次要矿物为黄铁矿、磁黄铁矿等；中矿段以石膏为主；东北矿段以硫铁矿为主，夹有磁铁矿体，主要矿石矿物为黄铁矿，次要矿物为磁黄铁矿、磁铁矿等.泥河铁矿由大型磁铁矿及大型硫铁矿、中型硬石膏矿、局部零星分布的铅锌矿组成的多矿种隐伏矿床，已探明磁铁矿1.87亿吨，硫铁矿1.39亿吨2).
2) 安徽五鑫矿业开发有限公司, 2010, 安徽省庐江县泥河铁矿勘探报告.
|
图 13 泥河铁矿布格重力异常Δg Fig. 13 Bouguer gravity anomaly map of the Nihe iron deposit |
图 13和图 14分别是泥河铁矿区1 : 10000布格重力异常图和1 : 5000地磁ΔT异常图(图中，x轴向北，y轴向东)；重力异常图西北部有北东向的罗河—霍家院子重力梯级带，Δg向北西降低，泥河矿区位于该重力梯级带上.泥河矿体的水平投影为北东向，与局部重力高不重合，且位于其北侧.地磁ΔT异常图西北部存在大规模的负异常，东南部存在跳跃和杂乱的正异常，难以识别矿异常(汪青松等，2012).
|
图 14 泥河铁矿地磁异常ΔT Fig. 14 Geomagnetic anomaly ΔT map of the Nihe iron deposit |
由物性资料(祁光等，2012)，西北部为杨湾组红层分布，地层密度(2.12~2.56 g·cm-3)小，为弱磁性或无磁性(磁化率：0~0.00023SI)；东南部为火山岩分布区，地层密度及磁化率(双庙组上段2.57~2.61 g·cm-3，0.02292~0.03383SI，双庙组下段2.51~2.60 g·cm-3，0.00507~0.01959SI，砖桥组2.61~2.69 g·cm-3，0.00134~0.00594SI)较大；闪长玢岩密度(2.85~2.95 g·cm-3)大，磁化率(0.00643~0.01590SI)较大；石膏矿密度(2.80~3.10 g·cm-3)大，磁性(0~0.00004SI)弱；黄铁矿(3.17~3.51 g·cm-3，0.00068~0.51000SI)和磁铁矿(3.34~3.61 g·cm-3，0.12250~0.20000SI)较围岩存在较大的密度和磁性差异.
用多项式拟合法、向上延拓法、插值切割法、匹配滤波法对研究区布格重力异常分离结果(刘彦等，2012)表明，这几种方法分离出的局部异常(或区域异常)数值范围不同，形态各异，就对泥河矿体的定性解释而言，多项式拟合法及插值切割法效果较好，向上延拓法次之，匹配滤波法失效.
本文对研究区布格重力异常Δg及地磁ΔT异常(Δx=Δy= 50 m)进行了异常分离，用等半径非线性迭代法和递减半径非线性迭代法分别计算，所取参数为：r= 30Δx，ner=6，n1=6，对重力异常q=0.8，对地磁异常q=0.3，异常分离结果为图 15—18.
|
图 15 等半径非线性迭代法重力异常分离 (a)区域异常; (b)局部异常; (c)去浅层后的局部异常. Fig. 15 Gravity anomalies separated by the equal radius nonlinear iterative method (a) R(x, y, r30); (b) L(x, y, r30); (c) L(x, y, r30)-L(x, y, r2). |
|
图 16 递减半径非线性迭代法重力异常分离 (a)区域异常; (b)局部异常; (c)去浅层后的局部异常. Fig. 16 Gravity anomalies separated by the decreasing radius nonlinear iterative method (a) R(x, y, r30); (b) L(x, y, r30); (c) L(x, y, r30)-L(x, y, r2). |
|
图 17 等半径非线性迭代法磁异常分离 (a)区域异常; (b)局部异常; (c)去浅层后的局部异常. Fig. 17 Geomagnetic anomalies separated by the equal radius nonlinear iterative method (a) R(x, y, r30); (b) L(x, y, r30); (c) L(x, y, r30)-L(x, y, r2). |
|
图 18 递减半径非线性迭代法磁异常分离 (a)区域异常; (b)局部异常; (c)去浅层后的局部异常. Fig. 18 Geomagnetic anomalies separated by the decreasing radius nonlinear iterative method (a) R(x, y, r30); (b) L(x, y, r30); (c) L(x, y, r30)-L(x, y, r2). |
递减半径非线性迭代法分离出的区域重力异常(图 16a)与区域磁异常(图 18a)总体呈北东走向，沿北西向数值递减，反映了地下砖桥组和双庙组火山碎屑岩及火山熔岩走向北东、平缓倾向北西的单斜地层特征；由于局部重力异常(图 16b)及局部磁异常(图 18b)得自观测场减去区域场，其自然包含观测场中的各种高频干扰，尤其是局部磁异常的高频干扰更多，难以识别矿异常，因此，将r=30Δx时的局部异常减去r=2Δx时的局部异常(图 16c和图 18c)，可见，重力局部异常(图 16c)形态、位置与矿体很好的对应，地磁局部异常(图 18c)中的正、负异常亦能反映斜磁化的泥河矿体.等半径非线性迭代法对重力异常的分离结果(图 15)中有假异常及明显的异常畸变，对地磁异常的分离结果(图 17)中有很多高频干扰，难以反映地下的地层和矿体.
5 结论本文通过导出一个新的八点圆周平均公式，由其传递函数构造了递减半径迭代传递函数，由后者给出计算位场区域异常的空间域递减半径线性迭代法，进一步给出计算位场区域异常的空间域递减半径非线性迭代法，并对该方法进行试验.具体结论如下：
(1) 导出了一个新的八点圆周平均公式及其传递函数，该传递函数以波数为自变量，以圆周半径为参变量，它具有以零波数点为中心的主瓣及其他规则分布点为中心的旁瓣，各负旁瓣极小值约为-0.6，主瓣及各正旁瓣极大值为1，规则分布的旁瓣使分离出的位场异常中出现规则分布的假异常及高频干扰，圆周半径越大，则旁瓣越多，使高频干扰越多.这揭示了圆周平均公式在圆周半径较大时分离效果差或失效的原因.
(2) 由固定半径的八点圆周平均公式的传递函数的乘幂构造了等半径迭代传递函数，若圆周半径大于一倍点距，其各正旁瓣依旧存在，只是变窄；这表明圆周平均公式的等半径多次迭代在圆周半径较大时分离效果亦差或失效.
(3) 由从大到小不同半径对应的八点圆周平均公式的传递函数的连乘积构造了递减半径迭代传递函数.该传递函数以大半径为参变量；该传递函数几乎无旁瓣，只存在主瓣，类似低通滤波器，截止波数与大半径成反比，与向上延拓滤波器相比，截止波数相近，而形态有差异。这表明圆周平均公式按递减半径迭代时异常分离效果得到显著的改善；还表明圆周半径较大，便得到较大深度的区域异常和较大深度上方各异常源的局部异常.
(4) 由递减半径迭代传递函数，给出了分离位场区域异常的空间域递减半径线性迭代法，其分离局部场的传递函数在零波数点为零，即分离局部场的积分为零。这表明空间域递减半径线性迭代法分离的重力局部场中真异常周围存在异号假异常及异常畸变.
(5) 通过构造非线性修正系数，把空间域递减半径线性迭代法中线性迭代公式变成非线性迭代公式，就得到空间域递减半径非线性迭代法，后者除具有前者的优点外，还能压制假异常及减少异常畸变.
(6) 空间域递减半径非线性迭代法(简称方法1) 中非线性迭代公式每次按相等半径迭代若干次，就是空间域等半径非线性迭代法(简称方法2).对理论模型和泥河矿区的重磁数据，分别应用这两种方法进行异常分离.对模型重磁异常，圆周半径较大(15倍点距)时，方法2分离结果中包含多个规则分布的假异常，异常畸变较大，不能压制高频干扰，方法1分离结果与模型理论异常接近，且有较强的抗高频干扰能力.对泥河矿区重磁实际资料，方法1分离的区域异常反映了地下单斜地层，分离的局部异常反映了泥河矿体；由于圆周半径较大(30倍点距)，方法2对重力异常的分离效果较差，对磁异常的分离失效.
致谢感谢安徽省勘查技术院汪青松及安徽省国土资源厅何义权两位教授级高工对本文工作的讨论和建议.
| Agocs W B. 1951. Least square residual anomaly determination. Geophysics, 16(4): 686-696. DOI:10.1190/1.1437720 | |
| Cheng F D, Liu D J, Yao R X. 1987. A study on the identification of regional and local gravity fields. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 9(1): 1-9. | |
| Fedi M, Quarta T. 1998. Wavelet analysis for the regional-residual and local separation of potential field anomalies. Geophysical Prospecting, 46(5): 507-525. DOI:10.1046/j.1365-2478.1998.00105.x | |
| Gonzalez R C, Woods R E. 2008. Digital Image Processing. 3rd ed. Upper Saddle River, NJ: Prentice Hall. | |
| Guo L H, Meng X H, Shi L, et al. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinnese continent. Chinese J. Geophys. (in Chinese), 55(12): 4078-4088. DOI:10.6038/j.issn.0001-5733.2012.12.020 | |
| Griffin W R. 1949. Residual gravity in theory and practice. Geophysics, 14(1): 39-56. DOI:10.1190/1.1437506 | |
| Hou Z Z, Yang W C. 1997. Wavelet transform and multi-scale analysis on gravity anomalies of China. Chinese J. Geophys. (in Chinese), 40(1): 85-95. | |
| Jacobsen B H. 1987. A case for upward continuation as a standard separation filter for potential-field maps. Geophysics, 52(8): 1138-1148. DOI:10.1190/1.1442378 | |
| Keating P, Pinet N. 2011. Use of non-linear filtering for the regional-residual separation of potential field data. Journal of Applied Geophysics, 73(4): 315-322. DOI:10.1016/j.jappgeo.2011.02.002 | |
| Li Q C, Pan Z S, Li J L. 1994. Order-variant sliding trend analysis and its application. Geophysical Prospecting for Petroleum (in Chinese), 33(3): 92-98, 114. | |
| Liu C Y, Yao C L, Zheng Y M. 2015. Preferential spatially varying filtering method in the wavelet domain for gravity anomaly separation. Chinese J. Geophys. (in Chinese), 58(12): 4740-4755. DOI:10.6038/cjg20151234 | |
| Liu D J. 1984. The regional-residual separation of the gravity anomaly. Changsha: The Central-South Institute of Mining and Metallurgy. | |
| Liu D J, Hong T Q, Jia Z H, et al. 2009. Wave number domain iteration method for downward continuation of potential fields and its convergence. Chinese J. Geophys. (in Chinese), 52(6): 1599-1605. DOI:10.3969/j.issn.0001-5733.2009.06.022 | |
| Liu D J, Hong T Q, Liao X T, et al. 2012. Iterative solution of integral equation for potential field continuation from an irregular surface to a horizontal plane. Chinese J. Geophys. (in Chinese), 55(10): 3467-3476. DOI:10.6038/j.issn.0001-5733.2012.10.030 | |
| Liu T Y, Wu Z C, Zhan Y L. 2007. Wavelet multi-scale decomposition of magnetic anomaly and its application in searching for deep-buried minerals in crisis mines: a case study from Daye iron mines. Earth Science—Journal of China University of Geosciences (in Chinese), 32(1): 135-140. | |
| Liu Y, Yan J Y, Wu M A, et al. 2012. Exploring deep concealed ore bodies based on gravity anomaly separation methods: A case study of the Nihe iron deposit. Chinese J. Geophys. (in Chinese), 55(12): 4181-4193. DOI:10.6038/j.issn.0001-5733.2012.12.030 | |
| Meng X H, Guo L H, Chen Z X, et al. 2009. A method for gravity anomaly separation based on preferential continuation and its application. Applied Geophysics, 6(3): 217-225. DOI:10.1007/s11770-009-0025-y | |
| Naudy H, et Dreyer H. 1968. Essai de filtrage non-linéaire appliqué aux profils aéromagnétiques. Geophysical Prospecting, 16(2): 171-178. DOI:10.1111/gpr.1968.16.issue-2 | |
| Oldham C H G, Sutherland D B. 1955. Orthogonal polynomials: their use in estimating the regional effect. Geophysics, 20(2): 295-306. DOI:10.1190/1.1438143 | |
| Pawlowski R S. 1995. Preferential continuation for potential-field anomaly enhancement. Geophysics, 60(2): 390-398. DOI:10.1190/1.1443775 | |
| Qi G, Lü Q T, Yan J Y, et al. 2012. Geologic constrained 3D gravity and magnetic modeling of Nihe deposit—A case study. Chinese J. Geophys. (in Chinese), 55(12): 4194-4206. DOI:10.6038/j.issn.0001-5733.2012.12.031 | |
| Rao B S R, Murthy I V, Rao C V. 1975. A successive approximation method of deriving residual gravity. Geoexploration, 13(2): 129-135. | |
| Spector A, Grant F S. 1970. Statistical models for interpreting aeromagnetic data. Geophysics, 35(2): 293-302. DOI:10.1190/1.1440092 | |
| Wang Q S, Wu M A, Yang P, et al. 2012. Characteristics of gravity and magnetic anomalies in the Nihe iron deposit of Lujiang County, Anhui Province. Geology and Exploration (in Chinese), 48(1): 148-154. | |
| Wang Y D, Li C, Wang Y H. 2015. Digital Image Processing (in Chinese). Xian: Xidian University Press. | |
| Wen B H, Cheng F D. 1990. A new interpolating cut method for identifying regional and local fields of magnetic anomaly. J. Cent. -South Inst. Min. Metall. (in Chinese), 21(3): 229-235. | |
| Wu M A, Wang Q S, Zheng G W, et al. 2011. Discovery of the Nihe iron deposit in Lujiang, Anhui, and its exploration significance. Acta Geologica Sinica (in Chinese), 85(5): 802-809. | |
| Xu D S, Zeng H L. 2000. Preferential continuation and its application to Bouguer gravity anomaly in China. Geoscience (in Chinese), 14(2): 215-222. | |
| Xu S Z, Yu H L, Li H X, et al. 2009. The inversion of apparent density based on the separation and continuation of potential field. Chinese J. Geophys. (in Chinese), 52(6): 1592-1598. DOI:10.3969/j.issn.0001-5733.2009.06.021 | |
| Yang W C, Shi Z Q, Hou Z Z, et al. 2001. Discrete wavelet transform for multiple decomposition of gravity anomalies. Chinese J. Geophys. (in Chinese), 44(4): 534-541. DOI:10.3321/j.issn:0001-5733.2001.04.012 | |
| Zeng H L, Xu D S, Tan H D. 2008. A model study for estimating optimum upward-continuation height for gravity separation with application to a Bouguer gravity anomaly over a mineral deposit, Jilin province, northeast China. Geophysics, 72(4): I45-I50. | |
| Zhao W G, Wu M A, Zhang Y Y, et al. 2011. Geological characteristics and genesis of the Nihe Fe-S deposit, Lujiang Country, Anhui Province. Acta Geologica Sinica (in Chinese), 85(5): 789-801. | |
| 程方道, 刘东甲, 姚汝信. 1987. 划分重力区域场与局部场的研究. 物探化探计算技术, 9(1): 1–9. | |
| 郭良辉, 孟小红, 石磊, 等. 2012. 优化滤波方法及其在中国大陆布格重力异常数据处理中的应用. 地球物理学报, 55(12): 4078–4088. DOI:10.6038/j.issn.0001-5733.2012.12.020 | |
| 侯遵泽, 杨文采. 1997. 中国重力异常的小波变换与多尺度分析. 地球物理学报, 40(1): 85–95. | |
| 李庆春, 潘作枢, 李九亮. 1994. 变阶数滑动趋势分析及其应用. 石油物探, 33(3): 92–98, 114. | |
| 刘彩云, 姚长利, 郑元满. 2015. 重力异常分离的小波域优化位变滤波方法. 地球物理学报, 58(12): 4740–4755. DOI:10.6038/cjg20151234 | |
| 刘东甲. 1984. 划分重力区域场与局部场[硕士学位论文]. 长沙: 中南矿冶学院. | |
| 刘东甲, 洪天求, 贾志海, 等. 2009. 位场向下延拓的波数域迭代法及其收敛性. 地球物理学报, 52(6): 1599–1605. DOI:10.3969/j.issn.0001-5733.2009.06.022 | |
| 刘东甲, 洪天求, 廖旭涛, 等. 2012. 位场曲化平积分方程的迭代解. 地球物理学报, 55(10): 3467–3476. DOI:10.6038/j.issn.0001-5733.2012.10.030 | |
| 刘天佑, 吴招才, 詹应林. 2007. 磁异常小波多尺度分解及危机矿山的深部找矿:以大冶铁矿为例. 地球科学—中国地质大学学报, 32(1): 135–140. | |
| 刘彦, 严家永, 吴明安, 等. 2012. 基于重力异常分离方法寻找深部隐伏铁矿——以安徽泥河铁矿为例. 地球物理学报, 55(12): 4181–4193. DOI:10.6038/j.issn.0001-5733.2012.12.030 | |
| 祁光, 吕庆田, 严家永, 等. 2012. 先验地质信息约束下的三维重磁反演建模研究——以安徽泥河铁矿为例. 地球物理学报, 55(12): 4194–4206. DOI:10.6038/j.issn.0001-5733.2012.12.031 | |
| 汪青松, 吴明安, 袁平, 等. 2012. 安徽省庐江县泥河铁矿重磁异常特征. 地质与勘探, 48(1): 148–154. | |
| 王一丁, 李琛, 王蕴红. 2015. 数字图像处理. 西安: 西安电子科技大学出版社. | |
| 文百红, 程方道. 1990. 用于划分磁异常的新方法—插值切割法. 中南矿冶学院学报, 21(3): 229–235. | |
| 吴明安, 汪青松, 郑光文, 等. 2011. 安徽庐江泥河铁矿的发现及意义. 地质学报, 85(5): 802–809. | |
| 许德树, 曾华霖. 2000. 优选延拓技术及其在中国布格重力异常图处理上的应用. 现代地质, 14(2): 215–222. | |
| 徐世浙, 余海龙, 李海侠, 等. 2009. 基于位场分离与延拓的视密度反演. 地球物理学报, 52(6): 1592–1598. DOI:10.3969/j.issn.0001-5733.2009.06.021 | |
| 杨文采, 施志群, 侯遵泽, 等. 2001. 离散小波变换与重力异常多重分解. 地球物理学报, 44(4): 534–541. DOI:10.3321/j.issn:0001-5733.2001.04.012 | |
| 赵文广, 吴明安, 张宜勇, 等. 2011. 安徽省庐江县泥河铁硫矿床地质特征及成因初步分析. 地质学报, 85(5): 789–801. | |