CN102496185A

CN102496185A - 基于多分辨率遥感影像离散点融合的dem构建方法

Info

Publication number: CN102496185A
Application number: CN2011104178457A
Authority: CN
Inventors: 刘永学; 李满春; 程亮; 李�真; 陈焱明; 胡炜; 童礼华; 蔡文婷; 杨康; 张雯
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2011-12-14
Filing date: 2011-12-14
Publication date: 2012-06-13
Anticipated expiration: 2031-12-14
Also published as: CN102496185B

Abstract

本发明公开了基于多分辨率遥感影像离散点融合的DEM构建方法，属于DEM构建领域。其步骤为：（1）将卫星成像时刻的水边线根据其遥感影像分辨率进行等间隔离散；（2）分别合并离散的等距潮位点，根据低空间分辨率遥感影像分辨率构建参考格网；（3）遍历低空间分辨率影像中的所有格网进行中值滤波；（4）以对应格网内中值滤波后的低空间分辨率影像离散点进行填充，最后在ArcGIS中实现融合；（5）利用线性内插构建最终的DEM结果。本发明提高了现有中低分辨率遥感影像离散点DEM构建的正确率，同时也提高了DEM构建的时间空间分辨率，能够快速、准确地从多期水边线数据中构建精确的DEM，降低了基础地理信息数据库的更新成本，提高了数据更新的效率。

Description

基于多分辨率遥感影像离散点融合的DEM构建方法

技术领域

本发明涉及一种基于离散点融合的DEM构建方法，更具体地说，是一种基于多分辨率遥感影像离散点融合的DEM构建方法。

背景技术

地表高程信息是地学研究的基础数据之一，快速获取、反演区域的高程信息一直是遥感信息反演的热点。一方面，对地表形态变化快速的区域而言，分析不同时间段内的地貌高程变化信息，可得出区域地貌的演化特征与规律，可为经济社会和自然环境的协调发展提供科学的决策依据；另一方面，由于滩面泥泞、潮沟密布、潮滩冲淤多变、海况复杂，使得海岸带潮滩高程测量成为地形测量中公认的困难区域。

目前通过航空、航天遥感可获取大范围地面高程信息，已形成较多的地表高程反演方法：张祖勋等人在数字摄影测量学中提出基于立体像对的航空影像高程信息反演；Wehr A，Axelsson P，杜国庆等人分别利用LiDAR（Light Detection And Ranging，机载激光探测和测距系统）进行不同地区的高程信息反演；Akira H，Stemphen G，Rabus B，金亚秋等人则利用立体像对的航天遥感影像高程信息反演；张鹰等人在《海岸带潮滩土壤含水量遥感测量》一文中提出了一种基于潮滩含水量的高程信息反演的思想，建立潮滩高程与滩面沉积物含水量的相关关系，反演滩面沉积物含水量，近似得出潮滩高程。

在数据源分辨率方面，党福星在《利用多波段卫星数据进行浅海水深反演方法研究》，田庆久在《江苏近海岸水深遥感研究》中以中空间分辨率遥感影像为数据源，从辐射传输理论方程、半理论半经验的水深遥感模型、实测水深与像元辐射值相关模型等角度出发反演水底地形，由于卫星重访周期长，不能获取连续的水边线信息，反演精度不高。刘永学在《基于MODIS遥感数据集的潮滩高程反演》中以低空间分辨率遥感影像MODIS为数据源，重访周期短，但是由于分辨率低，不能满足基础数据更新要求。

构建DEM需要足够的潮汐分布数据，但是中空间分辨率遥感影像不能满足长序列时间覆盖的要求，低空间分辨率遥感影像虽可以实现长序列时间覆盖，但是空间分辨率受到限制。因此，急需一种能够协同利用中空间分辨率遥感影像的空间覆盖和低空间分辨率遥感影像的时间覆盖，来改善潮滩高程反演精度的DEM构建方法。

发明内容

发明要解决的技术问题

本发明的目的在于克服现有技术中单一的中空间分辨率遥感影像不能满足长序列时间覆盖的要求，单一的低空间分辨率遥感影像的空间分辨率受到限制的不足，提供基于多分辨率遥感影像离散点融合的DEM构建方法，根据低空间分辨率影像的分辨率对所有遥感影像提取的水边线离散点进行格网划分，对于低空间分辨率遥感影像，根据该格网对离散点进行分区、标记、提取中值；对于中空间分辨率遥感影像，将没有离散点数据的网格用低空间分辨率中值滤波后的数据进行填充，再对融合后离散点进行内插，以达到构建DEM的目的。

技术方案

为达到上述目的，本发明提供的技术方案为：

本发明的原理：

根据遥感影像成像时刻潮位高程信息可知这一特点，结合多期遥感影像水边线数据，通过等距离散得到多期水边线离散点，根据最低遥感影像分辨率对离散点进行格网划分，以该格网对低空间分辨率遥感影像离散点进行分区、标记、提取中值，并对没有离散点数据的中空间分辨率影像格网填充中值滤波后的低空间分辨率遥感影像水边线离散点，最后通过内插得到DEM构建结果。

本发明的基于多分辨率遥感影像离散点融合的DEM构建方法，其步骤为：

步骤1：收集验潮站每日高低潮的潮时及潮高观测数据，以此内插得到卫星成像时刻潮位信息；结合离散卫星成像时刻提取的矢量水边线，将卫星成像时刻的水边线根据其遥感影像分辨率进行等间隔离散，得到一系列等间隔的离散点，并将上述计算所得的卫星成像时刻潮位信息赋予离散点，重复以上步骤，直到所有水边线被离散为带有潮位信息的等距离散点；

步骤2：分别合并离散的等距潮位点，根据低空间分辨率遥感影像分辨率构建参考格网，并对所有离散的等距潮位点进行标号，标明其所属格网ID；

步骤3：遍历低空间分辨率影像中的所有格网，针对各种不同情况采取下列操作：

（1）如果该格网内没有被标记的离散点，则不采取任何操作；

（2）如果该格网内被标记的离散点个数为奇数，设为2n+1，则将该格网内所有离散点潮位信息从小到大排序，保留潮位信息中数，即n+1点所对应的离散点，删除其他离散点；

（3）如果该格网内被标记的离散点个数为偶数，设为2n，则将该格网内所有离散点潮位信息从小到大排序，选取潮位信息最中间两个离散点，即n和n+1点，对其地理方位、潮位信息进行平均，生成一个新离散点，删除其他离散点；

其中，上述处理过的离散点写入一个新建的点shp文件中；

步骤4：对于中空间分辨率影像，将其中没有离散点的格网进行标记，并以对应格网内中值滤波后的低空间分辨率影像离散点进行填充，最后在ArcGIS中实现融合；

步骤5：在上述步骤的基础上利用线性内插构建最终的DEM结果，该线性内插过程在Erdas软件中利用反距离加权方法予以实现。

更进一步地，步骤1中水边线离散点生成实现步骤如下：遍历所有水边线shp文件，每遍历到一个新的水边线shp文件则生成一个投影信息与原shp文件一致的点文件，再遍历该水边线shp文件所有线性矢量记录，针对每条矢量记录，找到该矢量的起始节点A0作为线段一个端点，接着寻找矢量的下一个节点A1作为线段另一个端点，计算该线段长度，如果大于给定的间隔L，则生成一个距离A0点长度为L的离散点N0，并将N0作为起始端点，重复上述步骤直至线段长度小于或等于给定间隔时，将A0及所有生成的N0，N1，……Nn点写入遍历开始时的点shp文件，判断该矢量是否有下一个节点，如果存在，则将A1置为线段起点，A2置为线段终点，重复上述步骤，直至矢量文件遍历完成。

更进一步地，步骤2中构建中值滤波的参考格网及离散点标号，具体实现步骤如下：利用ArcGIS的Merge命令将所有离散点shp文件合并成一个shp文件，并求取该点shp文件所有离散点的极值X，Y坐标，将其作为网格的边界值，分析原始遥感影像空间分辨率S，并以其作为格网边长，遍历每一个格网，将落入其中的离散点标记该格网的ID号。

有益效果

相比现有技术，本发明实现了多分辨率遥感影像离散点融合的DEM构建的流程，引入了融合的方法，与已有的公知技术相比，具有如下显著效果：

（1）本发明提高了DEM构建的精度。由于多数基于离散点的DEM构建算法利用所有离散点信息，进行DEM构建时存在一定的随机取舍，而本方明通过对低分辨率遥感影像水边线离散点采用中值滤波的方法有效地利用了落在格网内所有离散点的高程信息，改善离散点空间分布不均衡的现象，提高了离散点高程精度；

（2）本发明提高了DEM构建的时间空间分辨率。由于中空间分辨率遥感影像不能满足时间覆盖，而低空间分辨率遥感影像虽可以满足时间覆盖，但是空间分辨率受到限制，本发明同时利用中空间分辨率遥感影像的空间覆盖和低空间分辨率遥感影像的时间覆盖，提高了DEM构建的时间空间分辨率；

（3）本发明实际应用表明，通过本发明能够快速、准确地从多期水边线数据中构建精确的DEM，能够为国家各级基础地理信息数据库DEM信息提供技术支撑，降低基础地理信息数据库的更新成本，提高数据更新的效率。

附图说明

图1为实施例中通过等距离散后的低分辨率遥感影像水边线离散点提取结果图；

图2为实施例中通过等距离散后的中分辨率遥感影像水边线离散点提取结果图；

图3为实施例中通过中值滤波后的低分辨率遥感影像水边线离散点提取结果图；

图4为实施例中的中低分辨率遥感影像水边线离散点融合结果图；

图5为实施例中的DEM构建结果图。

具体实施方式

下面结合附图和实施例对本发明作进一步的描述。

实施例

采用C#编程语言在Visual Studio 2008平台下实现本方法的5个步骤，矢量数据的读写操作通过开源地理数据格式转换类库GDAL 1.60以及ArcEngine实现。

步骤1：收集陈家坞验潮站每日高低潮的潮时及潮高观测数据，以此内插得到卫星成像时刻潮位信息；结合离散卫星成像时刻提取的矢量水边线，将卫星成像时刻的水边线根据其遥感影像分辨率进行等间隔离散，得到一系列等间隔的离散点，并将上述计算所得的卫星成像时刻潮位信息赋予离散点，重复以上步骤，直到所有水边线被离散为带有潮位信息的等距离散点。

本实施例中，采用附图1及附图2所示遥感影像水边线离散点数据作为构建DEM的参考数据，获取所有水边线数据均从Landsat TM、SPOT HRV、CBERS-2 CCD、IRS-P6、AWiFS/LISS、 Beijing-1 CCD 以及 ERS-2 SAR传感器所得的遥感影像中获取，时间为2006年全年。首先，利用IFeatureCursor进行各矢量线段的遍历，针对每个矢量线段新建长度为节点数的IPoint型数组用以存放每条矢量线段的节点，再新建一个长度为节点数-1的ILine型数组用于存放该矢量线段的各条直线段，根据最小划分阈值遍历直线段，划分线元，并将得到的离散点插入新的shapefile文件中，其实现的核心代码为：

/***********得到水边线和其各个点坐标 ***************/

//利用as实现多态IGeometry即shape向Ipolyline的转换

IPolyline pWaterline = pFeature.Shape as IPolyline;

//在Ipolyline工厂里利用IPolyline对象去访问IPointCollection4对象

IPointCollection4 pWaterline_points = pWaterline as IPointCollection4;

//建立一个存放点的数组

IPoint[] points = new IPoint[pWaterline_points.PointCount];

IGeometryBridge geometryBridge = new GeometryEnvironmentClass();

for (int i = 0; i < points.Length; i++)

{

points[i] = new PointClass();

｝

//将pointcollection中的点读到points数组中去

geometryBridge.QueryPoints(pWaterline_points, 0, ref points);

/********** 完毕 ***************/

/********* 得到每条水边线的直线段 **************/

//为每一条水边线段建立一个存放直线段的数组,并记录相邻两点的直线段

ILine[] Waterline_line_array = new ILine[pWaterline_points.PointCount - 1];

for (int i = 0; i < (pWaterline_points.PointCount - 1); i++)

{

Waterline_line_array[i] = new LineClass();

//用水边线的相邻两点构建直线段

Waterline_line_array[i].PutCoords(points[i], points[i + 1]);

}

//此处才是真正将第一个点写入shp中

if(Waterline_line_array[0].Length<2000)

{

OBJ = point_featureCursor.InsertFeature(point_featurebuffer);

point_featureCursor.Flush();

}

//用水边线的各个点来构建直线段

for (int j = 0; j < points.Length - 1; j++)

{

if (Waterline_line_array[j].Length < 2000)

{

//记录已经被划分的直线段的长度

double Waterline_line_length = 0;

//记录划分直线段线元的编号

int t = 0;

//遍历直线段，划分线元

while (Waterline_line_length < Waterline_line_array[j].Length)

{

//如果剩下的直线段已经不够再划分这个线元，那就用剩下的直线段来作为这个线元

if ((Waterline_line_array[j].Length - Waterline_line_length) < EL)

{

ILine pline = new LineClass();

IPoint from_point = new PointClass();

IPoint to_point = new PointClass();

//注意：每个线元的始末点划分方法，即用上一个线元划分的终点作为起点，再加上该线元的长度EL

Waterline_line_array[j].QueryPoint(esriSegmentExtension.esriNoExtension, Waterline_line_length, false, from_point);

to_point = Waterline_line_array[j].ToPoint;

TideDataset[index].MapToPixel(to_point.X, to_point.Y, out pcolunm, out prow);

obj = TideDataset[index].GetPixelValue(0, pcolunm, prow);

point_featurebuffer.set_Value(point_featurebuffer.Fields.FindFieldByAliasName(TideFieldName), obj);

point_featurebuffer.Shape = to_point;

OBJ = point_featureCursor.InsertFeature(point_featurebuffer);

point_featureCursor.Flush();

}

//否则，用计算出的该线元的长度划分线段

else

{

ILine pline = new LineClass();

IPoint from_point = new PointClass();

IPoint to_point = new PointClass();

Waterline_line_array[j].QueryPoint(esriSegmentExtension.esriNoExtension, Waterline_line_length + EL, false, to_point);

TideDataset[index].MapToPixel(to_point.X, to_point.Y, out pcolunm, out prow);

obj = TideDataset[index].GetPixelValue(0, pcolunm, prow);

point_featurebuffer.Shape = to_point;

OBJ = point_featureCursor.InsertFeature(point_featurebuffer);

point_featureCursor.Flush();

//用这两个点来构建线元

pline.PutCoords(from_point, to_point);

}

//递进，得到下一个直线段

t++;

Waterline_line_length += EL;

步骤2：分别合并离散的等距潮位点，根据低空间分辨率遥感影像分辨率构建参考格网，并对所有离散的等距潮位点进行标号，标明其所属格网ID。

本实施例中，将步骤1生成的低空间分辨率遥感影像水边线离散点文件和中空间分辨率遥感影像水边线离散点文件利用ArcGIS软件中Merge命令分别进行合并得到两个汇总离散点的shapefile文件，再利用ArcGIS软件中Create Fishnet功能进行格网构建，该功能中Template Extent选择汇总离散点的shapefile文件，Cell Size Width与Cell Size Height均为提取水边线的低空间分辨率遥感影像分辨率，再利用Feature to Polygon命令将上一步生成的Fishnet转换成polygon的格网，最后利用Spatial Join命令对汇总离散点进行标记，该功能中Target Features选择汇总离散点文件，Join Features选择由Fishnet转换成polygon的格网，其实现的核心代码为：

/// 确定格网内的点

private void PointsInGrids()

{

for (int i = 0; i < gridNumber; i++)

{

for (int j = 0; j < pointNumber; j++)

{

if (points[j].X > grids[i].X - gridSize / 2.0 && points[j].X < grids[i].X + gridSize / 2.0 &&

points[j].Y > grids[i].Y - gridSize / 2.0 && points[j].Y < grids[i].Y + gridSize / 2.0)

grids[i].Index.Add(j);

}

步骤3：对低空间分辨率遥感影像水边线离散点数据进行中值滤波，遍历低空间分辨率影像中的所有格网，针对各种不同情况采取下列操作：

（3）如果该格网内被标记的离散点个数为偶数，设为2n，则将该格网内所有离散点潮位信息从小到大排序，选取潮位信息最中间两个离散点，即n和n+1点，对其地理方位，潮位信息进行平均，生成一个新离散点，删除其他离散点；其中，上述处理过的离散点写入一个新建的点shp文件中；结果如图3所示。中值滤波代码如下：

/// 中值滤波

private void MedianFilter()

{

for (int i = 0; i < gridNumber; i++)

{

if (grids[i].Index.Count > 0)

{

List<PData> tempData = new List<PData>();

for (int j = 0; j < grids[i].Index.Count; j++)

tempData.Add(points[grids[i].Index[j]]);

for(int m = 0; m < tempData.Count; m++)

for (int n = m + 1; n < tempData.Count; n++)

{

if (tempData[m].TideValue > tempData[n].TideValue)

{

int tempOID = tempData[m].OID;

double tempX = tempData[m].X;

double tempY = tempData[m].Y;

double tempTideValue = tempData[m].TideValue;

tempData[m].OID = tempData[n].OID;

tempData[m].X = tempData[n].X;

tempData[m].Y = tempData[n].Y;

tempData[m].TideValue = tempData[n].TideValue;

tempData[n].OID = tempOID;

tempData[n].X = tempX;

tempData[n].Y = tempY;

tempData[n].TideValue = tempTideValue;

}

PData pData = new PData();

pData.OID = outPoints.Count;

if (tempData.Count % 2 == 0)

{

pData.X = (tempData[tempData.Count / 2 - 1].X + tempData[tempData.Count / 2].X) / 2.0;

pData.Y = (tempData[tempData.Count / 2 - 1].Y + tempData[tempData.Count / 2].Y) / 2.0;

pData.TideValue = (tempData[tempData.Count / 2 - 1].TideValue + tempData[tempData.Count / 2].TideValue) / 2.0;

}

else

{

pData.X = tempData[tempData.Count / 2].X;

pData.Y = tempData[tempData.Count / 2].Y;

pData.TideValue = tempData[tempData.Count / 2].TideValue;

}

outPoints.Add(pData);

}

步骤4：对于中空间分辨率影像，将其中没有离散点的格网进行标记，并以对应格网内中值滤波后的低空间分辨率影像离散点进行填充，最后在ArcGIS中实现融合，结果如图4所示，其实现的核心代码为：

private void SelectInsertPoints()

{

for (int i = 0; i < grids.Length; i++)

{

int flag = 0;

for (int j = 0; j < basepoints.Length; j++)

{

if (grids[i].X + gridSize / 2 > basepoints[j].X && grids[i].X - gridSize / 2 < basepoints[j].X

&& grids[i].Y + gridSize / 2 > basepoints[j].Y && grids[i].Y - gridSize / 2 < basepoints[j].Y)

{

flag = 1;

break;

}

if (flag == 0)

{

int index = SearchInsertPointsWithinGrid(i);

if (index != -1)

{

selectedinsertpointsindex.Add(index);

}

private int SearchInsertPointsWithinGrid(int gridindex)

{

int index = -1;

for (int i = 0; i < insertpoints.Length; i++)

{

if (grids[gridindex].X + gridSize / 2 > insertpoints[i].X && grids[gridindex].X - gridSize / 2 < insertpoints[i].X

&& grids[gridindex].Y + gridSize / 2 > insertpoints[i].Y && grids[gridindex].Y - gridSize / 2 < insertpoints[i].Y)

{

index = i;

break;

}

return index;

}

步骤5：在上述步骤的基础上利用线性内插构建最终的DEM结果，该线性内插过程在Erdas软件中利用反距离加权方法予以实现。

本实施例中，将步骤4中得到的离散点融合图利用Erdas软件中Data Prepration中的Create Surface功能实现。首先在对话框中点击file下的read功能，读取步骤4的离散点融合结果，点击Surface下的Surfacing，设置输出路径、输出文件名、输出结果分辨率等参数，在Surfacing Methods下选择Linear Rubber Shipping功能，点击运行，得到DEM结果，结果如图5所示。

本发明提高了现有中低分辨率遥感影像离散点DEM构建的正确率，同时也提高了DEM构建的时间空间分辨率。实际应用表明，通过本发明能够快速、准确地从多期水边线数据中构建精确的DEM，能够为国家各级基础地理信息数据库DEM信息提供技术支撑，降低基础地理信息数据库的更新成本，提高数据更新的效率。

Claims

1.基于多分辨率遥感影像离散点融合的DEM构建方法，其步骤为：

其中，上述处理过的离散点写入一个新建的点shp文件中；

2.根据权利要求1所述的基于多分辨率遥感影像离散点融合的DEM构建方法，其特征在于：步骤1中水边线离散点生成实现步骤如下：遍历所有水边线shp文件，每遍历到一个新的水边线shp文件则生成一个投影信息与原shp文件一致的点文件，再遍历该水边线shp文件所有线性矢量记录，针对每条矢量记录，找到该矢量的起始节点A0作为线段一个端点，接着寻找矢量的下一个节点A1作为线段另一个端点，计算该线段长度，如果大于给定的间隔L，则生成一个距离A0点长度为L的离散点N0，并将N0作为起始端点，重复上述步骤直至线段长度小于或等于给定间隔时，将A0及所有生成的N0，N1，……Nn点写入遍历开始时的点shp文件，判断该矢量是否有下一个节点，如果存在，则将A1置为线段起点，A2置为线段终点，重复上述步骤，直至矢量文件遍历完成。

3.根据权利要求1或2所述的基于多分辨率遥感影像离散点融合的DEM构建方法，其特征在于：步骤2中构建中值滤波的参考格网及离散点标号，具体实现步骤如下：利用ArcGIS的Merge命令将所有离散点shp文件合并成一个shp文件，并求取该点shp文件所有离散点的极值X，Y坐标，将其作为网格的边界值，分析原始遥感影像空间分辨率S，并以其作为格网边长，遍历每一个格网，将落入其中的离散点标记该格网的ID号。