来源：网络

编辑：好困 yaxin

【新智元导读】45亿年，从行星撞地球到生命的诞生，地球经过了怎么样的蜕变？一位小哥成功编写程序在GPU上实现了高精度的地球模拟，4分钟内演示了45亿年的历史，足以让人震撼！

45亿年前，地球是什么样子？

一块熔岩？

看过这张模拟图，是不是感觉自己亲眼目睹了45亿年前的地球。

一位小哥编写程序在GPU上实现了高精度的地球模拟，4分钟内演示了45亿年的历史，足以让人震撼！

项目地址: https://www.shadertoy.com/view/XttcWn

代码地址: https://github.com/DokimiCU/mg_tectonic

他表示，自己完全通过GLSL最终着色器编写，模拟刷新时间为60帧/秒。

别急，还有很多。

一览45亿年前地球的演进，原行星——构造板块——水力侵蚀——全球气候——生命诞生。

从原行星开始

原行星是在原行星盘内大小如同月球尺度的胚胎行星。

其中，原行星盘是在新形成的年轻恒星（如金牛T星）外围绕的浓密气体，因为气体会从盘的内侧落入恒星的表面，所以可以视为是一个吸积盘。

根据太阳星云形成的理论，原行星在轨道轻微的扰动下和因此导致的巨大撞击与碰撞下逐渐形成真正的行星。

早期的地球是一颗原行星，温度很高，被小行星撞击的痕迹很重。

由于模拟完全是程序化生成的，没有预先渲染的纹理，所以第一个任务是生成一张地形图。

为了计算出某一经纬度（lat ）和经度（lon）上的地形高度（hight），首先要将其转换为三维直角坐标：

vec3 p = 1.5 * vec3(  sin(lon*PI/180.) * cos(lat*PI/180.),  sin(lat*PI/180.),  cos(lon*PI/180.) * cos(lat*PI/180.));

由于撞击的小行星有各种不同的大小，因此产生的陨石坑也是如此。

为此着色器在五个级别的细节上进行迭代，将大小不一的陨石坑分层叠加。

为了使陨石坑有一个真实的凹凸不平的外观，它与一些分数布朗运动（fractional Brownian motion）的噪声混合。

同时，按比例对陨石坑的影响力进行调整，从而使越大的陨石坑对地形的影响也越大。

float height = 0.;for (float i = 0.; i < 5.; i++) {   float c = craters(0.4 * pow(2.2, i) * p);   float noise = 0.4 * exp(-3. * c) * FBM(10. * p);   float w = clamp(3. * pow(0.4, i), 0., 1.);   height += w * (c + noise);}height = pow(height, 3.);

陨石坑本身是在一个三维网格上生成的，从这个网格上刻出一个球体作为表面地形。

为了避免规律性的生成，陨石坑中心使用哈希函数（hash function）从网格点中得到一个伪随机的集合。

为了计算一个给定位置的陨石坑的影响，对属于附近网格点的陨石坑进行加权平均，其权重随着与中心的距离呈指数下降。

环形山的边缘是由一个简单的正弦曲线生成的：

float craters(vec3 x) {  vec3 p = floor(x);  vec3 f = fract(x);  float va = 0.;  float wt = 0.;  for (int i = -2; i <= 2; i++)    for (int j = -2; j <= 2; j++)      for (int k = -2; k <= 2; k++) {        vec3 g = vec3(i,j,k);        vec3 o = 0.8 * hash33(p + g);        float d = distance(f - g, o);        float w = exp(-4. * d);        va += w * sin(2.*PI * sqrt(d));        wt += w;      }  return abs(va / wt);}

生成的高度图

虽然看起来有些粗糙，但是在加入一定量的水之后，生成的地形与科学家公认的早期地球面貌相似。

构造板块

中学时期，地理课上的「板块构造说」让我们了解到地球大地构造运动和海陆分布规律是这样来的。

山脉、海沟和大陆地貌的形成需要一个地壳构造活动的模型。

模型将会随机生成板块的种子位置，并产生一个初始速度。

这些板块的大小随着时间的推移而增长，模型会随机选择相邻未分配的点，并将它们添加到板块中。

其中，板块内的所有像素都储存了该板块的运动速度。

板块将以一个像素为单位沿着水平或垂直方向来移动，也就是说以离散的时间步长进行。

每个板块的时间都是随机的，这样就能让其保持在设定的速度和方向上，并且使相邻的板块不会同时移动。

当一个板块的边界像素移动到另一个板块的像素所占据的位置时，就会发生板块碰撞。

这种情况就会产生隐没带（subduction zone），其建模方法是简单地增加碰撞地点的地形高度。

虽然这只发生在沿板块边界的像素上，但通过一个简单的热侵蚀模型，这种影响会逐渐扩散到邻近的像素上，并由此模拟山脉的形成。

水力侵蚀

自然界地形的崎岖不平在很大程度上是因河流而形成的，它们以一种分支模式侵蚀着地貌。

然而对于整个星球来说，地形图的分辨率是相当低的，因此，模型必须能够模拟宽度不超过一个像素的河流。

于是，程序会检查每个像素周围的8个像素，从而确定哪个方向有最大的海拔下降。

这个最大坡度的方向就是水从这个像素流出来的地方。

水最初通过降雨分布在各个单元中，然后在每个时间步长的相邻像素之间进行传输。

侵蚀将依据水流幂律（stream power law）进行：

elevation -= 0.05 * pow(water, 0.8) * pow(slope, 2.);

其中包含当前单元的海拔（elevation）和水量（water），以及水流动方向的坡度（slope）。

海拔的下降是有上限的，这样它就不会低于水的流向。

水流和侵蚀之间的相互作用导致了地形中河流盆地的自然形成。

通过给相连的水道着色（颜色由河口的位置决定），可以让人联想到真实的河流流域地图。

美国的河流流域

全球气候

模拟整个星球的气候系统是一项艰巨的任务，但事实证明，它可以由一个程序化生成的平均海平面压力（MSLP）地图进行近似。

根据「Geoff’s Climate Cookbook」，创建MSLP地图主要是依据地貌在海洋中的位置，以及纬度的影响。

从一个真实的地球MSLP地图上获取数据，并根据陆地或海洋划分位置。

通过绘制其MSLP与纬度的关系图，可以发现陆地和海洋有两条形状略有不同的正弦曲线。

因此，通过适当地调整参数，作者得出了一个粗略的年平均气压模型。

if (land) {  mslp = 1012.5 - 6. * cos(lat*PI/45.);} else { // ocean  mslp = 1014.5 - 20. * cos(lat*PI/30.);}

当然，这还不足以生成一个真实的MSLP地图，因为分别生成陆地和海洋的数值会导致边界出现明显的不连续。

在现实中，由于气体压力的局部扩散，MSLP会在海洋到陆地的过渡中产生平稳的变化。

这种扩散过程可以通过简单地在MSLP地图上应用高斯模糊（标准差为10-15度）而得到很好的近似。

为了让气候随着季节而变化，则需要对1月和7月之间的MSLP的差异也进行建模。

其中，陆地的数据再次表明，这遵循一个正弦模式。

通过确定参数和应用高斯模糊，并与每年的MSLP地图相结合之后，可以生成全年变化的动态气候模式：

if (land) {  delta = 15. * sin(lat*PI/90.);} else { // ocean  delta = 20. * sin(lat*PI/35.) * abs(lat)/90.;}

得到MSLP之后，就可以生成风向和温度了。

这需要更多的计算来产生现实的数值，其中，季节（season）在-1和1之间波动。

float temp = 40. * tanh(2.2 * exp(-0.5 * pow((lat + 5. * season)/30., 2.)))              - 15. - (mslp - 1012.) / 1.8 + 1.5 * land - 4. * elevation;

风会从高压向低压移动，但在全球范围内，还需要考虑科里奥利力（Coriolis force）。

它会使风在压力区周围循环，其中梯度是MSLP梯度向量。

vec2 coriolis = 15. * sin(lat*PI/180.) * vec2(-grad.y, grad.x); vec2 velocity = coriolis - grad;

尽管这是一个相对粗糙的模拟，但它产生了非常真实的风循环模式，例如在印度发生的风向逆转。

降水则可以利用水蒸气从海洋到陆地的风矢量场漂移来模拟。

生命诞生

气候影响着一个星球上生命的分布。

降雨模式和温度变化决定了植物的生长速度。

随着季节的变化，食草动物会迁移到有足够多植被的地区，以维持它们的生存。

而当食草动物在迁徙时，掠食者也会跟着它们。

这些动态都可以用Lotka-Volterra扩散模型来捕捉。

float dx = plant_growth - c.y; float dy = reproduction * c.x - predation * c.z - 1.; float dz = predation * c.y - 1.; float dt = 0.1; c.xyz += dt * c.xyz * vec3(dx, dy, dz);

其中，xyz元素分别代表植被、食草动物和食肉动物的种群。

在更大的尺度上，动物种群的动态会发生十分有趣的变化。