行星科学中的布格异常分析：理论推导与pyshtools实现

本文链接：https://blog.csdn.net/weixin_49036863/article/details/146512902

行星科学中的布格异常分析：理论推导与pyshtools实现

一、引言

布格异常是行星科学和地球物理学中的重要概念，它通过消除地形引起的重力效应，帮助我们探索行星内部结构。本文将从基础物理原理出发，逐步推导布格异常的计算公式，并结合pyshtools库展示实际的计算过程。

二、重力场理论基础

2.1 重力势能

任何质量体都会产生引力场。对于一个质点，引力势能为：

$\frac{Gm}{r}$

其中， $G$ 是万有引力常数， $m$ 是质量， $r$ 是距离。

对于形状不规则的行星，我们需要将其视为无数质点的集合，整体引力势能通过积分得到：

$G\iiint_{\text{body}} \frac{\rho(\vec{r'})}{|\vec{r}-\vec{r'}|} dV'$

其中， $\rho(\vec{r'})$ 是行星内部的密度分布函数。

2.2 球谐展开

因为上述积分难以直接求解，我们采用球谐函数展开法表示重力势能：

$V(r,\theta,\phi) = \frac{GM}{r}\left[1 + \sum_{l=1}^{\infty}\sum_{m=-l}^{l}\left(\frac{R_0}{r}\right)^l C_{lm}Y_{lm}(\theta,\phi)\right]$

其中：

$GM$ 是行星的引力参数（质量与引力常数的乘积）
$R_0$ 是参考半径
$C_{lm}$ 是球谐系数
$Y_{lm}$ 是球谐函数

推导说明：
球谐展开本质上是将势能函数在球面上分解为正交基函数的线性组合。这类似于傅里叶变换在一维情况下的应用，但扩展到了球面坐标系。这种分解使我们能够从重力场测量中推断行星内部的质量分布。

在pyshtools中，重力场模型通常以球谐系数的形式提供：

# 加载行星重力场模型
grav_model = pysh.datasets.Mars.MRO120F()

# 查看模型属性
print(f"GM = {grav_model.gm} m³/s²")
print(f"参考半径 R₀ = {grav_model.r0} m")
print(f"最大球谐阶数 = {grav_model.lmax}")

2.3 重力加速度

重力加速度是势能的梯度：

$\vec{g} = -\nabla V$

在球坐标系中，重力加速度有三个分量：

$g_r = -\frac{\partial V}{\partial r}$

$g_\theta = -\frac{1}{r}\frac{\partial V}{\partial \theta}$

$g_\phi = -\frac{1}{r\sin\theta}\frac{\partial V}{\partial \phi}$

将势能的球谐展开代入并求导，可以得到每个分量的具体表达式。例如，径向分量为：

$g_r = \frac{GM}{r^2}\left[1 + \sum_{l=1}^{\infty}\sum_{m=-l}^{l}(l+1)\left(\frac{R_0}{r}\right)^l C_{lm}Y_{lm}(\theta,\phi)\right]$

在pyshtools中，重力加速度计算如下：

# 计算重力加速度并展开为网格
grav_grid = grav_model.expand(a=a, f=f)

# grav_grid包含三个分量和总重力
# grav_grid.rad - 径向分量
# grav_grid.theta - 纬向分量
# grav_grid.phi - 经向分量
# grav_grid.total - 总重力

expand()方法内部执行了从球谐系数到空间域重力加速度的转换，实现了上述公式的计算。

三、重力异常概念

3.1 正常重力

要理解重力异常，首先需要定义"正常重力"。正常重力是参考椭球体上预期的重力加速度，通常用Somigliana公式计算：

$\gamma = \gamma_e\frac{1 + k\sin^2\phi}{\sqrt{1 - e^2\sin^2\phi}}$

其中：

$\gamma_e$ 是赤道处的正常重力
$k$ 是与椭球扁率和自转有关的常数
$e$ 是椭球的第一偏心率
$\phi$ 是大地纬度

参考椭球体的引力场通常只包含偶数阶、零阶的球谐系数（ $J_2$ , $J_4$ , $J_6$ 等）。

3.2 自由空气重力异常

自由空气重力异常是实际测量的重力与正常重力之间的差异：

$\Delta g_{fa} = g_{obs} - \gamma$

这种异常反映了地形高低和地下质量分布共同的影响。在pyshtools中：

# 计算自由空气重力异常
# normal_gravity=True表示从总重力中减去正常重力
grav_grid = grav_model.expand(a=a, f=f, normal_gravity=True)

# 显示自由空气重力异常图
plt.figure(figsize=(10, 6))
grav_grid.plot_total(cmap='RdBu_r', 
                     title='自由空气重力异常')
plt.show()

3.3 布格异常的概念

布格异常进一步消除了地形质量对重力的影响：

$\Delta g_b = \Delta g_{fa} - \Delta g_{topo}$

其中 $\Delta g_{topo}$ 是地形产生的重力效应。这一减法过程相当于"看穿"地表地形，直接观察深部结构。

四、布格异常计算方法

4.1 简单布格板模型

最简单的地形校正使用无限布格板模型：

$\Delta g_{topo} = 2\pi G \rho h$

其中：

$G$ 是引力常数
$\rho$ 是假设的地壳密度
$h$ 是地形高度

推导：
无限平板的引力可以通过积分得到：
$2\pi G \rho \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{h}{(x^2+y^2+h^2)^{3/2}}dxdy = 2\pi G \rho h$

这种模型虽然简单，但忽略了地形的曲率和有限范围，精度有限。

4.2 有限幅度方法

更精确的方法是有限幅度法，考虑地形的实际形状：

$\Delta V_{topo} = \frac{GM}{r}\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\left(\frac{R_0}{r}\right)^l\frac{\rho}{\bar{\rho}}\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{h}{r}\right)^n C_{lm}^{(n)}Y_{lm}(\theta,\phi)$

其中 $C_{lm}^{(n)}$ 是地形高度的n次方的球谐系数。

推导说明：
这个公式源于将引力势表达式中的 $(R + h) / r$ 项展开为泰勒级数，其中 $h$ 是地形相对于参考球面的高度。每一项 $h^n$ 都需要进行球谐展开，因此最终表达式包含多重求和。

在pyshtools中，有限幅度法通过from_shape方法实现：

# 计算地形产生的重力效应（有限幅度法）
topo_grav = pysh.SHGravCoeffs.from_shape(
    topo_model,  # 地形模型
    rho=2800.0,  # 假设密度为2800 kg/m³
    gm=grav_model.gm,  # GM值
    nmax=8  # 泰勒级数阶数，默认为8
)

# 调整参考半径以匹配重力场模型
topo_grav = topo_grav.change_ref(r0=grav_model.r0)

from_shape方法内部调用了CilmPlusDH函数，该函数实现了泰勒级数展开的有限幅度计算。nmax参数控制泰勒级数的阶数，值越大计算越精确，但计算量也越大。

4.3 布格异常计算步骤

完整的布格异常计算过程如下：

获取重力场球谐系数（通常来自卫星测量）
获取地形模型的球谐系数
计算自由空气重力异常
使用有限幅度法计算地形产生的重力效应
从自由空气重力异常中减去地形效应，得到布格异常

在pyshtools中实现如下：

# 从自由空气重力中减去地形效应，得到布格异常
bouguer_coeffs = grav_model.copy()  # 复制重力场系数
bouguer_coeffs.coeffs -= topo_grav.coeffs  # 减去地形效应

# 将布格异常展开为网格
bouguer_grid = bouguer_coeffs.expand(a=a, f=f)

# 显示布格异常图
plt.figure(figsize=(10, 6))
bouguer_grid.plot_total(cmap='RdBu_r', 
                        title='布格异常')
plt.show()

五、参考椭球体的处理

5.1 相对于球体的布格异常

最简单的布格异常计算是相对于球体的，但这会导致结果中包含行星扁率的影响。

5.2 相对于椭球体的布格异常

更精确的方法是计算相对于椭球体的布格异常，需要额外添加椭球体的引力效应：

# 创建参考椭球体
ellipsoid = pysh.SHGrid.from_ellipsoid(
    lmax=grav_model.lmax,
    a=a,  # 赤道半径
    b=a*(1-f)  # 极半径
)

# 计算椭球体产生的重力效应
ellipsoid_grav = pysh.SHGravCoeffs.from_shape(
    ellipsoid,
    rho=2800.0,  # 使用与地形相同的密度
    gm=grav_model.gm
)
ellipsoid_grav = ellipsoid_grav.change_ref(r0=grav_model.r0)

# 计算相对于椭球体的布格异常
bouguer_ellipsoid = bouguer_coeffs.copy()
bouguer_ellipsoid.coeffs += ellipsoid_grav.coeffs

# 展开为网格
bouguer_ellipsoid_grid = bouguer_ellipsoid.expand(a=a, f=f)

# 显示相对于椭球体的布格异常
plt.figure(figsize=(10, 6))
bouguer_ellipsoid_grid.plot_total(
    cmap='RdBu_r',
    title='相对于椭球体的布格异常'
)
plt.show()

六、地壳厚度估计

6.1 理论基础

布格异常与地壳厚度变化之间存在反相关关系。假设布格异常主要由地壳-地幔界面（莫霍面）的起伏引起，我们可以从布格异常推导地壳厚度：

$h(\theta,\phi) = h_0 + \frac{R^2}{4\pi G\Delta\rho}\Delta g_b(\theta,\phi)$

其中：

$h_0$ 是平均地壳厚度
$R$ 是行星半径
$\Delta\rho$ 是地壳与地幔的密度差
$\Delta g_b$ 是布格异常

推导说明：
这个公式基于阿里定理，假设地壳-地幔界面的起伏是引起布格异常的主要原因。通过球谐展开和势能表达式，可以推导出布格异常与界面起伏之间的线性关系。

6.2 有限幅度反演

更精确的方法是使用有限幅度反演。pyshtools提供了专门的函数BAtoHilmDH，实现从布格异常到界面起伏的反演：

# 创建初始地壳厚度模型
crustal_grid = topo_model.expand(grid='DH2')
crustal_grid.data = np.ones_like(crustal_grid.data) * 40000.0  # 40公里平均厚度

# 计算地壳厚度变化
nmax = 8  # 泰勒级数阶数
delta_rho = 2800.0 - 3300.0  # 地壳-地幔密度差

crust_coeffs = pysh.shtools.BAtoHilmDH(
    bouguer_coeffs.coeffs,  # 布格异常系数
    crustal_grid.data,      # 初始界面深度
    nmax,                   # 泰勒级数阶数
    grav_model.gm,          # GM
    grav_model.r0,          # 参考半径
    abs(delta_rho),         # 密度对比绝对值
    filter_type=1,          # 最小振幅滤波
    filter_deg=30,          # 滤波阶数
    lmax=60                 # 最大计算阶数
)

# 创建地壳厚度模型
crust_model = pysh.SHCoeffs.from_array(crust_coeffs)
crust_thickness = crust_model.expand()

# 显示地壳厚度图
plt.figure(figsize=(10, 6))
crust_thickness.plot(
    cmap='plasma',
    title='估计的地壳厚度'
)
plt.show()

BAtoHilmDH函数实现了一个迭代过程，考虑了有限幅度效应，使估计更加精确。函数的关键参数是：

nmax：泰勒级数阶数
delta_rho：密度对比
filter_type和filter_deg：用于稳定反演过程的滤波参数

七、完整代码示例

以下是一个完整的布格异常分析代码示例，从理论到实践，一步步实现布格异常计算和地壳厚度估计：

import numpy as np
import matplotlib.pyplot as plt
import pyshtools as pysh
from pyshtools import constants

# 设置图形风格
pysh.utils.figstyle(rel_width=0.75)

# 1. 加载数据
print("正在加载重力和地形数据...")
grav_model = pysh.datasets.Mars.MRO120F()  # 火星重力场模型
topo_model = pysh.datasets.Mars.MOLA_shape()  # 火星地形模型

# 2. 设置行星参数
print("设置火星物理参数...")
a = constants.Mars.a.value  # 赤道半径
f = constants.Mars.f.value  # 扁率
grav_model.set_omega(constants.Mars.omega.value)  # 自转角速度

# 3. 计算自由空气重力异常
print("计算自由空气重力异常...")
grav_grid = grav_model.expand(a=a, f=f, lmax=85)

# 4. 计算地形引起的重力效应
print("计算地形引起的重力效应...")
topo_grav = pysh.SHGravCoeffs.from_shape(
    topo_model,
    rho=2800.0,  # 假设密度为2800 kg/m³
    gm=grav_model.gm,
    lmax=85
)
topo_grav = topo_grav.change_ref(r0=grav_model.r0)

# 5. 计算布格异常
print("计算布格异常...")
bouguer_coeffs = grav_model.copy()
bouguer_coeffs.coeffs -= topo_grav.coeffs
bouguer_grid = bouguer_coeffs.expand(a=a, f=f)

# 6. 计算相对于椭球体的布格异常
print("计算相对于椭球体的布格异常...")
ellipsoid = pysh.SHGrid.from_ellipsoid(
    lmax=85,
    a=a,
    b=a*(1-f)
)
ellipsoid_grav = pysh.SHGravCoeffs.from_shape(
    ellipsoid,
    rho=2800.0,
    gm=grav_model.gm,
    lmax=85
)
ellipsoid_grav = ellipsoid_grav.change_ref(r0=grav_model.r0)
bouguer_ellipsoid = bouguer_coeffs.copy()
bouguer_ellipsoid.coeffs += ellipsoid_grav.coeffs
bouguer_ellipsoid_grid = bouguer_ellipsoid.expand(a=a, f=f)

# 7. 估计地壳厚度
print("估计地壳厚度...")
crustal_grid = topo_model.expand(grid='DH2')
crustal_grid.data = np.ones_like(crustal_grid.data) * 40000.0  # 40公里平均厚度

crust_coeffs = pysh.gravmag.BAtoHilmDH(
    bouguer_coeffs.coeffs,
    crustal_grid.data,
    8,  # nmax
    grav_model.gm,
    grav_model.r0,
    500.0,  # |Δρ|
    filter_type=1,
    filter_deg=30,
    lmax=60
)

crust_model = pysh.SHCoeffs.from_array(crust_coeffs)
crust_thickness = crust_model.expand()

# 8. 可视化结果
print("生成可视化结果...")
fig = plt.figure(figsize=(15, 12))

plt.subplot(2, 2, 1)
grav_grid.plot_total(
    cmap='RdBu_r',
    cmap_limits=[-400, 600],
    title='自由空气重力异常'
)

plt.subplot(2, 2, 2)
bouguer_grid.plot_total(
    cmap='RdBu_r',
    title='布格异常(相对于球体)'
)

plt.subplot(2, 2, 3)
bouguer_ellipsoid_grid.plot_total(
    cmap='RdBu_r',
    title='布格异常(相对于椭球体)'
)

plt.subplot(2, 2, 4)
crust_thickness.plot(
    cmap='plasma',
    title='估计的地壳厚度(米)'
)

plt.tight_layout()
plt.savefig('bouguer_analysis.png', dpi=300)
plt.show()

print("布格异常分析完成!")