5D Library › 复杂系统相变、临界慢化与生态韧性的非线性动力学 · English
第二章 临界慢化：恢复为什么越来越慢
变浅的山谷
上一章的折叠点有一个几何直观：系统所在的稳定状态像一个山谷，压力增加时山谷变浅、向旁边的山脊靠拢，到折叠点时两者合并、消失。山谷变浅的直接后果是：受到同样的推动之后，小球滚回谷底的速度变慢了。这就是临界慢化。
维塞尔在 1984 年把这一点写成一条「普适定律」：系统接近门槛时，受扰动后回到平衡的特征时间会增长，其增长方式可以用来从远离门槛处的测量外推门槛的位置；他用轮虫种群的初步收获实验检验了这一方法 (Wissel 1984)。范内斯与谢弗在六个生态模型中发现，受小扰动后恢复变慢在门槛之前很早就能看出来，并指出恢复速率（后文的「工程韧性」）可以作为吸引盆大小（「生态韧性」）的指标 (Nes and Scheffer 2007)。
平方根定律
推导
在平衡点 附近，令 ，一阶展开得 。 是恢复速率， 是特征恢复时间。在鞍结分岔处 ；在分岔附近
所以 随 近似线性地降到零 (Strogatz 2015)。维塞尔的外推正是利用这一点：在远离门槛的几个压力水平上测量恢复时间，画出 对压力的直线，它与横轴的交点就是门槛的估计。
下面的程序对放牧模型上方的稳定状态计算 、 与 （若平方根定律成立，这一列近似不变），做两种「脉冲实验」（把植被一次性减少 2% 或 10%，记录恢复 90% 的时间），并在三段压力范围内拟合 、外推门槛。
推导
# Recovery slows near a fold: return time of the upper (vegetated) state of the grazing
# model dx/dt = r x (1 - x/K) - c x^2/(x^2 + a^2) (r = 1, K = 10, a = 1) as c -> c2 = 2.604.
# Linearisation: lam = -f'(x*), return time T = 1/lam. Near a saddle-node lam ~ sqrt(c2 - c),
# so 1/T^2 falls roughly linearly to 0 at c2: the basis of Wissel's (1984) extrapolation.
import math
R, K, A, C2 = 1.0, 10.0, 1.0, 2.604365
def f(x, c):
return R * x * (1 - x / K) - c * x * x / (x * x + A * A)
def upper_eq(c):
lo, hi = 4.79, 10.0 # the upper equilibrium lies above the fold x = 4.78
for _ in range(80):
mid = (lo + hi) / 2
if f(mid, c) > 0: lo = mid
else: hi = mid
return lo
def lam(c, h=1e-6):
x = upper_eq(c)
return -(f(x + h, c) - f(x - h, c)) / (2 * h)
def pulse_recovery(c, kick, dt=0.001):
"""Push x down by a share kick; time until 90 % of the displacement has decayed."""
xs = upper_eq(c); x, t = xs * (1 - kick), 0.0
while xs - x > 0.1 * kick * xs:
x += f(x, c) * dt; t += dt
if x < 3.0: return None # fell into the lower basin: it tipped
return t
print(" c x* lam T = 1/lam T*sqrt(c2-c) pulse 2 %: 90 % back after pulse 10 %")
for c in (1.0, 1.5, 2.0, 2.3, 2.5, 2.58, 2.60):
l, big = lam(c), pulse_recovery(c, 0.10)
print(f" {c:4.2f} {upper_eq(c):4.2f} {l:5.3f} {1 / l:7.2f} {math.sqrt(C2 - c) / l:5.2f}"
f" {pulse_recovery(c, 0.02):7.2f} "
+ (f"{big:6.2f}" if big else "TIPPED"))
# Wissel-style forecast: fit 1/T^2 = lam^2 against c using only c <= cmax, extrapolate to 0
for cmax in (2.0, 2.3, 2.5):
top = int(cmax * 100)
pts = [(c / 100, lam(c / 100) ** 2) for c in range(top - 40, top + 1, 5)]
n = len(pts); mx = sum(p for p, _ in pts) / n; my = sum(q for _, q in pts) / n
slope = sum((p - mx) * (q - my) for p, q in pts) / sum((p - mx) ** 2 for p, _ in pts)
print(f"fit on c in [{cmax - 0.4:.1f}, {cmax:.1f}]: predicted threshold"
f" c = {mx - my / slope:.3f} (true {C2:.3f})")
# c x* lam T = 1/lam T*sqrt(c2-c) pulse 2 %: 90 % back after pulse 10 %
# 1.00 8.89 0.781 1.28 1.62 2.98 3.09
# 1.50 8.20 0.645 1.55 1.63 3.61 3.76
# 2.00 7.32 0.473 2.11 1.64 4.92 5.18
# 2.30 6.59 0.333 3.01 1.66 7.03 7.50
# 2.50 5.84 0.192 5.20 1.68 12.25 13.56
# 2.58 5.30 0.092 10.90 1.70 26.21 32.59
# 2.60 5.00 0.038 26.00 1.72 66.09 TIPPED
# fit on c in [1.6, 2.0]: predicted threshold c = 2.584 (true 2.604)
# fit on c in [1.9, 2.3]: predicted threshold c = 2.592 (true 2.604)
# fit on c in [2.1, 2.5]: predicted threshold c = 2.598 (true 2.604)从 1.62 缓慢变到 1.72，平方根定律在很大范围内近似成立。压力从 1.0 增加到 2.5 时，恢复时间从 1.28 增加到 5.20，到 2.60 时为 26.0；2% 的扰动恢复 90% 所需的时间从 2.98 增加到 66.1。一个细节值得注意： 时，10% 的扰动已经把系统推过了不稳定平衡点，系统翻转了。离门槛越近，吸引盆越窄，脉冲实验本身可能就是把系统推过门槛的那一下。
维塞尔外推在这个没有噪声的模型中效果很好：用 1.6—2.0 的数据外推得到 2.584，用 2.1—2.5 的数据得到 2.598，真值是 2.604。但这里的 是精确计算的；真实测量中，恢复速率要从带噪声的数据中估计。
噪声中的慢化
更常见的情形是：系统不断受到小的随机扰动，我们只能观察它的波动。
推导
推导。 平衡点附近的线性化系统加上白噪声，就是奥恩斯坦—乌伦贝克过程 ，它的解是
第二项与 独立、均值为零。两边乘以 取期望，得到间隔 的一阶自相关 ；由伊藤等距公式，平稳方差为
时，，方差 ：系统「记得」越来越久以前的扰动，波动也越来越大 (Scheffer et al. 2009)。反过来，。程序用欧拉—丸山法（步长 0.01）模拟，每 1 个时间单位采样一次，共 20,000 个样本。
推导
# Critical slowing down in its simplest form: the Ornstein-Uhlenbeck process
# dx = -lam x dt + sigma dW (lam = recovery rate, the linearised restoring force)
# Sampled every dt_s: lag-1 autocorrelation AC(1) = exp(-lam dt_s), variance = sigma^2/(2 lam).
# Check both formulas by Euler-Maruyama simulation (step 0.01, sampled every dt_s = 1).
import math, random
rng = random.Random(3)
def simulate(lam, sigma, n=20000, h=0.01, every=100):
x, out, s = 0.0, [], sigma * math.sqrt(h)
for _ in range(2000): # burn-in
x += -lam * x * h + s * rng.gauss(0, 1)
for i in range(n * every):
x += -lam * x * h + s * rng.gauss(0, 1)
if i % every == every - 1:
out.append(x)
return out
def ac1_var(y):
m = sum(y) / len(y)
v = sum((u - m) ** 2 for u in y) / len(y)
c = sum((y[i] - m) * (y[i + 1] - m) for i in range(len(y) - 1)) / (len(y) - 1)
return c / v, v
print(" lam sigma AC(1) theory simulated Var theory simulated lam = -ln AC(1)")
for lam, sigma in ((1.0, 0.5), (0.5, 0.5), (0.2, 0.5), (0.05, 0.5), (0.5, 1.0)):
a, v = ac1_var(simulate(lam, sigma))
print(f" {lam:4.2f} {sigma:4.2f} {math.exp(-lam):.3f} {a:.3f}"
f" {sigma ** 2 / (2 * lam):7.3f} {v:7.3f} {-math.log(a):.3f}")
print("halving lam: AC(1) rises and Var doubles; doubling sigma: Var x4, AC(1) unchanged")
# lam sigma AC(1) theory simulated Var theory simulated lam = -ln AC(1)
# 1.00 0.50 0.368 0.351 0.125 0.125 1.046
# 0.50 0.50 0.607 0.609 0.250 0.252 0.495
# 0.20 0.50 0.819 0.820 0.625 0.631 0.199
# 0.05 0.50 0.951 0.954 2.500 2.636 0.048
# 0.50 1.00 0.607 0.612 1.000 1.004 0.491
# halving lam: AC(1) rises and Var doubles; doubling sigma: Var x4, AC(1) unchanged公式与模拟吻合： 从 1 降到 0.05 时，自相关从 0.37 升到 0.95，方差从 0.125 升到 2.5。 时模拟的方差比理论高约 5%（2.636 对 2.500）：恢复越慢，相邻样本越相关，20,000 个样本中有效独立的观测只有几百个。越接近门槛，指标越大，估计它所需的数据也越多。 最后一行说明另一件事：把噪声强度加倍，方差变为四倍（0.252 → 1.004），自相关几乎不变（0.609 → 0.612）。单是方差上升，可能只是外界扰动变大了。 迪特列夫森与约翰森据此强调，预警需要方差与自相关同时上升，只依据其中一个信号得出的结论是无效的 (Ditlevsen and Johnsen 2010)。
一条接近门槛的轨迹
推导
把噪声加进放牧模型（），让压力在 1,000 个时间单位里从 1.0 升到 2.8，每个时间单位记录一次。找到崩溃的时刻（），去掉最后 10 个点，用高斯核平滑去趋势（带宽 25），在长度为记录一半的移动窗口中计算残差的一阶自相关与方差，用肯德尔 衡量上升趋势（参照达科斯等人 2008、2012 年的做法）。
推导
# Early-warning indicators before a fold: dx = [x (1 - x/10) - c x^2/(x^2+1)] dt + 0.15 dW,
# c rising from 1.0 to 2.8 over 1,000 time units (fold at 2.604), one sample per unit.
# Before the shift: Gaussian-kernel detrending (bandwidth 25); window = half the record;
# AC(1) and variance of the residuals; trend by Kendall's tau (after Dakos et al. 2008, 2012).
import math, random
def simulate(seed, T=1000, sigma=0.15, c0=1.0, c1=2.8, dt=0.05):
rng = random.Random(seed); x, xs, cs = 8.9, [], []
for t in range(T):
c = c0 + (c1 - c0) * t / T
for _ in range(int(1 / dt)):
x += ((x * (1 - x / 10) - c * x * x / (x * x + 1)) * dt
+ sigma * math.sqrt(dt) * rng.gauss(0, 1))
xs.append(x); cs.append(c)
return xs, cs
def detrend(y, bw=25):
out = []
for i in range(len(y)):
js = range(max(0, i - 3 * bw), min(len(y), i + 3 * bw + 1))
w = [math.exp(-0.5 * ((j - i) / bw) ** 2) for j in js]
out.append(y[i] - sum(a * y[j] for a, j in zip(w, js)) / sum(w))
return out
def stats(w): # lag-1 autocorrelation, variance
m = sum(w) / len(w); v = sum((u - m) ** 2 for u in w)
return sum((w[i] - m) * (w[i + 1] - m) for i in range(len(w) - 1)) / v, v / len(w)
def kendall(y):
s = sum((y[j] > y[i]) - (y[j] < y[i]) for i in range(len(y)) for j in range(i + 1, len(y)))
return s / (len(y) * (len(y) - 1) / 2)
def indicators(pre, step=5):
res = detrend(pre); win = len(res) // 2
idx = list(range(win, len(res) + 1, step))
a, v = zip(*(stats(res[i - win:i]) for i in idx))
return idx, a, v
if __name__ == "__main__":
xs, cs = simulate(11)
shift = next(t for t, x in enumerate(xs) if x < 3.0)
print(f"shift (x < 3) at t = {shift}, c = {cs[shift]:.3f}; deterministic fold at c = 2.604")
idx, a, v = indicators(xs[:shift - 10]) # drop the last 10 points: the fall itself
print(" window ends at t c AC(1) variance")
for k in sorted({round(i * (len(idx) - 1) / 6) for i in range(7)}):
print(f" {idx[k]:4d} {cs[idx[k] - 1]:.2f} {a[k]:.3f} {v[k]:.4f}")
print(f"Kendall tau: AC(1) {kendall(a):+.2f}, variance {kendall(v):+.2f}")
# shift (x < 3) at t = 913, c = 2.643; deterministic fold at c = 2.604
# window ends at t c AC(1) variance
# 451 1.81 0.479 0.0179
# 526 1.94 0.469 0.0185
# 601 2.08 0.513 0.0197
# 676 2.21 0.556 0.0212
# 751 2.35 0.571 0.0225
# 826 2.48 0.596 0.0232
# 901 2.62 0.737 0.0352
# Kendall tau: AC(1) +0.91, variance +0.92
自相关从 0.48 升到 0.74，方差从 0.018 升到 0.035，肯德尔 分别为 +0.91 与 +0.92；上升在后半段明显加速，正是平方根定律的样子。最后一个窗口结束时 ，已越过确定性的折叠点，系统却还在上方：折叠点留下的「瓶颈」使系统要过一段时间才离开（第 13 本第四章）。这是一个理想化的单次演示：系统确实在缓慢接近门槛，参数已知，噪声恒定，记录很长。下一章系统地评估在什么条件下这种预警可信。
证据
证据
综述与古气候。 谢弗等人 2009 年的综述把恢复变慢、自相关与方差上升、偏度变化与闪烁归纳为可能普遍适用于生态系统、气候与金融市场的「早期预警信号」(Scheffer et al. 2009)；卡彭特与布罗克在湖泊富营养化模型中说明了方差上升可作为前导指标 (Carpenter and Brock 2006)。达科斯等人分析了八次重建的突然气候转变（约 3,400 万年前温室地球的结束、新仙女木期的结束、博林—阿勒罗德转变、约 5,000 年前北非的沙漠化与四次冰期的结束），八次之前自相关都有所上升，替代数据检验的联合概率 ；但作者也说明，自助法检验显示约一半的正趋势可能来自偶然，三条最强的记录各有 450 个以上的点，其余不到 150 个点，点数少的结果可靠性较低，去趋势方式会影响结果，而且慢化只在系统逐渐接近门槛时才出现 (Dakos et al. 2008)。
证据
实验。 卡彭特等人在美国威斯康星州的彼得湖于 2008—2009 年分批投放大口黑鲈，逐步改变食物网，邻近的保罗湖作为参照；在食物网完成转变之前一年多，彼得湖叶绿素的方差、自相关等指标就已异常，参照湖没有。这是有对照的整体生态系统实验，但只有一个被操纵的湖，没有真正的重复，确切的非线性机制也不清楚 (Carpenter et al. 2011)。戴等人用以蔗糖为食的出芽酵母做实验：合作分解蔗糖造成强烈的阿利效应，每日稀释倍数相当于死亡率，在 500 到约 1,600 倍之间系统有存活与灭绝两个稳态。接近折叠点时，变异系数增加四倍以上，自相关时间从不到 2 天增加到一周以上；偏度没有系统变化 (Dai et al. 2012)。
这些正面证据有一个共同点：要么系统是被人为地、缓慢地推向门槛的（彼得湖、酵母），要么分析是在转变已经发生之后回头做的（古气候）。后一种情形有一个容易被忽略的统计问题，下一章讨论。
本章小结
接近折叠点时，恢复速率按平方根趋于零；在放牧模型中，用远离门槛处的恢复时间可以较准确地外推门槛的位置。在随机扰动下，恢复变慢表现为一阶自相关 与方差 的上升；只有方差上升可能只是噪声变大。实验湖泊、酵母与古气候记录中观察到了类似的信号。但离门槛越近，测量本身越危险，估计所需的数据也越多。
参考文献
Fan, Y. (2026). Phase Transitions, Critical Slowing Down and Resilience in Complex Dynamical Systems: From Saddle-Node Bifurcations and Early-Warning Signals to Cascading Failures, Read Beside Classical Chinese Thought on Incipient Change [复杂系统相变、临界慢化与生态韧性的非线性动力学] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23770082 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).