尧图网页设计 
简介这是一份面向地球物理勘探、航空磁测及地磁场建模方向学习者的技术文档原题为《利用IGRF模型计算全张量地磁梯度》以PDF单文件形式提供压缩包约2.03MB共1个文件适合需要查阅球谐分析、勒让德多项式与地磁张量计算细节的研究人员与研究生。文档梳理了国际地磁参考场IGRF的计算原理推导出球谐展开的全张量磁梯度表达式给出任意点位地磁场七要素与全张量的计算流程并用NOAA公开数据对七要素结果进行对照验证说明算法准确可靠。文中还绘制了某地区地磁全张量磁梯度等值线图结果满足Laplace方程可为航空全张量磁梯度测量中选取学习飞行工区和飞行高度提供理论依据。目前已有313人学习适合作为算法复现、公式推导和工程应用参考。1. 从地磁七要素到全张量梯度航空磁测到底缺哪一块数据航空全张量磁梯度测量FTMG在选学习飞行工区时最先要回答的问题不是仪器噪声多少而是这片空域的主磁场本身有多陡。IGRF 这类国际地磁参考场模型一直只解决七要素北向、东向、垂向三分量再加水平强度、总场、偏角、倾角。要拿梯度过去只能对相邻网格点做有限差分0.1° 的网格间距下差分截断误差和边界点的处理方式直接决定梯度图有没有物理意义。这篇工作做的是另一条路从球谐磁位的二阶偏导数出发把全张量梯度的六个独立分量写成解析表达式再经过两级坐标旋转落到沿椭球面法线向下的北东下局部坐标系。4°×4° 工区、0.1° 网格、大地高 1 km 的算例里六个分量之和落在 −0.00110.0011 nT/km等价于数值上满足 Laplace 方程四个城市点与 NOAA 在线结果的偏差最大 0.1 nT 和 0.0001°。做航磁测量设计、地磁数据处理的同行照着这套推导加递推实现就能自己算出任意点位的张量值。2. 磁位球谐展开到张量分量的推导与坐标系变换2.1 内源场磁位与拉普拉斯方程的解球谐函数法把地球当成均匀球体在球坐标系 (r, θ, φ) 下场源区外任意一点的磁位满足拉普拉斯方程。用分离变量法解这个方程的 Dirichlet 问题得到磁位的级数形式其中含1/r^(n1)因子的那部分是内源场含r^n因子的是外源场。地磁场总强度里 99% 以上来自内源场所以工程实现只保留第一项就够了磁位写成V(r, θ, φ) a · Σ(n1..k) (a/r)^(n1) · Σ(m0..n) (g_n^m·cos(mφ) h_n^m·sin(mφ)) · P_n^m(cosθ)这里 a 是参考球半径取 6371.2 kmk 是球谐展开的最大阶数g_n^m 和 h_n^m 就是 IAGA 每五年发布一次的高斯球谐系数也是 IGRF 模型的本体。P_n^m(cosθ) 是施密特半规格化伴随勒让德多项式不是普通勒让德多项式规格化方式选错后面所有分量的量级都会整体偏掉。2.2 一阶偏导给出地磁三分量磁场强度定义为磁位的负梯度在球坐标系下分别对 r、θ、φ 求偏导再取负就得到三个分量。径向分量每项多出一个 (n1) 因子θ 方向分量带着勒让德多项式对 θ 的一阶导数φ 方向分量则把 cos(mφ)、sin(mφ) 换成 −sin(mφ)、cos(mφ) 并乘上 m。这一步的实现难度全在勒让德多项式导数的递推上公式本身不长但 θ 接近 0 或 π 时 1/sinθ 会炸是后面坐标变换里最容易出 NaN 的地方。2.3 二阶偏导与 Laplace 约束对三个一阶分量继续求 r、θ、φ 的偏导得到九个二阶偏导量Uθθ、Uφθ、Urθ、Uφφ、Urφ、Urr加上径向和 θ 方向的一阶项参与组合。物理上磁张量是对称的且无源区里迹为零所以九个分量里只有六个独立。这个约束非常好用它同时是正确性的自检条件算完之后把主对角线三项 Bxx、Byy、Bzz 相加如果结果远大于计算量级说明递推、规格化或坐标变换里至少有一处错了。二阶偏导量是否含勒让德二阶导主要量级来源Uθθ是一阶导数项主导低阶系数占大头Uφθ否只含一阶导带 m 因子赤道附近最强Urθ否只含一阶导带 (n1)/r 因子Uφφ否只含零阶带 m² 因子高阶衰减慢Urφ否只含零阶带 m(n1)/r 因子Urr否只含零阶带 (n1)(n2)/r 因子2.4 从地心球坐标到北东下局部坐标的三级变换算出来的量都在地心球坐标系里实际测量要的是局部直角坐标。变换分三级先用 Balmino 给出的关系把球坐标分量换成北东上坐标系的 Bx、By、Bz 和六个张量分量再把 z 轴反向得到右手定则下更常用的北东下坐标系最后绕 y 轴旋转一个角度 d把 z 轴从指向地心转到沿旋转椭球面法线向下。d 是地心余纬度 θ 与地理余纬度 θ′ 的差这一步不做中高纬度点位的水平分量会出现零点几度的方向偏差。大地坐标到地心坐标的换算原文给的是式 (18)(20) 的解析式。工程上我一般用卯酉圈曲率半径 N 走一遍直角坐标再反算写法短、不容易抄错import numpy as np def geodetic_to_geocentric(lat_deg, h_km): 大地纬度(度)、大地高(km) - 地心距 r、地心余纬度 theta、旋转角 d(弧度) a 6378.137 # 椭球长半轴 km f 1.0 / 298.257223563 # 地球扁率 e2 f * (2.0 - f) # 第一偏心率平方 lat np.radians(lat_deg) N a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2) # 卯酉圈曲率半径 x (N h_km) * np.cos(lat) # 赤道面内投影 z (N * (1.0 - e2) h_km) * np.sin(lat) r np.hypot(x, z) theta np.arctan2(x, z) # 地心余纬度 d theta - (np.pi / 2.0 - lat) # 旋转角 return r, theta, d def rotate_to_ellipsoid_normal(bx, by, bz, d): 把北东下分量绕 y 轴旋转 d使 z 轴沿椭球面法线向下 cd, sd np.cos(d), np.sin(d) bx2 cd * bx sd * bz by2 by bz2 -sd * bx cd * bz return bx2, by2, bz2geodetic_to_geocentric里 a 取 6378.137 km、f 取 1/298.257223563与原文一致返回的 θ 是地心余纬度不是纬度后面所有勒让德多项式都吃这个角。rotate_to_ellipsoid_normal只处理三分量张量要走 R·T·Rᵀ 的两侧乘法写成T2 R T R.T更省事。注意旋转角 d 在高纬地区会明显变大如果只做单侧旋转只乘 R 不乘 R张量会出现不对称用np.allclose(T2, T2.T)一查就露馅。3. 施密特半规格化勒让德多项式的一阶二阶递推实现3.1 为什么不直接用 scipy 的伴随勒让德函数scipy.special.lpmn算的是未经规格化的伴随勒让德多项式而且 m 从 0 开始、含 Condon-Shortley 相位约定与 IGRF 用的施密特半规格化差着一整套归一化因子。直接拿来乘高斯系数低阶项还能凑合看到 n8 以上数值会迅速失衡因为未规格化的 P_n^m 在高阶时量级跨度极大求和的相对重要性完全被打乱。更麻烦的是全张量梯度需要 P 对 θ 的一阶和二阶导数scipy 不直接给用差分近似会把 Laplace 自检毁掉。自己写递推还有个工程理由递推出来的 P 表在同一纬度圈上是复用的。0.1° 网格、4° 工区只有 41 个不同的纬度值经度方向 41 个点可以共用同一张 P 表和三张导数表这比逐点重算快一个数量级。3.2 对角项与非对角项的递推关系施密特半规格化下对角项递推是 P_n^n √(1 − 1/(2n)) · sinθ · P_(n−1)^(n−1)非对角项用三项递推P_n^m a_nm · cosθ · P_(n-1)^m - b_nm · P_(n-2)^m a_nm sqrt( ((2n-1)^2 - m^2) / (n^2 - m^2) ) b_nm sqrt( ((n-1)^2 - m^2) / (n^2 - m^2) )一阶导数只要对 θ 用乘积法则展开二阶导数再展开一次。原文式 (37)、(38) 给的就是这两个展开结果。实现时要注意 n1、m0 这种起点位置b_nm 的分子里会出现 (n−1)²−m² 为负的情况必须显式置零否则开根号直接出 NaN。阶数 k待求项数量含 m≤n二阶导项数量典型用途62828快速估算场形态106666一般区域磁场建模13105105与 NOAA 在线结果对齐16 以上153153岩石圈异常细分需谨慎3.3 Python 实现与数值稳定性处理下面这段是核心三个数组一次填满供后面的分量求和直接索引import numpy as np def schmidt_legendre(theta, K): 施密特半规格化伴随勒让德多项式及对 theta 的一阶、二阶导数 theta: 地心余纬度(弧度) K: 最大阶数 返回 P, dP, dP2形状均为 (K1, K1)索引 [n][m] st, ct np.sin(theta), np.cos(theta) P np.zeros((K 1, K 1)) dP np.zeros((K 1, K 1)) dP2 np.zeros((K 1, K 1)) P[0, 0] 1.0 # 递推起点 for n in range(1, K 1): if n 1: P[1, 1] st dP[1, 1] ct dP2[1, 1] -st else: c np.sqrt(1.0 - 1.0 / (2.0 * n)) P[n, n] c * st * P[n - 1, n - 1] dP[n, n] c * (st * dP[n - 1, n - 1] ct * P[n - 1, n - 1]) dP2[n, n] c * (-st * P[n - 1, n - 1] 2.0 * ct * dP[n - 1, n - 1] st * dP2[n - 1, n - 1]) for m in range(0, n): a np.sqrt(((2 * n - 1) ** 2 - m ** 2) / (n ** 2 - m ** 2)) # m n-1 时 P[n-2][m] 不存在系数置 0 b 0.0 if m n - 2 else np.sqrt(((n - 1) ** 2 - m ** 2) / (n ** 2 - m ** 2)) p1, p0 P[n - 1, m], P[n - 2, m] if m n - 2 else 0.0 d1, d0 dP[n - 1, m], dP[n - 2, m] if m n - 2 else 0.0 e1, e0 dP2[n - 1, m], dP2[n - 2, m] if m n - 2 else 0.0 P[n, m] a * ct * p1 - b * p0 dP[n, m] a * (ct * d1 - st * p1) - b * d0 dP2[n, m] a * (-2.0 * st * d1 ct * e1 - ct * p1) - b * e0 return P, dP, dP2theta必须是弧度从上一节的geodetic_to_geocentric拿。K是截断阶数与高斯系数的最大阶一致。二阶导数组 dP2 只在 Bxx 和 Bzz 的求和里用到如果只做三分量可以传一个开关跳过这段计算省时间。m n−2时把 b 置零是必须的否则P[n-2, m]会取到零行甚至负索引结果看着像对的其实是脏数据。3.4 高斯系数的读取与时间归算IGRF 系数文件按行给出 n、m、g、h 以及各年代的值和末段的年变率。给定日期 t 时落在两个publihed epoch 之间就线性插值落在最后一个 epoch 之后就用年变率外推常见写法是 g(t) g(t0) (t − t0)·。解析这类文本文件时注意三点h_n^0 恒为零、m0 的项不参与 φ 方向的求和、单位是 nT 而不是 T。系数表拿到手第一件事是抽查 n1 那一行的量级g_1^0 在 −29000 附近才是正常值差三个数量级说明单位或者符号读反了。4. 4°×4° 工区网格计算与 NOAA 逐点校验4.1 工区、网格与高度参数原文算例选的是经度 103.3056°107.3056°、纬度 27.3056°31.3056° 的 4°×4° 区域网格步长 0.1°大地高 1 km时间取 2019 年 4 月 7 日。这几个参数不是随便定的网格步长要小到能分辨张量等值线的形态又不能小到让等值线图被计算噪声糊住高度取 1 km 是因为航空测量本身有一个巡航高度同时主磁场梯度随高度衰减拿 1 km 的结果做学习飞行的高度设计比较贴近实际。4.2 批量计算的组织方式按纬度分组是提高效率的关键。同一个纬度上geodetic_to_geocentric的 r、θ、d 只算一次schmidt_legendre的三张表也只算一次经度方向的循环只更新 cos(mφ)、sin(mφ) 和 (a/r)^(n1) 这类标量。def compute_tensor_at_point(lat, lon, h, t, coeffs, K): 返回单点的七要素与六个独立张量分量(北东下) r, theta, d geodetic_to_geocentric(lat, h) phi np.radians(lon) P, dP, dP2 schmidt_legendre(theta, K) Ur Uth Uph 0.0 Uthth Uphth Urth Uphph Urph Urr 0.0 for n in range(1, K 1): ratio (6371.2 / r) ** (n 1) for m in range(0, n 1): g, hc coeffs[n][m] cm, sm np.cos(m * phi), np.sin(m * phi) gc_hs g * cm hc * sm # g·cos h·sin mgh m * (-g * sm hc * cm) # m·(-g·sin h·cos) Ur (n 1) * ratio / r * gc_hs * P[n, m] Uth - 6371.2 * ratio * gc_hs * dP[n, m] Uph - 6371.2 * ratio * mgh * P[n, m] Uthth - 6371.2 * ratio * gc_hs * dP2[n, m] Uphth - 6371.2 * ratio * mgh * dP[n, m] Urth (n 1) * ratio / r * gc_hs * dP[n, m] Uphph 6371.2 * ratio * m * mgh * P[n, m] / m if m else 0.0 Urph (n 1) * ratio / r * mgh * P[n, m] Urr - (n 2) * (n 1) * ratio / (r * r) * gc_hs * P[n, m] st, ct np.sin(theta), np.cos(theta) bx, by, bz -Uth / r, Uph / (r * st), Ur bxx Ur / r Uthth / r ** 2 bxy ct / (r ** 2 * st ** 2) * Uph - Uphth / (r ** 2 * st) bxz Uth / r ** 2 - Urth / r byy Ur / r Uth / (r ** 2 * np.tan(theta)) Uphph / (r ** 2 * st ** 2) byz Urph / (r * st) - Uph / (r ** 2 * st) bzz Urr return (bx, by, bz), (bxx, bxy, bxz, byy, byz, bzz)注意Uphph那一行在 m0 时不成立代码里用if m else 0.0挡掉实际写的时候把 m0 单独拎出去更清楚。奇点保护也要加θ 小于 1e-6 或者 sinθ 接近 0 时直接跳过该点或者把该点标记后单独用邻域值插补硬算出来的是无意义的巨大值。4.3 与 NOAA 七要素对照选取工区内四个点做对照本文算法与 NOAA 在线计算的结果差异如下Bx、By、Bz、Bh、B 单位 nTD、I 单位度点位来源BxByBzBhBDI成都 (30.67N,104.07E)本文计算33866.2−1213.937855.433887.950807.6−2.052848.1652成都NOAA33866.2−1213.937855.433887.950807.7−2.052948.1652自贡 (29.35N,104.78E)本文计算34644.5−1268.736051.734667.650015.7−2.097246.1212自贡NOAA34644.5−1268.736051.734667.750015.8−2.097346.1212泸州 (28.87N,105.43E)本文计算34909.9−1334.735359.934935.349707.1−2.189545.3460泸州NOAA34909.9−1334.735359.934935.449707.2−2.189545.3460德阳 (31.13N,104.38E)本文计算33580.9−1266.438443.833604.751060.8−2.159748.8424德阳NOAA33580.9−1266.438443.833604.851060.8−2.159748.8424三分量在保留一位小数的情况下完全一致水平强度、总场和偏角最大绝对误差分别是 0.1 nT、0.1 nT 和 0.0001°倾角无差异。这个量级的残差主要来自 NOAA 结果本身的有效位数以及大地坐标到地心坐标换算时的舍入不是推导错误。做校验时要注意统一时间拿 2019 年的算例去对 2020 年的在线结果残差会跳到几十 nT。4.4 Laplace 方程自检与张量对称性检查没有第二来源的全张量数据可比所以正确性只能靠物理约束兜底。主对角线三项相加应当恒为零原文算例里这个和落在 −0.00110.0011 nT/km相对张量分量本身的量级nT/km 到几十 nT/km已经算数值零。实际编码时按下面三件事查一遍张量矩阵对称性np.allclose(T, T.T, atol1e-6)不对称说明坐标旋转只做了一半。迹的分布图画成等值线如果出现明显的条纹或边界跳变多半是 θ 接近极点时的除零污染。阶数收敛k 取 8、10、13 各跑一遍张量分量的变化应该随 k 单调减小并趋于平稳如果 k 增大后反而发散递推的规格化因子写错了。4.5 容易翻车的几个点纬度单位弄错是最高频的错误。所有勒让德多项式的自变量是地心余纬度 θ等于 90° 减地心纬度而输入的大地纬度还要先转成地心纬度中间差一个旋转角 d。漏掉这一步在中纬度地区会带来零点几度的方向偏差映射到张量上就是百分之几的量级错误图看着还挺像回事。第二高频的是规格化方式混用。高斯规格化、施密特半规格化、完全规格化三套因子互不相同混着用会让高阶项和低阶项的相对权重彻底跑偏。判断方法很直接算 n1、m0 这一项施密特半规格化下 P_1^0 cosθ如果代码输出的是别的值说明归一化因子挂错了。第三是角度与弧度的边界。递推表内部的 θ、φ 都是弧度但网格生成、绘图和输入参数多半用度转换点最好只留在入口和出口各一处中间层全部用弧度避免在循环里反复np.radians。5. 大批量网格的向量化、阶数收敛与云端批处理5.1 递推表复用与数组化单点跑通之后真正的开销在网格上。4°×4°、0.1° 网格是 41×41 共 1681 个点单点循环在笔记本上还能忍换成全国范围 0.05° 网格就是几十万个点。两个优化收益最大按纬度分组的递推表缓存把 P、dP、dP2 的计算次数从点数降到纬度数以及把经度方向的求和数组化用np.einsum或者预先把 cos(mφ) 矩阵乘上去把内层 m 循环消掉。前者是纯算法层面的复用后者靠 NumPy 的向量化两者叠加通常能压掉一个数量级以上的时间。5.2 阶数收敛验证IGRF 的球谐系数阶数有限梯度对高阶项的敏感度比场值本身高因为每一项都带着 (n1)(n2)/r 这类放大因子。验证方式是把 k 从 6 递增到模型最大阶逐档记录工区内 Bxx、Byy、Bzz 的最大变化量截断阶数 kBxx 最大变化 (nT/km)说明6基准偶极子主导形态大体正确10较 k6 有可见下降四极、八极项贡献开始收敛13继续下降但幅度变小接近模型上限超过模型阶数不变系数为零加了也没用判断标准落在仪器噪声上如果 k 从 13 加到 16 带来的张量变化已经低于航空测量系统的噪声底继续加阶数只是白烧算力。反过来如果 k6 和 k13 的梯度图形态差异明显说明测试工区落在高阶异常较强的区域用低阶结果做飞行设计会低估梯度。5.3 参数固化成技术文档与云上批量跑整个计算链没有 IO 依赖输入是纬度、经度、高度、时间和系数表输出是六个数天然适合打包成容器放进云上按网格分块跑。我的做法是把椭球参数、参考球半径、最大阶数、时间归算方式写在配置里随代码一起提交同时把每个参数的含义、单位和典型取值范围整理成一页技术文档谁接手都能改对地方。需要注意的是容器里浮点行为要一致np.float64和某些加速库的默认精度混用会让 Laplace 自检的残留量从 1e-3 跳到 1e-1一批结果直接作废。跑批量任务前先在本地用四个对照点做一次回归通过以后再上云能省掉大量返工。本文还有配套的精品资源点击获取 IGRF球谐模型二阶偏导与航空全张量磁梯度计算
简介这是一份面向地球物理勘探、航空磁测及地磁场建模方向学习者的技术文档原题为《利用IGRF模型计算全张量地磁梯度》以PDF单文件形式提供压缩包约2.03MB共1个文件适合需要查阅球谐分析、勒让德多项式与地磁张量计算细节的研究人员与研究生。文档梳理了国际地磁参考场IGRF的计算原理推导出球谐展开的全张量磁梯度表达式给出任意点位地磁场七要素与全张量的计算流程并用NOAA公开数据对七要素结果进行对照验证说明算法准确可靠。文中还绘制了某地区地磁全张量磁梯度等值线图结果满足Laplace方程可为航空全张量磁梯度测量中选取学习飞行工区和飞行高度提供理论依据。目前已有313人学习适合作为算法复现、公式推导和工程应用参考。1. 从地磁七要素到全张量梯度航空磁测到底缺哪一块数据航空全张量磁梯度测量FTMG在选学习飞行工区时最先要回答的问题不是仪器噪声多少而是这片空域的主磁场本身有多陡。IGRF 这类国际地磁参考场模型一直只解决七要素北向、东向、垂向三分量再加水平强度、总场、偏角、倾角。要拿梯度过去只能对相邻网格点做有限差分0.1° 的网格间距下差分截断误差和边界点的处理方式直接决定梯度图有没有物理意义。这篇工作做的是另一条路从球谐磁位的二阶偏导数出发把全张量梯度的六个独立分量写成解析表达式再经过两级坐标旋转落到沿椭球面法线向下的北东下局部坐标系。4°×4° 工区、0.1° 网格、大地高 1 km 的算例里六个分量之和落在 −0.00110.0011 nT/km等价于数值上满足 Laplace 方程四个城市点与 NOAA 在线结果的偏差最大 0.1 nT 和 0.0001°。做航磁测量设计、地磁数据处理的同行照着这套推导加递推实现就能自己算出任意点位的张量值。2. 磁位球谐展开到张量分量的推导与坐标系变换2.1 内源场磁位与拉普拉斯方程的解球谐函数法把地球当成均匀球体在球坐标系 (r, θ, φ) 下场源区外任意一点的磁位满足拉普拉斯方程。用分离变量法解这个方程的 Dirichlet 问题得到磁位的级数形式其中含1/r^(n1)因子的那部分是内源场含r^n因子的是外源场。地磁场总强度里 99% 以上来自内源场所以工程实现只保留第一项就够了磁位写成V(r, θ, φ) a · Σ(n1..k) (a/r)^(n1) · Σ(m0..n) (g_n^m·cos(mφ) h_n^m·sin(mφ)) · P_n^m(cosθ)这里 a 是参考球半径取 6371.2 kmk 是球谐展开的最大阶数g_n^m 和 h_n^m 就是 IAGA 每五年发布一次的高斯球谐系数也是 IGRF 模型的本体。P_n^m(cosθ) 是施密特半规格化伴随勒让德多项式不是普通勒让德多项式规格化方式选错后面所有分量的量级都会整体偏掉。2.2 一阶偏导给出地磁三分量磁场强度定义为磁位的负梯度在球坐标系下分别对 r、θ、φ 求偏导再取负就得到三个分量。径向分量每项多出一个 (n1) 因子θ 方向分量带着勒让德多项式对 θ 的一阶导数φ 方向分量则把 cos(mφ)、sin(mφ) 换成 −sin(mφ)、cos(mφ) 并乘上 m。这一步的实现难度全在勒让德多项式导数的递推上公式本身不长但 θ 接近 0 或 π 时 1/sinθ 会炸是后面坐标变换里最容易出 NaN 的地方。2.3 二阶偏导与 Laplace 约束对三个一阶分量继续求 r、θ、φ 的偏导得到九个二阶偏导量Uθθ、Uφθ、Urθ、Uφφ、Urφ、Urr加上径向和 θ 方向的一阶项参与组合。物理上磁张量是对称的且无源区里迹为零所以九个分量里只有六个独立。这个约束非常好用它同时是正确性的自检条件算完之后把主对角线三项 Bxx、Byy、Bzz 相加如果结果远大于计算量级说明递推、规格化或坐标变换里至少有一处错了。二阶偏导量是否含勒让德二阶导主要量级来源Uθθ是一阶导数项主导低阶系数占大头Uφθ否只含一阶导带 m 因子赤道附近最强Urθ否只含一阶导带 (n1)/r 因子Uφφ否只含零阶带 m² 因子高阶衰减慢Urφ否只含零阶带 m(n1)/r 因子Urr否只含零阶带 (n1)(n2)/r 因子2.4 从地心球坐标到北东下局部坐标的三级变换算出来的量都在地心球坐标系里实际测量要的是局部直角坐标。变换分三级先用 Balmino 给出的关系把球坐标分量换成北东上坐标系的 Bx、By、Bz 和六个张量分量再把 z 轴反向得到右手定则下更常用的北东下坐标系最后绕 y 轴旋转一个角度 d把 z 轴从指向地心转到沿旋转椭球面法线向下。d 是地心余纬度 θ 与地理余纬度 θ′ 的差这一步不做中高纬度点位的水平分量会出现零点几度的方向偏差。大地坐标到地心坐标的换算原文给的是式 (18)(20) 的解析式。工程上我一般用卯酉圈曲率半径 N 走一遍直角坐标再反算写法短、不容易抄错import numpy as np def geodetic_to_geocentric(lat_deg, h_km): 大地纬度(度)、大地高(km) - 地心距 r、地心余纬度 theta、旋转角 d(弧度) a 6378.137 # 椭球长半轴 km f 1.0 / 298.257223563 # 地球扁率 e2 f * (2.0 - f) # 第一偏心率平方 lat np.radians(lat_deg) N a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2) # 卯酉圈曲率半径 x (N h_km) * np.cos(lat) # 赤道面内投影 z (N * (1.0 - e2) h_km) * np.sin(lat) r np.hypot(x, z) theta np.arctan2(x, z) # 地心余纬度 d theta - (np.pi / 2.0 - lat) # 旋转角 return r, theta, d def rotate_to_ellipsoid_normal(bx, by, bz, d): 把北东下分量绕 y 轴旋转 d使 z 轴沿椭球面法线向下 cd, sd np.cos(d), np.sin(d) bx2 cd * bx sd * bz by2 by bz2 -sd * bx cd * bz return bx2, by2, bz2geodetic_to_geocentric里 a 取 6378.137 km、f 取 1/298.257223563与原文一致返回的 θ 是地心余纬度不是纬度后面所有勒让德多项式都吃这个角。rotate_to_ellipsoid_normal只处理三分量张量要走 R·T·Rᵀ 的两侧乘法写成T2 R T R.T更省事。注意旋转角 d 在高纬地区会明显变大如果只做单侧旋转只乘 R 不乘 R张量会出现不对称用np.allclose(T2, T2.T)一查就露馅。3. 施密特半规格化勒让德多项式的一阶二阶递推实现3.1 为什么不直接用 scipy 的伴随勒让德函数scipy.special.lpmn算的是未经规格化的伴随勒让德多项式而且 m 从 0 开始、含 Condon-Shortley 相位约定与 IGRF 用的施密特半规格化差着一整套归一化因子。直接拿来乘高斯系数低阶项还能凑合看到 n8 以上数值会迅速失衡因为未规格化的 P_n^m 在高阶时量级跨度极大求和的相对重要性完全被打乱。更麻烦的是全张量梯度需要 P 对 θ 的一阶和二阶导数scipy 不直接给用差分近似会把 Laplace 自检毁掉。自己写递推还有个工程理由递推出来的 P 表在同一纬度圈上是复用的。0.1° 网格、4° 工区只有 41 个不同的纬度值经度方向 41 个点可以共用同一张 P 表和三张导数表这比逐点重算快一个数量级。3.2 对角项与非对角项的递推关系施密特半规格化下对角项递推是 P_n^n √(1 − 1/(2n)) · sinθ · P_(n−1)^(n−1)非对角项用三项递推P_n^m a_nm · cosθ · P_(n-1)^m - b_nm · P_(n-2)^m a_nm sqrt( ((2n-1)^2 - m^2) / (n^2 - m^2) ) b_nm sqrt( ((n-1)^2 - m^2) / (n^2 - m^2) )一阶导数只要对 θ 用乘积法则展开二阶导数再展开一次。原文式 (37)、(38) 给的就是这两个展开结果。实现时要注意 n1、m0 这种起点位置b_nm 的分子里会出现 (n−1)²−m² 为负的情况必须显式置零否则开根号直接出 NaN。阶数 k待求项数量含 m≤n二阶导项数量典型用途62828快速估算场形态106666一般区域磁场建模13105105与 NOAA 在线结果对齐16 以上153153岩石圈异常细分需谨慎3.3 Python 实现与数值稳定性处理下面这段是核心三个数组一次填满供后面的分量求和直接索引import numpy as np def schmidt_legendre(theta, K): 施密特半规格化伴随勒让德多项式及对 theta 的一阶、二阶导数 theta: 地心余纬度(弧度) K: 最大阶数 返回 P, dP, dP2形状均为 (K1, K1)索引 [n][m] st, ct np.sin(theta), np.cos(theta) P np.zeros((K 1, K 1)) dP np.zeros((K 1, K 1)) dP2 np.zeros((K 1, K 1)) P[0, 0] 1.0 # 递推起点 for n in range(1, K 1): if n 1: P[1, 1] st dP[1, 1] ct dP2[1, 1] -st else: c np.sqrt(1.0 - 1.0 / (2.0 * n)) P[n, n] c * st * P[n - 1, n - 1] dP[n, n] c * (st * dP[n - 1, n - 1] ct * P[n - 1, n - 1]) dP2[n, n] c * (-st * P[n - 1, n - 1] 2.0 * ct * dP[n - 1, n - 1] st * dP2[n - 1, n - 1]) for m in range(0, n): a np.sqrt(((2 * n - 1) ** 2 - m ** 2) / (n ** 2 - m ** 2)) # m n-1 时 P[n-2][m] 不存在系数置 0 b 0.0 if m n - 2 else np.sqrt(((n - 1) ** 2 - m ** 2) / (n ** 2 - m ** 2)) p1, p0 P[n - 1, m], P[n - 2, m] if m n - 2 else 0.0 d1, d0 dP[n - 1, m], dP[n - 2, m] if m n - 2 else 0.0 e1, e0 dP2[n - 1, m], dP2[n - 2, m] if m n - 2 else 0.0 P[n, m] a * ct * p1 - b * p0 dP[n, m] a * (ct * d1 - st * p1) - b * d0 dP2[n, m] a * (-2.0 * st * d1 ct * e1 - ct * p1) - b * e0 return P, dP, dP2theta必须是弧度从上一节的geodetic_to_geocentric拿。K是截断阶数与高斯系数的最大阶一致。二阶导数组 dP2 只在 Bxx 和 Bzz 的求和里用到如果只做三分量可以传一个开关跳过这段计算省时间。m n−2时把 b 置零是必须的否则P[n-2, m]会取到零行甚至负索引结果看着像对的其实是脏数据。3.4 高斯系数的读取与时间归算IGRF 系数文件按行给出 n、m、g、h 以及各年代的值和末段的年变率。给定日期 t 时落在两个publihed epoch 之间就线性插值落在最后一个 epoch 之后就用年变率外推常见写法是 g(t) g(t0) (t − t0)·。解析这类文本文件时注意三点h_n^0 恒为零、m0 的项不参与 φ 方向的求和、单位是 nT 而不是 T。系数表拿到手第一件事是抽查 n1 那一行的量级g_1^0 在 −29000 附近才是正常值差三个数量级说明单位或者符号读反了。4. 4°×4° 工区网格计算与 NOAA 逐点校验4.1 工区、网格与高度参数原文算例选的是经度 103.3056°107.3056°、纬度 27.3056°31.3056° 的 4°×4° 区域网格步长 0.1°大地高 1 km时间取 2019 年 4 月 7 日。这几个参数不是随便定的网格步长要小到能分辨张量等值线的形态又不能小到让等值线图被计算噪声糊住高度取 1 km 是因为航空测量本身有一个巡航高度同时主磁场梯度随高度衰减拿 1 km 的结果做学习飞行的高度设计比较贴近实际。4.2 批量计算的组织方式按纬度分组是提高效率的关键。同一个纬度上geodetic_to_geocentric的 r、θ、d 只算一次schmidt_legendre的三张表也只算一次经度方向的循环只更新 cos(mφ)、sin(mφ) 和 (a/r)^(n1) 这类标量。def compute_tensor_at_point(lat, lon, h, t, coeffs, K): 返回单点的七要素与六个独立张量分量(北东下) r, theta, d geodetic_to_geocentric(lat, h) phi np.radians(lon) P, dP, dP2 schmidt_legendre(theta, K) Ur Uth Uph 0.0 Uthth Uphth Urth Uphph Urph Urr 0.0 for n in range(1, K 1): ratio (6371.2 / r) ** (n 1) for m in range(0, n 1): g, hc coeffs[n][m] cm, sm np.cos(m * phi), np.sin(m * phi) gc_hs g * cm hc * sm # g·cos h·sin mgh m * (-g * sm hc * cm) # m·(-g·sin h·cos) Ur (n 1) * ratio / r * gc_hs * P[n, m] Uth - 6371.2 * ratio * gc_hs * dP[n, m] Uph - 6371.2 * ratio * mgh * P[n, m] Uthth - 6371.2 * ratio * gc_hs * dP2[n, m] Uphth - 6371.2 * ratio * mgh * dP[n, m] Urth (n 1) * ratio / r * gc_hs * dP[n, m] Uphph 6371.2 * ratio * m * mgh * P[n, m] / m if m else 0.0 Urph (n 1) * ratio / r * mgh * P[n, m] Urr - (n 2) * (n 1) * ratio / (r * r) * gc_hs * P[n, m] st, ct np.sin(theta), np.cos(theta) bx, by, bz -Uth / r, Uph / (r * st), Ur bxx Ur / r Uthth / r ** 2 bxy ct / (r ** 2 * st ** 2) * Uph - Uphth / (r ** 2 * st) bxz Uth / r ** 2 - Urth / r byy Ur / r Uth / (r ** 2 * np.tan(theta)) Uphph / (r ** 2 * st ** 2) byz Urph / (r * st) - Uph / (r ** 2 * st) bzz Urr return (bx, by, bz), (bxx, bxy, bxz, byy, byz, bzz)注意Uphph那一行在 m0 时不成立代码里用if m else 0.0挡掉实际写的时候把 m0 单独拎出去更清楚。奇点保护也要加θ 小于 1e-6 或者 sinθ 接近 0 时直接跳过该点或者把该点标记后单独用邻域值插补硬算出来的是无意义的巨大值。4.3 与 NOAA 七要素对照选取工区内四个点做对照本文算法与 NOAA 在线计算的结果差异如下Bx、By、Bz、Bh、B 单位 nTD、I 单位度点位来源BxByBzBhBDI成都 (30.67N,104.07E)本文计算33866.2−1213.937855.433887.950807.6−2.052848.1652成都NOAA33866.2−1213.937855.433887.950807.7−2.052948.1652自贡 (29.35N,104.78E)本文计算34644.5−1268.736051.734667.650015.7−2.097246.1212自贡NOAA34644.5−1268.736051.734667.750015.8−2.097346.1212泸州 (28.87N,105.43E)本文计算34909.9−1334.735359.934935.349707.1−2.189545.3460泸州NOAA34909.9−1334.735359.934935.449707.2−2.189545.3460德阳 (31.13N,104.38E)本文计算33580.9−1266.438443.833604.751060.8−2.159748.8424德阳NOAA33580.9−1266.438443.833604.851060.8−2.159748.8424三分量在保留一位小数的情况下完全一致水平强度、总场和偏角最大绝对误差分别是 0.1 nT、0.1 nT 和 0.0001°倾角无差异。这个量级的残差主要来自 NOAA 结果本身的有效位数以及大地坐标到地心坐标换算时的舍入不是推导错误。做校验时要注意统一时间拿 2019 年的算例去对 2020 年的在线结果残差会跳到几十 nT。4.4 Laplace 方程自检与张量对称性检查没有第二来源的全张量数据可比所以正确性只能靠物理约束兜底。主对角线三项相加应当恒为零原文算例里这个和落在 −0.00110.0011 nT/km相对张量分量本身的量级nT/km 到几十 nT/km已经算数值零。实际编码时按下面三件事查一遍张量矩阵对称性np.allclose(T, T.T, atol1e-6)不对称说明坐标旋转只做了一半。迹的分布图画成等值线如果出现明显的条纹或边界跳变多半是 θ 接近极点时的除零污染。阶数收敛k 取 8、10、13 各跑一遍张量分量的变化应该随 k 单调减小并趋于平稳如果 k 增大后反而发散递推的规格化因子写错了。4.5 容易翻车的几个点纬度单位弄错是最高频的错误。所有勒让德多项式的自变量是地心余纬度 θ等于 90° 减地心纬度而输入的大地纬度还要先转成地心纬度中间差一个旋转角 d。漏掉这一步在中纬度地区会带来零点几度的方向偏差映射到张量上就是百分之几的量级错误图看着还挺像回事。第二高频的是规格化方式混用。高斯规格化、施密特半规格化、完全规格化三套因子互不相同混着用会让高阶项和低阶项的相对权重彻底跑偏。判断方法很直接算 n1、m0 这一项施密特半规格化下 P_1^0 cosθ如果代码输出的是别的值说明归一化因子挂错了。第三是角度与弧度的边界。递推表内部的 θ、φ 都是弧度但网格生成、绘图和输入参数多半用度转换点最好只留在入口和出口各一处中间层全部用弧度避免在循环里反复np.radians。5. 大批量网格的向量化、阶数收敛与云端批处理5.1 递推表复用与数组化单点跑通之后真正的开销在网格上。4°×4°、0.1° 网格是 41×41 共 1681 个点单点循环在笔记本上还能忍换成全国范围 0.05° 网格就是几十万个点。两个优化收益最大按纬度分组的递推表缓存把 P、dP、dP2 的计算次数从点数降到纬度数以及把经度方向的求和数组化用np.einsum或者预先把 cos(mφ) 矩阵乘上去把内层 m 循环消掉。前者是纯算法层面的复用后者靠 NumPy 的向量化两者叠加通常能压掉一个数量级以上的时间。5.2 阶数收敛验证IGRF 的球谐系数阶数有限梯度对高阶项的敏感度比场值本身高因为每一项都带着 (n1)(n2)/r 这类放大因子。验证方式是把 k 从 6 递增到模型最大阶逐档记录工区内 Bxx、Byy、Bzz 的最大变化量截断阶数 kBxx 最大变化 (nT/km)说明6基准偶极子主导形态大体正确10较 k6 有可见下降四极、八极项贡献开始收敛13继续下降但幅度变小接近模型上限超过模型阶数不变系数为零加了也没用判断标准落在仪器噪声上如果 k 从 13 加到 16 带来的张量变化已经低于航空测量系统的噪声底继续加阶数只是白烧算力。反过来如果 k6 和 k13 的梯度图形态差异明显说明测试工区落在高阶异常较强的区域用低阶结果做飞行设计会低估梯度。5.3 参数固化成技术文档与云上批量跑整个计算链没有 IO 依赖输入是纬度、经度、高度、时间和系数表输出是六个数天然适合打包成容器放进云上按网格分块跑。我的做法是把椭球参数、参考球半径、最大阶数、时间归算方式写在配置里随代码一起提交同时把每个参数的含义、单位和典型取值范围整理成一页技术文档谁接手都能改对地方。需要注意的是容器里浮点行为要一致np.float64和某些加速库的默认精度混用会让 Laplace 自检的残留量从 1e-3 跳到 1e-1一批结果直接作废。跑批量任务前先在本地用四个对照点做一次回归通过以后再上云能省掉大量返工。本文还有配套的精品资源点击获取