DEM随机噪声误差去除方法研究--以黄土沟沿线地形剖面为例.docx

1、地理与地理信息科学GeographyandGeo-InformationSciencedoi：10.3969/j.issn.1672-0504.2017.02.003DEM随机噪声误差去除方法研究以黄土沟沿线地形剖面为例曹建军",那嘉明】，汤国安七杨昕】（1.南京师范大学地理科学学院，虚拟地理环境教育部重点实验室，江苏南京210023；2.南京晓庄学院环境科学学院，江苏南京211171）摘要:在黄土高原数字地形分析中，基于DEM提取的地形参数和地貌对象是研究黄土地形地貌特征的基础。但在使用这些对象和参数进行小流域尺度的数字地形分析时,DEM数据中隐含的噪声在一定程度上会影响分析结果的

2、可靠性。该文以分别位于陕北黄土丘陵沟壑区（绥德）和黄土堀区（淳化）的两个小流域为试验区，以国家基础地理数据库5m分辨率DEM为实验数据,以黄土地貌中最具典型特征的沟沿线为研究对象，对其地形剖面所隐含噪声的特征、提取及去除方法进行了研究。实验中采用EMI>MFDFA和EEMI>MFDFA两种组合方法对上述两个小流域黄土沟沿线地形剖面数据进行了分析。实验结果显示，沟沿线地形剖面数据噪声与地形复杂程度有很大关系洞EMD-MFDFA方法相比,EEMDMFDFA对去除噪声成分具有更好的效果；EEMDMFDFA方法在上述两个样区的沟沿线原始剖面与降噪后重构地形剖面的RMSE分别为2.053m和

3、3.188m0关键词:DEM；沟沿线；噪声；EEMD；MF-DFA中图分类号:P931文献标识码：A文章编号：1672-0504（2017）02-0012-070引言在基于DEM进行黄土高原地貌形态的分析与研究中,DEM误差与不确定性极大地影响分析结果的可信性。DEM误差类型可分为人为误差、系统误差、随机误差等。对于前两种误差，很多学者从误差的来源、分布规律2F和解决方法等做了深入的研究和分析。但对于DEM随机噪声的来源、影响因素和去除方法的研究还很有限。通常，地表可通过公式表达，其中Z表示点危以）的高程值,从信息论的角度考虑，地表可被认定为一个二维的非线性、非平稳的信号。这些信号中的细节信号

4、越多，与局部初始地形就越接近，但加入过多的细节信号可能会带来随机噪声问题。由于在信号数据中噪声普遍存在，因此在很多研究领域，如GPS坐标时间序列数据、反映地表地物变化的遥感数据、水文时间序列数据等分析中，均需要先去除噪声。在基于要素采样的地形分析中，地形随机噪声信号通常与地形分析的尺度有关9】。在大尺度的地形分析中，这种随机噪声可以被忽略，但在小尺度（如小流域）的地形分析中，这种影响局部地形细节变化的随机噪声就不能忽略，否则地形的局部细节就不能反映真实的地形变化，对揭示流域的宏观地形变化趋势造成困难。对于如何消除或减小数据噪声，前人提出了很多方法，如小波变换皿、集合经验模态分解、希尔伯特一黄变

5、换方、马尔科夫变换14、低秩矩阵重建口刃等方法，且均取得了很好的效果。其中，集合经验模态分解（EnsembleEmpiricalModeDecomposition,EEMD）方法作为经验模态分解（EmpiricalModeDecomposition,EMD）方法算的改进，在数据噪声分解中具有局部细节自适应性的优势，在地震信号17,、SAR信号、电磁信号2。等方面的噪声分解中均可取得很好的效果，但无法实现对噪声的识别。而多重分形消除趋势分析（MultifractalDtrendedFluctuationAnalysis,MF-DFA）方法提出的广义Hurst指数对细节随机信号噪声具有很好的量化识

6、别能力，将其与EEMD方法相结合，可有效解决噪声去除问题2习。但以上方法在复杂地形地貌的DEM噪声去除中的应用尚待探索。剖面分析是基于DEM的数字地形分析方法的重要方法S,前人已对地形剖面25、流域边界剖面等进行了研究。沟沿线是黄土地貌正负地形的分界线，地形在沟沿线处存在较大的高程变化和坡收稿日期：2016-09-08；修回日期：2017-02-07基金项目：国家自然科学基金项目（41671389.41471316.41401441）作者简介：曹建军（1976）,男，博士后，副教授，研究方向为黄土高原数字地形分析。*通讯作者E-mail：tangguoan度转折27.28,以沟沿线处的地形高程

7、采样数据所形成的地形剖面为对象进行分析，对研究黄土高原小流域局部地形变化规律具有指示意义。但这种沟沿线采样数据往往包含噪声，因此在分析时必须考虑噪声去除问题。本文使用EEMD和MF-DFA、EMD与MF-DFA相结合的方法，基于黄土沟沿线5mDEM高程采样数据，对采样数据噪声进行分析，得到数据噪声成分的识别和去除方法，并对比两种方法的分析结果和准确程度。1方法与原理1.1研究区和数据本文选择分别位于地形极为破碎的黄土丘陵沟壑区（绥德）和相对平缓的黄土堀区（淳化）的两个典型小流域为研究样区（表1）,采用国家基础地理数据库5m分辨率DEM为实验数据，对于已提取的经过校正的两个小流域的沟沿线数据（图

8、1）,从流域的出口开始，基于DEM沿沟沿线顺时针走向对高程数据进行依次采样，得到高程随采样距离变化的二维信号数据，即沟沿线地形剖面。其中，绥德样区的采样点有16000多个,淳化样区的采样点有7000多个。1.2研究方法经验模态分解（EMD）方法K依据数据自身的时间尺度特征进行信号分解，无须预先设定任何基函数，与传统的建立在先验性的谐波基函数和小波基函数上的傅里叶分解与小波分解方法具有本质性的差别。所以，EMD方法在理论上可以应用于任何类型信号的分解，在处理非平稳及非线性数据上具有明显优势，适合于分析非线性、非平稳信号序列,具有很高的信噪比2刃。（b）淳化样区表1研究区描述Table1Geogr

9、aphicdescriptionofstudyareas地貌类型地理位置描述黄土弗状丘陵沟壑黄士墟110°15'00"110°22'30"E37°32'30"37°37'30"N108°22'30"108°30'00"E34°50'00”34°55'OO"N位于绥德.海拔814-1188m,平均坡度约29°沟壑密度9.39km/km2；该区域土壤侵蚀严重，沟壑纵横各部分与整

10、体的形态自相似性特征较强位于淳化,海拔768-1188m,地面平均坡度12°.沟粮密度1.79km/km2面地形平缓侵蚀作用较弱，周围往往被沟谷环绕图1研究样区数据及沟沿线地形剖面Fig,1Iliedataofstudyareaandprofileofshoulderline而集合经验模态分解（EEMD）方法是Wu等针对EMD方法混频的不足，提出的一种局部自适应的时间序列分析技术，其本质是一种噪声辅助数据分析方法。该方法通过将白噪声加入待分解信号,将信号和噪声的组合看成一个总体，进而得到模态一致的内禀函数（IntrinsicModeFunctions,IMFs）分量。当把信号加在整个

11、时频空间分布一致的白噪声背景上时，由于白噪声频谱均匀分布，不同尺度的信号会自动分布到合适的参考尺度上。由于噪声的零均值特性，经过多次平均运算处理后，噪声会相互抵消，集成均值的结果就可直接作为最终分解结果。多重分形消除趋势分析(MF-DFA)方法是一种十分有效的非平稳信号序列多重分形特征分析方法，该方法首先通过去趋势过程去除信号序列中的非平稳因素影响，再使用不同阶数的波动函数计算信号序列不同尺度层次上的标度行为，通过对各信号序列的多重分形谱参数特征的分析和对比，最后选择多重分形谱参数作为噪声信号的量化特征，可以充分揭示隐藏在非平稳信号序列中的噪声信号。本文分别使用两种不同的方法组合进行去噪对比实

12、验。首先对沟沿线地形剖面进行EMD和EEMD分解，获得高程信号数据的分解分量;然后使用MF-DFA方法对分解结果分量进行判定，以区分哪些分量为噪声，哪些分量为具有物理意义的高程信号数据;最后,选取具有物理意义的高程信号数据对沟沿线地形进行重构，得到去噪后的沟沿线地形剖面，对于沟沿线地形剖面中隐含噪声大小的判定，通过原始数据与重构数据的均方根误差(RMSE)来表征。2结果与讨论2.1数据分解与Hurst指数计算图2、图3依次为绥德样区、淳化样区沟沿线地形剖面的EMD和EEMD分解结果。受样本点数影响,绥德样区分解结果有13个IMFs,其中IMF13是一个残差趋势项(RES),而淳化样区分解结果有

13、12个IMFs,其中IMF12是一个残差趋势项(RES)0从两个样区的分解结果看，随着分解分量周期的增大,频率越来越低，因此最后的残差项应视为研究样区的地形变化总趋势。此外，由于EMD分解存在分解模式混叠的缺陷5,对于同一样区的EMD和EEMD分解结果具有一定差异。为了对两个研究样区沟沿线地形剖面数据IMFs分量中的噪声成分进行识别，本文采用了MF-DFA方法对所有IMFs分量进行分析，计算每个IMFs的Hurst指数(Hq)。将IMFs划分成不重叠的数据单元，并计算其多重分形消除趋势指数F(q)。5005050t.sol二1050fal000土950200040006000800010000

14、120001400016000样本点数(b)EEMD分解结果020004000600080001000()12000!400016000样本点数(a)EMD分解结果0500009502宜图2绥德样区沟沿线地形剖面分解结果对比Fig.2ResultcomparisonofEMDandEEMDmethodsinSuide200-20+'一，1AdAVJ人50j.,p-r11-T)f.i>yi,£SI2SI&W10-JWJ-5011*-ZJWl6U.N2U.ZXI样本点数(a)EMD分解结果样本点数(b)EEMD分解结果图3淳化样区沟沿线地形剖面分解结果对比Fig.3

15、ResultcomparisonofEMDandEEMDmethodsinChunhua在具体计算过程中，输入的参数(如数据单元的大小、波动阶数q和趋势阶数玖等)均采用了文献21关于MF-DFA算法的最优化结果。.图4、图5(彩图见封3)为两个样区两种分解方法结果的IMFs分量Hurst指数计算结果。从中可以看出，对于每个IMFs,Hurst指数(Hq)总体呈现随波动阶数q的增大而减小的趋势。当q=5时，所有IMFs分量的Hq均达到最小，因此在对各分量进行噪声判别时，可选取9=5作为Hq计算的统一标准。(a)EMD分解的Hurstlt数(b)EEMD分解的Hurst指数图4绥德样区Hurst指

16、数计算结果对比Fig.4HurstexponentresultsoftwodecompositionmethodsinSuide(a)EMD分解的Hurst指数£23450(b)EEMD分解的Hurst指数图5淳化样区Hurst指数计算结果对比Fig.5HurstexponentresultsoftwodecompositionmethodsinChunhua2.2沟沿线地形剖面噪声去除形数据变化的复杂性、非线性和非平稳性等特征。当Hq>0.5时,数据序列具有强的长程相关性,所以Hq<o.5的IMFs分量都可认为是噪声成分。因此在本研究中，对于EEMD分解结果,g=5时=

17、0.5可作为判断是否为噪声成分的阈值;但山于EMD分解时存在模式混叠i3。,在噪声判别时Hq阈值范围应取稍大一些，这里参考文献22取值为0.7。图6为q=5时两种分解结果的】MFs的Hurst指数变化趋势。从图6可以看出，对于两个样区的EMT)分解内禀函数中,IMF1.1MF2的Hq均小于0.7,可视为噪声成分;而对EEMD分解的内禀函数,IMF1JMF2JMF3和IMF4的他均小于0.5,可视为噪声成分。IMFsIMFs(a)绥德样区(b)淳化样区图6波动阶数q=5时EMD和EEMD分解的内禀函数(IMFs)的Hurst指数变化Fig.6ThevarietyofHurstexponentof

18、IMFsfromEMDandEEMDdecompositionatfluctuationorderq=5基于此，对绥德和淳化两个研究样区沟沿线地形剖面的噪声成分进行了提取(图7)0从图7中可以看出，对于两个研究样区,EMT>MFDFA方法提取的噪声成分的振幅较大，周期也较长，频率较低;而EEMI>MFDFA方法提取结果的振幅较小，周期短,但频率较高。这两种方法所提取的地形噪声具有完全不同的形态特征和物理含义，体现了黄土高原地样本点数(a)绥德样区20r样本点数(b)淳化样区图7EMD-MFDFA与EEMD-MFDFA消除地形数据的噪声.对比Fig.7Thenoiseremovalo

19、fterraindatabyusingEMD-MFDFAandEEMD-MFDFA图8为采用MF-DFA方法分别对EMD.EEMD分解结果噪声进行去除后，重构剖面与原始剖面的比较。对于绥德样区，两种方法去除噪声后的地形与原始地形的总体变化趋势一致，但在局部变化上具有差异，说明原始地形数据在局部变化尺度上存在噪声。而EMD-MFDFA和EEMI>MFDFA对局部噪声的去除效果也有不同，两种方法得到的重构地形与原始地形之间的均方根误差(RMSE)分别为3.704m和2.053mo对于淳化样区，原始地形与经过噪声去除后的地形之间的差别要比绥德样区大,说明淳化样区地形数据中所隐含的噪声成分更多,

20、其RMSE分别为7.686m和3.188mo3结论由于多种原因，DEM及基于DEM提取的沟沿线地形剖面均存在一定的数据噪声。对于不同类型EMD-MFDFA去噪地形，EMD-MFDFA去噪地形，0200040006000800010000120001400016000EEMD-MFDFA去噪地形1000OQA,二，q4.一0200040006000800010000120001400016000样本点数(a)绥德样区图8EMD-MFDFA和EEMD-MFDFA去除噪声的地形与原始地形比较ComparisonoftheterrainafternoiseremovalusingEMI>MFDF

21、AandEEMD-MFDFAwiththeoriginalterrainFig.8的地形地貌而言，其沟沿线地形剖面的数据噪声特征是不同的，地形变化复杂的黄土昴沟沿线数据噪声以低频噪声为主，地形变化相对单一的黄土堀沟沿线以高频噪声为主，因此基于一定的科学方法对数据噪声的去除是必要的。鉴于EMD、EEMD在噪声信号分解中具有局部自适应性，以及MF-DFA在局部噪声识别中具有很好的量化能力，两者相结合的方法对于沟沿线地形剖面噪声去除是有效的。实验结果表明，同EMD-MFDFA方法相比，EEMDMFDFA方法能够分离和去除更多的数据噪声成分，其沟沿线原始剖面与降噪后重构地形剖面的均方根更小，对沟沿线地

22、形剖面数据噪声去除具有更好的效果。实验结果显示，沟沿线地形剖面数据的噪声与地貌类型、地形复杂程度具有很大关系。对于黄土地貌数字地形分析，数据噪声与地貌类型、地形复杂度之间到底存在怎样的内在机理与规律，有待下一步进行深入的探索。参考文献：1 汤国安，赵牡丹,李天文，等.DEM提取黄土高原地面坡度的不确定性J.地理学报，2003,58(6)：824-830.2 杨灿灿，王春，江岭，等.DEM地形描述误差的概率分布特征研究口.测绘工程,2016,25(8):19-23.3 张勇，汤国安，彭酎.数字高程模型地形描述误差的量化模拟以黄土丘陵沟壑区的实验为例J.山地学报，2003,21(2):252-25

23、6.4 赵尚民，何维灿，王莉.DEM数据在黄土高原典型地貌区的误差分布J.测绘科学，2016,41(2)：67705J王春,顾留碗，陶昉，等.DEM地形描述误差(Et)计算模型研究EJ1地球信息科学学报,2014,18(5):699-706.6姜卫平，周晓慧.澳大利亚GPS坐标时间序列跨度对噪声模型建立的影响分析J.中国科学(地球科学)，2014,44(11):2461-2478.7J李喜佳，肖志强,王锦地，等.双集合卡尔曼滤波估算时间序列LAIEJ1遥感学报,2014,18(1):27-44.8 桑秀丽，苏俞真，肖汉杰，等.基于01测试方法的含噪声降雨一径流时间序列混沌特征分析口.云南大学学

24、报(自然科学版)，2014,36(2)：233240.9 孙嘉骏，黄天勇，陈科基于DEM地形要素提取方法的实验JJ.地理空间信息,2014,12(2):136-138.10 秦承志，呼雪梅栅格数字地形分析中的尺度问题研究方法EJ1地理研究,2014,33(2):270283.口1李耀恩，吴小平.小波变换在双频激电数据处理中的应用J.地球物理学进展,2015,30(1)：460465.12 WUZHfHUANGNE.Ensembleempiricalmodedecomposition：Anoiseassisteddataanalysismethod-J.AdvancesinAdaptiveDat

25、aAnalysis,2009,1(1)：141.13 时丕旭，杨志强,计国锋.希尔伯特一黄变换下的磁悬浮陀螺仪数据处理方法研究J.测绘通报,2015(Sl)：119-123.14 许春冬，战鸽，应冬文，等.基于隐马尔可夫模型的非监督噪声功率谱估计J.数据采集与处理,2015,30(2)：359364.口5孟樊，杨晓梅,周成虎.一种以低秩矩阵重建的图像混合噪声去除算法J.地球信息科学学报,2015,17(3)：344-352.16 HUANGNE,SHENZ,LONGSR,etaLTheEmpiricalModeDecompositionandtheHilbertSpectrumforNonli

26、nearandNonStationaryTimeSeriesAnalysisR.TheRoyalSociety,1998.903-995.17 颜中辉，王赞，李攀峰，等.基于希尔伯特一黄变换的多分量地震去麟物理学版,2015,30(6):28462854.口8贾瑞生，赵同彬,孙红梅，等.基于经验模态分解及独立成分分析的微震信号降噪方法口.地球物理学报，2015,58(3)：1013-1023.19郭际明，黄长军，喻小东，等.一种改进的二维经验模态分解SAR图像噪声抑制方法口.武汉大学学报(信息科学版)，2013,38(12):1445-1450.20J李成武，黄利辉,王启飞，等.煤岩微弱电磁信

27、号的噪声源识别及去噪方法J.煤炭学报,2016,41(8):1933-1940.21 KANTELHARDTJW,ZSCHIEGNERSA,KOSCIELNY-BUNDEE,etal.MultifractaldetrendedfluctuationanalysisofnonstationarytimeseriesJ.PhysicaA：StatisticalMechanicsandItsApplications,2002,316(1):87114.22 MIRIYALAS,KOPPIREDDIPR,CHANAMALLUSR.Robustdetectionofionosphericscintill

28、ationsusingMF-DFAtechnique.EarthsPlanetsandSpace,2015(67)：15.23 杨昕，汤国安,刘学军，等.数字地形分析的理论、方法与应用J.地理学报，2009,64(9)：10581070.24 汤国安.我国数字高程模型与数字地形分析研究进展J.地理学报，2014,69(9)：13051325.25 吴瑞，王兰辉，汤国安.中国陆地边界地形剖面谱研究初探口.地理与地理信息科学，2012,28(5)：5154.26 贾旖旎，基于DEM的黄土高原流域边界剖面谱研究D.南京:南京师范大学,2010.27 刘玮，李发源,熊礼阳，等.基于区域生K的黄土地貌沟

29、沿线提取方法与实验J.地球信息科学学报,2016,18(2)：220226.28 张磊，汤国安,李发源，等.黄土地貌沟沿线研究综述J,地理与地理信息科学,2012,28(6):44-48.29 戴吾蛟，丁晓利,朱建军，等.基于经验模式分解的滤波去噪法及其在GPS多路径效应中的应用J.测绘学报,2006,35(4)：321-327.30 AHMETM,AYDINA.DetrendedfluctuationthresholdingforempiricalmodedecompositionbaseddenoisingJ,DigitalSignalProcessing,2014(32)：4856.Me

30、thodforNoiseRemovalofDEMs：ACaseStudyonTerrainProfileofLoessShoulderLineCAOJianjun1*2,NAJiaming1,TANGGuoan1,YANGXin1(1.SchoolofGeographyScience,KeyLaboratoryofVirtualGeographicEnvironmentofMinistryofEducation,NanjingNormalUniversity9Nanjing210023;2.SchoolofEnvironmentScience»NanjingXiaozhuangUni

31、versityNanjing211171,China)Abstract：TerrainparametersandobjectsderivedfromtheDEMsarethebasisofgeomorphologicstudied,especiallyinloessplateauofChina.However,whenanalyzinginthescaleofwatershedthetopographicnoisesofDEMsareoftenignored.Inthispaper,takingSuideandChunhuaasthestudyareas,andtheloessshoulder

32、-linesasresearchobjects,thecharacteristics,extraction,andremovalmethodofthetopographicnoisesarestudied.Firstly»terrainelevationissampledbyloessshouldeilineextractedbyDEMs.Then,EMDandEEMDisperformedtoobtainthecomponentsofthesamplingdata.Finally,theHurstindex(Hq)inMF-DFAisanalyzedforthespatialdis

33、tributionofthetopographicnoises.Theresultsshowthat：1)thenoiseoftheterrainprofilealongtheshoulder-linehasagreatrelationshipwiththeterraincomplexity；2)comparedtotheEMD-MFDFAmethod,theEEMT>MFDFAmethcxlhasabettereffectinmitigatingthetopographicnoiseofloessshoulder-line；3)theRMSEbetweentheoriginalterr

34、ainandtheterrainafternoiseremovalusingtheEEMI>MFDFAmethodis2.053minSuideand3.188minChunhua.Keywords：DEM；shouldcilinc；noise；EEMD；MFDFA(上接第11页)ResearchonaTwo-LevelRetrievalMethodofDistributedSpatialDataServiceLIANGQi-yang1,WUBao-guo1,ZHAOChun-jiang2，ha。Xing-yao2*3*4,PANYu-chun2-3*4(1SchoolofJnforma

35、tionBeijingForestryUniversity»Beijing100083;2.NationalEngineeringResearchCenterforInformationTechnologyinAgricultureBeijing100097；3.KeyLaboratoryofAgri-informaticsMinistryofAgriculture,Beijing100097；4.BeijingEngineeringResearchCenterofAgriculturalInternetofThingsBeijing1000979China)Abstract：Atp

36、resent,theorganizationmethodofspatialdataservicebasedonserviceregistryandmetadatamanagementisoflowefficiencytoretrievedistributedmulti-scalespatialdataservice,whichcan'tmeetthedemandofintegratedapplicatiorxTosolvethisproblem,adistributedmulti-scalespatialdataindexingmethodbasedontwo-levelspatial

37、indexstructureisproposed.Basedonspatialgridindex,combinedwithspatialdataservicemetadatamanagement,indexingandretrievalofspatialdataserviceanditsspatialobjectarerealized.Inthedataretrievalprocess>accordingtothedatacontent,scopeandotherinformation,firstly,thefirstlevelindexisretrievedtofindouttargetspatialdataservice,thenthesecondlevelindexisusedtoretrievetargetdataobjectintheservice.Inourexperiments,theapplicationofindexconstructionandretrievali