一种基于矢量网络的地形湿度指数模拟方法与流程

本发明属于数字地形分析技术领域，特别涉及一种基于矢量网络的地形湿度指数模拟方法。

背景技术：

beven和kirkby(文献1)提出了地形湿度指数，它是一种常用的通过量化汇水累积量和局部地形的坡度来描述土壤湿度和表面饱和度空间分布特征的地形指标。自从地形湿度指数被提出，广泛应用于水文(文献2、3)、农业(文献4、5)、生态(文献6、7)等其他领域(文献8、9)。地形湿度指数的模拟精度在实际应用中是至关重要的且许多研究者仍致力于地形湿度指数模拟精度的提高。

地形湿度指数的模拟主要包括单位汇水面积(sca)和坡度的模拟。目前用于模拟地形湿度指数中的sca方法主要包括单流向算法和多流向算法。在单流向算法中，d8算法(文献10)是最常用的算法，该算法是将像元中的所有流量全部流向其周围八个像元中具有最大坡降的一个像元中。dinf算法(文献11)也是比较常用的算法，它将网格划分成三角面，并将其八个三角面中具有最大坡度的三角面的水流方向视为像元的水流方向。随机八方向算法(rho8)(文献12)有时也被用于模拟地形湿度指数，该算法是通过随机参数的引入，在周围的八个网格中随机选择水流汇聚的像元点。这些方法简单、高效并且在凹坡面上具有较强的处理优势，但是它们过于简化水流在坡面上的流动情况且未考虑到水流漫散。为了解决这些问题，多流向算法被用于模拟地形湿度指数。最早的fd8算法(文献13)和fmfd算法(文献14)都是根据各自的计算公式分配该网格中的流量到其周围所有的下游网格中的流量。tmfd算法(文献15)是将fd8和dinf相结合得到的，该算法就是一个考虑了多流向的dinf算法。相对于单流向算法，这些算法在处理平面或者凸坡面时具有优势，更适用于模拟sca，尤其是在水流漫散区域，但是流量的精确分配仍然是难题且在平坦地区表现较差。

对于地形湿度指数中的坡度计算方法而言，三阶反距离平方权差分算法(文献16)是最常用的算法，该算法通过计算网格单元到其八个相邻单元在x和y方向上的变化率来获得坡度。该方法简单高效但可能不能较好地描述水流，因为它只考虑了一个单元距离内的下坡地形。为了更好的描述水流情况，hjerdt等(文献17)提出了通过计算预设的垂直距离与相应的水平移动距离的比例来模拟坡度。但该方法依赖于预设的垂直距离。

另外，还有很多研究者提出了改进的地形湿度指数模拟算法。例如,yong等(文献18)提出了nmfd算法，该算法改进了计算汇水面积的数学方程并且能更加精确地计算有效等高线长度。另外，该算法还利用tfd算法计算平坦区域的坡度。qin等(文献19)提出了一种新的地形湿度指数模拟方法，该算法分别利用自适应多流向算法和最大坡降算法模拟sca和坡度。ma等(文献20)将地形位置考虑进fd8得到了改进的fd8算法。尽管现在有许多的模拟地形湿度指数方法，但是它们对水流路径确定算法和坡度计算方法的依赖仍是一个关键问题。

现有的发明多数是关于地形湿度指数的应用。例如，史舟等(文献21)提出了一种基于低分辨率卫星遥感数据的环境变量动态筛选建模降尺度方法。它先将由1km分辨率的地形湿度指数、植被指数、数字高程模型、白天地表温、晚上地表温、坡度、坡向、坡长坡度聚合得到的25km分辨率的这8个环境变量因子作为自变量，25km的降水数据作为因变量，然后建立分区的多元回归模型，最后将1km的8个环境变量因子带入回归模型，得到1km的降水数据。最终能够比较准确地预测复杂地区的降水数据。董张玉等(文献22)提出了一种基于ndvi与lswi相结合的湿度信息提取方法。该发明首先利用经过预处理的遥感数据计算得到ndvi和lswi(湿度指数)，然后建立湿度信息提取的决策模型，最终不受季节限制地实现了湿地信息提取，从而为生态环境的规划和管理等提供支撑依据。冯光胜等(文献23)提出了一种冲洪积扇信息的自动提取方法。该发明通过对遥感提取指数、植被指数和地形湿度指数经过彩色合成的图像进行主成分变换之后得到第一主成分灰度图像并对其进行阈值分割，最终能够更加精确地获得冲洪积扇信息。

通过对以上内容的了解，发现应用中地形湿度指数模拟精度的提高仍至关重要。因此，本发明旨在提出一种基于矢量网络的地形湿度指数模拟方法，提高地形湿度指数的模拟精度，满足实际应用中精度的需求。

文献1.k.j.bevenandm.j.kirkby,1979.aphysicallybased,variablecontributingareamodelofbasinhydrology[j].hydrologicalsciencesbulletin,24(1),43–69；

文献2.patriciarull,2016.hydrologicalmodelingenhancementusinggis:animprovedtopographicwetnessindexforwetlanddetection[d].sweden,stockholmuniversity；

文献3.peitao,qincheng-zhi,zhua-xing,yanglin,luoming,libaolinandzhouchenghu,2010.mappingsoilorganicmatterusingthetopographicwetnessindex:acomparativestudybasedondifferentflow-directionalgorithmsandkrigingmethods[j].ecologicalindicators,10,610-619；

文献4.b.p.buchanan,m.fleming,r.l.schneider,b.k.richards,j.archibald,z.qiu,andm.t.walter,2014.evaluatingtopographicwetnessindicesacrosscentralnewyorkagriculturallandscapes[j].hydrologyandearthsystemsciences18,3279-3299；

文献5.a.posluschnyandk.lambers,2007.topographicwetnessindexandprehistoriclanduse[c].caa2007-layersofperception.proceedingsofthe35thinternationalconferenceoncomputerapplicationsandquantitativemethodsinarchaeology(caa),berlin,germany,april2–6,2007；

文献6.a.petroselli,f.vessella,l.cavagnuolo,g.piovesanandb.schirone,2013.ecologicalbehaviorofquercussuberandquercusilexinferredbytopographicwetnessindex(twi)[j].trees,27,1201-1215；

文献7.martinand2010.usingtopographicwetnessindexinvegetationecology:doesthealgorithmmatter？[j].appliedvegetationscience,13,450-459；

文献8.m.sujit,2015.upslopecontributingarea,topographicwetnessandlandsliding:acasestudyoftheshivkholawatershed,darjilinghimalaya[j].internationalresearchjournalofearthsciences3,23-29；

文献9.r.d.marjerison,h.dahlke,z.m.easton,s.seifertandm.t.walter,2011.aphosphorusindextransportfactorbasedonvariablesourceareahydrologyfornewyorkstate.journalofsoil&waterconservation,66(3),149-157；

文献10.j.f.o’callaghanandd.m.mark,1984.theextractionofdrainagenetworksfromdigitalelevationdata.computervisiongraphicsimageprocess,28,323–344；

文献11.davidg.tarboton,1997.anewmethodforthedeterminationofflowdirectionsandupslopeareasingriddigitalelevationmodels.waterresourcesresearch,33(2),662-670；

文献12.johnfairfieldandpierreleymarie,1991.drainagenetworksfromgriddigitalelevationmodels.waterresourcesresearch,27(5),709–717；

文献13.p.quinn,k.beven,p.chevalierando.planchon,1991.thepredictionofhillslopeflowpathsfordistributedhydrologicalmodelingusingdigitalterrainmodels[j].hydrologicalprocesses,5,59–79；

文献14.t.g.freeman,1991.calculatingcatchmentareawithdivergentflowbasedonaregulargrid[j].computers&geosciences,17(3),413-422；

文献15.janseibertandbrainl.mcglynn,2007.anewtriangularmultipleflowdirectionalgorithmforcomputingupslopeareasfromgriddeddigitalelevationmodels[j].waterresourcesresearch,43(4),306-320；

文献16.p.a.burroughandr.a.mcdonnell,1998.principlesofgeographicinformationsystems(p.333).newyork:oxforduniversitypress；

文献17.k.n.hjerdt,j.j.mcdonnell,j.seibertanda.rodhe,2004.anewtopographicindextoquantifydownslopecontrolsonlocaldrainage[j].waterresourcesresearch40,w05602；

文献18.binyong,li-liangren,yanghong,jonathanj.gourley,xichen,you-jingzhang,xiao-liyang,zeng-xinzhangandwei-guangwang,2012.anovelmultipleflowdirectionalgorithmforcomputingthetopographicwetnessindex[j].hydrologyresearch,43,135-145；

文献19.cheng-zhiqin,a-xingzhu,taopei,bao-linli,thomasscholten,thorstenbehrensandcheng-huzhou,2011.anapproachtocomputingtopographicwetnessindexbasedonmaximumdownslopegradient[j].precisionagriculture,12,32-43；

文献20.majianchao,linguanfa,chenjunmingandyangliping,2010.animprovedtopographicwetnessindexconsideringtopographicposition[c].internationalconferenceongeoinformatics,1-4；

文献21.史舟,马自强,吕志强,刘用.一种基于低分辨率卫星遥感数据的环境变量动态筛选建模降尺度方法:中国,201610307333.8[p].2016-10-12；

文献22.董张玉,王宗明,刘殿伟,任春颖,贾明明,丁智.一种基于ndvi与lswi相结合的湿地信息提取方法:中国,201310520474.4[p].2014-05-07；

文献23.冯光胜,杨树文,高山,顾湘生,李小和,赵新益,曹柏树.冲洪积扇信息的自动提取方法:中国,201110191094.1[p].2011-11-02.

技术实现要素：

为了解决现有的地形湿度指数模拟算法无法满足实际应用需求中的精度需求问题，本发明提供了一种基于矢量网络的地形湿度指数模拟方法以实现高精度的地形湿度指数模拟。

本发明所采用的技术方案是：一种基于矢量网络的地形湿度指数模拟方法，其特征在于，包括以下步骤：

步骤1：利用dem生成规则三角面的矢量网络；

步骤2：利用三角面网络tfn算法从矢量网络上获取汇水累积量，计算步骤1中所得到的矢量网络上每个三角面的水流方向，并选取每个三角面的重心为降雨源点，利用tfn算法追踪得到流水线路径，统计每个像元上通过的流水线数量，即是每个像元的汇水累积量；

步骤3：从步骤1中得到的矢量网络上获取坡向和坡度；

步骤4：利用步骤2中得到的汇水累积量和步骤3中得到的坡向计算单位汇水面积sca；

步骤5：将步骤4中的sca和步骤3中得到的坡度相结合模拟地形湿度指数。

依照本发明所提供基于矢量网络的地形湿度指数模拟方法，可以高精度地进行地形湿度指数模拟,从而解决现有算法无法满足实际应用中的精度需求问题。因为地形湿度指数是一个重要的复合地形属性，因此本发明特别应用于数字地形分析领域。

附图说明

图1为本发明实施例的流程图；

图2为本发明实施例的利用dem生成规则三角面的矢量网络原理示意图；

图3为本发明实施例的利用三角面网络(tfn)算法从矢量网络上获取汇水累积量原理示意图；

图4为本发明实施例的利用新算法计算坡向和坡度的原理示意图；

图5为本发明实施例的基于矢量网络的地形湿度指数模拟方法精度评估示意图。

具体实施方法

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明要解决的核心问题是：利用矢量网络上每个三角面恒定的坡度和坡向解决了水流路径的不确定性和计算坡度坡向模型本身的误差问题，从而得到一个高精度的地形湿度指数模拟方法。

请见图1，本发明提供的一种基于矢量网络的地形湿度指数模拟方法，包括以下步骤：

步骤1：利用dem生成规则三角面的矢量网络；

请见图2，具体实施时，用2x2的窗口将dem剖分成规则三角面，将4个像元的中心点(假设为①②③④)连接起来形成三角面的公共边。利用双变量立方样条曲面来模拟p点的高程值(假设为a)，计算①和③高程的平均值(假设为a1)，计算②和④高程的平均值(假设为a2)。如果a1更接近于a，则连接①和③进行剖分，否则连接②和④进行剖分。

步骤2：利用三角面网络(tfn)算法从矢量网络上获取汇水累积量，计算步骤1中所得到的矢量网络上每个三角面的水流方向，并选取每个三角面的重心为降雨源点，利用tfn算法跟踪得到流水线路径，统计每个像元上通过的流水线数量，即是每个像元的汇水累积量；

请见图3，具体实施时，根据步骤1中得到的矢量网络上每个三角面三顶点的坐标值计算每个三角面的坡度和坡向，从而得到每个三角面的水流方向，每个三角面的坡向是三角面的水流方向，坡度表示的是水流长度。在此基础上，选取每个三角面的重心为降雨源点，结合每个三角面的水流方向，利用tfn算法追踪得到从每个降雨源点出发的所有水流路径，统计每个像元上通过的流水线数量。每个像元上经过的水流路径的数量即是每个像元的汇水累积量。更加详细的流水线追踪方法可参见q.zhou,p.andy.chen,2011.estimatingsurfaceflowpathsonadigitalelevationmodelusingatriangularfacetnetwork[j].waterresourcesresearch,47(7),1-12。例如，图3中的以p为中心点的像元的汇水累积量即是5。

步骤3：利用新算法从步骤1中得到的矢量网络上获取坡向和坡度；

从矢量网络上获取坡向和坡度的新算法是通过遍历矢量网络上的三角面获得以待求坡向和坡度的像元中心点(假设为s)为顶点的所有三角面(假设为t1,t2,…,tn,0＜n≤8)，这些三角面的坡向(假设为α1,α2,...,αn,0＜n≤8)和坡度(假设为β1,β2,...,βn,0＜n≤8)的平均值即是该像元的坡向(坡度),坡向坡度具体实施步骤如下：

步骤3.1：基于新算法从步骤1中得到的矢量网络上获取坡向；

请见图4，具体实施时，遍历步骤1中的矢量网络上每个三角面，寻找以像元中心点p为顶点的所有三角面，假设是如图4所示的6个三角面(①②③④⑤⑥)，根据这6个三角面三个顶点的坐标值计算得到坡向(假设为α1,α2,α3,α4,α5,α6)，那么p点坡向

步骤3.2：基于新算法从步骤1中得到的矢量网络上获取坡度。

请见图4，具体实施时，遍历步骤1中的矢量网络上每个三角面，寻找以像元中心点p为顶点的所有三角面，假设是如图4所示的6个三角面(①②③④⑤⑥)，根据这6个三角面三个顶点的坐标值计算得到坡度(假设为β1,β2,β3,β4,β5,β6)，那么p点坡度

步骤4：利用步骤2中得到的汇水累积量和步骤3中得到的坡向计算单位汇水面积(sca)；

具体实施时，利用步骤2计算得到的汇水累积量(假设为flowaccum)和步骤3.1计算得到的坡向(假设为α)，结合costa-cabral和burges(1994)的公式sca＝flowaccum*(g/2)/((fabs(cosα)+fabs(sinα)))计算得到sca，其中g为dem栅格单元大小。

步骤5：将步骤4中的sca和步骤3中得到的坡度相结合模拟地形湿度指数。

具体实施时，将步骤4中计算得到的sca和步骤3.2计算得到的坡度(假设为β)，结合beven和kirkby(1979)的公式twi＝ln(sca/tan(β))模拟得到地形湿度指数，其中twi为地形湿度指数。

请见图5，本实施例对基于矢量网络的地形湿度指数模拟方法的精度进行评估，包括定量评价和定性评价：

定量评价，选择均方根误差(rmse)为评价指标，与其他常用的地形湿度指数模拟方法进行对比验证基于矢量网络的地形湿度指数模拟方法的精度；

具体实施时，利用zhou和liu(2002)提出的四种数学曲面(凸椭球面、凹椭球面、平面和马鞍面)得到不同分辨率(例如1m，5m，10m，20m，…)的dem，基于数学曲面的表达式分别计算得到不同数学曲面下不同分辨率下的理论地形湿度指数(理论值)。其他常用的地形湿度指数模拟方法：第一种是利用d8算法模拟sca，利用三阶反距离平方权差分算法计算坡度，从而得到地形湿度指数；第二种方法是利用dinf算法模拟sca，利用三阶反距离平方权差分算法计算坡度，从而得到地形湿度指数；第三种方法是利用fd8算法模拟sca，利用三阶反距离平方权差分算法计算坡度，从而得到地形湿度指数；第四种方法是利用自适应多流向算法模拟sca，利用最大坡降算法计算坡度，从而得到地形湿度指数。利用这四种常用的地形湿度指数模拟方法和基于矢量网络的地形湿度指数模拟方法分别计算不同数学曲面下不同分辨率下的地形湿度指数(模拟值)，计算每种方法在不同数学曲面下不同分辨率下的rmse进行精度的定量评价。

在实例中，rmse的计算公式如下：

其中，t′i是理论值，ti是模拟值，n表示像元总数，下标i表示第i单元，该值越小表示模拟精度越高。

定性评价，地形湿度指数空间分布特征作为定性评价参考。

具体实施时，将基于矢量网络的地形湿度指数模拟方法和步骤6.1中四种常用的地形湿度指数模拟方法分别应用于真实的dem数据，对计算得到的地形湿度指数空间分布特征(例如在河流的干流上会因为水流的汇集地形湿度指数偏高和在干流的周围因为水流的漫散地形湿度指数偏低等其他特征)进行分析从而进行精度的定性评价；如果基于矢量网络的地形湿度指数模拟方法所得到的地形湿度指数在河流干流和支流上相对偏高，在河流干流和支流周围相对偏低，在平坦地区基本相等，那表示基于矢量网络的地形湿度指数模拟方法的模拟精度更高。那是因为水流在河流干流和支流上容易汇聚使得地形湿度指数相当偏高；水流在河流干流和支流周围容易散漫使得地形湿度指数相对偏低；水流在平坦区域分布比较均匀使得地形湿度指数基本相等。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。