径直切入正题吧，此前于7月27日，我设想了关于模拟地球生物量的数学模型，并尝试将其与行星的物理性质简单挂钩。 基于该随笔，尽管我并不得知诸多参数的实际数值范围，但我试图使用模糊处理的方式，将一系列参数看作整体，通过调整数值来逼近现实中可能的生物量增长曲线。 不过值得注意的是，截至撰写本篇随笔时，我暂时还没有确认任何已经考证的生物量增长曲线，所以这种逼近是带有主观性的，但也值得进行一个先验。

依照Python我简单搭建了一段代码，来模拟可能的生物量增长情况：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import truncnorm
K = 6.875 * 10**26
initial_lambda = 1.5 * 10**(-4) 
mu = 0 
sigma = 9 * 10**(-6) / 16 
lower_bound = -7.5 * 10**(4)
upper_bound = 7.5 * 10**(-4) 
a, b = (lower_bound - mu) / sigma, (upper_bound - mu) / sigma
def f(t, lambda_val):
    return (1 / K + np.exp(-lambda_val * t)) ** (-1)
t_values = np.linspace(0, 380000, 380000) 
delta_t = t_values[1] - t_values[0] 
lambda_values = [initial_lambda]
for i in range(1, len(t_values)):
    if i % 100 == 0:
        mu_sample = truncnorm.rvs(a, b, loc=mu, scale=sigma)
        new_lambda = lambda_values[-1] + mu_sample
    else:
        new_lambda = lambda_values[-1]
    lambda_values.append(new_lambda)
f_values = []
extincted = False
for t, lambda_val in zip(t_values, lambda_values):
    f_val = f(t, lambda_val)
    if f_val < 1 and not extincted:
        f_values.append(0)
        extincted = True
    else:
        f_values.append(f_val)
plt.figure(figsize=(10, 6))
plt.plot(t_values, f_values, label=f'μ ~ N({mu}, {sigma**2:.2e})')
plt.xlabel('t')
plt.ylabel('M(t)')
plt.title('Biomass M(t)')
if extincted:
    plt.text(0.5, 0.5, 'extincted', transform=plt.gca().transAxes,
             fontsize=20, color='red', ha='center', va='center')
plt.legend()
plt.grid(True)
f_max = max(f_values)
plt.xlim(250000, 400000)
yticks = np.linspace(0, f_max, 5)
plt.yticks(yticks, [f'{y:.2e}' for y in yticks])
plt.show()
                

这段代码旨在建立一个模糊模型，对应公式：\[M(t)=\frac{1}{K^{-1}+e^{-(\lambda + \sum \mu) t}}\]

其中环境容纳量\(K\)粗略地定在\(6.875\times 10^{26} \space bm\)，这是使用地球人口数粗略定下的。依照 Bar-On (2018), 认为目前的生物量为\(550 \space Gt\)，对应同时期人口约80亿，换算到人口极限100亿之下，得到粗略的生物量上限并换算至\(bm\)下得到\(K\)的数值；而每一个单位时间\(t\)对应1万地球年。

此外统合non-Markov Chain的整体效果，将其纳入参数\(\lambda+\mu\)的变动之中，基础速率\(\lambda=1.5 \times 10^{-4}\)， 而漂变参数\(\mu \sim N(0,(\frac{9}{16}\times 10^{-6})^2)\)，且该正态分布被\(\pm 7.5 \times 10^{-4}\)所截断，模拟随机的突变+环境影响，并且限制了影响程度。 \(\mu\)值按照\(\Delta t=100\)的间隔进行抽取，步长对应于每100万年产生一次影响整体的重大突变。

此外追加判定，若\(M(t)\)出现小于\(1\space bm\)的情况，则群落演化被判定为 Extincted (灭绝)。

基于以上设定，在GPT-4o的辅助下，经过几次修改得到之前的代码并成功运行 (这不得不提我这个编程白痴了，没有系统学习过代码知识，全靠询问AI起手)，得到两批各100张模拟图像。 所幸\(\lambda\)和\(\mu\)的数值在事先计算后的投放表现非常不错，模拟曲线展现了相对合理的波动，并且没有出现被判断灭绝的情况。不过，这极有可能只是“较为合理”的众多数组之一， 其他的组合也能达到甚至超过本批次的模拟结果。

对于本批次图像，大致分为五类：J型、类J型、S型、类S型和不规则型，但尚不具备严格的量化，曲线的平滑程度由人为划分认定，我从分类结果中抽取一些经典的图像作为展示： (值得注意的是，这几张图像由于y轴代码编写的疏忽，\(M(t)\)误写为\(f(t)\)。)

Fig.1   J型模拟曲线示意图，几乎是在模拟最末期才出现生命大爆发，说明前期的增长状态并不顺利。 这可能对应于一个环境剧变的地球规模的行星，生命在中早期经历了抑制，而在生命投放后约37亿年时，环境改良加上进化适应，产生了生命超级爆发的现象 (这种爆发规模比地球上的生命大爆发更加庞大)。
Fig.2   类J型模拟曲线示意图，在投放第35亿年开始出现可观的波动，并且出现几个较大的峰谷交替，最后达到较高的水平。 这种演化方式就可能与地球最为相似，每个高峰对应生命大爆发，而低谷对应大灭绝。
Fig.3   S型模拟曲线示意图，在投放的第27亿年出现大规模增长，领先前两幅图约10亿年，可能对应于一个比较舒适的行星环境。 在历经大约2亿年的持续增长后，生物量稳定在了环境容纳量附近，并且长期没有衰退现象，这是非常符合理想情况的一种随机组合。
Fig.4   类S型模拟曲线示意图，在投放的第31亿年开始大幅波动，虽然也是峰谷交替，但整体有如\(y=x+sinx\)那样的上升趋势， 最终在第37亿年时达到环境容纳量并保持平稳。这类图像所代表的行星性质可能与地球类同，但环境容纳量相对较小，或者前期的宜居度更高。
Fig.4   不规则模拟曲线示意图，在第35亿年前可以看作类J型增长，但生物量并未保持在一个相对稳定的水平，而是出现了快速的衰减。 此后则是频繁的峰谷交替，这可能代表一个环境持续剧烈变化的行星，但每次变化带来的灭绝效果又不足以灭绝所有物种。

以上五类随机图像的比例为：

J型：\(\frac{6}{200}=3\%\)

S型：\(\frac{21}{200}=10.5\%\)

类J型：\(\frac{24}{200}=12\%\)

类S型：\(\frac{41}{200}=20.5\%\)

不规则型：\(\frac{108}{200}=54\%\)

不规则的情况明显居多，可能它才是宇宙中生命演化最普遍的一种模式，只是由于人择原理 (Cirkovic, 2002)，我们“恰好”处于一个类J型的演化模型中。 但我在编写程序是遗漏了一个重要的环节，即后台备份每一次图像生成过程中\(\mu\)的随机数组，尽管目前我们仍然未知这些数值大小对应的生物物理含义，但它们仍然对演化事件的对应具有极其重要的启示作用。

还有一处未知点，即突变影响和环境影响的相对比例，此前谈到过\(\mu\)的取值是环境+突变的结果，而这个加法是带有随时间变化的权重的。 至少我现在知道，随着时间推演，单次突变对整体生物量的影响会是越来越小的，也就是说，环境影响的权重随时间推移而增大。 具体的数值探讨将会是一件非常困难的事情，希望之后能够早日解决。