洪涝有源淹没算法及淹没结果分析

约定不等于承诺〃 2022-03-30 05:39 172阅读 0赞

洪涝模拟仿真的实现方法主要有两种：一种是基于水动力学的洪水演进模型；还有一种是基于DEM的洪水淹没分析。详细分析例如以下：

![20141119210749844][]![20141119210753057][]

我是GIS从业者，从我们的专业角度出发，选择基于DEM的洪水淹没分析来做洪涝的模拟仿真。而基于DEM的洪水淹没分析方法主要分为有源淹没和无源淹没。

![20141119211933467][]

本篇博客採用有源淹没算法实现洪涝的模拟，算法为八领域种子扩散算法。

採用C\#版本号GDAL编写了FloodSimulation类，以下给出所有源码：

class FloodSimulation
        {
            #region 类成员变量
    
            //点结构体
            public struct Point
            {
                public int X;          //行号
                public int Y;          //列号
                public int Elevation;  //像素值(高程值)
                public bool IsFlooded; //淹没标记
    
            };
            private bool[,] IsFlood;                //淹没区域标记二维数组，用于标记淹没栅格
            private List<Point> m_FloodBufferList;  //淹没缓冲区堆栈
          
            public Dataset m_DEMDataSet;            //DEM数据集
            public Dataset m_FloodSimulatedDataSet; //洪涝淹没范围数据集
            public int m_XSize;                     //数据X方向栅格个数
            public int m_YSize;                     //数据Y方向栅格个数
            public OSGeo.GDAL.Driver driver;        //影像格式驱动
            public int[] m_FloodBuffer;            //填充缓冲区（洪涝淹没范围）
            public int[] m_DEMdataBuffer;          //DEM数据（存储高程值） 
    
            public double m_AreaFlooded;            //水面面积
            public double m_WaterVolume;            //淹没水体体积
            /* 这里的GeoTransform(影像坐标变换參数)的定义是：通过像素所在的行列值得到其左上角点空间坐标的运算參数
                比如:某图像上(P,L)点左上角的实际空间坐标为：
                Xp = GeoTransform[0] + P * GeoTransform[1] + L * GeoTransform[2];
                Yp = GeoTransform[3] + P * GeoTransform[4] + L * GeoTransform[5];                                                                     */
            public double[] m_adfGeoTransform;   
    
            #endregion
            
            //构造函数
            public FloodSimulation()
            {
                m_adfGeoTransform = new double[6];
                m_FloodBufferList = new List<Point>();
                
            }
    
            /// <summary>
            /// 载入淹没区DEM,并创建淹没范围影像
            /// </summary>
            /// <param name="m_DEMFilePath">DEM文件路径</param>
            /// <returns></returns>
            public void loadDataSet(string m_DEMFilePath)
            {
                //读取DEM数据
                m_DEMDataSet = Gdal.Open(m_DEMFilePath, Access.GA_ReadOnly);
                m_XSize = m_DEMDataSet.RasterXSize;
                m_YSize = m_DEMDataSet.RasterYSize;
                
                //获取影像坐标转换參数
                m_DEMDataSet.GetGeoTransform(m_adfGeoTransform); 
                //读取DEM数据到内存中
                Band m_DEMBand = m_DEMDataSet.GetRasterBand(1); //获取第一个波段
                m_DEMdataBuffer = new int [m_XSize * m_YSize];
                m_DEMBand.ReadRaster(0, 0, m_XSize, m_YSize, m_DEMdataBuffer, m_XSize, m_YSize, 0, 0);
    
                //淹没范围填充缓冲区
                m_FloodBuffer = new int[m_XSize * m_YSize];
                IsFlood=new bool[m_XSize,m_YSize];
    
                //初始化二维淹没区bool数组
                for (int i = 0; i < m_XSize; i++)
                {
                    for (int j = 0; j < m_YSize; j++)
                    {
                        IsFlood[i, j] = false;
                    }
                }
    
                //创建洪涝淹没范围影像
                string m_FloodImagePath = System.IO.Path.GetDirectoryName(System.Windows.Forms.Application.ExecutablePath) + "\\FloodSimulation\\FloodedRegion.tif";
                if (System.IO.File.Exists(m_FloodImagePath))
                {
                    System.IO.File.Delete(m_FloodImagePath);
                }
    
                //在GDAL中创建影像,先须要明白待创建影像的格式,并获取到该影像格式的驱动
                driver = Gdal.GetDriverByName("GTiff");
                //调用Creat函数创建影像
                // m_FloodSimulatedDataSet = driver.CreateCopy(m_FloodImagePath, m_DEMDataSet, 1, null, null, null);
                m_FloodSimulatedDataSet = driver.Create(m_FloodImagePath, m_XSize, m_YSize, 1, DataType.GDT_Float32, null);
                //设置影像属性
                m_FloodSimulatedDataSet.SetGeoTransform(m_adfGeoTransform); //影像转换參数
                m_FloodSimulatedDataSet.SetProjection(m_DEMDataSet.GetProjectionRef()); //投影
                //m_FloodSimulatedDataSet.WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 1, null, 0, 0, 0);
                
                //输出影像
                m_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m_XSize, m_YSize, m_FloodBuffer, m_XSize, m_YSize, 0, 0);
                m_FloodSimulatedDataSet.GetRasterBand(1).FlushCache();
                m_FloodSimulatedDataSet.FlushCache();
    
            }
    
            /// <summary>
            /// 种子扩散算法淹没分析
            /// </summary>
            /// <param name="m_Lat">种子点纬度</param>
            /// <param name="m_Log">种子点经度</param>
            /// <param name="m_FloodLevel">淹没水位</param>
            public void FloodFill8Direct(double m_Lat,double m_Log,double m_FloodLevel)
            {
                //首先依据种子点经纬度获取其所在行列号
                Point pFloodSourcePoint = new Point();
                int x, y;
                geoToImageSpace(m_adfGeoTransform, m_Log, m_Lat, out x, out y);
                pFloodSourcePoint.X = x;
                pFloodSourcePoint.Y = y;
                
                //获取种子点高程值
                pFloodSourcePoint.Elevation = GetElevation(pFloodSourcePoint);
                m_FloodBufferList.Add(pFloodSourcePoint);
    
                //推断堆栈中时候还有未被淹没点，如有继续搜索。没有则淹没分析结束。

while (m\_FloodBufferList.Count!=0) \{ Point pFloodSourcePoint\_temp = m\_FloodBufferList\[0\]; int rowX = pFloodSourcePoint\_temp.X; int colmY = pFloodSourcePoint\_temp.Y; //标记可淹没,并从淹没堆栈中移出 IsFlood\[rowX, colmY\] = true; m\_FloodBuffer\[getIndex(pFloodSourcePoint\_temp)\] = 1; m\_FloodBufferList.RemoveAt(0); //向中心栅格单元的8个临近方向搜索连通域 for (int i = rowX - 1; i < rowX + 2; i++) \{ for (int j = colmY - 1; j < colmY + 2; j++) \{ //推断是否到达栅格边界 if (i<=m\_XSize&&j<=m\_YSize) \{ Point temp\_point = new Point(); temp\_point.X = i; temp\_point.Y = j; temp\_point.Elevation = GetElevation(temp\_point); //搜索能够淹没且未被标记的栅格单元 if ((temp\_point.Elevation<m\_FloodLevel||temp\_point.Elevation <= pFloodSourcePoint\_temp.Elevation) && IsFlood\[temp\_point.X, temp\_point.Y\] == false) \{ //将符合条件的栅格单元增加堆栈，标记为淹没，避免反复运算 m\_FloodBufferList.Add(temp\_point); IsFlood\[temp\_point.X, temp\_point.Y\] = true; m\_FloodBuffer\[getIndex(temp\_point)\] = 1; \} \} \} \} \} //统计淹没网格数 int count = 0; for (int i = 0; i < m\_XSize; i++) \{ for (int j = 0; j < m\_YSize; j++) \{ if (IsFlood\[i,j\]==true) \{ count++; \} \} \} \} /// <summary> /// 输出洪涝淹没范围图 /// </summary> public void OutPutFloodRegion() \{ m\_FloodSimulatedDataSet.GetRasterBand(1).WriteRaster(0, 0, m\_XSize, m\_YSize, m\_FloodBuffer, m\_XSize, m\_YSize, 0, 0); // m\_FloodSimulatedDataSet.WriteRaster(0, 0, m\_XSize, m\_YSize, m\_FloodBuffer, m\_XSize, m\_YSize, 1, null, 0, 0, 0); m\_FloodSimulatedDataSet.GetRasterBand(1).FlushCache(); m\_FloodSimulatedDataSet.FlushCache(); \} // public void OutPutFloodedInfo() // \{ // \} /// <summary> /// 获取第x行第y列栅格索引 /// </summary> /// <param name="point">栅格点</param> /// <returns>该点在DEM内存数组中的索引</returns> private int getIndex(Point point) \{ return point.Y\* m\_XSize + point.X; \} /// <summary> /// 获取高程值 /// </summary> /// <param name="m\_point">栅格点</param> /// <returns>高程值</returns> private int GetElevation(Point m\_point) \{ return m\_DEMdataBuffer\[getIndex(m\_point)\]; \} /// <summary> /// 从像素空间转换到地理空间 /// </summary> /// <param name="adfGeoTransform">影像坐标变换參数</param> /// <param name="pixel">像素所在行</param> /// <param name="line">像素所在列</param> /// <param name="x">X</param> /// <param name="y">Y</param> public void imageToGeoSpace(double\[\] m\_GeoTransform, int pixel, int line, out double X, out double Y) \{ X = m\_GeoTransform\[0\] + pixel \* m\_GeoTransform\[1\] + line \* m\_GeoTransform\[2\]; Y = m\_GeoTransform\[3\] + pixel \* m\_GeoTransform\[4\] + line \* m\_GeoTransform\[5\]; \} /// <summary> /// 从地理空间转换到像素空间 /// </summary> /// <param name="adfGeoTransform">影像坐标变化參数</param> /// <param name="x">X(经度)</param> /// <param name="y">Y（纬度）</param> /// <param name="pixel">像素所在行</param> /// <param name="line">像素所在列</param> public void geoToImageSpace(double\[\] m\_GeoTransform, double x, double y, out int pixel, out int line) \{ line = (int)((y \* m\_GeoTransform\[1\] - x \* m\_GeoTransform\[4\] + m\_GeoTransform\[0\] \* m\_GeoTransform\[4\] - m\_GeoTransform\[3\] \* m\_GeoTransform\[1\]) / (m\_GeoTransform\[5\] \* m\_GeoTransform\[1\] - m\_GeoTransform\[2\] \* m\_GeoTransform\[4\])); pixel = (int)((x - m\_GeoTransform\[0\] - line \* m\_GeoTransform\[2\]) / m\_GeoTransform\[1\]); \} \}

模拟结果在ArcGlobe中的显示效果图：

![Center][]

欢迎大家留言交流。

[20141119210749844]: https://img-blog.csdn.net/20141119210749844?
[20141119210753057]: /images/20220330/b1b5f97903464fbcb39256269e3d1457.png?
[20141119211933467]: /images/20220330/1f2fb383538d477888ff273b518c12ea.png?
[Center]: /images/20220330/e7eacd606c1048399b971c601a5fe6c2.png