第6章 冰川中心线与几何参数化

centerlines.py 约2454行，是OGGM中算法密度最高的模块。它的任务是将一个二维冰川掩码（第5章的产物）转化为一组合适的流线（flowlines）——每条流线都携带几何属性（宽度、面积-海拔分布），并组织为有向无环图（DAG）来模拟支流汇入主干的拓扑关系。本章沿着数据转换的完整链路，深入Kienholz et al. (2014) 最短路径算法的每个细节。

6.1 Centerline类：几何与拓扑的容器

Centerline 类（约第69行）封装了一条流线的所有信息。流线在代码中表示为沿中心线的等距坐标点序列：

class Centerline():
    def __init__(self, line, dx=None, surface_h=None):
        self.line = line       # shapely LineString: N个等距点
        self.dis_on_line = ... # 距终点的沿程距离 [m]
        self.nx = len(line.coords)
        self.order = ...       # Strahler-like序列号
        self.flows_to = ...    # 下游汇入的另一条Centerline对象
        self.inflows = []      # 汇入本线的上游Centerline对象列表

关键设计是 flows_to 和 inflows 构建的有向无环图（DAG）。通过inflows可以遍历所有上游支流，通过line_order()可以建立Strahler排水次序式的拓扑排序——从末梢支流开始在依次向上游推进到主干。

入门

6.2 Voronoi成本网格（_make_costgrid）

_make_costgrid（约第684-714行）是中心线提取的核心，将冰川掩码转换为"通过成本最小的路径从末端走到头部"的成本图。成本函数源自Kienholz et al. (2014) 的方法：

cost = ((dmax - dis) / dmax * f1)^a  +  ((z - zmin) / (zmax - zmin) * f2)^b

成本由两项加权和构成：

• 距离项：(dmax-dis)/dmax——距冰川轮廓边界越远的像素成本越低（中央像素受到偏好）。f1 默认0.3，控制距离项的贡献。
• 高程项：(z-zmin)/(zmax-zmin)——高程越低的像素成本越低（模拟冰的河流式流动特征）。f2 默认0.7，使得高程项权重更高，确保流线向下行走。

指数a和b（默认均为1）控制各项的敏感度。此外，对冰川掩码外部和边界像素施加50倍惩罚因子，有效阻止路径"逃出"冰川区域。成本网格的可视化通常呈现为从末端到头部的平滑成本梯度。

6.3 冰川头部检测（_get_centerlines_heads）

冰川"头部"是中心线的起点，位于冰川支流的最上游端。_get_centerlines_heads（约第775-817行）使用局部高程极大值检测和自适应过滤来识别它们。

6.3.1 局部极大值检测

对地形数据应用 scipy.signal.argrelmax，扫描高程场的每个像素，当中心像素值高于其邻域所有像素时标记为局部极大值。参数order=5 控制比较范围的窗口大小。

6.3.2 自适应过滤

Kienholz et al. (2014) 中定义了一个与冰川面积相关的缓冲区半径：

radius = q1 * area + q2

其中q1和q2是经验系数（默认q1=0.15, q2=0）。_filter_heads 函数在此半径内仅保留高程最高的那个局部极大值作为头部候选点。候选点还需要高于高程分布的第33百分位数，过滤掉地形噪声中的微小峰值。

参数rmax 限制了能检测到的最大头部数量。对于喜马拉雅等山地冰川，通常每个流域产生2-5个头部。

6.4 最优路径路由（compute_centerlines）

compute_centerlines（约第874-977行）是中心线提取的主控函数。对每个检测到的头部，从该像素到冰川末端计算一条成本最低的路径。

6.4.1 最短路径算法

路径搜索通过 skimage.graph.route_through_array 实现，它在成本网格上直接执行Dijkstra-like搜索。指定起点（头部）和终点（末端）后，返回一条成本积分最小的路径。

6.4.2 末端识别

末端的确定取决于冰川类型。对于入海型冰川（tidewater），末端为冰川高程分布的某百分位数处（默认第5百分位，确保位于冰舌断裂线附近）。对于陆地终止冰川，末端为冰川范围内的最低高程像素。

6.4.3 线的过滤与连接

中心线后处理包含四个级联步骤：

• _filter_lines：保留最长的一条线，然后逐步缓冲剩余线——短于主线的重叠线被移除，避免多条线堆积在同一主干上。
• _filter_lines_slope：移除沿高程上升行走的线（物理上不合理的"上坡中心线"），通过检查沿线高程的趋势符号实现。
• _join_lines：将支流通过几何投影连接到主干的对应交点位置。支流的终点不是主干的终点，而是在主干的side pixel处终止。
• _line_order：对整个DAG进行拓扑排序，给每条线赋予order值，类似于河流的Strahler排序，主干最高、小支流最低。

def _join_lines(cls, line, main_line, gdir):
    # 1. 获取line的终点像素坐标
    # 2. 在main_line上找到距离最近的像素（几何投影）
    # 3. 沿main_line回溯到line的终点高程之上2dx，作为交汇点
    # 4. 截断line使其在交汇点结束
    # 5. 设置 line.flows_to = main_line

6.5 宽度计算（_point_width）

_point_width（约第1375-1444行）计算中心线每个点的冰川宽度。算法如下：

1. 在离散点处沿法线方向取长度为 kbuffer * map_dx 的线段。
2. 该线段与冰川集水区多边形求交集。
3. 交集的长度即为该点处的冰川宽度。
4. 如果该点在nunatak（岩壁）内部，从宽度中减去nunatak覆盖的部分。

参数kbuffer（默认3.0）确保法线段足够长以跨越整个冰川宽度区域。Nunatak的减法处理是OGGM区别于许多其他模型的重要特性，保证了宽度参数的真实性。

6.6 集水区面积与宽度校正

6.6.1 集水区面积

catchment_area（约第1549-1662行）通过 salem.Grid 的分水岭算法将冰川划分为多边形的子流域——每个多边形代表一条中心线的汇水区域。集水区边界遵循冰面地形的水文线路，类似于流域边界的概念。

6.6.2 宽度校正（catchment_width_correction）

catchment_width_correction（约第1935-2081行）是一个关键的质量保证步骤。它确保各中心线的面积-高度分布与RGI数据库中该冰川的已知面积相匹配：

1. 将每条中心线的宽度-高程积分为面积-高度分布（分箱为50m高程间隔）。
2. 将各中心线的曲线累加，得到总分布。
3. 对每个高程箱进行缩放，使总分布与RGI面积一致（bin-based校正）。
4. 最终对每条线的所有点的宽度做全球缩放，使总面积精确匹配RGI面积。

这个校正步骤确保了OGGM中的几何参数化与RGI的全球冰川清单保持一致，使得区域到全球尺度的体积估算具有可比性。未校正的宽度积分面积与该冰川RGI面积之间的偏差通常小于5%，但仍需精确校正。

6.7 海拔带流线（elevation_band_flowline）

elevation_band_flowline（约第2181-2352行）是OGGM 1.4+ 中引入的"折叠冰川"概念的实现（Huss 2012, Werder 2019）。它将每条中心线压缩为一维表示：

# 沿中心线将每个海拔段表示为一个"盒状"单元
dx = bsize / tan(slope)   # 盒子的长度：高程步长 / 坡度
width = area / dx         # 盒子的宽度：面积 / 长度

其中bsize是固定的海拔步长（通常200m），slope是该海拔段顶部高程处的当地表面坡度，area是该海拔段内中心线集水区的总面积。这个方案将二维冰川问题映射为一维，极大地加速了冰流模型的数值求解（第10章），同时保留了冰川的整体几何特征。

进阶

6.8 下游线（Downstream Line）

每条中心线都附带一条下游线（downstream_line）——从冰川末端沿地形向更低海拔延伸的线。它用于：

• 物质平衡模型中确定冰川末端以下的高程外推点
• 冰流模型中定义冰川前缘之外区域的等温条件
• 为入海型冰川计算海底线（submarine bed）以估算崩解通量

下游线由 _downstream_line 函数通过沿地形梯度步进生成，步进终点为指定的最低高程或最大长度。

6.9 关键参数速查

参数名	所属函数	默认值	描述
flowline_dx	Centerline构造	map_dx * 2.0	中心线离散点的间距
q1	_get_centerlines_heads	0.15	头部缓冲半径的面积系数
q2	_get_centerlines_heads	0.0	头部缓冲半径的偏移量
rmax	_get_centerlines_heads	10	最大头部数量限制
f1	_make_costgrid	0.3	距离项在成本中的权重
f2	_make_costgrid	0.7	高程项在成本中的权重
a	_make_costgrid	1.0	距离项指数（灵敏度）
b	_make_costgrid	1.0	高程项指数（灵敏度）
kbuffer	_point_width	3.0	法线缓冲宽度倍数
tidewater_percentile	compute_centerlines	5	入海型冰川末端高程百分位数

6.10 总结

centerlines.py 体现了OGGM工程设计的核心思想：将复杂的冰川几何问题转换为图论中的最短路径问题。Voronoi成本网格、自适应头部检测、流线合并和宽度校正这四个模块相互衔接，每一步都有清晰的物理直觉作支撑。中心线的质量直接决定了后续物质平衡、冰厚反演和动力学演化的精度——一条偏心的中心线会导致错误的宽度、错误的面积-高度分布，并最终影响冰量变化的所有预测。

在下一章中，将看到气候输入数据如何通过climate.py和shop/模块变换为标准化的climate_historical.nc文件，为物质平衡模型提供气象驱动。