SRTM与ASTER数据精度分析和比较本科毕业论文.doc

商丘师范学院2013届本科毕业论文（设计）SRTM与ASTER数据的精度分析和比较引言由于可获得的SRTM每个90m的数据点是由9个30m的数据点算术平均得来的，我们可以了解SRTM数据所记录的并不是真正的地面点数据。ASTER以地面获得的实际数据点的值为矩阵的元素，它能反映真实的地貌特征。通过选取同一个地区的不同数据形式的数字地面模型2，在Global Mapper中绘制等高线，分析等高线的特点，观察其3维模型的特点，再在Matlab中进行矩阵的运算，求的此地区两数据的差值，来分析SRTM DEM与ASTER DEM的精度，从而确定哪种数据精度更高，在什么地区运用什么种类的DEM更好、更合适。本文在介绍了SRTM数据、ASTER数据、Global Mapper、Matlab软件之后，运用Global Mapper软件进行等高线的绘制，等高线的叠加，发现两种数据的等高线并不一致，等高线的走向以及弯曲程度不一致。再观察其3维模型，基于ASTER的3维模型更细致的表达了地面的起伏变化。从而得出ASTER的精度高于SRTM的精度。利用Matlab软件进行矩阵的运算，计算差值矩阵的最大值，最小值，平均值以及方差等，分析两数据的在平原和山区的差异，得出在平原地区无论是30m的数据好似是90m的数据均能很好地表达地貌，在实际生活和普通科研中，二者均能满足需要。而在山区二者的差异很大，平均值在SRTM所描绘的地貌比ASTER所描绘的地貌平均高出6m的高程，而且其差值不稳定，上可达270m左右，下可达-160m，方差达到15m，所以两种数据所描绘的山区地貌有很大差异，在表现地面起伏形态时，在山区益选择ASTER数据。而在垂直度上需要借助其他数据进行比较。1 SRTM与ASTER 数据10及Global Mapper和Matlab简介SRTM数据与ASTER数据都是通过遥感传感器对地面上的信息进行采集，获得的地面上的地物地貌三维数据，即数字高程模型4，因此很多时候又称为SRTM DEM和ASTER DEM5。在大面积范围内，两者都近似的表现了地面的高低起伏变化，为研究地球，商业发展，国防军事，科学研究等提供了方便。但二者又有一些不同之处，下面就对其分别介绍，并简要介绍一下Global Mapper 和Matlab软件。1.1 SRTM 数据8航天地形测绘是指以人造地球卫星、宇宙飞船、航天飞机等航天器为工作平台，对地球表面所进行的遥感测量。以往的航天测绘由于其精度有限，一般只能制作中小比例尺的地图。2000年2月11日上午11时44分，美国“奋进”号航天飞机在佛罗里达州卡那维拉尔角的航天发射中心发射升空，执行耗资3.64亿美元，称之为“航天飞机雷达地形测绘使命（Shuttle Radar Topography Mission ,简称SRTM）”的空间飞行任务。“奋进”号上搭载的SRTM系统共计进行了3222小时23分钟的数据采集工作，获取北纬60至南纬56之间面积超过1.19亿平方公里9.8万亿字节的影像数据，覆盖全球陆地表面的80%以上，覆盖了中国全境。SRTM数据之大，精度之高在测绘史上前所未有，全部数据的处理约需两年的时间，最终获取平面精度为20m,高程精度为16m的全球数字高程模型DEM，可以将美军现有的全球DEM精度提高约30倍。借助如此精准的地图，美国国防部将能更好地选择部署部队，甚至发射导弹的地点，因此这项计划所得的最佳品质的图像被列为机密文件，供美国军方和情报局使用，只将一些分辨率较低的地图向公众公开，用于普通科研和商业。SRTM数据每经纬度方格提供一个文件，精度有1arcsecond和3 arcseconds两种，称作SRTM1和SRTM3或者30m和90m数据，SRTM1的文件里面包含36013601个采样点的高程数据，SRTM3的文件里面包含12011201个采样点的高程数据。目前，能够免费获取的中国境内的SRTM3文件是90m的数据，每个90m的数据点是由9个30m的数据点算术平均得来的。用户可以通过中国科学院计算机网络信息中心国际科学服务器平台免费获得（/），本文选择网址：/SELECTION/inputCoord.asp。下载的数据，其产品类型是SRTM 90m DEM version 4，文件名为srtm_59_06.zip（东经E110E115，北纬N30N35的区域），包括山区和平原，便于后续工作的进展。下载完成后解压数据，一次解压后不修改文件名，放在指定位置。再次解压得到文件srtm_59_06.hdr(头文件)，srtm_59_06.tfw（tiff影像坐标信息的文本文件），srtm_59_06.tif（二进制的数字高程模型）。1.2 ASTER数据ASTER是在2009年6月30日，美国航天局NASA与日本经济产业省（METI）共同推出的最新的地球电子地形数据ASTER GDEM(先进星载热发射和反射辐射仪全球数字高程模型)，该数据是根据NASA的新一代对地观测卫星TERA的详细观测结果制作完成的。这一全新的地球数字高程模型包含了先进星载热发射和反辐射及（ASTER）搜集的130万个立体影像。ASTER测绘数据覆盖范围为北纬83到南纬83之间的所有陆地区域，比以往任何地形图都要广得多，达到了地球陆地表面的99%。此前，最完整的地形数据是NASA的航天飞机雷达地形测绘任务(SRTM)提供的，此项任务对位于北纬60和南纬56间地球80%的陆地进行了测绘。ASTER GDEM是采用全自动化的方法对150万景的ASTER存档数据进行处理生成的，其中包括通过立体相关生成的1264118个基于独立场景的ASTER GDEN对应区域的最后像素值。纠正剩余的异常数据，再按11分片，生成全球ASTER GDEM数据。目前，ASTER GDEM数据可以在网上免费获取。本文所使用数据的下载网址为 /xiazai/5241.html，从中选择河南省的网址：http:/115.com/file/cl326i69。下载并保存数据，解压后选择处于平原地区的压缩包N34E114.zip和处于高山地区的压缩包N33E110.zip作为分析数据。再次解压得到平原区的ASTGTM2_N34E114_dem.tif（DEM文件）和ASTGTM2_N34E114_num.tif (质量评估文件)以及高山区的ASTGTM2_N33E110_dem.tif（DEM文件）和ASTGTM2_N33E110_num.tif(质量评估文件)，放在指定位置。1.3 Global Mapper简介Global Mapper 是一款地图绘制软件，不仅能够将数据（例如：SRTM数据）显示为光栅地图、高程地图、矢量地图，还可以对地图作编辑、转换、打印、记录GPS及利用数据的GIS(地理信息系统)功能。它可以转换数据集的投影方式以符合你的项目的坐标系统，并可以同时对数据集的范围进行裁剪。它还提供距离和面积计算，光栅混合、对比度调节、海拔高度查询、视线计算，以及一些高级功能，如图像校正、通过地表数据进行轮廓生成、通过地表数据观察分水岭、对3Dpoint数据转换为三角多边形和网格化等。通过内建的脚本语言或众多的批处理转换选项能够高效地完成重复性任务标题栏菜单栏工具栏图1-1 Global Mapper界面1.4 Matlab1简介Matlab是由美国MathWorks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。Matlab一年有两次更新上半年为a，下半年为b，本文所用的是Matlab 2012b。菜单栏获得帮助工作间浏览器输入指令或函数工具栏工作间浏览器在历史窗中查看或再执行以前的指令图1-2 Matlab 2012b新界面2 等高线、三维模型的生成和比较72.1等高线的绘制、叠加和比较由于下载的srtm_59_06.tif和ASTGTM2_N33E110_num.tif数据量太大，在生成等高线时速度缓慢，为清楚显示等高线细部也不需要如此大的数据，并且，正射影像图的精确显示，也要求选择部分数据。所以首先选择部分数据，本文选择北纬33.633.7，东经110.6110.7的区域。保存成.dem的格式，保存名字为srtmshangqu.dem和astershangqu.dem。绘制等高距为25m等高线，并保存成srtmshanqudenggaoxian.KMZ和astershanqudenggaoxian.KMZ。再在Global Mapper中同时打开srtmshanqudenggaoxian.KMZ和astershanqudenggaoxian.KMZ实现等高线的叠加，查看高程为750m的一条等高线的走向。图2-1 高程为750m的等高线由图2-1可以看出两种数据的等高线走向大体一致，个别地区会出现偏差,在至高点数量上和等高线长度上不相等，在细部表现上也有所不同，其中astershanqudenggaoxian.KMZ的等高线弯弯曲曲，比srtmshanqudenggaoxian.KMZ的等高线更细致平滑、更显得准确。下面观看其折射影像图，看有何特点。2.2 正射影像图的比较9在Global Mapper中打开srtmshangqu.dem和astershangqu.dem可看到其正射影像，观察其特征。图2-2 正射影像图 两种数据的三维模型有明显的不同点，astershangqu.dem比srtmshanqu.dem棱角分明，在表现地貌方面更棱角分明、更清晰、更细致。然而在垂直方向上，两种数据是否有差距，差距为多少，在山区和平原区又有什么样的特征，以下将重点介绍。3 山区与平原区在垂直方向上的比较6已知srtm_59_06文件的第12012401行，11201列数据与文件ASTGTM2_N33E110_dem相对应，且是属于山区，又因为srtm_59_06文件中的一个方格对应地面90m的距离，而ASTGTM2_N33E110_dem文件中的一个方格对应地面30m的距离，所以必须使矩阵中代表同一地面点的数据参加运算才能得到两数字高程模型的差值矩阵，否则发生错位，得不到想要的结果。以此为原则，可采用数据提取法和数据内插法3等来实现运算，本文提供了三种算法：（1）数据提取法：在ASTER数据中直接提取与SRTM数据对应的3点数据进行差值分析。（2）均值赋值法：利用ASTER数据中3及其周围的八个点求平均值，然后把均值赋给3点。再进行两种数据的差值分析。（3）移动曲面拟合赋值法：利用ASTER数据中3点周围的八个点拟合内插出3点得高程，然后进行两种数据的差值分析。在Matlab 2012b软件中载入srtm_59_06.tif、ASTGTM2_N33E110_dem.tif和ASTGTM2_N34E114_dem.tif 分别利用以上几种方法在“Command window“中编写编程，求在山区和平原区的差值矩阵，以及两差值矩阵的最大值、最小值、平均值和方差。并分析插值前后插值矩阵的特点，具体方法和步骤见附录1。结果无论是均值法还是拟合法赋值，前后的差值矩阵表明，ASTER在山区的插值后的数据精度稍微低于平原区，在平原地区均值后的数据与原始数据相差不大，差值矩阵平均值都非常小，接近与零，并且差值的最大值、最小值和方差都不大，所以，我们认为赋值前后数据非常相近，可以将赋值前后的数据等同起来看待。ASTER数据取自地面的实际数据，如果SRTM3数据也达到同样的精度，均值或拟合内插的ASTER数据也有同样的精度，则赋值前后ASTER数据的差值矩阵与赋值后两数据的差值矩阵相近，而结果却是相差很大。ASTER赋值后与SRTM3数据在平原地区的差值超过了与原始数据在山区的差值，更远远大于在平原区差值，可见ASTER数据与SRTM3数据在山区具有较大的差，所描绘的地貌形态在垂直方向上并不一致。不管用何种方法，插值后两数据在山区差值均值都高达6m之多，最大值在270m左右，最小值-160m左右，方差在15m左右，所以在山区两数据在垂直方向上的差值大，且不稳定性大，两种数据所描绘的地貌在垂直方向上并不重合。而在平原区，两数据差值矩阵的均值接近于0m，最大值15m左右、最小值-103m左右，且方差在4m和5m之间，说明在平原区，二者的相差不是很大，两种数据所描绘的地貌均能正确表达地面的形态，因此在平原地区，两种数据的精度相当，在实际生活中二者均可用。而在山区，若要得到更为精确的地貌特征，根据等高线和正射影像图的特点，宜选择ASTER数据。结语利用Global Mapper软件分析SRTM和ASTER数据的等高线和3维模型，根据等高线的走向和弯曲程度，平滑程度，分析得ASTER数据更能细致的表现地貌的特征，更能准确的表现地形的特征。利用Matlab软件对矩阵进行提值、插值、求差矩阵，对其元素进行分析，根据计算所得结果得知，无论是用提取法还是用插值法，这些数据充分证明了二者在山区的差值较大，且非常不稳定，此种结果造成，在山区SRTM数据所描绘的地形与ASTER数据所描绘的地形在垂直方向上有很大的差别，在某些地区可能给商业经济等带来损失。由于ASTER1数据是采取的真实的30m的地面点，在表现地貌、地形方面，其结果无论是在山区还是在平原都比SRTM数据（每个90m的数据点是由9个30m的数据点算术平均得来的）更为细致、精确、有优势，所以在山区宜选择ASTER数据。在平原区二者相差不大，都可以为人们应用。但在（山区）垂直方向上究竟哪种数据更贴近于实际的地面，还需要第三种数据进行比较。就目前而言，追求SRTM、 ASTER两种数据和其他数据的结合已成为发展趋势，用以获得更为准确和完备的全球地形图。参考文献1 王家文Matlab7.6图形图像处理M国防工业大学出版社，20092 张剑清摄影测量学M武汉大学出版社，20033 王佩军，徐亚明摄影测量学M武汉大学出版社，20104 汤国安，刘学军，闾国年数字高程模型及地学分析的原理与方法M科学出版社，20065 汤国安，龚建雅，陈正江数字高程模型地形描述精度量化模拟研究(a)J测绘学报，2001，4(18)：80816 胡德勇，李京，陈云浩不同模式Radarsat影像DEM提取及精度比较J遥感技术与应用，2006，6(15)：68707 闫业超，张树文，岳书平东北川岗地区SRTM数据质量评价J中国科学院研究生院学报，2008，1(28)：1041058 惠凤鸣，田庆久，李应成ASTER数据的DEM生产及精度评价J遥感信息，2004， 1(25)：1001029 夏涛，杨武年，刘汉湖利用ASTER立体像对提取相对DEM及正射影像地图制作J测绘科学，2007，3（32）：14414510 赵国松，杜耘，凌峰，李晓冬ASTER GDEM与SRTM3高程差异影响因素分析J测绘科学，2012，4(37)：156158致 谢在此次毕业论文写作过程中，非常感谢老师热情的指点和帮助，当我有什么不懂的问题，老师就一遍遍详细地讲解，有时电话解答，有时当面解决。让我在写作过程中即巩固了专业知识，也学到了很多其他的东西，懂得了科学研究就是一个严谨的创作过程，来不得半点的马虎和虚假。老师们知识渊博，积极乐观，热真负责，鼓励着我们要时刻充实自己。通过这次论文写作，让我了解到，我国在测绘方面所做出的努力和所取得的成就，以及与世界各国的差距，让我们看到世界未来的发展方向，中国未来的发展任务和趋势，我们新一代的年轻人在建设中国，创造美好世界的过程中，更应该脚踏实地的做好该做之事，能做之事。再次感谢老师！附录1“Command window”中的编程及结果如下： whos %显示文件信息 Name Size Bytes Class Attributes ASTGTM2_N33E110_dem 3601x3601 25934402 int16 ASTGTM2_N34E114_dem 3601x3601 25934402 int16 srtm_59_06 6001x6001 72024002 int16 %数据提取的方法 srtmshanqu=srtm_59_06(1201:2401,1:1201); % srtm山区数据3 astershanqu=ASTGTM2_N33E110_dem(1:3:3601,1:3:3601); % aster山区( 1)数据提取(3) srtmpingyuan=srtm_59_06(1:1201,4801:6001); %srtm平原数据3 asterpingyuan=ASTGTM2_N34E114_dem(1:3:3601,1:3:3601); % aster平原区( 1)数据提取(3) %处理过程 shanqucha1=srtmshanqu-astershanqu; %山区对应3数据点矩阵的差值矩阵 pingyuancha1=srtmpingyuan-asterpingyuan; %平原区对应3数据点矩阵的差值矩阵 sqcmax1=double(max(max(shanqucha1); %山区差值矩阵的最大值 sqcmin1=double(min(min(shanqucha1); %山区差值矩阵的最小值 pycmax1=double(max(max(pingyuancha1) ; %平原差值矩阵的最大值 pycmin1=double(min(min(pingyuancha1) ; %平原差值矩阵的最小值 sqcmean1=double(mean2(shanqucha1); %山区差值矩阵的平均值 pycmean1=double(mean2(pingyuancha1); %平原差值矩阵的平均值 sqcstd1=double(std2(shanqucha1) ; %山区差值矩阵的方差 pycstd1=double(std2(pingyuancha1) ; %平原差值矩阵的方差 jieguo1= sqcmax1 sqcmin1 sqcmean1 sqcstd1;pycmax1 pycmin1 pycmean1 pycstd1jieguo1 = 288.0000 -161.0000 -6.1272 15.5155 114.0000 -103.0000 0.1884 5.1480 clear astershanqu asterpingyuan shanqucha1 pingyuancha1 sqcmax1 sqcmin1 pycmax1 pycmin1 sqcmean1 pycmean1 sqcstd1 pycstd1 %进行山区(aster)3数据点均值赋值 k= ASTGTM2_N33E110_dem ; k2=k; sizk=size(k2); for i=4:3:sizk(1)-1 for j=4:3:sizk(2)-1 %从第四行第四列开始算起，避免计算最后一行最后一列数据 m=k(i-1,j-1),k(i-1,j),k(i-1,j+1),k(i,j-1),k(i,j),k(i,j+1),k(i+1,j-1),k(i+1,j),k(i+1,j+1); %提取3数据点及周围的八个数据点 k2(i,j)=mean(m); %九个数据点求平均值并赋值给3数据点 end end %进行平原区3数据点均值赋值 g = ASTGTM2_N34E114_dem; g2=g; sizg=size(g2); for i=4:3:sizg(1)-1 for j=4:3:sizg(2)-1 %从第四行第四列开始算起，避免计算最后一行最后一列数据 m=g(i-1,j-1),g(i-1,j),g(i-1,j+1),g(i,j-1),g(i,j),g(i,j+1),g(i+1,j-1),g(i+1,j),g(i+1,j+1); g2(i,j)=mean(m); %九个数据点求平均值并赋值给3数据点 end end clear i j m sizk sizg shanqucha2=k(1:3:3601,1:3:3601)-k2(1:3:3601,1:3:3601); %对应3数据点相减后所得的矩阵 pingyuancha2=g(1:3:3601,1:3:3601)-g2(1:3:3601,1:3:3601); %对应3数据点相减后所得的矩阵 %山区和平原区(插值前后的)差值矩阵的比较 sqcmax2=double(max(max(shanqucha2); pycmax2= double(max(max(pingyuancha2); sqcmin2= double(min(min(shanqucha2); pycmin2= double(min(min(pingyuancha2); sqcmean2= double(mean2(shanqucha2); pycmean2= double(mean2(pingyuancha2); sqcstd2= double(std2 (shanqucha2); pycstd2= double(std2 (pingyuancha2); jieguo2=sqcmax2 sqcmin2 sqcmean2 sqcstd2;pycmax2 pycmin2 pycmean2 pycstd2jieguo2 = 34.0000 -42.0000 0.0020 3.9706 15.0000 -15.0000 -0.0013 1.9601 clear shanqucha2 pingyuancha2 sqcmax2 pycmax2 sqcmin2 pycmin2 sqcmean2 pycmean2 sqcstd2 pycstd2 %srtm数据与插值后的aster数据分析 shanqucha3=srtmshanqu(2:1200,2:1200)-k2(4:3:3598,4:3:3598); %对应3数据点相减后所得的矩阵 pingyuancha3=srtmpingyuan(2:1200,2:1200)-g2(4:3:3598,4:3:3598);%对应3数据点相减后所得的矩阵 sqcmax3= double(max(max(shanqucha3); pycmax3= double(max(max(pingyuancha3); sqcmin3= double(min(min(shanqucha3); pycmin3= double(min(min(pingyuancha3); sqcmean3= double(mean2(shanqucha3); pycmean3= double(mean2(pingyuancha3); sqcstd3= double(std2 (shanqucha3); pycstd3= double(std2 (pingyuancha3); jieguo3=sqcmax3 sqcmin3 sqcmean3 sqcstd3;pycmax3 pycmin3 pycmean3 pycstd3jieguo3 = 262.0000 -155.0000 -6.1303 14.3963 112.0000 -101.0000 0.1857 4.3857 clear shanqucha3 pingyuancha3 sqcmax3 pycmax3 sqcmin3 pycmin3 sqcmean3 pycmean3 sqcstd3 pycstd3 k2 g2 %移动曲面拟合法内插3点坐标 zs=k; zp=g; dx=1 1 1 0 0 -1 -1 -1; %待插值点距离周围八个数据的X坐标差 dy=1 0 -1 1 -1 1 0 -1; %待插值点距离周围八个数据的Y坐标差 m1=dx.2; m2=dx.*dy; m3= dy.2; m4=dx; m5=dy; m6=ones(8,1); %生成一个8*1元素为1的矩阵 M=m1 m2 m3 m4 m5 m6; %系数矩阵 clear dx dy m1 m2 m3 m4 m5 m6 %清除垃圾变量 P=diag(0.5 1 0.5 1 1 0.5 1 0.5); %权阵(待插点与周围数据点的距离的平方的倒数) for i=4:3:3598 for j=4:3:3598 vk=double(k(i-1,j-1), k(i-1,j), k(i-1,j+1), k(i,j-1), k(i,j+1), k(i+1,j-1), k(i+1,j), k(i+1,j+1); %下面的*运算要求数据为双精度型 xk=inv(M*P*M)*M*P*vk; zs(i,j)=xk(length(xk); %xk中最后一个数据便为3插值，将其赋值到zs中 vg=double(g(i-1,j-1), g(i-1,j), g(i-1,j+1), g(i,j-1),g(i,j+1), g(i+1,j-1), g(i+1,j), g(i+1,j+1); %下面的*运算要求数据为双精度型 xg= inv(M*P*M)*M*P*vg; zp(i,j)=xg(length(xg); %xg中最后一个数据便为3插值，将其赋值到zp中 end end clear i j vk vg xk xg M P %清除垃圾变量 %比较(aster)拟合前后数据之间的差值 shanqucha4=zs(1:3:3601,1:3:3601)-k(1:3:3601,1:3:3601); %对应3数据点相减后所得的矩阵 pingyuancha4=zp(1:3:3601,1:3:3601)-g(1:3:3601,1:3:3601); %对应3数据点相减后所得的矩阵 sqcmax4= double(max(max(shanqucha4); pycmax4= double(max(max(pingyuancha4); sqcm