0 引 言

化极是将总磁异常化向地磁极，以消除斜磁化的影响，是磁异常处理、解释的基础，也是磁法勘探的研究热点问题.在区域矿产地质调查中，由于工作区往往包含多个图幅，纬度变化范围大，斜磁化的影响程度随纬度不同而变化，为了更好地发挥磁测资料在综合地质、地球物理解释中的作用，将总磁异常进行化极，就显得尤为重要.国内外许多学者都研究过磁异常化极方法(王纪恒，1993；柴玉璞等，1998; Li et al., 2000; 许东青等，2006; Maus et al., 2009; 骆遥等，2010; 彭利丽等，2012; 汪青松等，2012)，并取得了一定成果.其中Baranov(1957)最早提出化极的概念，并在空间域进行化极；Arkan等(1988)提出了一种区域差值化极方法，这种分区计算存在计算速度慢、不适用于大数据量计算等问题；姚长利等(2004，2012)；林晓星等(2012)对直接阻尼法、迭代化极方法进行了改进，使其能适用于低纬度地区，并且提高了迭代的效率；石磊等(2012)采用伪倾角化极因子，对南海海域总场异常进行化极，使得化极更为稳定；彭利丽等(2012)在频率域对多种化极因子处理方法进行对比，并对部分方法作了改进，在南海地区取得了明显的应用效果；张培琴等(1996)、方迎尧等(2006)和刘振军等(2010)提出用滑动窗口分带变倾角磁方向转换方法进行化极来解决采用单一倾角对大区域进行化极处理带来的误差在一定程度上可提高计算效率，但在全区处理中保证计算精度方面需要进一步优化.

本文针对在高、中磁纬度区大面积磁法测量及编图工作的磁异常化极处理中，可能存在的测区内磁纬度的变化范围很大的问题进行讨论，不再采用单一的磁倾角作化极处理，试验研究在频率域进行大区域逐点变磁倾角实现化极的方法，即对测区内每一个点用实际对应的磁倾角和磁偏角逐一进行化极，并取该点的化极结果，逐点进行循环处理.处理过程采用快速傅里叶变换，达到减少计算时间，提高计算效率的目的，并避免了窗口分区拼接产生的误差.理论模型试算和实际资料计算结果表明，该方法适应快速高精度处理多图幅磁法联测的整体处理与解释，能明显改善目前大面积矿产资源磁法勘探资料连片处理中的化极效果.

1 基本原理

1.1 频率域化极的基本公式

磁异常化极处理既可在空间域进行，也可在频率域进行，空间域表现为较为复杂的褶积关系，而频率域表现为乘积形式，即由实测异常Δ T(x，y)的傅里叶变换频谱ST(u，v)乘以相应的转换因子，便可得到变换后的异常频谱S(u，v)，再经傅里叶反变换，就得到转换的结果.

频率域化极中的转换因子的一般形式为

进行化极时，I₂=90°(垂直的测量方向)，I₃=90°(垂直的磁化强度方向)，则q₂=q₃= u²+v²，此时测量总场磁异常ΔT所对应的测量方向就是地磁场方向.

假设磁化强度方向与地磁场方向一致，就有q₀=q₁，I₀，D₀为正常地磁场方向(磁化方向)倾角和偏角，此时化极因子可简化为

1.2 球谐分析计算磁倾角和磁偏角

地面高精度磁测时，对每个测点均进行了高精度GPS坐标定位，利用测点平面坐标信息可以转换为大地经纬度坐标，再借助地磁场球谐分析方法便可计算出理论地磁场磁化方向倾角和偏角.球谐分析法于1838年由高斯首先提出，他将球谐分析法用于分析地磁场，成功地区分了内外源场，找到了模拟地磁场的数学表达式，该方法是表示全球范围地磁场分布及其长期变化的一种数学方法.

如果已知球谐系数和某点地理坐标经纬度，便可计算地球表面和它外部任意一点的地磁要素值，并进一步根据式(2)得到化极因子.

2 变倾角化极实现步骤

变倾角化极主要针对测区含多个图幅且数据点较多的情况，实现在频率域化极处理，其计算过程和主要步骤如下：

(1)输入全区实测总场磁异常ΔT(x_i，y_i)，利用快速傅里叶变换计算其频谱S_T(u_i，v_i)；

(2)根据球谐分析法计算出每一个实测点(x_i，y_i)的地磁要素——磁倾角I_i和磁偏角D_i；

(3)利用公式(2)计算每一个点的化极因子H(u_i，v_i)；

(4)由(x_i，y_i)点的化极因子H(u_i，v_i)，根据公式(3)计算全区所有点磁异常化极后的频谱S(u，v)；

(5)将求得的频谱S(u，v)进行反傅里叶变换，得到全区所有点化极结果Z_a⊥(x，y)；只取Z_a⊥(x_i，y_i)作为待计算点(x_i，y_i)的化极结果；

(6)循环(4)~(5)，得到整个工区所有点的变磁倾角化极结果.

为提高计算效率，在计算过程中采取如下优化方案：将(5)求得的Z_a⊥(x_i，y_i)及时存入硬盘，既节省内存空间，又防止因程序运行忽然中断等原因导致数据丢失，实现动态监控；对特别大的数据量，设置子块计算参数，分别由不同计算机同时计算，之后再合并结果，实现伪并行计算，减少整体运算时间.

3 理论模型试验

为验证上述方法的有效性，设计了以下变磁倾角理论模型进行试算：异常由两个棱柱体产生，其磁参数及模型位置参数见表 1，按140 m×200 m网格计算的总磁异常场如图 1a所示，理论化极结果如图 1b所示，图中实线绘制的四边形表示棱柱体模型边界的水平投影.

表 1(Table 1)

表 1 理论模型参数 Table 1 Theoretical model parameters 棱柱体编号 磁化强度10-3A/m 磁倾角 磁偏角 X起止边界(m) Y起止边界(m) Z起止边界(m) A 2000 40° 10° 620~770 420~570 100~180 B 2000 60° 5° 620~770 1480~1630 100~180

表 1 理论模型参数 Table 1 Theoretical model parameters

图 1

Fig. 1

图 1 磁异常理论模型及化极结果 (a)理论总磁异常场；(b)理论化极异常. 等值线单位nT，坐标刻度单位m，四边形表示棱柱体边界的水平投影. Fig. 1 Theoretical total filed anomalies and RTP results (a)Theoretical total filed anomalies;(b)RTP results in theory.

分别对模型中的棱柱体采用单一磁倾角进行化极，化极结果分别如图 2，3，4所示，图中实线表示理论化极结果，总体异常幅值变化范围在-5~+175 nT；虚线表示单一磁倾角化极结果.

图 2

Fig. 2

图 2 采用棱柱体A的磁倾角的化极结果与理论化极结果对比 Fig. 2 Comparison of the RTP results with inclination from rectangle A and those in theory

其中图 2为对棱柱体A和B同时采用棱柱体A的磁倾(偏)进行化处理得到的结果，其异常幅值变化范围在-10~+195 nT，可以看出，棱柱体A化极后的异常极大值为175 nT，与理论结果一致，正负异常的中心位置及异常整体形态也与理论异常基本吻合，棱柱体B的化极结果与理论相比出现较大偏差，异常极大值相差约20 nT，正异常中心位置向北移动幅度偏大，负异常出现明显畸变.图 3是采用棱柱体B的磁倾角和磁偏角进行化极得到的结果，其异常幅值在-18~+175 nT，由图可知，棱柱体B化极后的异常幅值变化范围为-5~+175 nT，其异常值大小、中心位置及整体异常形态均与理论结果基本吻合，棱柱体A化极后异常极大值与理论结果相差较大，差约25 nT，负异常形态出现明显畸变，极小值相差约5 nT，正负异常中心位置相对理论结果偏南，这是因为化极所用的磁倾角大于理论值.由以上分析可知，采用某一个异常体的理论磁倾角和磁偏角进行化极时，只能使其中一个异常体得到较好的恢复，而另一个异常体的极值大小、中心位置及异常总体形态与理论相比均出现较大偏差.

图 3

Fig. 3

图 3 采用棱柱体B的磁倾角的化极结果与理论化极结果对比 Fig. 3 Comparison of the RTP results with inclination from rectangle B and those in theory

图 4是对棱柱体A和B同时采用工区中心位置磁倾角和磁偏角进行化极得到的结果，整体异常幅值范围在-12~+184 nT，其中，棱柱体A异常极大值160 nT，极小值 -12 nT，棱柱体B异常极大值184 nT，极小值-7 nT，与理论值对比可知，两个异常体的极值均存在一定偏差，正负异常的中心位置有明显移动，异常整体形态出现明显畸变，尤其负异常畸变更严重.通过分析表明，当实测工区中不同位置的磁倾角变化范围较大，而采用单一磁倾角进行化极时，化极结果与理论相比出现较大的偏差，从而进一步导致解释误差.

图 4

Fig. 4

图 4 采用工区中间位置磁倾角的化极结果与理论化极结果对比 Fig. 4 Comparison of the RTP results with inclination from center of surveying area and those in theory

图 5为采用本文提出的全区变磁倾角方法得到的化极结果，其中实线表示理论化极结果，虚线表示全区变磁倾角化极结果，其异常幅值变化范围在-5~+175 nT.由理论化极结果与实际化极结果对比可见，全区变磁倾角的化极结果的极值大小、正负异常中心位置及异常整体性形态均与理论值相吻合，虽然负异常极值中心位置与理论有微小偏差，但误差在允许范围内，相对于采用单一磁倾角来说，精度最高，表明用该方法处理磁倾角(偏)角变化较大的磁测数据是可行的.

图 5

Fig. 5

图 5 全区变磁倾角化极结果与理论化极结果对比 Fig. 5 Comparison of the RTP results with varying inclinations and those in theory

通过理论模型分析表明，采用单一磁倾角进行化极时，若使用的磁倾角小于理论值，则导致正异常中心位置向北移动幅度比理论移动幅度大；若使用的磁倾角大于理论值，则导致正异常的中心位置向北移动幅度比理论移动幅度小；只有当使用变磁倾角进行化极时，才能得到正异常中心位置移动幅度适中的化极结果，这为变磁倾角处理实际资料的正确性、可行性提供了依据. 4 实际资料处理结果对比

本文在理论模型试验的基础上，进行了实际资料试处理.以我国某地1 5万高精度地面磁测实测磁异常资料为例.该区由12个1 5万标准图幅组成，经度跨越1°15′，纬度跨越40′，理论地磁倾角变化53.5′，总面积达5000多平方公里，测点总数98000多个；单一磁倾角化极所用磁倾角、偏角位置位于工区中心，选取的典型异常区位于工区右上方.此典型异常区共有3528个测点，程序运行用时95 s.

图 6为单一磁倾角化极与变磁倾角化极效果对比.图 6a是实测磁异常色块图，图 6b是单一磁倾角化极色块图，图 6c是变磁倾角化极色块图.整体上，单一磁倾角与变磁倾角化极后，异常整体形态均与化极前形似，正异常的中心位置均明显向北移动；由于单一磁倾角化极所用磁倾、偏角小于典型异常区各测点的实际值，变磁倾角化极后正异常中心位置不仅应向北移动，且移动幅度应小于单一磁倾角化极.

图 6

Fig. 6

图 6 单一磁倾角化极与变磁倾角化极效果对比 (a)实测磁异常；(b)单一磁倾角化极；(c)变磁倾角化极. Fig. 6 Comparison of the RTP results with single inclination and those with varying inclinations (a)Magnetic total field anomaly;(b)RTP results with single inclination;(c)RTP results with varying inclinations.

选择3个典型异常体A、B、C进行具体分析.异常体A：化极前正负异常伴生，异常幅值范围在-430~+800 nT间跳跃，正异常中心位置坐标约为(2124 m，7655 m)；单一磁倾角化极后，异常体A正异常极值增加到900 nT，负异常极值减小到-160 nT，正异常中心位置坐标约为(2124 m，8088 m)，向北明显移动约433 m；变磁倾角化极后，异常体A范围缩小，形态简化，正负异常伴生明显，异常幅值范围为-270~+550 nT，正异常中心位置坐标约为(2124 m，8054 m)，向北移动约399 m，移动幅度小于单一磁倾角化极结果，符合理论模型的验证结果.异常体B：化极前正负异常伴生，负异常展布范围比伴生的正异常稍大，异常幅值范围为-200~+600 nT，正异常中心位置坐标约为(4815 m，2725 m)；单一磁倾角化极后，异常体B以正异常为主，周围区域为负异常，正异常极值增加到650 nT，中心位置坐标约为(4815 m，3520 m)，向北移动约795 m；变磁倾角化极后，异常体B异常幅值范围为-100~+300 nT，正负异常伴生及相互分离更明显，异常范围缩小，正异常中心位置坐标约为(4815 m，2960 m)，相对化极前向北移动约560 m，移动幅度小于单一磁倾角化极结果，符合理论模型的验证结果.异常体C：化极前以正异常为主，周围区域为负异常，异常幅值范围在-300~+500 nT间跳跃，正异常中心位置坐标约是(3804 m，555 m)；单一磁倾角化极后，异常体C幅值范围为-350~+460 nT，正异常范围缩小，中心位置坐标约为(3804 m，747 m)，向北移动约192 m；变磁倾角化极后，异常体C正负异常伴生明显，异常幅值范围为-275~+500 nT，正异常中心位置坐标约为(3804 m，738 m)，相对化极前向北移动约183 m，移动幅度小于单一磁倾角化极结果，符合理论模型的验证结果.

由实际资料分析可知，变磁倾角化极结果与理论模型得出的结论相一致，相比单一磁倾角化极方法而言，变磁倾角化极结果减少了有效信号的化极特征的损失，更符合实际情况，具有实际意义.

5 结 论

对于大面积磁测数据的化极若采用单一磁倾角化极会使化极后的大部分地区产生超过误差允许范围的不可忽视的误差，难以取得理想的化极效果.本文针对大面积磁测数据处理采用单一磁倾角化极不尽合理这一问题展开研究，改进了以往的化极实现过程，尝试了频率域进行大区域逐点变磁倾角实现化极的方案，首先对全区总磁异常通过快速傅里叶变换得到频谱，然后借助球谐分析法确定每一个测点的地磁要素，并进一步得到相应的化极因子和化极后的频谱，最后通过反傅里叶变换得到化极结果，对于测区内的每个点进行循环计算，实现每个测点的高精度化极；计算时在普通单机上采用伪并行计算，计算时间不超过100 s，在保证化极效果的基础上提高了化极计算的效率；化极原理简单且方便实用.理论模型和实际资料均表明，变倾角磁异常化极方法算法稳定，收敛速度快，避免了窗口分区拼接产生的误差，提高了转换计算精度，化极结果可靠，能适用于测区跨度大磁异常化极问题的解释，对目前大面积矿产资源磁法勘探资料连片处理也是非常有效的.