时序数值型遥感专题数据变化过程的地图可视化表达方法与流程

本发明涉及一种时序数值型遥感专题数据变化过程的地图可视化表达方法。

背景技术：

针对数值型时序数据的变化分析，现有技术主要采用趋势分析的方法，在地图上通过色彩的不同反映变化的趋势，只能反映整个时段的总体变化趋势，无法在地图上反映具体的变化过程，无法同时在一张地图上表示突变产生的时间、不同阶段的变化趋势等信息。或者只能针对个别的点进行断点检测和分段趋势分析，然后分别针对变化类型、主要变化的持续时间、主要变化时段的变化趋势等进行地图可视化，而无法在一张地图上进行多种信息的同时表达。

技术实现要素：

本发明的目的在于提供一种时序数值型遥感专题数据变化过程的地图可视化表达方法，以克服现有技术中存在的缺陷。

为实现上述目的，本发明的技术方案是：一种时序数值型遥感专题数据变化过程的地图可视化表达方法，包括如下步骤：

步骤s1：构建时间序列数值型遥感专题数据集；

步骤s2：对象单元划分；

步骤s3：进行总体趋势估计运算；

步骤s4：进行趋势分割；

步骤s5：进行断点判别；

步骤s6：根据断点判别的结果，进行分段趋势估计运算；

步骤s7：进行地图可视化表达。

相较于现有技术，本发明具有以下有益效果：本发明所提出的一种时序数值型遥感专题数据变化过程的地图可视化表达方法，把时序数值型数据的变化过程在一张地图上表达出来，可以同时在一张地图上反映突变点、持续时间、变化方向和强度等多种信息，为时序数据变化分析提供了一种新的地图可视化方法。可以应用于年际尺度的变化分析，如植被变化、温度变化、径流变化等。

附图说明

图1为本发明一种时序数值型遥感专题数据变化过程的地图可视化表达方法的流程图。

图2为本发明一实施例中a县2000-2015年植被覆盖度变化趋势示意图。

图3为本发明一实施例中a县2000-2015年植被覆盖度断点数目及趋势示意图。

图4为本发明一实施例中基于多尺度分割对象单元的a县2000-2015年植被覆盖度fvc断点检测与分段变化趋势示意图。

图5(a)为本发明一实施例中基于栅格的a县2000-2015年植被覆盖度fvc断点检测与分段趋势鸟瞰图。

图5(b)为本发明一实施例中基于栅格的a县2000-2015年植被覆盖度fvc断点检测与分段趋势由西部至东部方向的侧视图。

图5(c)为本发明一实施例中基于栅格的a县2000-2015年植被覆盖度fvc断点检测与分段趋势由东部至西部方向的侧视图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

本发明通过一种时序数值型遥感专题数据变化过程的地图可视化表达方法，通过准备时间序列数值型遥感专题数据集，进行对象单元划分，按单元进行断点检测，并分别计算不同时间段的变化趋势。按分段情况建立柱状图，柱子的分段数目代表变化的不同阶段，柱子标注变化时间，柱子高度代表绝对变化强度的大小，柱子颜色代表正向变化、负向变化或基本不变，对象单元可以根据不同的颜色代表整个时段的总体变化趋势。

为了让本领域技术人员进一步了解本发明所提出的方法，下面结合具体技术内容进行说明。还包括如下步骤：

步骤s1：数据集构建。

进一步的，在本实施例中，构建时间序列数值型栅格数据集。

步骤s11：对于低空间分辨率遥感专题数据，如modis的归一化植被指数ndvi、植被覆盖度fvc、叶面积指数lai、地表温度lst，avhrrgimmsndvi等，一般以旬、月等为间隔，年内具有周期性变化，时间序列几十年，进行s-g滤波处理；

步骤s12：对于中等空间分辨率遥感数据，一般以年度为间隔，属于非周期性数据。如landsat系列数据，选择云雨比较少的秋季作为预设时间段影像，并对影像进行几何校正和大气校正；

步骤s13：以年季为对象，如果中等空间分辨率遥感数据不完整，如缺乏某个年季，则通过如下时空融合方法生成缺失年份的数据：

式中，m为中等分辨率反射率数据，l为低空间分辨率反射率数据，ω为窗口大小；(xω/2,yω/2)为参与融合的中心像元，wij为归一化权重系数，由光谱、时间与空间距离三个维度构成；x、y分别代表空间位置横坐标和纵坐标，其下标分别代表不同的位置；t代表时间，t0为已有的数据获取时间，tp为拟融合预测的数据时间；

步骤s14：对中等空间分辨率遥感数据(如landsat数据)计算生成中等空间分辨率遥感专题数据，包括：归一化植被指数ndvi、植被覆盖度fvc、叶面积指数lai、地表温度lst等。

步骤s2：对象单元划分。

进一步的，在本实施例中，根据研究目的，按一定的规则划分对象单元，可以是地貌单元、地块、小流域、林业小班；或者按行政单元；或者直接按栅格数据集本身特征进行聚类来划定等，也可以直接以栅格作为对象单元。

步骤s3：总体趋势分析。

进一步的，在本实施例中，采用最小二乘回归分析法或sen+mann-kendall趋势分析法进行总体趋势估计。

步骤s4：趋势分割。

进一步的，在本实施例中，根据时序数据的特点和实际需要，获取各点水平移动阈值。

步骤s41：对于周期性数据而言，利用基于loess的周期-趋势分解(stl)方法将原始时间序列分解为周期项、趋势项和余项，并用趋势项参与趋势分割计算；非周期性数据可直接转至步骤s42进行断点检测计算；

步骤s42：先利用f函数获取波峰、波谷点

δvi(i-1，i)＝vi(i)-vi(i-1)

δvi(i，i+1)＝vi(i+1)-vi(i)

式中，sign为符号函数，vi表示某一数值型遥感专题数据，下标i、i-1或i+1表示不同的年份或时间；所有f＝1的点就是时间序列中的波峰、波谷点，而f＝0的点包括无明显变化的点和可能存在大幅变动的转点。

步骤s43：使用g函数获取转点

式中，g函数的作用是标记f＝0的点中大于距离阈值条件的转点，d(i)表示相邻波峰、波谷之间的连线和相距最远数据点间的最小垂直距离，ε是距离阈值，可以根据实际需要自定义。

步骤s44：使用h函数获取局部趋势变化(trendlocalchange，tlc)

式中，局部趋势变化，即tlc是相邻两个转点i和z(即g＝1的点)之间的差值，非转点的tlc值被赋为0，全部计算完成后将tlc序列按降序排列，tlc值与变化幅度成正比。

步骤s5：断点判别。

进一步的，在本实施例中，用贝叶斯信息准则(bayesianinformationcriterion,bic)推测断点数量n。bic准则能够推测合适的断点数量，可显著降低时间序列最小二乘拟合的误差，而又不会造成过拟合现象。

步骤s6：分段趋势估计。

进一步的，在本实施例中，根据bic得到的断点数量n，把各对象单元时间序列数据分为n+1个阶段，对每个阶段进行趋势分析，通过拟合出变化斜率作为各分阶段的变化强度。采用最小二乘回归分析法或sen+mann-kendall趋势分析法分别对n+1个阶段进行分段趋势估算。

步骤s7：地图可视化表达。

进一步的，在本实施例中，引入分层柱状和颜色，对每个对象单元进行地图可视化表达。每个单元的底色代表总体变化趋势，采用渐变色的形式；每个单元柱状分为n+1段：每段变化趋势以颜色表示，一般正向变化趋势用蓝色，负向变化趋势用红色，基本不变用黄色(颜色可以根据个人喜好进行调整)；每段柱子两端标注突变时间；柱子高度用趋势分析的斜率绝对值，其中基本不变化的采用固定长度表示。每个对象单元的柱状图包括了柱状的段数和时间(代表发表变化时间及整个时段变化分段情况)、颜色(代表变化趋势)、高度(代表变化强度)，具体表达的详细程度可以根据实际情况选用。

为了让本领域技术人员进一步了解本发明所提出的方法，下面结合具体实例进行说明。

进一步的，在本实施例中，传统植被变化趋势表达为：(1)总体趋势，如图2所示；2、断点数目及某个断点分段趋势，如图3所示。采用本实施例提出的方法的表达为：如图4以及图5所示。

如图4所示，为利用2000-2015年空间分辨率为250m，时间分辨率为16天的modisndvi产品生成的时间序列植被覆盖度fvc数据集。首先进行相似对象单元划分，计算各对象单元总体趋势，底图颜色从浅绿到深绿表示不同对象在2000-2015年间的总体变化趋势，颜色越深表示植被覆盖度增加幅度越大。柱子旁边的数值表示断点开始的时间，柱子颜色代表正向或负向变化，柱子高度表示变化的强度，其中平稳变化的高度不代表变化强度。

如图5所示，为基于栅格的a县2000-2015年植被覆盖度fvc断点检测与分段趋势图，以栅格作为对象单元，柱子的颜色代表植被覆盖度上升、下降或平稳，柱子高低代表变化幅度，分段的柱子高度代表某一时间节点到下一时间节点之间植被的累积变化量，越高表明变化程度越剧烈。由于以栅格为单元，如果对a县每个像元逐一进行断点时间标记，则会造成柱子排列过密、相互拥挤的情况。断点时间可以通过对每个栅格查询的方式获得，如图5(c)中查询的部分表示2011年开始增长，植被累积增长变化为0.21。右下角是控制旋转和缩放的按钮，地图的漫游可直接通过鼠标拖曳控制。图5(a)是a县植被变化三维时空表达的鸟瞰图，图中左下角的图例表示堆叠柱形图的元素含义，蓝色表示有植被增加的趋势，黄色表示基本保持稳定，红色表示有下降的趋势，绿色的底图表示基本没有发生变化的区域。图5(b)是从a县西部至东部方向的侧视图，可以看到无论是断点数量还是变化级别都少于东部；图5(c)是从a县东部至西部方向的侧视图，从图中能够直观地发现a县中东部的植被变化由早期的降低到后来的升高状态。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。