harmonica

Harmonica - 重力与磁法

快速参考

import harmonica as hm
import numpy as np

# Forward model - prism gravity
prism = [-500, 500, -500, 500, -2000, -500]  # (west, east, south, north, bottom, top)
gravity = hm.prism_gravity(coordinates, prism, density=500, field='g_z')

# Terrain correction
layer = hm.prism_layer((easting, northing), surface=topo, reference=0,
                        properties={'density': 2670})
terrain_effect = layer.gravity(coordinates, field='g_z')

# Equivalent source gridding
eqs = hm.EquivalentSources(depth=10000, damping=10)
eqs.fit(coordinates, gravity_data)
grid = eqs.grid(spacing=5000, data_names=['gravity'])

# Upward continuation (requires gridded xarray)
upward = hm.upward_continuation(gravity_grid, height_displacement=1000)

主要函数

函数	用途
point_gravity	计算点质量产生的重力
prism_gravity	计算长方体棱柱产生的重力
tesseroid_gravity	计算球面棱柱产生的重力（区域/全球）
prism_magnetic	计算棱柱产生的磁异常
prism_layer	根据地形创建棱柱层
EquivalentSources	使用等效源对离散数据进行网格化
upward_continuation	基于 FFT 的向上延拓
bouguer_correction	简单布格板改正

基本操作

正演模型 - 长方体棱柱

# Define prism: (west, east, south, north, bottom, top) in meters
prism = [-500, 500, -500, 500, -2000, -500]
density = 500  # kg/m3 density contrast

# Observation grid
x_obs, y_obs = np.meshgrid(np.linspace(-5000, 5000, 100), np.linspace(-5000, 5000, 100))
z_obs = np.zeros_like(x_obs)

# Calculate gravity (mGal). Fields: 'g_z', 'g_north', 'g_east', 'potential'
gravity = hm.prism_gravity((x_obs.ravel(), y_obs.ravel(), z_obs.ravel()),
                           prism, density, field='g_z')

地形改正

import xarray as xr

topo = xr.open_dataarray('dem.nc')
layer = hm.prism_layer((topo.easting.values, topo.northing.values),
                       surface=topo.values, reference=0,
                       properties={'density': 2670})
terrain_effect = layer.gravity((obs_easting, obs_northing, obs_height), field='g_z')
bouguer_anomaly = free_air_anomaly - terrain_effect

等效源网格化

import verde as vd

# Project to Cartesian
projection = vd.get_projection(longitude, latitude)
easting, northing = projection(longitude, latitude)

eqs = hm.EquivalentSources(depth=10000, damping=10)
eqs.fit((easting, northing, altitude), gravity_mgal)
grid = eqs.grid(spacing=5000, data_names=['gravity'])

磁法正演模型

prism = [-500, 500, -500, 500, -2000, -500]
magnetization = hm.magnetic_vector(intensity=5.0, inclination=60, declination=10)
b_total = hm.prism_magnetic(coordinates, prism, magnetization, field='b_total')

导数滤波器

dx = hm.derivative_easting(gravity_grid)
dy = hm.derivative_northing(gravity_grid)
dz = hm.derivative_upward(gravity_grid)

thg = np.sqrt(dx**2 + dy**2)  # Total horizontal gradient
tilt = np.arctan2(dz, thg)     # Tilt angle

坐标系

Harmonica 使用右手坐标系：

• 东向（x）：向东为正
• 北向（y）：向北为正
• 向上（z）：向上为正（高度为正，深度为负）

单位采用 SI 制：距离单位为米，密度单位为 kg/m3，重力单位为 mGal。

适用场景与替代方案

使用场景	工具	原因
重力/磁法正演建模	Harmonica	专为此用途打造，属于 Fatiando 生态系统
位场反演	SimPEG	具备正则化功能的完整反演框架
商业重力数据处理	Oasis Montaj	行业标准图形界面，使用专有格式
仅需简单布格改正	Custom numpy	一次性计算所需的依赖更少
等效源网格化	Harmonica	位场领域的最佳开源选项
区域/全球尺度	Harmonica（tesseroids）	原生支持球面几何
磁数据化极	Harmonica	针对网格数据的基于 FFT 的滤波器
教学/原型开发	Harmonica	API 简洁，文档完善

选择 Harmonica 的情况：需要使用开源工具处理重力/磁法数据，包括正演建模、地形改正或等效源网格化。它与 Verde 配合使用，可进行投影和网格化。Harmonica 属于 Fatiando a Terra 生态系统。

选择 SimPEG 的情况：需要对位场数据进行反演，以求得地下物性分布（密度或磁化率模型）。

选择 Oasis Montaj 的情况：在行业环境中工作，需要使用专有格式、商业支持或基于图形界面的交互式处理。

常见工作流

进行重力测量数据处理、地形改正与网格化

• 加载原始重力观测数据和测站坐标
• 进行纬度、自由空气和潮汐改正
• 加载 DEM，并使用 hm.prism_layer() 构建棱柱层
• 计算观测点处的地形效应
• 从自由空气异常中减去地形效应，得到布格异常
• 使用 verde.get_projection() 将坐标投影到笛卡尔坐标系
• 如果测站密度不均，进行分块降采样
• 使用 hm.EquivalentSources() 拟合等效源
• 将布格异常网格化到规则网格上
• 应用 vd.distance_mask() 屏蔽远离数据的区域
• 应用导数滤波器（水平梯度、倾斜角）进行解释
• 进行向上延拓以突出区域特征
• 将网格数据导出为 NetCDF

常见问题

问题	解决方案
重力符号错误	检查 z 轴约定（向上为正）
等效源拟合效果不佳	调整 depth 和 damping 参数
地形改正速度慢	降低 DEM 分辨率或使用更大的棱柱
FFT 滤波器出现边缘效应	应用 upward_continuation 前先对网格进行填充
坐标不匹配	确保统一使用投影坐标或地理坐标

提示

1. 使用投影坐标（米）进行局部测量
2. 使用曲球体进行区域／全球尺度建模
3. 等效源能很好地处理间距不规则的数据
4. 选择合适的密度（上地壳的典型值为 2670 kg/m3）
5. 检查符号约定——深度为负的 z 值

参考资料

• 重力／磁力改正 - 标准改正和异常
• 正演建模 - 详细的正演建模方法

脚本

• scripts/gravity_processing.py - 处理重力测量数据