2. 河流网络拓扑构建与Strahler分级

任务边界与数据前提

本案例的输入是例1已经提取并经过基本清理的河网矢量（线图层），属性中至少包含唯一河段编号 reach_id 与汇流累积量 grid_code（或 accum）。输出需要三样东西：一张带拓扑关系的河段表、每个河段的 Strahler 级数、以及按级别汇总的数量与长度统计。核心难点不在分级算法本身，而在于先把“谁接谁”的拓扑关系建对——拓扑错了，分级必然错。

拓扑构建的关键判断是：一条河段的下游端点若与另一条河段的上游端点重合（容差内），则前者是后者的上游支流。GIS 中常用“端点捕捉 + 空间连接”实现，但更稳妥的做法是先把线拆成单段（在节点处打断），再按端点坐标匹配。下面给出可直接在 QGIS 字段计算器 / PyQGIS 或 ArcGIS Pro 的 Python 窗口中运行的实现思路，核心逻辑与平台无关。

拓扑关系与 Strahler 分级的可视化

下图展示了一个典型树状河网的拓扑结构与分级结果。注意 Strahler 规则的两个要点：同级支流汇合升一级，不同级汇合取高级；河源（无上游的河段）恒为 1 级。图中用颜色区分级别，箭头表示水流方向（由上游指向下游）。

拓扑构建与分级实现

以下代码以 PyQGIS 为例（ArcGIS Pro 可用 arcpy 对应改写），假设河网图层名为 river，已按节点打断。核心分三步：端点匹配建拓扑、按拓扑排序自上游向下游递推分级、汇总统计。

from qgis.core import QgsProject, QgsFeature, QgsGeometry
from collections import defaultdict, deque

layer = QgsProject.instance().mapLayersByName("river")[0]
feats = {f["reach_id"]: f for f in layer.getFeatures()}
TOL = 1e-6  # 坐标容差，按数据坐标系单位调整

# 1) 端点坐标 -> 河段映射（起点=上游端，终点=下游端）
def pt_key(pt):
    return (round(pt.x()/TOL), round(pt.y()/TOL))

start_at, end_at = defaultdict(list), defaultdict(list)
for rid, f in feats.items():
    g = f.geometry().asPolyline()
    start_at[pt_key(g[0])].append(rid)
    end_at[pt_key(g[-1])].append(rid)

# 2) 建拓扑：downstream[rid] = 下游河段；upstream[rid] = 上游河段列表
downstream, upstream = {}, defaultdict(list)
for rid, f in feats.items():
    end_pt = pt_key(f.geometry().asPolyline()[-1])
    for cand in start_at.get(end_pt, []):
        if cand != rid:
            downstream[rid] = cand
            upstream[cand].append(rid)
            break

# 3) 拓扑排序（Kahn），保证先算上游再算下游
indeg = {rid: len(upstream[rid]) for rid in feats}
q = deque([r for r in feats if indeg[r] == 0])
order = []
while q:
    r = q.popleft(); order.append(r)
    if r in downstream:
        d = downstream[r]; indeg[d] -= 1
        if indeg[d] == 0: q.append(d)

# 4) Strahler 分级
strahler = {}
for rid in order:
    ups = upstream[rid]
    if not ups:
        strahler[rid] = 1
    else:
        levels = [strahler[u] for u in ups]
        mx = max(levels)
        strahler[rid] = mx + 1 if levels.count(mx) >= 2 else mx

# 5) 写回属性并统计
layer.startEditing()
layer.dataProvider().addAttributes([QgsField("strahler", QVariant.Int)])
layer.updateFields()
idx = layer.fields().indexOf("strahler")
for f in layer.getFeatures():
    f["strahler"] = strahler[f["reach_id"]]
    layer.updateFeature(f)
layer.commitChanges()

stat = defaultdict(lambda: {"n": 0, "len": 0.0})
for rid, f in feats.items():
    s = stat[strahler[rid]]
    s["n"] += 1
    s["len"] += f.geometry().length()
for lv in sorted(stat):
    print(f"{lv}级: 数量={stat[lv]['n']}, 总长={stat[lv]['len']:.1f}")

结果判读与合理性检查

分级完成后，用三条经验规则快速验证拓扑与分级是否正确：

• 级别单调性：沿任一从河源到出口的路径，级别不应下降。若出现下降，说明拓扑连接方向反了（把下游当成了上游）。
• 数量金字塔：1 级河段数量应最多，级别越高数量越少；若 1 级数量异常少，可能是河网未在节点处打断，导致多条支流被合并成一条。
• 长度特征：高级别河段通常单段更长、总长占比更大；若 4 级干流总长反而小于 3 级，需检查干流是否被错误拆分。

下表是一个示例流域的统计输出样式（数值为教学示例，非实测）：