用程序推演一个架空世界的地形、气候与河流

先上成品：World Map（右上角可以切 3D 地球，左侧可以切换地形 / 气候 / 洋流 / 水系等视图）

技术栈 Python（numpy / scipy / numba / OpenCV），部分数据源：NOAA ETOPO 2022 高程数据。

感谢Claude模型提供的大量帮助

前言

事情始于高二时和好友 @电子鱼一起画过的一张地图。

那时的我们并不满足于只是绘制地图的轮廓，还尝试真正推演其中的地理过程，让地形、气候和水系彼此关联，并尽可能遵循真实世界的地理规律。

但是比较遗憾的是，在当时受限于技术力，虽然有推理分析，但没有一个完整的数学建模的过程，大多都是定性的，定量计算少之又少（不过这是不是也是高中地理的“特色”呢）。

为了给当年的野心交一份答卷，便有了这个项目：把当年那些定性的推理一环一环换成可以计算的方程，去推演这个世界的地形、风向、洋流、降水、气候、河流、水系等。

做下来发现，整个过程一直在和同一个矛盾打交道：手绘图是我们想要的那个世界，可只要把手绘的线条直接拿来计算，结果就会很假。所以每一步要做的，都是让手绘退到约束的位置，具体的形状交给地理过程去生成。

地形：真实高程数据的尺度分解与移植

之前的手绘图已经有了大陆的轮廓和山脉的位置，缺的是高程。第一阶段的目标是一张覆盖陆地和海底的全球高程图，大格局要和手绘一致，放大之后的细节又要像真实的地形。

原始手绘图层处理

对 PS 手绘底图进行预处理与矢量化提取后，得到以下基础输入数据：

海陆掩膜：陆地轮廓的二值化 Mask；
山脉骨架：经过骨架化算法提取出的山脊折线；
地形区面域：高原、盆地、丘陵的多边形边界；
河流线：提取水系走向，并将近海端点标记为河口；
板块边界：标记张裂（离散）与消亡（汇聚）边界。并结合两侧海陆属性，细分为陆陆碰撞、洋陆俯冲、岛弧三类构造。

图 3 · 解析后的结果：红线是汇聚边界，紫线是离散边界，彩色块是高原、盆地、丘陵区域

第一次尝试: 距离场加分形噪声

从线条得到高程，最直接的办法是程序化地形的常规流程：以到海岸的距离作为基底高程，沿山脉线叠加高斯剖面的山脊，高原和盆地的范围整体抬升或压低，再叠六个倍频的分形噪声（fBm），最后做 80 步 stream-power 河流侵蚀。

结果如上图所示，大体尺度基本OK，但是放大以后还是有很多细节瑕疵：

首先形状上，距离场和高斯核把线条直接扩散成高程，完全继承了手绘形状。而PS里面的手绘线缺乏真实山脉的曲折，放大会暴露出严重的等宽山体与同心圆伪影。

其次细节上，fBm 生成的起伏缺少真实地形的结构。它的功率谱可以调到和真实地形一致（谱斜率约为 2），但各频率分量的相位是随机的。真实地形里，水系的树状组织、构造线的走向、山脊与谷地的不对称，都体现为相位之间的相关性。

虽然侵蚀模拟可以逐步建立这种相关性，但需要的迭代量远远超过几十步。说白了，fBm 的起伏在统计上像山，但不会自己长出水系和山脊。

所以下一步的目标就是：细节得换一个来源，手绘的线条也得换一种方式进入计算（至少不能是直接提取轮廓来参与计算了）。

第二次尝试: 尺度分解

对于细节部分，我们直接使用实测高程来进行偷懒（bushi），这里用的是 ETOPO 2022数据源。但真实地形不能整块搬来，我们只需要地形相对微观的细节，这可以通过按空间尺度拆分高程场做到：

z = \underbrace{G_\sigma * z}_{\text{低频}} + \underbrace{\left(z - G_\sigma * z\right)}_{\text{高频}}, \qquad G_\sigma\ \text{为高斯核}

低频决定区域基底高度，高频负责沟谷、山脊等纹理细节，两者的输入完全独立：

低频（宏观格局）：高原、盆地等区域赋以目标高程，交界处设 300–400 km 软过渡带并叠加密度扰动，消除台阶状的痕迹。同时，手绘线预先经过随机扰动（海岸附近不扭曲以保证海岸线不动），扭曲之后的线仍然在手绘的位置附近，但形状就会从我鼠标画出来的那种粗糙简陋的圆润曲线变成有曲折扭曲的更真实的山脉线。
高频（微观细节）：高频的细节部分来自 ETOPO的地球高程数据。对真实地形做高通滤波，保留 1200 km 以内所有尺度的起伏，只丢掉整体高度。滤波之前，源区域里的海域先用最近的陆地高程填补，来避免海岸线处出现很大的假边缘。

这样拆开之后两部分互不干扰（这个按尺度拆开的办法，后面处理河谷和气候降尺度时还会再用两次），那么剩下的问题就是高频部分从地球上的哪里取、放到图上的哪里去。

山脉和成片的平原情况不同，先介绍山脉。

山脉：曲线坐标系下的条带移植

山脉是线状构造，主脊、支脉、山前丘陵都沿走向组织。

如果只是随机取几块山地纹理拼在一起，得不到一条连贯的山脉。

所以山脉要整条移植：选一条真实的山脉，让它沿着目标线的走向弯过去。

弯过去之前先要拉直。在地球上沿真实山脉的走向建立曲线坐标 (s, t)，走向先做 60 km 平滑再求法向，按 2 km 步长采样，得到一个矩形数组，相当于一条拉直了的山脉。

拉直得到的条带同样按尺度拆成两层：120 km 低通的部分是山体，残差是山体上的沟谷。两层用不同的羽化宽度贴回去，山体 500 km，细节 300 km，让山前因此有一条逐渐抬升的丘陵带。

这样得到的山脉，主脊连续，两侧是真实的支脉和冰川谷，宽度和高度沿线变化，拐弯处没有接缝。它的位置和手绘线大致一致，但不完全重合。

陆地背景：重叠相加平铺

山脉的办法用不到平原上。山脉以外的陆地是成片的区域，没有固定的走向，面积又远大于任何一块可用的真实地区，只能用小块纹理去铺。这里的要求是块与块之间看不出接缝，整体看不出重复。

纹理的来源先按手绘的气候分区选取（毕竟程序化的气候推理要等到地形模拟出来之后，但好在我自己推理的气候分区还是大差不差的），让地貌的质感和气候相称：沙漠区取撒哈拉和阿拉伯，雨林区取亚马逊和刚果，温带取西欧，冻土带取泰梅尔和西伯利亚。

接缝的问题借用信号处理里的重叠相加（overlap-add）解决。窗函数取 Hann 窗 w(u) = sin²(πu)，相邻两块错开半个窗宽，此时任意位置上各窗的权重之和恒为 1，加权平均不会带来明暗起伏。

重复的问题靠随机化：每一块从源区域随机取窗口，随机旋转任意角度。另外东西方向按 cos φ_目标 / cos φ_源校正比例，否则高纬的地貌会被横向拉长。

不同气候类型的权重场先做平滑再归一化，类型之间是渐变过渡。

再把铺出来的平原和丘陵带着真实的水系切割纹理，这样子让干旱区和湿润区的质感也有区别。

板块边界和海洋

陆地铺完，还剩海底。海底的形态主要由板块构造决定。第一版的大陆架和大陆坡只是离岸距离的函数，海岸线外面因此套着一圈等宽的浅色带，看上去很假，claude提出的新思路是改为按板块边界的类型布置。

图 6 · 第一版的海洋：大陆架只是离岸距离的函数，所有海岸线外面都套着一圈同样宽的浅色带

汇聚边界在陆地一侧、没有手绘山脉的地方补一条山脉，在海的一侧贴马里亚纳式的海沟和岛弧。离散边界在海里贴洋中脊，在陆上贴东非裂谷式的地形。做法和山脉的条带移植相同，来源也是 ETOPO。

大陆架的宽度取对数正态分布，中位数 110 km，范围 30 到 600 km，俯冲一侧偏窄，被动大陆边缘偏宽。

深海的深度用洋壳年龄估算：把到洋中脊的距离按每百万年 45 km 的半扩张速率折算成年龄，代入 Parsons–Sclater 的年龄-深度关系， $\text{深度} = 2600 + 320\sqrt{\text{年龄}}$ （米），再叠加真实的深海纹理。最后结果看起来就真实多了。

合成的结果

几项统计量和地球对照：陆地平均海拔 550 m（地球约 800 m），海拔 1000 m 以上的面积占 12.6%，3000 m 以上占 3% 到 4%，最高 6760 m，最深 -10.4 km，量级相符。

残存问题：非常平坦的平原上偶尔能看出平铺块之间的淡接缝，高原块体的轮廓偏明显，部分边缘还是有明显的分界线，山脉的一些走向还是不够自然。先留个坑吧后面再修。

气候：从气压场到柯本分类

地形定下来，海陆分布和山脉的位置就都有了，可以往下推气候。手绘图里原本也有一张气候分区图，是参照地球类比画的，推演的结果最后要和它对照：

要得到气候类型，需要气温和降水。气温和水汽由风从海上带到陆地，海面的温度受洋流影响，洋流由风驱动，风又由气压梯度驱动。追到头是气压，所以我们先从气压算起。

模型做了很大的简化。它不是大气环流模式，没有时间积分，也不解原始方程组，而是每个环节各用一个可以直接求稳态解的方程。

网格约 52 km 一格，经度方向周期，计算 1 月、春秋分、7 月三个时相，整套运行约 5 分钟。

气压

需要得到冬夏两季海平面气压的分布，尤其是高低压中心的位置，盛行风向和季风都由它们决定。气压场由三项叠加：

p = \underbrace{p_z(\varphi - 6^\circ c)}_{\text{纬向气压带}} - \underbrace{\alpha\, G_{650\,\text{km}}\big[(T - \overline{T}^{\lambda})\,\mathrm{eff}(\varphi)\big]}_{\text{海陆热力差异}} + \underbrace{A\, B(\varphi)\, G_{300\,\text{km}}[E]}_{\text{副高的东西不对称}}

其中 $c \in \{-1, 0, 1\}$ 为季节相位， $\alpha = 1.7\ \text{hPa}/^\circ\text{C}$ ， $\overline{T}^{\lambda}$ 是同纬度的平均气温。

第一项是三圈环流对应的纬向气压廓线：赤道低压、副热带高压、副极地低压、极地高压，整体随季节南北移动 6°。廓线用 PCHIP（保形分段三次插值）通过 13 个节点，普通三次样条会在节点之间过冲，产生多余的极值。

只有第一项的话，气压只随纬度变化。实际的气压中心来自海陆差异，这是第二项。它的依据是流体静力学：气柱偏暖则密度小，地面气压偏低。夏季陆地比同纬度的海洋暖，形成热低压，冬季相反。系数 eff(φ) 在低纬减小，因为科氏力弱的地方气压梯度很快被环流抹平；高纬也减小，因为冷气团浅薄。

只有前两项时，副热带高压的中心会偏向大洋西部，与观测相反，第三项用来修正这一点。它来自 Rodwell 和 Hoskins 的解释：夏季大陆上的加热激发 Rossby 波，在其西侧也就是大洋东部产生下沉。E 是一个指示场，标记东边 12 到 24 个经度内有陆地的洋面。

三项叠加之后，1 月北半球几块大陆上各出现一个 1029 到 1032 hPa 的冷高压，副极地洋面上是 986 hPa 左右的低压；7 月大陆上转为 1001 到 1008 hPa 的热低压，副热带高压的中心落在大洋上，这些中心的位置都由海陆分布去模拟算出来的。

风

气压场决定了风。中高纬的风接近地转风 $\mathbf{v} = \mathbf{k}\times\nabla p/(\rho f)$ ，但这个公式在赤道附近发散，因为科氏参数 f 趋于 0。而信风带和赤道辐合带恰好在低纬，不能回避。这里用带线性摩擦的定常动量方程，它有解析解：

f\,\mathbf{k}\times\mathbf{v} = -\frac{\nabla p}{\rho} - r\,\mathbf{v}

u = -\frac{r\,\partial_x p + f\,\partial_y p}{\rho\,(r^2 + f^2)}, \qquad v = \frac{f\,\partial_x p - r\,\partial_y p}{\rho\,(r^2 + f^2)}

f 远大于 r 时，解接近地转风，沿等压线吹，略偏向低压。f 趋于 0 时，解变成 −∇p/(ρr)，风顺着气压梯度直接吹向低压。两种情形之间平滑过渡，分母不会为零。

偏转角为 arctan(r/f)，陆地的摩擦系数取海洋的两倍，所以陆上风速更小，也更偏向低压。

越过赤道时 f 变号，偏转方向随之改变，南半球的东南信风越过赤道后转为西南风，不需要另外处理。

后面几步对风的要求不一样，摩擦系数因此取了三套。

地面风用于出图和计算辐合。
摩擦较小的自由低层风用来驱动洋流和输送热量。
摩擦更小、接近 850 hPa 的风用来输送水汽。对于由地面风输送的水汽，后面我们在模拟降水分布的时候会让他进不了内陆，确保严谨。

得到的风场里，低纬是信风，中纬是西风，极地是东风，风带的位置随季节南北移动。大陆东岸和南岸冬夏风向相反的区域，是后面判定季风气候的依据。

洋流

风吹在海面上，驱动洋流。洋流改变海面的温度，进而影响沿岸的气候：寒流经过的大陆西岸干燥，暖流经过的地方温和湿润。所以在求气温和降水之前，得先知道洋流。

洋流的主体是大洋尺度的风生环流，用 Stommel 模型[4]描述。模型是正压、线性的，带底摩擦，未知量是流函数 ψ：

r\nabla^2\psi + \beta\,\frac{\partial \psi}{\partial x} = \frac{\operatorname{curl}\boldsymbol{\tau}}{\rho}, \qquad \beta = \frac{\mathrm{d}f}{\mathrm{d}y}

右端的风应力就来自刚才的风场。大洋内区摩擦项很小，方程近似为 Sverdrup 平衡 $\beta\,\partial_x\psi = \operatorname{curl}\boldsymbol{\tau}/\rho$ ，风应力旋度决定经向的输送。

摩擦项只在西边界附近起作用，在那里形成宽度约 r/β 的窄而强的边界流，对应黑潮和湾流。

摩擦系数的取值使边界层宽约 150 km，大致三个网格，模拟的时候试过再窄网格就分辨不了了，解会振荡。

边界条件需要特别处理。这个世界有 42 块大陆和较大的岛，海洋是多连通区域。ψ 在每块陆地的海岸上是常数，但不同陆地的常数并不相同，两块陆地的 ψ 之差就是它们之间海峡的体积输送。如果把所有海岸都设成 ψ = 0，等于禁止海峡通流，绕极的西风漂流也就不存在。

处理办法是把每块大陆整体当作一个未知数节点，和海洋格点一起做有限体积离散。大陆节点对应的方程是摩擦通量沿整条海岸线的环路积分为零，这是 Godfrey 岛屿法则在 Stommel 模型里的离散形式。最后得到一个约 20 万未知数的稀疏线性方程组，直接求解。

在环流之上，表层再叠加风海流。Ekman 输送的辐散给出上升流，只保留近岸和大洋东部赤道的部分，对应沿岸寒流和赤道冷舌。

解出来的流函数在 −34 到 23 Sv 之间，最大流速约 0.5 m/s。副热带洋面上形成环流圈，西侧是向极的暖流，东侧是向赤道的寒流；中高纬有西风漂流；海峡里的流向和两侧的环流一致。

海温、气温、水汽：同一个输送方程求解

风和洋流都有了，可以求气温和降水了。海温、气温、水汽这三个场的物理过程相近：被流场携带，同时扩散，并且向某个平衡状态调整或者有源有汇。它们可以写成同一个方程，用同一个求解器：

\nabla\cdot(\mathbf{U} q) - \nabla\cdot(\kappa \nabla q) + \lambda q = S

离散用有限体积加一阶迎风。这样得到的系数矩阵是 M-矩阵，离散极值原理成立，解不会振荡，也不会出现负的水汽。代价是有一定的数值扩散，在 52 km 的网格上和物理扩散同量级，可以接受。方程不做时间积分，稀疏矩阵一次 LU 分解直接得到稳态。

先解海温。流场是上面的洋流，温度向随纬度变化的平衡值松弛，时间尺度 220 天，上升流区再降温。暖流把低纬的热量带向高纬，寒流相反。

有了海温再解气温。海上直接取海温；陆上沿风向从海洋的值向大陆平衡温度松弛，时间尺度 2.6 天，最后按每千米 6 ℃ 的递减率订正到地面。风速 8 m/s 时，松弛的特征距离约 1800 km。西风带里，空气登陆后一路向东，大陆西岸受海洋影响大，越往东大陆性越强，同纬度东西两岸的温度差异由此产生。

最后是水汽，方程是柱水汽的收支，海面蒸发，陆面有一部分降水重新蒸发。

降水里最重要的是地形雨，速率取 $\max(\mathbf{U}\cdot\nabla h,\ 0)/H_q$ ， $H_q = 2600$ m 是水汽标高，气流每爬升一个标高凝结掉 1/e 的水汽；迎风坡消耗了水汽，背风坡水汽少，又处在下沉气流里，雨影区由此形成。

其余的动力降水由低压、地面风辐合、暖季对流、温带风暴轴几项相加。降水速率依赖水汽本身，方程是非线性的，用 Picard 迭代求解。

到这一步，海陆位置、洋流、地形对气候的影响都已经体现在气温和降水里：暖流沿岸的气温高于同纬度，寒流沿岸偏低；山脉迎风一侧多雨，背风一侧少雨；大陆内部的降水随离海距离减少。

气候分类

气温和降水是连续的场。要和手绘的分区图对照，还需要归成离散的气候类型。

分类判据用的是逐月数据，而模型只算了三个时相，所以先要插值。三个时相用一次谐波加二次谐波插成 12 个月： $X(m) = a_0 + a_1\cos\theta + a_2\cos 2\theta$ ，三个系数由 1 月、7 月、春秋分三个值解出。

二次谐波表达的是春秋季和冬夏平均值的差别，赤道附近的双雨季特点的气候判定需要这一项。

之后按柯本分类[6]的判据分类，再归并为中学教科书上的 12 种气候类型。

粗网格的结果还要降尺度到约 13 km，和地形的分辨率相称。气温按细分辨率的地形做递减率订正。这里又是按尺度拆开：降水乘一个地形因子，这个因子只用粗网格分辨不了的那部分地形计算，避免地形雨被重复计入。

最后是清理分类图上的碎斑和锯齿。把每一类拆成一张 0/1 图，分别做高斯模糊，再逐像素取最大的那一类，相当于高斯加权的众数滤波。

图 11 · 左：逐像素分类的原始结果，边界呈锯齿状，夹着很多碎斑。右：清理之后，边界平滑了许多

陆地上各类型的面积占比：热带雨林 7.1%，热带草原 11.0%，热带季风 5.3%，热带沙漠 15.7%，亚热带季风和湿润 4.7%，地中海 4.8%，温带海洋性 5.1%，温带季风 0.8%，温带大陆性 13.5%，亚寒带针叶林 3.5%，苔原及冰原 24.1%，高原山地 4.4%。

回头和高二手绘的分区图对照，大格局是对得上的：赤道附近的雨林和草原、副热带的沙漠带、北方大陆上的温带大陆性气候和亚寒带针叶林，位置都和手绘差不多。

对不上的主要有三处。地中海气候比手绘多出不少，而且偏北；温带季风反而少了，只占陆地的 0.8%；右下角那块大陆的南缘，手绘画了一整条亚寒带针叶林，推演出来大多是温带大陆性气候和高原山地，针叶林只剩零星几块。

调整过程和现存问题

模型里的经验系数不少，需要对照结果调整。每次运行输出一张诊断总览图，前后调整了 9 轮。下面是第一轮和最后一轮：

图 13 · 左边第 1 轮，右边第 9 轮。从上到下依次是气压与风、气温、降水、季风指数、气候分类

其中几处主要的调整：

现象	原因	处理
赤道两侧各有一条很窄的雨带，内陆极干	气压廓线在赤道两侧各有一个槽；水汽由高摩擦的地面风输送	改为单一赤道槽；水汽改用低摩擦的风输送
洋流几乎都被判为寒流	只按海温距平判断	先看流向，向极为暖流，向赤道为寒流，其次才看海温
季风区成了横贯大洋的条带	只看冬夏风向反转，赤道辐合带的季节移动在洋面上也造成反转	再乘上夏雨占优、近海陆地两个条件
副高中心位于大洋西部	只有热力项时的固有偏差	加入东西不对称项
赤道 1 月气温只有 20 ℃	赤道上升流作用于整个洋盆	上升流限制在大洋东部

现存的问题：信风带的大陆东岸偏干；部分大陆西岸的地中海气候延伸到北纬 47 到 49 度，偏北；西北大陆内陆的年降水接近 0；全球陆地平均降水约 1030 mm，高于地球的约 800 mm。这也导致了后面河流的流量也整体偏大。

河网：最小代价汇流与洼地水量平衡

地形决定水往哪里流，气候决定有多少水，现在两样都有了，可以继续推理河网。

做地形的时候其实顺带出过一版河网，用的是最常见的三步：填洼，D8 最陡下降，汇水面积超过阈值的像元算作河道。把手绘的大河叠上去对比：

对比下来，还是问题不少的：沙漠和雨林的河网一样密；平原上有大量笔直、互相平行的河道；手绘的 12 条大河基本没有被还原；没有湖泊，也没有内流区。河网于是整个重做。

产流计算

先说疏密。旧河网只看汇水面积，面积够大就算河，所以沙漠和雨林没有区别，所以这一版要根据气候来修正。

河流的大小其实取决于流域里有多少水可以流走，也就是降水减去蒸发。上一章已经给出了降水和气温，可以先求每个像元的年径流深。用 Budyko 曲线[7]。先由月均温计算 Holdridge 生物温度，换算成潜在蒸发 PET，令 φ = PET/P：

\frac{E}{P} = \sqrt{\varphi\,\tanh\frac{1}{\varphi}\,\left(1 - e^{-\varphi}\right)}, \qquad R = P - E

曲线的两端对应两种极限。φ 很小时气候湿润，蒸发受能量限制，E 接近 PET。φ 很大时气候干旱，蒸发受水分限制，E 接近 P，径流趋于 0。

陆地（不含冰盖）平均径流深 532 mm。湿润区的径流深大，干旱区接近 0，河网的疏密由此拉开。像元面积乘 cos φ，高纬的汇水面积不再被高估。

汇流：不填洼的最小代价搜索

每个像元产多少水有了，接下来是这些水流向哪里。地形是移植的真实地形，洼地很多，而我们的汇流算法怎样处理洼地，直接决定了平原上的河道形态。

旧河网里那些笔直的河道，就出在填洼这一步。洼地被填成平面以后没有坡度，流向只能由算法的遍历顺序决定，于是成片地指向同一个方向。

于是我们改用最小代价搜索，思路来自 GRASS GIS 的 r.watershed。从海岸出发，用优先队列向内陆生长，每个像元流向最先到达它的那个邻居。

换用这个算法后，平原上的直线河道消失了，河道沿真实的谷地延伸。让Claude自检一下，河道像元里下游比上游高出 5 m 以上的只占 0.18%，还算不错。

手绘大河：走廊内的最小代价路径

直线河道没有了，但纯粹由地形决定的河网，不会恰好经过手绘的那 12 条大河的位置。这些大河关系到小说里的世界设定，需要保留。

另一方面，手绘线只是示意，水不一定流得过去，照着描，手绘线又成了直接写进结果的数据。

折中的办法是把手绘线当作一条走廊的中线，实际的河道在走廊里顺着地形选取，和山脉线扭曲那里差不多。

每条河在走廊内用 Dijkstra 搜索一条路径，每一步的代价为：

\text{代价} = \ell\left(W(\mathbf{x}) + \frac{d^2}{\sigma^2}\right) + w_{\text{up}}\,\Delta h^{+}

其中 $\ell$ 是这一步的步长， $\Delta h^{+}$ 是这一步的爬升高度。d 是到手绘线的距离，σ 控制走廊的宽度，这一项把路径留在手绘线附近。爬升项让路径沿真实的谷地走。

光有这两项，山区没有问题，平原上却不行：平原起伏小，爬升项几乎不起作用，只剩下偏离惩罚，最优路径就是手绘线本身，等于照描。W(x) = exp(1.5·fbm) 是为此加的，它是一个随机的通行难度场，有了它，路径在平原上会自然摆动。

路径找到以后，还要让前面的汇流服从它。用的是 stream burning：把河道像元的搜索键设为一个很小的负数，并且向上游递增。汇流搜索从河口开始，会先沿整条大河上溯到源头，两岸的像元随后汇入。

河谷：下切量的尺度分解

河道定下来以后，地形要反过来和它一致，因为我们之前的算法带来的小问题就是河流会在中间穿越一些小坎（见下图），那么真实的河流经过的地方应该是谷地，而且越往下游海拔越低才对。

目前的困难在于手绘的河和生成的地形并不配合：山脉线被随机扭曲过，有的河要翻过三四千米的山脊。要让河流得过去，沿河的高程必须单调下降，需要下切的量很大：

如果直接沿河下切，平均深度有 300 到 500 米，结果是高原上一条很深的窄缝，不像河谷。

真实的大河流经的地方，往往整片地势就是低的，只有穿过局部山梁的地方才是峡谷。这和地形那里的情形相似，办法也一样，对下切量做一次尺度分解：

把沿河的下切量沿河平滑，平滑的部分说明这一带的地势整体应该更低，把它向两侧展宽成约 150 km 的平滑场，从低频地势里减掉，地形细节保持不变。展宽时只有河道像元上有值，用归一化卷积处理，即分别对值和掩膜做高斯卷积再相除。

剩下的残差才是局部的山梁，切成峡谷，两岸按下切深度放坡。

处理之后峡谷的平均下切降到 85 米。有两条河的源头位于山脉的另一侧，沿途要翻过 3000 米以上的山脊，这种情况不硬切，适当地把河源移到山脊的近侧，各舍去 300 到 400 km 的上游。

湖泊和内流区

还剩湖泊。前面的汇流让水穿过洼地流走，相当于默认每个洼地都会被灌满并溢出。实际并非如此。地形上有一千多个较大的洼地，一个洼地是积水成湖再溢出，还是成为内流区的终点，不由地形单独决定，要看来水和蒸发哪个大。

对每个面积和深度足够大的洼地做水量平衡。湖面越高，湖面积越大，湖面的净蒸发也越大。预先算出蒸发量随水位变化的累积曲线，再和入湖径流比较：

水位达到溢出口之前，蒸发就和入流平衡：湖面停在平衡水位。这是内流湖，咸水，它的整个集水区是内流区。
灌到溢出口仍有余量：水经溢出口流出。这是吞吐湖，淡水。

判定为内流的洼地，把湖心设为汇流的终点，汇流重新做一遍。

最后保留了 69 个较大的湖，内流区占陆地面积的 23%（地球约为 18%），主要出现在干旱地区。手绘图上的北丰盆地没有做任何特殊处理，结果是一个内流盆地，盆底有咸水湖。

流量、河道和流域

汇流定了，剩下的是把它变成地图上的线和面。

流量沿汇流方向累积。干旱区的河段按干燥度指数沿程扣减，对应渗漏和蒸发损失，源于湿润山区、流经沙漠的河会越流越小。多年平均流量超过 160 m³/s、集水面积超过 1 万 km² 的像元算常年河，干旱区按集水面积另外补充季节河。

河道按主支流关系拆成干流和各级支流，线宽取 log₂(Q/Q₀)，向下游逐渐加粗。流域由汇流搜索时记录的根节点直接得到。

最终有干支流约 2800 条，大流域 214 个。最大的一条河干流长约 1 万 km，流域面积 1280 万 km²，河口流量约 6.6 万 m³/s。

极地冰盖：物质平衡判据与滞后阈值

到这里，已经完整地算过了一遍因果链。回头检查，两极还有一处不合理：冰盖的边界几乎是一条纬线。

原因在手绘图上。陆地图层里，两极各是一整条矩形的陆地带，做地形时把整行都是陆地的部分直接当作了冰盖。南极冰盖的边界因此就是南纬 75 度，地形上是一道平行于纬线的陡坎。

这又是手绘的形状直接变成了地形（后悔当初的手绘版本绘制的如此粗糙）。冰盖的范围本来应该由气候决定，而气温和降水现在都已经有了，那么问题不大，我们直接开干。

物质平衡判据

哪里一年积下的雪多于夏天化掉的雪，哪里才有冰。把这个条件写成一个标量场，也就是物质平衡判据：

S = T_c(P) + A\,N(\mathbf{x}) + \Delta T_0 - \left(T_w - \Gamma\, z_{\text{bed}}\right)

先看括号里的部分。T_w 是最热月的海平面气温，Γ = 6 ℃/km，z_bed 是冰下地形的高程，合起来是冰下地面的最热月气温。用最热月而不用年平均，是因为冰盖能否维持主要取决于消融季，消融量近似正比于正积温，冬天再冷也只是不融化。

Tc(P) 是允许的夏季气温上限，随降水增加： $T_c = 1.0 + 1.1\log_2(P/150\ \text{mm})$ ，限制在 −1.5 到 5 ℃。降水多则积累多，可以容忍更暖的夏天，海洋性冰川和大陆性冰川的差别就在这里。

S > 0 的地方积累大于消融，可以形成冰。式子里还有 N 和 ΔT₀ 两项，稍后再说。

从判据到冰盖边界

S 只描述每个点自身的条件。而冰盖是一个整体：冰会从积累区向外流动；冰面抬高以后，上面的气温也随之下降。需要一条规则把这种非局部的效应考虑进去。

规则是滞后阈值：冰盖范围取与 {S > 0} 连通的 {S > −1.5 ℃}。这是 Canny 边缘检测[12]里用的同一个方法，本质上是形态学重建[13]。

1.5 ℃ 大约相当于 250 m 的冰面抬升，是冰缘附近冰厚的量级。自身条件略差、但与积累区相连的地方，算作冰流能够到达的范围。

不与积累区相连的孤立区域，即使 S 略高于低阈值也不保留。冰盖内部不与海岸连通的无冰区也填为冰，理由同样是四周冰面抬高之后的降温。

如果只用单一阈值，会留下大量孤立的小斑块，冰盖主体的边缘又偏保守。

按上面的判据和规则算出来的边界，仍然接近纬线，因为模型给出的极区温度场几乎是纬向对称的。实际冰盖的非纬向形态来自模型没有分辨的过程，比如云量、风吹雪、海冰、冰流动力学。判据里的 N 用来代表这些因素，它是一个幅度约 3 ℃ 的相关噪声，对应冰缘约 ±4 个纬度的进退。

这个噪声必须在球面上各向同性。平面 fBm 直接按经纬度取值是不行的：纬度 75° 处一个像素东西向的实际长度只有赤道的四分之一，噪声贴到球上会成为南北向的细条，在 ±180° 经线处还不连续。

做法是在三维空间里生成噪声，再取球面上各点的值。每个倍频在 [−1,1]³ 上放一个随机格点，对每个像素的单位球坐标做三次样条插值。三维各向同性的场限制到球面上仍然各向同性，也没有接缝。

剩下 ΔT₀。判据里的阈值是物理量，对气温的绝对值很敏感。而模型的绝对温度有一两度的系统偏差，直接计算，冰盖面积可能相差一倍。

因此引入全球统一的偏移量 ΔT₀，使冰盖总面积（按球面面积计算，不按像素数）等于设定值，面积是 ΔT₀ 的单调函数，用二分法求解一定收敛。这样处理，偏移量只影响总量，冰盖的形状仍然完全由 S 的空间分布决定。

以上步骤在 1⁄2 分辨率上完成，确定冰盖的大形。再到全分辨率，只在边界附近的窄带里重新判定，这时加入细尺度的地形和更短波长的噪声。谷地比两侧低，气温高一些，冰缘在谷地后退，在山脊上前伸。

图 19 · 冰盖范围的变化。白色保留，浅蓝新增，橙色退去。上半是北极，下半是南极

替换冰穹

冰盖范围改变后，地形也要跟着改：新结冰的地方要有冰穹，退冰的地方要露出冰下的陆地，而且旧边界的位置不能留下痕迹。

当初合成地形时，冰面是这样得到的：

h = h_{\text{lf}} + \max\left(\text{dome} - h_{\text{lf}},\ 0\right)\cdot s + \text{细节}

hlf 是低频地势；dome 是冰穹剖面，随到冰缘的距离 d 增高；s 是冰盖的软掩膜。替换的办法是减去旧的第二项，加上新的第二项。旧的 dome 和 s 可以按原来的随机种子精确重算。未知的只有 hlf，它在冰穹内部被 max 运算覆盖了，只能估计：用无冰陆地的低频高程向冰盖内部做最近邻外推。

估计不准会有多大影响，可以直接算出来。设 hlf 的估计误差为 δ，结果的误差是：

\delta\,\left(s_{\text{old}} - s_{\text{new}}\right)

冰盖状态不变的地方 s_old = s_new，误差为零，与 δ 的大小无关。误差只出现在退冰和新结冰的过渡带，并且是一个低频量乘以一个光滑的权重，不会形成陡坎。这个结论的前提是 s_old 精确，所以旧冰盖必须按原来的种子重算，不能用估计值代替。

地形上的痕迹可以这样避免，冰盖范围本身却出过同样的问题：新的冰缘总是落回旧的纬线位置。原因在输入数据。判据用到的冰下高程在旧边界处有几百米的台阶，降尺度降水里也带有旧冰穹造成的地形雨台阶，任何一项都足以把新边界固定在原处。

排查的办法是把判据的每一项沿垂直于旧边界的方向输出剖面，看哪一项有阶跃，再分别处理。

改动之后

南极冰缘所在的纬度由全线 75 度变为在 63 到 78 度之间变化，有的地方冰舌伸到海岸，有的地方退到很高的纬度。冰盖总面积和原来相近，退冰区成为苔原，有河流发育。

冰盖和地形改过之后，气候和河网按前面的方法依次重算了一遍。判据用的是海平面气温，对冰穹的高度不敏感，所以这个循环走一遍就够，不需要迭代。

冰缘变弯之后，±180° 经线附近第一次出现了无冰的陆地，暴露出几处没有跨经线环绕的运算。连通域标记用并查集把首末两列相邻的区域合并，距离变换和平滑改为经向环绕。洼地填充函数原先也没有环绕，接缝旁的一个盆地因此被填到了 1962 m，修正后恢复正常。

河流的水文特征：逐月水量平衡

推河网时只用了年平均的降水和气温，每条河因此只有多年平均流量一个属性。而描述一条河，通常还要说它的汛期、补给类型、结冰期，以及是不是季节河。这些都取决于径流在一年之内怎样分配，需要月尺度的模型。逐月的气温和降水在推气候时已经算过，可以直接用。

径流模型

每个格点设三个状态量：积雪 W、土壤水 M、地下水 G。按月更新，运行三年，取最后一年：

\begin{aligned} \text{降雪：}\quad & P_s = P_m \cdot \operatorname{clip}\!\left(\frac{1.5 - T}{3},\ 0,\ 1\right), \qquad P_r = P_m - P_s \\ \text{融雪：}\quad & M_e = \min\big(W + P_s,\ k_d \max(T, 0)\big), \qquad k_d = 100\ \text{mm}/(^\circ\text{C}\cdot\text{月}) \\ \text{土壤：}\quad & M' = \max\left(M + P_r + M_e - \mathrm{PET}_m,\ 0\right), \quad X = \max\left(M' - 150\ \text{mm},\ 0\right), \quad M \leftarrow M' - X \\ \text{出流：}\quad & q = 0.5X + G/3, \qquad G \leftarrow G + 0.5X - G/3 \end{aligned}

其中 $P_r$ 是降雨， $X$ 是当月的产流。

这是 Thornthwaite–Mather[14-15] 月水量平衡加一个线性地下水库，参数取文献里的典型值。三个状态量各对应一个决定年内分配的过程。

积雪是季节性的存储。冬季的降水以雪的形式留在地面，气温回到 0 ℃ 以上后的一两个月内集中融化，形成春汛。

土壤要蓄满之后才产流。旱季过后，雨季最初一两个月的降水先补足土壤的亏缺，河流的汛期因此比雨季滞后。地中海气候区和热带草原气候区的河流是这种情况。

地下水库的出流与蓄量成正比，退水呈指数衰减，时间常数 3 个月。停止补给半年以后，基流还剩约 14%，一部分河流枯季不断流靠的就是它。

冰盖上不走这套模型，产流按 max(T + 2, 0) 为权重分配到各月，集中在夏季。

与河网的衔接

月模型给出的年径流总量，与前面 Budyko 曲线的结果并不相等。河网的阈值和分级都是按 Budyko 的年径流调定的，所以这里只取月模型的年内分配比例，年总量保持不变。已有河流的位置和年均流量因此都不受影响。

至于各月的流量怎么算，则是利用了汇流的一个性质：汇流对输入是线性的，Q = A·q，矩阵 A 只由汇流的拓扑关系和沿程损失决定，与 q 无关。12 个月的流量只需要在同一张汇流图上累积 12 次，每次是一趟 O(N) 的扫描。

融雪补给和冰盖补给的占比各再累积一次，得到的占比按流量加权，上游大支流的补给类型自然占主导。

特征与结果

由逐月流量得到各河段的特征：

特征	判据
汛期	月流量超过年平均 1.25 倍的月份，按 12 个月首尾相接合并成区间
补给类型	冰盖融水或积雪融水的份额超过 20% 即列出，其余为雨水
结冰期	河段所在位置月均温低于 −3 ℃ 的月数
季节河	最枯月流量低于 8 m³/s

季节河此前只按年均流量判断，现在旱季断流的河段会画成虚线，约 6700 个河道像元因此改判。

结果和各河所处的气候对得上。中纬度的几条大河是积雪融水和雨水混合补给，汛期在春季到夏季（4 到 6 月、4 到 7 月、6 到 8 月不等）。低纬的河是雨水补给，汛期跟随雨季，北半球 6 到 9 月，南半球 1 到 4 月。

有一条河的上游有高山冰帽，冰川融水的份额超过了 20%。6 号大河的河口位于北纬 74 度，年均流量约 23343 m³/s，最丰月为年均值的 4.1 倍，汛期 6 到 8 月，结冰期约 5 个月。

不过模型忽略了水在河道里的传播时间，也没有考虑湖泊的调蓄，大河下游的汛期可能比实际偏早一个月左右。

小结

行文至此，现在可以回收开头的那个遗憾了。高二的时候设计这张地图，凭借的仅仅是我高中的地理知识。

我知道迎风坡多雨、背风坡少雨，但说不清雨影区到哪里为止；知道寒流会让沿岸变干，但说不清它能影响到内陆多远。

而这次程序的推演在每一环都换成了可以计算的方程。以前只能说偏高、偏多的地方，现在有了数字：在我们这个世界里，冬季几块大陆上的冷高压在 1030 hPa 上下，最大的那条河河口流量有六万多立方米每秒，北纬 74 度入海的那条河每年要封冻五个月左右，内流区占了陆地面积的 23%。

当年的手绘图也没有白画。山脉、大河、盆地大多还在当年画的位置附近，只是从直接写进结果的数据，变成了计算要遵守的约束。

终于可以痛快地写小说了，哈哈哈哈哈。

最后需要说明的是，这套模型仍然很粗糙。它是一个稳态的简化模型，里面有不少经验系数，对着诊断图调了九轮才到现在的样子，并不是严格的模拟演化推理。现在依然存在信风带东岸偏干、全球降水偏多这些问题。

地图地址：World Map

附录：参数表

正文里略去的实现参数。

环节	参数	取值
地形 · 陆地背景平铺	平铺块的边长上限	源区域的 0.85 倍
地形 · 陆地背景平铺	气候类型权重场的平滑尺度	60 km
气候 · 网格	网格尺寸	776×400（约 52 km 一格）
气候 · 网格	微分算子与平滑	带球面度量；纬向平滑在频域实现，宽度随纬度调整，保证各纬度是同一个物理尺度
气候 · 洋流	风海流大小	风速的 0.9%
气候 · 洋流	风海流偏转角	45°·tanh(φ/10°)
气候 · 洋流	保留上升流的近岸范围	350 km
河网 · 手绘大河走廊	偏离惩罚的 σ	142 km
河网 · 手绘大河走廊	走廊上限	480 km
河网 · 手绘大河走廊	通行难度场 W(x) 的波长	约 180 km
冰盖 · 球面噪声	倍频数	5
冰盖 · 球面噪声	波长范围	800 km 到 50 km
冰盖 · 冰穹	冰穹剖面 dome	200 + 3000·√(d/1200 km)，d 为到冰缘的距离

参考文献

地形数据：[1] NOAA National Centers for Environmental Information. ETOPO 2022 15 arc-second global relief model[DB/OL]. NOAA National Centers for Environmental Information, 2022[2026-09-19].
海底深度：[2] Parsons B, Sclater J G. An analysis of the variation of ocean floor bathymetry and heat flow with age[J]. Journal of Geophysical Research, 1977, 82(5): 803-827.
副热带高压的东西不对称：[3] Rodwell M J, Hoskins B J. Subtropical anticyclones and summer monsoons[J]. Journal of Climate, 2001, 14(15): 3192-3211.
洋流：[4] Stommel H. The westward intensification of wind-driven ocean currents[J]. Transactions, American Geophysical Union, 1948, 29(2): 202-206.
洋流：[5] Godfrey J S. A Sverdrup model of the depth-integrated flow for the world ocean allowing for island circulations[J]. Geophysical & Astrophysical Fluid Dynamics, 1989, 45(1-2): 89-112.
气候分类：[6] Peel M C, Finlayson B L, McMahon T A. Updated world map of the Köppen-Geiger climate classification[J]. Hydrology and Earth System Sciences, 2007, 11(5): 1633-1644.
径流：[7] Budyko M I. Climate and life[M]. New York: Academic Press, 1974.
汇流：[8] Ehlschlaeger C. Using the A^T search algorithm to develop hydrologic models from digital elevation data[C]//Proceedings of International Geographic Information Systems (IGIS) Symposium '89. Baltimore, MD, 1989: 275-281.
汇流：[9] Metz M, Mitasova H, Harmon R S. Efficient extraction of drainage networks from massive, radar-based elevation models with least cost path search[J]. Hydrology and Earth System Sciences, 2011, 15(2): 667-678.
汇流：[10] GRASS Development Team. r.watershed: GRASS GIS manual[EB/OL]. [2026-09-19].
填洼：[11] Barnes R, Lehman C, Mulla D. Priority-flood: an optimal depression-filling and watershed-labeling algorithm for digital elevation models[J]. Computers & Geosciences, 2014, 62: 117-127.
滞后阈值：[12] Canny J. A computational approach to edge detection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1986, PAMI-8(6): 679-698.
形态学重建：[13] Vincent L. Morphological grayscale reconstruction in image analysis: applications and efficient algorithms[J]. IEEE Transactions on Image Processing, 1993, 2(2): 176-201.
月水量平衡：[14] Thornthwaite C W, Mather J R. The water balance[M]. Centerton, NJ: Drexel Institute of Technology, Laboratory of Climatology, 1955. (Publications in Climatology, 8(1)).
月水量平衡：[15] McCabe G J, Markstrom S L. A monthly water-balance model driven by a graphical user interface[R/OL]. Reston, VA: U.S. Geological Survey, 2007[2026-09-19]. (Open-File Report 2007-1088).

前言