内容摘要：磁法勘探作为重力勘探的姊妹手段，常用于金属矿藏探测，对铁磁性物质敏感。磁化率和密度作为地球介质的两个物性重要参数，在地下结构模型构建中具有重要的应用与研究价值。地磁和重力场常统称为位场，对于同一个场源体，两者之间可以通过泊松公式实现直接转换（图1），从位场反演的核函数构建来讲，两种方法之间并无差别。但地磁异常的处理分析上与重力异常还略有差别，今天我们一起讨论几个必须知道的关键概念。
1、地磁模型正演
（1）地磁异常的概念
地磁场的理论分布是有变化的。而实际上测得的地球磁场强度和理论磁场强度是有区别的，这种区别称地磁异常。在实际工作中，会发现某地区实测地磁场要素的数据与正常值有显著的差别，这种现象称作地磁异常。
一般野外实测的地磁数据，通常包括三部分：主磁场、外源场、岩石圈磁场。通常磁法勘探中，关注的主要是岩石圈磁场。主磁场一般根据IGRF或WMM模型计算得到即可，外源场可以在勘探区位置架设一个日变站，然后把每测量结果与其观测结果相减即可。
异常是磁法勘探常提到的，那这个量如何定义的呢？
如图所示，表示测量到的总磁场强度，
为主磁场强度，
为磁异常强度。
但在实际计算的过程中，为了简化计算令，
。
由于测量点的总场方向与主磁场方向总有一个差异角度，直接减会引起误差
，有多大呢？如下面表所示。
（2）地磁异常与居里面深度
19世纪末，著名物理学家皮埃尔·居里（居里夫人的丈夫）在自己的实验室里发现磁石的一个物理特性，就是当磁石加热到一定温度时，原来的磁性就会消失。后来，人们把这个温度叫“居里点“。
根据这个理论在地壳内部温度发生变化时，岩石会失去磁性，因此，通过磁异常图可以研究岩石圈的温度状态，居里面就是这样一个温度界面。
居里面以下的物质一般认为没有磁性，因此这个界面通常是地磁异常场源解释的底界面。
（3）程序实现
看过前面重力模型正演方法介绍后，其实磁异常正演也非常类似，直接看代码吧！
inc, dec = 30, -15
bounds = [-5000, 5000, -5000, 5000, 0, 5000]
model = [
geometry.Prism(-4000, -3000, -4000, -3000, 0, 2000,
{'magnetization': giutils.ang2vec(1, inc, dec)}),
geometry.Prism(-1000, 1000, -1000, 1000, 0, 2000,
{'magnetization': giutils.ang2vec(1, inc, dec)}),
# This prism will have magnetization in a different direction
geometry.Prism(2000, 4000, 3000, 4000, 0, 2000,
{'magnetization': giutils.ang2vec(3, -10, 45)})]
# Create a regular grid at 100m height
shape = (200, 200)
area = bounds[:4]
xp, yp, zp = gridder.regular(area, shape, z=-500)
# Calculate the anomaly for a given regional field
tf = prism.tf(xp, yp, zp, model, inc, dec)
# Plot
plt.figure()
plt.title("Total-field anomaly (nT)")
plt.axis('scaled')
giplt.contourf(yp, xp, tf, shape, 15)
plt.colorbar()
plt.xlabel('East y (km)')
plt.ylabel('North x (km)')
giplt.m2km()
plt.show()
2、地磁分量转换
地磁场的七个分量，一般知道三个可以计算剩余四个。图3是地磁场分量示意图。
在任一点P上，磁场矢量（红线）通常用其方向、在该方向上的总值F以及H和Z，F的局部水平分量和垂直分量表示。角D和I表示磁场矢量的定向。磁偏角D为水平面上H和地理北之间的夹角。磁倾角I为磁场矢量和H所在水平面之间的夹角。
from geoist.pfm.magdir import DipoleMagDir
# Give the centers of the dipoles
centers = [[3000, 3000, 1000], [7000, 7000, 1000]]
# Estimate the magnetization vectors
solver = DipoleMagDir(x, y, z, tf, inc, dec, centers).fit()
# Print the estimated and true dipole monents, inclinations and declinations
print('Estimated magnetization (intensity, inclination, declination)')
for e in solver.estimate_:
print(e)
估计结果：
Estimated magnetization (intensity, inclination, declination)
[2514.452175001507, -20.01579001683044, -10.036363448371244]
[4181.363570901708, 3.0373253097068478, -66.96549442175747]
3、地磁数据化极
化到地磁极（RTP，reduced to the pole），是—种将斜磁化或
磁异常换算为各种垂直磁化磁异常的磁场换算方法，因为地球磁场在每个位置方向都不相同，在两极近似垂直磁化，赤道位置近似为水平磁化（图5），故称化到地磁极，简称：化极。
一般来讲低纬度地区磁异常比较复杂，采用化到地磁极，异常直观简单，较易于解释。化极后的磁异常一般就可以直接进行解释了，这也是地磁异常解释相比重力异常略复杂的一点。
程序实现：
from geoist.pfm import prism, pftrans
# Reduce to the pole using FFT. Since there is only induced magnetization, the
# magnetization direction (sinc and sdec) is the same as the geomagnetic field
pole = pftrans.reduce_to_pole(x, y, tf, shape, inc, dec, sinc=inc, sdec=dec)
# Calculate the true value at the pole for comparison
true = prism.tf(x, y, z, model, 90, 0, pmag=giutils.ang2vec(10, 90, 0))
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax in axes:
ax.set_aspect('equal')
plt.sca(axes[0])
plt.title("Original total field anomaly")
giplt.contourf(y, x, tf, shape, 30, cmap=plt.cm.RdBu_r)
plt.colorbar(pad=0).set_label('nT')
giplt.m2km()
plt.sca(axes[1])
plt.title("True value at pole")
giplt.contourf(y, x, true, shape, 30, cmap=plt.cm.RdBu_r)
plt.colorbar(pad=0).set_label('nT')
giplt.m2km()
plt.sca(axes[2])
plt.title("Reduced to the pole")
giplt.contourf(y, x, pole, shape, 30, cmap=plt.cm.RdBu_r)
plt.colorbar(pad=0).set_label('nT')
giplt.m2km()
plt.tight_layout()
plt.show()
一句话总结：最好的学习方式就是看代码。从代码中学习，感受编码人的思路和逻辑，不看文档就能猜出来函数如何调用，这时候其实你就在不知不觉之间进步了。
