第9章 冰厚反演

inversion.py 约1180行，实现了OGGM的质量守恒冰厚反演方法（Farinotti et al. 2009; Huss and Farinotti 2012）。反演把表观物质平衡（apparent mass balance）积分得到的冰通量与SIA通量关系联系起来，从而求解沿流线的冰厚。没有冰厚信息，后续就无法可靠计算冰流速度、崩解通量或冰川对气候变化的动力响应。

9.1 物理基础：质量守恒与SIA

冰厚反演的物理核心来源于两个约束：(1) 质量守恒——通过给定横截面的冰通量（ice flux）必须等于上游集水区总物质平衡的积分；(2) 浅冰近似（Shallow Ice Approximation, SIA）——给定表面坡度和冰厚，可以计算竖向积分的冰通量。

关键公式（Huss and Farinotti 2012）：通过横截面的冰通量 q 为：

q = int_0^H u(z) dz = f_d * (2*A/(n+1)) * (rho*g*alpha)^n * H^(n+2)
                         + f_s * (2*A) * (rho*g*alpha)^n * H^n

其中：

• A —— Glen流动律参数（默认2.4e-24 Pa^-3 s^-1）
• n —— Glen指数，= 3
• rho —— 冰密度（900 kg m^-3）
• g —— 重力加速度（9.81 m s^-2）
• alpha —— 表面坡度（来自DEM）
• H —— 冰厚（未知量，反演目标）
• f_d —— 形变流动的形状因子（shape factor；取决于冰床横截面形状 bed cross-section geometry：抛物线/矩形/梯形）
• f_s —— 滑动的形状因子（shape factor；通常 f_s=0，即假设基底冻结）

通量q来自上游表观物质平衡（apparent mass balance）的积分；在离散流线中，它本质上是上游净质量收支换算得到的冰体积流率。因此反演任务变为：已知 q 和 alpha，求解 H。

进阶

9.2 准备阶段（prepare_for_inversion）

prepare_for_inversion（约第53行）为每条流线的每个横截面位置预先计算冰厚求解所需的所有辅助量：

def prepare_for_inversion(gdir, inversion_input=None):
    # 1. 计算表面坡度：负梯度（高程向下降方向）
    slope = -np.gradient(hgt, dx)

    # 2. 从表观物质平衡（apparent mass balance）计算质量通量
    flux = fl.flux * dx**2 / SEC_IN_YEAR / rho   # m3 s-1

    # 3. 预计算通量/宽度比
    flux_a0 = flux / width * shape_factor           # 用于SIA多项式

表观物质平衡（apparent mass balance, apparent MB）考虑了冰川几何相对于参考几何的变化——不同于第8章的固定几何MB，它需要调用 apparent_mb_from_any_mb 来产生实际的物质收支。

坡度计算中的关键细节：对坡度进行截断，alpha = max(alpha, min_slope)，其中 min_slope = tan(1.5°)。当实际坡度低于此阈值时（如在冰斗或宽阔河谷中），强制取这个最小值——因为零坡度会导致alpha^n = 0，使得冰厚方程退化。

9.2.1 形状因子（shape factor）

形状因子（shape factor）f_d取决于假定的冰床横截面形状（bed cross-section geometry）：

• 半圆/抛物线（f_d = 0.5）：经典值，适用于山谷冰川
• 矩形（f_d = 1.0）：无侧向约束（流速仅受厚度和坡度控制）
• 梯形：介于两者之间，取决于形状参数trapezoid_lambdas

为了兼容默认抛物线方案和物理上更合理的梯形方案，OGGM使用 mixed_min_shape 阈值进行平滑过渡。

9.3 厚度求解器

对于质量通量 q、坡度 alpha、形状因子（shape factor） fs 和 f_d，SIA方程化为一个关于 H 的五次多项式：

H^5 + a3 * H^3 + a0 = 0

9.3.1 多项式求解（_inversion_poly）

当fs = 0（无基底滑动）时，多项式退化为 H^5 + a0 = 0，解为 H = (-a0)^(1/5)（对应_inversion_simple，约第153行）。

当fs != 0时（有基底滑动项），_inversion_poly（约第146行）使用np.roots()求解五次多项式的所有根，并选择物理上可接受的正实根作为厚度。多项式系数a3和a0由SIA方程系数组装。

9.3.2 遍历求解（_compute_thick）

_compute_thick（约第159行）对方程网格中所有flux > 0的横截面逐点求解。循环结构为：

for i in range(nx):
    if flux[i] > 0:
        thick[i] = _inversion_poly(a3[i], a0[i])
    else:
        thick[i] = 0     # 无物质积累的位置 = 零冰厚

9.3.3 Brent优化求解（sia_thickness_via_optim）

sia_thickness_via_optim（约第194行）是对_inversion_poly的替代方案，使用brentq数值求解非线性方程：

u = H^(n+1) * fd * (rho*g*alpha)^n + H^(n-1) * fs * (rho*g*alpha)^n
sect * u - flux = 0

其中sect表示由冰厚和横截面形状确定的截面积。Brent方法比多项式求根更稳健，尤其适合梯形横截面等无法直接化成简单解析形式的情况，但每点需要多次函数求值。

9.3.4 sia_thickness 解析函数

sia_thickness（约第264行）是对外的公共接口，接受表面斜率、通量、宽度等参数，在内部调用_compute_thick进行解算。这是外部模块（如动力学、崩解等）使用SIA厚度估计的标准入口。

9.4 核心任务：mass_conservation_inversion

mass_conservation_inversion（约第359行）是完整的反演任务函数。其算法流程为：

1. 读取输入：从gdir加载流线、质量通量、DEM高程、宽度、坡度等信息。
2. 裁剪坡度：alpha = max(slope, min_slope)，防止数值问题。
3. 计算多项式系数：对每个横截面组装 a0（来自flux_a0）和 a3（来自fs项）。
4. 求解厚度：调用 _compute_thick 获取每个横截面的冰厚。
5. 计算体积：V = sum(thick * width * dx) 得到冰川总体积。
6. 梯形后备方案：对于bed_shape低于mixed_min_shape的横截面，切换到梯形SIA解算器（sia_thickness_via_optim），以更好地表示宽谷或平坦冰床横截面（bed cross-section）中的冰流。

反演的三个关键输入参数是：

• glen_a（默认2.4e-24）：Glen流动律参数，控制冰的柔软度。越大的值代表越软的冰（更快的流动），从而获得更薄的厚度。
• fs（默认0）：基底滑动因子。非零值时允许冰床滑动，降低形变厚度。
• trapezoid_lambdas：梯形形状参数。在梯形横截面方案中控制谷底的平坦程度。

进阶

9.5 后处理

9.5.1 过滤与平滑（filter_inversion_output）

filter_inversion_output（约第505行）对反演的原始厚度场进行两种后处理：

• 末端凹槽平滑：冰川末端的反演厚度常出现非物理的"末端凹槽"（terminal depression）——因为通量在末端前急剧下降为0，导致局部厚度过薄。函数沿前缘方向应用低通滤波来平滑这些凹槽。
• 转换为梯形表示：将抛物线-基底表示转换为更物理的梯形形状，以便后处理算法（如2D厚度场重建）使用。

9.5.2 反演速度计算（compute_inversion_velocities）

compute_inversion_velocities（约第632行）从反演的厚度场计算表面速度：

u_deformation = (2*A/(n+1)) * tau^n * H     # 内部形变贡献
u_basal = A * tau^n / lambda * ...               # 基底滑动贡献（可选）
u_surface = u_deformation + u_basal

当fs=0时，深度平均速度与表面速度之间满足 u_mean/u_surface = (n+1)/(n+2)；因此在n=3时，表面速度约为深度平均速度的1.25倍。速度和应力可用于验证反演的一致性，也可作为后续诊断量。

9.5.3 体积输出（get_inversion_volume）

get_inversion_volume（约第625行）计算反演的总冰体积：V = sum(thick * width * dx)，以m³为单位。这是在区域到全球尺度上汇总冰储存量的核心指标。

9.6 厚度场扩展（2D可视化）

9.6.1 按海拔分配厚度（distribute_thickness_per_altitude）

distribute_thickness_per_altitude（约第728行）使用反距离加权（IDW）方法将一维沿中心线的厚度扩展到二维网格上。核心思想是：

1. 对每个海拔带，收集该带内所有中心线点的厚度值。
2. 在该海拔带内的所有冰川掩码像素上对厚度进行IDW插值——离中心线越近的像素获得更接近中心线点的厚度权重。
3. 不同海拔带之间独立处理——避免了横谷的厚度平滑伪影。

这种"按海拔分带"的IDW方法比全局IDW插值更能保持冰川厚度场中海坡梯度的物理合理性。

9.6.2 立方插值（distribute_thickness_interp）

distribute_thickness_interp（约第856行）是替代的二维插值方案，使用scipy.interpolate.griddata的cubic模式进行三阶插值。对于光滑的大型冰川（如阿拉斯加沿海冰川），这个方案能产生视觉上更平滑的厚度图，但可能模糊冰川边缘的厚度梯度。

9.7 崩解通量的厚度依赖

9.7.1 calving_flux_from_depth

calving_flux_from_depth（约第960行）使用Oerlemans and Nick (2005)的崩解速率公式：

q_calving = k * water_depth * ice_thickness * width

其中k是崩解常数；OGGM v1.6.3中反演使用的默认inversion_calving_k为0.6 yr^-1。water_depth是末端处的水深，可由冰厚、出水高度（freeboard）和水位约束估算。这个公式表达的是经验关系：更深水位通常对应更高的崩解通量。

9.7.2 迭代一致性求解（find_inversion_calving_from_any_mb）

find_inversion_calving_from_any_mb（约第1033行）解决了一个"鸡和蛋"问题：冰厚依赖于通量，而通量依赖于崩解损失，崩解又依赖于冰厚。通过迭代求解实现SIA反演与崩解通量之间的自洽：

1. 用当前冰厚估计崩解通量。
2. 将崩解通量加入物质平衡通量。
3. 用修正的flux求解新冰厚。
4. 重复直至收敛。

一般在5-10次迭代内收敛。这是处理入海型冰川的关键技术，因为崩解是这类冰川质量损失的主要形式。

9.8 关键参数速查

参数名	默认值	作用	单元
glen_a	2.4e-24	Glen流动律参数	Pa-3 s-1
fs	0.0	基底滑动因子	无量纲
min_slope	tan(1.5°)	最小坡度阈值	m/m
mixed_min_shape	0.001	梯形切换阈值	无量纲
trapezoid_lambdas	变化	梯形形状参数	无量纲
inversion_calving_k	0.6	反演阶段使用的崩解速率常数	yr-1
use_kcalving_for_inversion	False	是否在反演中启用k-calving校正	-

9.9 总结

inversion.py 将第5章到第8章的成果——地形、几何、气候、物质平衡——融合为一个统一的冰厚度反演框架。物理上，它基于质量守恒和浅冰近似这两个冰川学的基石原理；数值上，它将问题转化为五次多项式的求根或Brent数值优化；工程上，它用迭代方法和自适应形状因子（shape factor）处理从山谷冰川到入海型冰川的多样化冰川形态。

反演的冰厚场连接了预处理和正演模拟：在此之前是GIS处理、几何参数化和物质平衡建模，在此之后是冰流正演、冰川演化和海平面贡献估计。反演精度会影响后续动力学模拟：若反演偏薄，冰川响应通常偏快；若反演偏厚，响应通常偏慢。因此，冰厚反演的参数选择和区域校准（特别是inversion_glen_a）是OGGM应用中必须认真检查的环节。