10.2 分形与地形
站内 /graphics-terrain-large-world/
那篇关心的是怎么把一座大地形高效画出来：LOD 层级、四叉树切分、GPU 上的三角形预算。本节关心的是它的上游一步——这座地形的高度值是怎么算出来的。卷三的场景《烛影山谷》只有一份地形，定义在 config/scene_spec.json 里：gridN = 16、extentQ16 = 1441792，也就是 289 个顶点、512 个三角形。本节不另立一份地形，只把它「怎么来」这件事讲清楚，并逐顶点验证算出来的就是规格里那一份。
一个 octave 太光滑，所以要叠
上一节的噪声在 4.0 世界单位的格距上平滑变化——直接当高度场用，地形会像一块被揉过的布：有起伏，但没有细节层次。真实的侵蚀地形在不同尺度上都该有起伏，这正是分形的用武之地。
fBm（fractal Brownian motion）的写法只有一行：
h(x, z) = Σ_k gain^k · noise(lacunarity^k · x, lacunarity^k · z, seed_k) / Σ_k gain^k
三个参数各管一件事：octaves 管叠几层，lacunarity 管每层频率翻几倍，gain 管每层振幅缩多少。本节取 lacunarity = 2、gain = 0.5、octaves = 4，这是地形生成里最常见的一组。
三个整数化决定
lacunarity = 2用左移实现。fx << k就是把采样坐标乘2^k，等价于格点间距减半、频率翻倍。整数左移没有精度损失，也不需要pow。gain用 Q16.16 连乘。gain = 32768就是0.5，第 k 层的权重是gpow，每层结束后gpow = fmul(gpow, gain)。权重序列是ONE, 32768, 16384, 8192, …，全整数。- 归一化用
fdiv而不是//。分母是增益和gsum（4 层时是ONE + 32768 + 16384 + 8192 = 122880），分子分母都非负，但fdiv是显式的向零截断，与 TS 侧BigInt的/语义一致；Python 的//是向下取整，负数时会对不上。
每层的哈希种子由卷三底座的 SplitMix64（cgtrace.rng_next）从场景 seed 逐层派生，取 z % 2^31 后作为 noiseSeed。两侧跑的是同一串状态推进，所以种子序列逐位相同——这是 10.3「同一 seed 同一场景」的起点。
{
"name": "lantern-fbm-terrain",
"seed": 1337,
"octaves": 4,
"referenceOctaves": 8,
"gainQ16": 32768,
"lacunarityShift": 1,
"latticeShift": 18,
"profileSamples": 512,
"profileStepQ16": 65536,
"roughnessLags": [1, 2, 4, 8],
"statsGridN": 64,
"expectedTerrainVertices": 289,
"expectedTerrainTriangles": 512,
"expectedSceneVertices": 487,
"expectedSceneTriangles": 896
}
expected* 四个字段是场景规格一致性的断言口径：地形 289 顶点 / 512 三角形，加上 3 个球共 487 顶点 / 896 三角形。脚本从规格读参数重建一遍，与这四个数逐项比对——这样「本章的地形就是规格里那一份」不是一句话，而是一次检查。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_10_2_fbm_terrain.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_10_2_fbm_terrain.py
"""10.2《分形与地形》：fBm 叠加、高度场构造与「与场景规格一致」的实测。
一个 octave 的噪声太平滑，撑不起地形。fBm 把它叠起来：第 k 层用 `lacunarity^k` 的频率、
`gain^k` 的振幅，再整体归一化回 [0, ONE)。本节量四件事：
1. **叠加**：`lacunarity = 2` 用**左移采样坐标**（`fx << k`）实现，`gain` 用 Q16.16 连乘，
全程整数、无浮点；每层哈希种子由卷三底座的 SplitMix64 从场景 seed 派生。
2. **与场景规格一致**：地形必须就是 `config/scene_spec.json` 那一份——`gridN = 16`、
`extentQ16 = 1441792`，289 个地形顶点 / 512 个地形三角形，加 3 个球共 487 / 896。
3. **代理指标**：高度场统计、与单层参照和与 8 层参照解的最大 / 平均偏差、结构函数斜率。
这些只用来做**相对比较**（层数多寡、参数改动），不冒充「地形好不好看」。
4. **1 层 == 底座值噪声**：`octaves = 1` 时归一化因子恰为 ONE，fBm 必须**逐位等于**
`cglib.terrain_height`——这条恒等式是「叠加没有引入额外误差」的判据。
产物：out/cg_trace_10_2_fbm.json、out/cg_trace_10_2_profile.csv、out/cg_trace_10_2_height.png
"""
import json
import math
import os
import sys
import numpy as np
HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
sys.path.insert(0, HERE)
import cglib # noqa: E402
import cgtrace # noqa: E402
from cg_trace_10_noise import (autocorr, digest_i32, fbm, loglog_slope, octave_seeds, # noqa: E402
roughness, stats)
ONE = cglib.ONE
fmul = cglib.fmul
fdiv = cglib.fdiv
CONFIG = os.path.join(ROOT, "config", "cg_trace_10_2_fbm.json")
def fbm_height(fx, fz, seeds, gain, shift, height_scale):
"""归一化后的 fBm（[0, ONE)）乘上高度尺度，得到世界单位下的高度。"""
return fmul(fbm(fx, fz, seeds, gain, shift), height_scale)
def main():
cfg = json.load(open(CONFIG, encoding="utf-8"))
spec = cglib.load_spec()
terr = spec["terrain"]
hs, shift, gain = terr["heightScaleQ16"], cfg["latticeShift"], cfg["gainQ16"]
extent, grid_n, oct_k = terr["extentQ16"], terr["gridN"], cfg["octaves"]
seeds = octave_seeds(cfg["seed"], cfg["referenceOctaves"])
rng_state = cfg["seed"] & cgtrace.RNG_MASK
for _ in range(cfg["referenceOctaves"]):
rng_state, _z = cgtrace.rng_next(rng_state)
out_dir = os.path.join(ROOT, "out")
os.makedirs(out_dir, exist_ok=True)
# ---- 1. 与场景规格一致：顶点数 / 三角形数 / 逐顶点高度 ----
tv, tt = cglib.build_terrain(grid_n, extent, hs, terr["noiseSeed"])
tris, _ = cgtrace.build_scene_tris(spec)
scene_verts = len(tv) + sum(len(cglib.build_sphere(s["centerQ16"], s["radiusQ16"], s["subdiv"])[0])
for s in spec["spheres"])
counts = {"terrainVertices": len(tv), "terrainTriangles": len(tt),
"sceneVertices": scene_verts, "sceneTriangles": len(tris),
"extentQ16": extent, "gridN": grid_n}
spec_ok = (len(tv) == cfg["expectedTerrainVertices"] and len(tt) == cfg["expectedTerrainTriangles"]
and scene_verts == cfg["expectedSceneVertices"]
and len(tris) == cfg["expectedSceneTriangles"] and extent == 1441792)
# ---- 2. 1 层 == 底座值噪声（逐位恒等式） ----
one_layer_max = 0
for (x, _y, z) in tv:
one_layer_max = max(one_layer_max,
abs(fbm_height(x, z, seeds[:1], gain, shift, hs)
- cglib.terrain_height(x, z, seeds[0], hs)))
# ---- 3. 高度场：4 层 vs 单层 vs 8 层参照 ----
k4 = [fbm_height(x, z, seeds[:oct_k], gain, shift, hs) for (x, _y, z) in tv]
# 单层参照用**同一层 0 种子**：偏差只反映「多叠了三层」，不掺 seed 之差
k1 = [cglib.terrain_height(x, z, seeds[0], hs) for (x, _y, z) in tv]
k8 = [fbm_height(x, z, seeds, gain, shift, hs) for (x, _y, z) in tv]
dev_k1 = [abs(k4[i] - k1[i]) for i in range(len(tv))]
dev_k8 = [abs(k4[i] - k8[i]) for i in range(len(tv))]
# ---- 4. 剖面、结构函数与自相似性指数 ----
n, step = cfg["profileSamples"], cfg["profileStepQ16"]
prof4 = [fbm_height(i * step, 0, seeds[:oct_k], gain, shift, hs) for i in range(n)]
prof8 = [fbm_height(i * step, 0, seeds, gain, shift, hs) for i in range(n)]
lags = cfg["roughnessLags"]
rough = {str(L): roughness(prof4, L) for L in lags}
local = {str(lags[i]): round(math.log(rough[str(lags[i + 1])] / rough[str(lags[i])])
/ math.log(2.0), 6) for i in range(len(lags) - 1)}
slopes = {str(K): round(loglog_slope(
[fbm_height(i * step, 0, seeds[:K], gain, shift, hs) for i in range(n)], [1, 2, 4]), 6)
for K in range(1, cfg["referenceOctaves"] + 1)}
# 口径敏感性：同一个斜率对「lag 窗口」和「采样点数」都敏感，故它只能同口径相比
sens = {f"lags{lab}": round(loglog_slope(prof4, lg), 6) for lab, lg in
(("[1,2,4]", [1, 2, 4]), ("[1,2,4,8]", [1, 2, 4, 8]),
("[1..32]", [1, 2, 4, 8, 16, 32]), ("[1..128]", [1, 2, 4, 8, 16, 32, 64, 128]))}
sens["4096pts lags[1,2,4]"] = round(loglog_slope(
[fbm_height(i * step, 0, seeds[:oct_k], gain, shift, hs) for i in range(4096)], [1, 2, 4]), 6)
# ---- 5. 高度场图（64x64，覆盖整幅 extent） ----
g = cfg["statsGridN"]
step2 = fdiv(2 * extent, g << 16)
field = [fbm_height(-extent + fmul(step2, ix << 16), -extent + fmul(step2, iz << 16),
seeds[:oct_k], gain, shift, hs)
for iz in range(g) for ix in range(g)]
cglib.save_gray_png(os.path.join(out_dir, "cg_trace_10_2_height.png"),
np.array(field, dtype=np.int64).reshape(g, g))
cglib.write_csv(os.path.join(out_dir, "cg_trace_10_2_profile.csv"),
["i", "xQ16", "fbm4Q16", "fbm8Q16"],
[[i, i * step, prof4[i], prof8[i]] for i in range(n)])
st4, st1, st8 = stats(k4), stats(k1), stats(k8)
report = {
"seed": cfg["seed"], "octaves": oct_k, "referenceOctaves": cfg["referenceOctaves"],
"gainQ16": gain, "lacunarityShift": cfg["lacunarityShift"], "latticeShift": shift,
"heightScaleQ16": hs, "octaveSeeds": seeds, "rngFinalStateHex": hex(rng_state),
"specConsistency": {**counts, "ok": spec_ok},
"oneLayerIdentityMaxQ16": one_layer_max,
"heightField": {"fbm": st4, "singleOctave": st1, "reference": st8},
"deviationVsSingleOctave": {"maxQ16": max(dev_k1), "maeQ16": round(sum(dev_k1) / len(tv), 4)},
"deviationVsReference": {"maxQ16": max(dev_k8), "maeQ16": round(sum(dev_k8) / len(tv), 4)},
"roughnessQ16": rough, "roughnessLocalSlope": local,
"profileAutocorrPermille": {str(L): autocorr(prof4, L) for L in [1, 2, 4]},
"selfSimilarityByOctaves": slopes, "selfSimilaritySensitivity": sens,
"profileDigest4": digest_i32(prof4), "profileDigest8": digest_i32(prof8),
"heightFieldDigest4": digest_i32(k4), "heightFieldDigest8": digest_i32(k8),
}
with open(os.path.join(out_dir, "cg_trace_10_2_fbm.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
fh.write("\n")
print("== 10.2 分形与地形 ==")
print(f"octave 种子 : {seeds[:oct_k]}")
print(f"场景规格一致 : 地形 {counts['terrainVertices']} 顶点 / {counts['terrainTriangles']} 三角形，"
f"全场景 {counts['sceneVertices']} / {counts['sceneTriangles']}，extent {extent} -> {spec_ok}")
print(f"1 层 == 值噪声 : 最大偏差 {one_layer_max} Q16.16")
print(f"4 层高度场 : mean {st4['meanQ16']} std {st4['stdQ16']} "
f"range [{st4['minQ16']}, {st4['maxQ16']}] 幅 {st4['rangeQ16']}")
print(f"单层 / 8 层参照 : 幅 {st1['rangeQ16']} / {st8['rangeQ16']}")
print(f"偏差 vs 单层 : max {max(dev_k1)} MAE {report['deviationVsSingleOctave']['maeQ16']} Q16.16")
print(f"偏差 vs 8 层参照 : max {max(dev_k8)} MAE {report['deviationVsReference']['maeQ16']} Q16.16")
print("结构函数 Q16.16 : " + " ".join(f"L{k}={v}" for k, v in rough.items()))
print("局部斜率 : " + " ".join(f"{k}->{int(k) * 2}={v}" for k, v in local.items()))
print("自相似性 vs 层数 : " + " ".join(f"K{k}={v}" for k, v in slopes.items()))
print("口径敏感性 : " + " ".join(f"{k}={v}" for k, v in sens.items()))
return 0
if __name__ == "__main__":
sys.exit(main())
第 3 步里有一个容易被忽略的细节：单层参照必须用同一个第 0 层种子，而不是规格里的 noiseSeed。第一版写成了 cglib.terrain_height(x, z, terr["noiseSeed"], hs)，结果测出来的偏差高达六位数——那个偏差量的是两个不同 seed 的场之差，跟「叠了几层」没关系。
# snippet: 本机实跑（工作区 lantern/）
/opt/homebrew/bin/python3 sim/cg_trace_10_2_fbm_terrain.py
== 10.2 分形与地形 ==
octave 种子 : [332046347, 1742079433, 1409821530, 2129642082]
场景规格一致 : 地形 289 顶点 / 512 三角形，全场景 487 / 896，extent 1441792 -> True
1 层 == 值噪声 : 最大偏差 0 Q16.16
4 层高度场 : mean 97483 std 26415 range [35421, 168075] 幅 132654
单层 / 8 层参照 : 幅 184836 / 119754
偏差 vs 单层 : max 67551 MAE 18108.9654 Q16.16
偏差 vs 8 层参照 : max 7551 MAE 1956.2699 Q16.16
结构函数 Q16.16 : L1=16506 L2=24229 L4=34818 L8=38430
局部斜率 : 1->2=0.553744 2->4=0.523098 4->8=0.1424
自相似性 vs 层数 : K1=0.801562 K2=0.698804 K3=0.579878 K4=0.538421 K5=0.53584 K6=0.533568 K7=0.533106 K8=0.53347
口径敏感性 : lags[1,2,4]=0.538421 lags[1,2,4,8]=0.418083 lags[1..32]=0.23029 lags[1..128]=0.122529 4096pts lags[1,2,4]=0.505766
这些数字怎么读
1 层 == 值噪声 的最大偏差是 0，这是本节最强的一条结论：octaves = 1 时归一化因子 gsum = ONE，fdiv(acc, ONE) == acc，所以 fBm 逐位退化成 cglib.terrain_height。它证明叠加这件事本身没有引入额外误差——后面所有偏差都只可能来自「多叠了几层」。
4 层高度场的幅是 132654（mean 97483、std 26415），而单层的幅是 184836。幅反而变小，是因为单层噪声是纯 [0, ONE) 的均匀插值，端点能取到接近 0 和 1 的极值；叠了四层之后，高层的细节把极值「拉平」了，分布更集中。8 层参照解的幅是 119754，比 4 层更小——继续加层，幅会继续收窄并收敛。
偏差是有方向的：与单层参照的最大偏差 67551、MAE 18109（这是换掉了场，量级自然是整个高度尺度 196608 的几分之一）；与 8 层参照的最大偏差只有 7551、MAE 1956（这是同一族场只差层数，量级小一个数量级）。两个数字放在一起才说明问题：4 层与 8 层已经足够接近，说明层数已经接近收敛；4 层与单层差得远，说明「叠」这件事确实在起作用。
自相似性：代理指标，只做相对比较
地形「像不像分形」不能靠眼睛判断。本节用结构函数：相隔 L 个采样点的平均高度差 E|h(x+L) − h(x)|。对自相似的场，这个量随 L 幂律增长，双对数坐标下的斜率就是自相似性指数 H。
实测结构函数是 L1 = 16506 → L2 = 24229 → L4 = 34818 → L8 = 38430，相邻两点的局部斜率分别是 0.553744、0.523098、0.1424——注意最后一段（L4→L8）已经明显塌下来，因为 4 层 fBm 在格距 4.0 世界单位处截止了，比这更粗的尺度上不再有新的细节。
把层数从 1 加到 8，log-log 最小二乘斜率从 0.801562 单调下降到 0.533470（K4 时 0.538421，K6 时 0.533568），K≥4 之后变化已经很小。
但这个数不是 Hurst 指数的估计，它强烈依赖量它的口径。脚本里专门加了一组敏感性检查：同一份 4 层剖面，把 lag 窗口从 [1,2,4] 放宽到 [1,2,4,8]、[1..32]、[1..128]，斜率从 0.538421 掉到 0.418083、0.23029、0.122529；把采样点从 512 加到 4096（仍是 [1,2,4]），斜率从 0.538421 变成 0.505766。
为什么会掉这么多？因为这份地形只有 4 层 octave，最粗一层在 32 个世界单位处就截止了：lag 一旦接近这个尺度，E|h(x+L) − h(x)| 就不再增长（已经饱和），斜率自然往下掉。局部斜率那一行已经露了头——4->8 只有 0.1424，而 1->2 是 0.553744。所以这个指标只在同口径下才有意义。
口径必须写进结论里：本章所有自相似性数字都限定在「512 点剖面 + lag [1,2,4] + 4 层」这一个口径内。在这个口径下，它能说的是「4 层已经到 0.538、8 层到 0.533，层数接近饱和」；它不能说的是「地形是分形」，也不能说的是「比另一个 0.533 的地形好看」。这就是把质量写成代理指标的意思——数字只有带着口径才能读。
P 层把高度场钉死
P 层断言三件事：octave 种子序列与 S 层逐位相同、289 个顶点的高度序列 sha256 与 S 层相同、octaves = 1 的恒等式在每个顶点上成立。
# snippet: core/test/cg_trace_10_2_fbm_terrain.test.ts（节选：恒等式与负例）
# verify: node --test 'core/test/cg_trace_10_2_fbm_terrain.test.ts'
test("恒等式：octaves = 1 时 fBm 逐位等于底座值噪声", () => {
for (const [x, z] of XZ) {
assert.equal(fbmHeight(x, z, [seeds[0]]), terrainHeight(x, z, seeds[0], HS));
}
assert.equal(golden.oneLayerIdentityMaxQ16, 0);
});
test("负例：lacunarity 换成 4（fx << 2k）高度场必须变", () => {
assert.notEqual(digest(wrong), golden.heightFieldDigest4, "换 lacunarity 必须换场");
});
负例把 fx << k 改成 fx << 2k（等价于 lacunarity = 4），摘要立刻分叉——这条负例守的是「lacunarity 真的参与了计算」，而不只是配置里写了个数字。
小结
fBm 把单个 octave 的噪声叠成有细节层次的高度场：lacunarity = 2 由左移采样坐标实现，gain = 0.5 由 Q16.16 连乘实现，归一化用向零截断的 fdiv，全程整数。本节的地形就是 scene_spec.json 那一份（289 顶点 / 512 三角形 / 全场景 896），octaves = 1 时逐位等于底座值噪声（最大偏差 0），4 层与 8 层参照解的最大偏差只有 7551，结构函数斜率随层数收敛到 0.533——作为相对比较的代理指标，它说明 4 层已经接近饱和。下一节把这份高度场推进成完整的程序化场景：加位移、分材质、渲染，并用 sha256 把「同一 seed 同一场景」钉住。
阅读导航：上一节：10.1 噪声的构造 · 下一节：10.3 程序化场景的生成 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。