CN102436679A

CN102436679A - 一种基于中值滤波的中分辨率遥感影像离散点dem构建方法

Info

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

Abstract

本发明属于中分辨率遥感影像离散点DEM构建方法领域，公开了一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法。它包括以下步骤：步骤1：计算卫星成像时刻潮高信息；步骤2：离散卫星成像时刻提取的矢量水边线，步骤3：合并上述步骤离散的等距潮位点，根据遥感影像分辨率构建中值滤波的参考格网，并对所有离散的等距潮位点进行标号，标明其所属格网ID；步骤4：遍历所有格网；步骤5：在中值滤波结果的基础上利用线性内插构建最终的DEM结果。本发明的方法通过中值滤波的方法有效地利用了落在格网内所有离散点的高程信息，提高了离散点高程精度。

Description

一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法

技术领域

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

背景技术

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

目前通过航空、航天遥感可获取大范围地面高程信息，已形成较多的地表高程反演方法：①Wehr A，Axelsson P，杜国庆等人分别在《Airborne laser scanning：An introduction and overview. ISPRS Journal of Photogrammetry and Remote Sensing》，《Airborne laser scanning：An introduction and overview. ISPRS Journal of Photogrammetry and Remote Sensing》以及《LIDAR技术在江苏沿海滩涂测绘中的应用》文章中利用LiDAR（Light Detection And Ranging，机载激光探测和测距系统）进行不同地区的高程信息反演，但由于LiDAR数据在浅水地区的测量精度不高，因此需要后期纠正；②Akira H，Stemphen G，Rabus B，金亚秋等人分别在《Mapping from ASTER stereo image data: DEM validation and accuracy assessment》、《Extraction of topography from side-looking satellite systems: A case study with SPOT simulation data》、《The shuttle radar topography mission-a new class of digital elevation models acquired by spaceborne radar》以及《用合成孔径雷达(SAR)单次飞行全极化图像数据进行地面数字高程(DEM)反演》等文章中，利用立体像对的航天遥感影像高程信息反演，但是数据源获取难度大，成本高；③张鹰等人在《海岸带潮滩土壤含水量遥感测量》一文中提出了一种基于潮滩含水量的高程信息反演的思想，建立潮滩高程与滩面沉积物含水量的相关关系，反演滩面沉积物含水量，近似得出潮滩高程，由于该方法需要实测卫星同步的地面潮滩土壤信息，不易实现业务化操作；④党福星在《利用多波段卫星数据进行浅海水深反演方法研究》，田庆久在《江苏近海岸水深遥感研究》中从辐射传输理论方程、半理论半经验的水深遥感模型、实测水深与像元辐射值相关模型等角度出发反演水底地形，由于卫星成像机理复杂，反演精度不高。

由于易变海岸带的潮汐过程提供了卫星成像时刻海岸带潮滩高度的相关“隐性信息”：即每日潮涨潮落揭示了遥感影像成像时刻潮滩水边线的高度信息，如能在易变海岸带潮滩高程信息反演过程中，挖掘卫星成像时刻的隐性潮滩高度信息，将极大改善潮滩高程的反演精度，增强潮滩高程反演结果的细节信息。

发明内容

发明要解决的技术问题

针对现有技术的不足，本发明提供了一种基于中值滤波的中分辨率遥感影像离散点DEM提取方法，它通过中值滤波的方法有效地利用了落在格网内所有离散点的高程信息，提高了离散点高程精度。

本发明的技术方案如下：

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

一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法，包括以下步骤：

步骤1：计算卫星成像时刻潮高信息，可以根据如下公式进行遥感影像成像时刻潮位信息的计算：

（1） H=H_low+A×(1-COS(t/T×180o)/2

（2） H=H_hi-A×(1-COS(t/T×180o)/2

式中，H为待插值卫星成像时刻潮位，H_hi为对应潮周期最高潮潮高，H_low为对应潮周期最低潮潮高，T为该潮周期涨潮历时，A为该潮周期潮差，t为待插值时刻与低潮时或高潮时的时间间隔。

所述步骤1具体实现方法如下：建立一个Access数据库用以存储各验潮站潮位信息，其存储字段为ID、年、月、日、时、分、潮高、纬度以及经度，提取遥感影像矢量水边线，并以_YYYYMMDDHHMM.shp的方式命名，其中YYYY为年，MM为月，DD为天，HH为时，MM为分，进行遥感影像成像时刻潮位高度计算时，先提取矢量水边线名称中的成像时间，并将GMT时间换算到东八区时间，按照步骤1中的公式进行计算。

步骤2：离散卫星成像时刻提取的矢量水边线，将卫星成像时刻的水边线根据其遥感影像分辨率进行等间隔离散，得到一系列等间隔的离散点，对该离散点shapefile文件新建一个字段，并将上述步骤计算所得的卫星成像时刻潮位信息写入离散点该文件字段，重复以上步骤直到所有水边线被离散为带有潮位信息的等距离散点，即等距潮位点。

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

步骤3：合并上述步骤离散的等间隔潮位点，根据遥感影像分辨率构建中值滤波的参考格网，并对所有离散的等距潮位点进行标号，标明其所属格网ID。

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

步骤4：遍历所有格网，针对各种不同情况采取下列操作：

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

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

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

所述步骤4中处理过的离散点需要写入一个新建的点shp文件中；

步骤5：在中值滤波结果的基础上利用Erdas软件的线性内插命令构建最终的DEM结果，所述步骤5中空间内插是在Erdas软件中利用反距离加权方法实现。

有益效果

相比现有技术，本发明实现了基于中值滤波的中分辨率遥感影像离散点DEM构建的流程，引入了中值滤波的方法，提高了DEM的构建精度。具体有益效果如下：

第一，本发明提高了DEM构建的精度。由于多数基于离散点的DEM构建算法利用所有离散点信息，进行DEM构建时存在一定的随机取舍；而本方法通过中值滤波的方法有效地利用了落在格网内所有离散点的高程信息，提高了离散点高程精度；

第二，本发明降低了DEM构建的时间复杂度。由于多数基于离散点的DEM构建算法利用所有离散点信息，数据量巨大，在进行DEM构建时需要进行频繁的数据空间查询，并在构建DEM时存在一定程度的随机取舍，因此时间效率地下；而本方法通过中值滤波的方法有效地利用了落在格网内所有离散点的高程信息，大大降低了离散点的数据量，提高了DEM构建的时间效率；

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

附图说明

图1为待反演DEM的水边线矢量数据图；

图2为通过等距离散后的水边线离散点提取结果图；

图3为通过中值滤波后的水边线离散点提取结果图；

图4为DEM构建结果图。

具体实施方式

下面结合附图对本发明做进一步描述。

实施例1

本实施例以江苏省中部辐射沙洲为例，由于江苏中部沿海大多属于正规性半日潮，即在一个太阴日(24小时25分)中有两次低潮和两次高潮，相邻的低潮或相邻的高潮的潮高大体相等的潮汐现象，因此可收集陈家坞验潮站每日的高低潮潮时、潮高观测数据。

采用附图1所示遥感影像水边线数据作为所有构建DEM的参考数据，共有34个水边线矢量文件，获取时间为2006年全年，所有水边线数据均从Landsat TM、SPOT HRV、CBERS-2 CCD、IRS-P6、AWiFS/LISS、 Beijing-1 CCD 以及 ERS-2 SAR传感器所得的遥感影像中获取。采用C#编程语言在Visual Studio 2008平台下实现本方法的5个步骤，矢量数据的读写操作通过开源地理数据格式转换类库GDAL 1.60以及ArcEngine实现。

步骤1：利用sql语句查询得到与卫星成像时刻最相邻的两个高低潮时刻，并计算高低潮时间差以及卫星成像时刻与卫星成像时刻最相近的后一时刻的时间差，根据步骤1中的公式进行计算，如果没有计算结果，则赋值为-9999：

//计算潮汐表中与卫星成像时刻最相近的后一时刻

DateTime pTideStationTimeNow = DateTime.Parse(pDataTable[ii].Rows[Now_Index]["年"].ToString() + "-" + pDataTable[ii].Rows[Now_Index]["月"].ToString() + "-" + pDataTable[ii].Rows[Now_Index]["日"].ToString() + " " + pDataTable[ii].Rows[Now_Index]["时"] + ":" + pDataTable[ii].Rows[Now_Index]["分"] + ":00");

//计算潮汐表中与卫星成像时刻最相近的前一时刻

DateTime pTideStationTimeLast = DateTime.Parse(pDataTable[ii].Rows[Now_Index - 1]["年"].ToString() + "-" + pDataTable[ii].Rows[Now_Index - 1]["月"].ToString() + "-" + pDataTable[ii].Rows[Now_Index - 1]["日"].ToString() + " " + pDataTable[ii].Rows[Now_Index - 1]["时"] + ":" + pDataTable[ii].Rows[Now_Index - 1]["分"] + ":00");

//计算卫星成像时刻

DateTime pSatelliteTime = DateTime.Parse(year + "-" + month + "-" + day + " " + hour + ":" + minute + ":00");

pSatelliteTime = pSatelliteTime.AddHours(8);

//计算高低潮时间差

System.TimeSpan NM = pTideStationTimeNow - pTideStationTimeLast;

double T = Math.Abs(NM.Minutes) + Math.Abs(NM.Hours) * 60;

//计算卫星成像时刻与卫星成像时刻最相近的后一时刻的时间差

NM = pTideStationTimeNow - pSatelliteTime;

double t = Math.Abs(NM.Minutes) + Math.Abs(NM.Hours) * 60;

double TempTide;

//根据公式计算卫星成像时刻的潮高

if (t == 0)

TempTide = Now_Tide;

else

{

if (Last_Tide > Now_Tide)

TempTide = Last_Tide - (Last_Tide - Now_Tide) * (1 - Math.Cos((T - t) / T * 3.1416)) / 2;

else

TempTide = Now_Tide - (Now_Tide - Last_Tide) * (1 - Math.Cos(t / T * 3.1416)) / 2;

}

Tide.Add(TempTide);

}

else

Tide.Add(-9999);

}

步骤2：首先，利用IFeatureCursor进行各矢量线段的遍历，针对每个矢量线段新建长度为节点数的IPoint型数组用以存放每条矢量线段的节点，再新建一个长度为节点数-1的ILine型数组用于存放该矢量线段的各条直线段，根据最小划分阈值遍历直线段，划分线元，并将得到的离散点插入新的shapefile文件中，所得结果如图2所示，其实现的核心代码为：

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

//利用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;

步骤3：

将步骤2生成的所有离散点文件利用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的格网；

步骤4，结果如图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);

}

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

Claims

1.一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法，其特征在于，它包括以下步骤：

（1） H=H_low+A×(1-COS(t/T×180o)/2；

（2） H=H_hi-A×(1-COS(t/T×180o)/2；

式中，H为待插值卫星成像时刻潮位，H_hi为对应潮周期最高潮潮高，H_low为对应潮周期最低潮潮高，T为该潮周期涨潮历时，A为该潮周期潮差，t为待插值时刻与低潮时或高潮时的时间间隔；

步骤2：离散卫星成像时刻提取的矢量水边线，将卫星成像时刻的水边线根据其遥感影像分辨率进行等间隔离散，得到一系列等间隔的离散点，对该离散点shapefile文件新建一个字段，并将上述步骤计算所得的卫星成像时刻潮位信息写入离散点该文件字段，重复以上步骤直到所有水边线被离散为带有潮位信息的等距离散点，即等距潮位点；

步骤3：合并上述步骤离散的等距潮位点，根据遥感影像分辨率构建中值滤波的参考格网，并对所有离散的等距潮位点进行标号，标明其所属格网ID；

步骤4：遍历所有格网，针对各种不同情况采取下列操作：

步骤5：在中值滤波结果的基础上利用线性内插构建最终的DEM结果。

2.根据权利要求1所述的一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法，其特征在于，所述步骤1具体实现方法如下：建立一个Access数据库用以存储各验潮站潮位信息，其存储字段为ID、年、月、日、时、分、潮高、纬度以及经度，提取遥感影像矢量水边线，并以_YYYYMMDDHHMM.shp的方式命名，其中YYYY为年，MM为月，DD为天，HH为时，MM为分，进行遥感影像成像时刻潮位高度计算时，先提取矢量水边线名称中的成像时间，并将GMT时间换算到东八区时间，按照步骤1中的公式进行计算。

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

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

5.根据权利要求1所述的一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法，其特征在于步骤4中处理过的离散点需要写入一个新建的点shp文件中。

6.根据权利要求1所述的一种基于中值滤波的中分辨率遥感影像离散点DEM构建方法，其特征在于步骤5中空间内插是在Erdas软件中利用反距离加权方法实现。