程序水文学：河流和湖泊的动态模拟

注意：该项目的完整源代码发布在Github [ 这里 ]。该存储库还包含有关如何阅读和使用代码的详细信息。

说明：我不是地质学家，所以我根据自己的知识创建了系统。

水文学概念

• 河流和溪流迁移
• 天然瀑布
• 峡谷形成
• 土壤膨胀和洪泛区

• 浮雕影响水的流动
• 侵蚀和沉积影响地形

• 在流动的溪流中，水的行为有所不同。
• 在立式游泳池中，水的行为有所不同

简化的水文模型

注意：不要忘记将它们建模为2D系统。

流和池地图

注意：流图和池图是与高度图相同大小的数组。

水为颗粒

水文循环与地图互动

• 在地形上创建粒子
• (.. ).
• , .
• , .
• , ( ) .
• .

实作

粒子类

struct Drop{
  
  //... constructors

  int index;                         //Flat Array Index
  glm::vec2 pos;                     //2D Position
  glm::vec2 speed = glm::vec2(0.0);
  double volume = 1.0;
  double sediment = 0.0;

  //... parameters
  const double volumeFactor = 100.0; //"Water Deposition Rate"

  //Hydrological Cycle Functions
  void descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale);
  void flood(double* h, double* pool, glm::ivec2 dim);
};

下降算法

void Drop::descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale){

  glm::ivec2 ipos; 

  while(volume > minVol){

    ipos = pos; //Initial Position
    int ind = ipos.x*dim.y+ipos.y; //Flat Array Index

    //Register Position
    track[ind] = true;

    //...
  }
};

//...
  //Effective Parameter Set
  double effD = depositionRate;
  double effF = friction*(1.0-0.5*stream[ind]);
  double effR = evapRate*(1.0-0.2*stream[ind]);
//...

注意：我发现这样的参数更改效果很好。

//... nind is the next position after moving the particle
  
  //Out-Of-Bounds
  if(!glm::all(glm::greaterThanEqual(pos, glm::vec2(0))) ||
     !glm::all(glm::lessThan((glm::ivec2)pos, dim))){
       volume = 0.0;
       break;
  }

  //Slow-Down
  if(stream[nind] > 0.5 && length(acc) < 0.01)
    break;

  //Enter Pool
  if(pool[nind] > 0.0)
    break;

//...

注意：溢出条件还会重置音量，以使粒子不会移入泛洪算法。

泛洪算法

void Drop::flood(double* height, double* pool, glm::ivec2 dim){

  index = (int)pos.x*dim.y + (int)pos.y;
  double plane = height[index] + pool[index];  //Testing Plane
  double initialplane = plane;                 //Water Level

  //Flood Set
  std::vector<int> set;
  int fail = 10; //Just in case...

  //Iterate while particle still has volume
  while(volume > minVol && fail){

    set.clear();
    bool tried[dim.x*dim.y] = {false};

    //Lowest Drain
    int drain;
    bool drainfound = false;

    //Recursive Flood-Fill Function
    std::function<void(int)> fill = [&](int i){

      //Out of Bounds
      if(i/dim.y >= dim.x || i/dim.y < 0) return;
      if(i%dim.y >= dim.y || i%dim.y < 0) return;

      //Position has been tried
      if(tried[i]) return;
      tried[i] = true;

      //Wall / Boundary of the Pool
      if(plane < height[i] + pool[i]) return;

      //Drainage Point
      if(initialplane > height[i] + pool[i]){

        //No Drain yet
        if(!drainfound)
          drain = i;

        //Lower Drain
        else if( pool[drain] + height[drain] < pool[i] + height[i] )
          drain = i;

        drainfound = true;
        return; //No need to flood from here
      }

      //Part of the Pool
      set.push_back(i);
      fill(i+dim.y);    //Fill Neighbors
      fill(i-dim.y);
      fill(i+1);
      fill(i-1);
      fill(i+dim.y+1);  //Diagonals (Improves Drainage)
      fill(i-dim.y-1);
      fill(i+dim.y-1);
      fill(i-dim.y+1);
    };

    //Perform Flood
    fill(index);

    //...

注意：为简单起见，此处使用八向填充算法。将来，将有可能更有效地实施它。

    //...

    //Drainage Point
    if(drainfound){

      //Set the Particle Position
      pos = glm::vec2(drain/dim.y, drain%dim.y);

      //Set the New Waterlevel (Slowly)
      double drainage = 0.001;
      plane = (1.0-drainage)*initialplane + drainage*(height[drain] + pool[drain]);

      //Compute the New Height
      for(auto& s: set) //Iterate over Set
        pool[s] = (plane > height[s])?(plane-height[s]):0.0;

      //Remove some sediment
      sediment *= 0.1;
      break;
    }

    //...

注意：当水位因泄漏而下降时，我发现这在泄漏率低的情况下效果最佳。另外，去除部分沉积岩有助于实施。

    //...

    //Get Volume under Plane
    double tVol = 0.0;
    for(auto& s: set)
      tVol += volumeFactor*(plane - (height[s]+pool[s]));

    //We can partially fill this volume
    if(tVol <= volume && initialplane < plane){

      //Raise water level to plane height
      for(auto& s: set)
        pool[s] = plane - height[s];

      //Adjust Drop Volume
      volume -= tVol;
      tVol = 0.0;
    }

    //Plane was too high and we couldn't fill it
    else fail--;

    //Adjust Planes
    float approach = 0.5;
    initialplane = (plane > initialplane)?plane:initialplane;
    plane += approach*(volume-tVol)/(double)set.size()/volumeFactor;
  }

  //Couldn't place the volume (for some reason)- so ignore this drop.
  if(fail == 0)
    volume = 0.0;

} //End of Flood Algorithm

侵蚀膜

class World {

public:
  void generate();            //Initialize
  void erode(int cycles);     //Erode with N Particles

  //...

  double heightmap[256*256] = {0.0};
  double waterstream[256*256] = {0.0};
  double waterpool[256*256] = {0.0};

};

注意：高度图是使用Perlin噪声初始化的。

//...

//Spawn Particle
glm::vec2 newpos = glm::vec2(rand()%(int)dim.x, rand()%(int)dim.y);
Drop drop(newpos);

int spill = 5;
while(drop.volume > drop.minVol && spill != 0){

  drop.descend(heightmap, waterstream, waterpool, track, dim, scale);

  if(drop.volume > drop.minVol)
    drop.flood(heightmap, waterpool, dim);

  spill--;
}

//...

注意：颗粒很少进入池中，并在下降阶段完全蒸发之前留出一到两次以上，但为防万一，我添加了这一点。

void World::erode(int N){

  //Track the Movement of all Particles
  bool track[dim.x*dim.y] = {false};

  //Simulate N Particles
  for(int i = 0; i < N; i++){
   
    //... simulate individual particle

  }

  //Update Path
  double lrate = 0.01;  //Adaptation Rate
  for(int i = 0; i < dim.x*dim.y; i++)
    waterstream[i] = (1.0-lrate)*waterstream[i] + lrate*((track[i])?1.0:0.0);

}

树木

//... descend function
double effD = depositionRate*max(0.0, 1.0-treedensity[ind]);
//...

注意：可在文件“ vegetation.h”和函数“ World :: grow（）”中找到树模型。

其他详情

结果

比较缩小流量的效果

• 基于粒子的侵蚀（基本侵蚀）接收流量和池图。池仍会影响生成
• 基本侵蚀，其参数由水文地图更改（结合侵蚀）
• 结合侵蚀与水文地图更改的参数并影响树木的侵蚀

池影响

注意：我得到了几粒种子，其中三个湖泊依次流入，但是我不想花很多时间为本文找到合适的种子。我已经生成了太多的图片和视频。

仿真速度

注意：在本文发表后不久，我对计算表面法线的方法进行了小改动，从而提高了仿真速度。在仿真过程中，这种效果非常明显，但是由于启动和泛洪所需的时间差异很大，因此很难对其进行基准测试。据我估计，仿真速度提高了一倍。

更漂亮的视频

结论与未来工作