我国国土资源丰富, 单草地面积就高达约392万km², 居世界第二位. 国土资源生态覆被不同类型之间的相互转换一方面取决于自然环境的变化, 另一方面, 国家的相关政策也对生态系统具有重要的影响, 比如退耕还林政策、退牧还草工程^[1]等都是对不同生态类型的再定义和重新分配, 以满足我国经济发展的需要和人民生产生活的需要. 生态系统的研究主要集中在不同类型之间的面积变化及其比例的对比, 从而分析变化趋势, 判断政策的执行力度和效果, 为下一步发展做好决策, 生态变化数据的研究是生态安全的基础技术支撑, 为生态决策提供支持.

当前我国生态覆被分析具有以下3点迫切需求:

需求1. 对多个类型生态覆被转移变化进行整体和局部的对比分析;

需求2. 对生态覆被时序变化的分析. 在考虑整体变化的影响下, 分析某一区域生态覆被变化;

需求3. 生态覆被空间划分分析. 对具有相同变化趋势的区域进行空间聚类, 研究不同空间的相似性.

本文针对以上需求研发了ECOVIS系统. 该系统利用可视化交互设计实现了对生态数据的可视分析. 解决了生态数据的空间划分和时序分析需求. 本文的主要贡献包括:

(1)设计了一种有效的转移矩阵可视化视图, 可以从全局和局部交互式地对生态变化进行可视化;

(2)实现了一套针对生态变化数据的可视分析系统ECOVIS, 可以在考虑整体变化背景下分析某一生态类型的时序变化;

(3)实现了基于T-SNE降维算法的散点图, 可以对生态数据进行空间划分.

本文第1节介绍了生态系统数据可视化的研究现状; 第2节介绍了基于转移矩阵设计的数据模型; 第3节介绍了针对数据模型的可视化方法; 第4节介绍了可视化交互方法; 第5节利用用户使用案例来证明本文方法的有效性; 最后一节对本文系统进行了总结和展望.

1 相关工作

随着遥感技术的发展, 以遥感数据作为生态系统监测与评价的基础已成为宏观生态学研究的重要手段^[2], 其中用到的研究数据特征中最为重要的是生态系统转移矩阵^[3]. 借助生态系统的转移矩阵特征, 可以较好地分析区域生态系统变化的特征, 研究各类型的转移转化关系^[4]. 矩阵数据也被越来越多的应用在生态系统研究的各个方面^[5-7]. 转移矩阵中包含多种特征, 关系错综复杂, 如何充分有效地利用生态系统转移矩阵深度挖掘生态系统的详细特征是当前生态系统研究领域的一个难题^[8], 以往该领域的研究主要利用人工统计的方法得出若干结论, 或通过曲线图对特定维度的数据进行可视化. 这些方式较难体现出生态数据内部包含的逻辑关系^[9]. 转移矩阵来源于系统分析中对系统状态与状态转移的定量描述. 一般用二维表来表示, 从中可以查看各个生态类型间相互转化的情况. 比如某一类别的土地转化成其他类型土地的比例等, 描述了前后两个时刻之间覆被发生转变的对比信息. 转移矩阵可以从栅格图像中计算得到, 也可以从两个矢量文件中计算获得. 以上研究主要通过表格和热力图的方式加以呈现. 这些方法较难给用户以直观的数据呈现, 并且只能刻画某一个时刻生态覆被类型的分布, 无法体现出生态转移的整体状态, 也无法较好地让用户对感兴趣数据做更进一步的可视化交互分析. 本文针对生态系统格点数据中隐含的转移特征信息, 设计有效的可视化算法和交互方案, 为生态覆被类型的转化特征从全部到局部提供可视化分析方法, 为分析生态系统中各覆被类型的变化趋势提供技术手段.

2 数据模型 2.1 数据介绍

本文实验数据是我国2000年和2010年的生态系统格点数据, 数据编码采用LUCC分类土地利用体系. 每一个数据格点代表不同的生态覆被类型, 包括耕地、林地、草地、水域、城乡、未利用等, 其中每一类包含多项子类, 比如未利用土地包括沙地、戈壁、盐碱地、沼泽地、裸土等多项子类. 根据格点数据代表的生态覆被类型, 对比两个年份的变化, 统计不同转换发生的总量累计和, 可以得到不同覆被类型之间转化的矩阵数据, 即转移矩阵. 矩阵对角线表示没有状态变化的覆被类型占比数量, 转移矩阵主要反映了区域内生态变化情况.

2.2 数据处理

如表1所示, 即为一个典型的转移矩阵, 转移矩阵统计局部区域格点的变化值, 该值基于格点的数量计算. 转移矩阵是一个非对称的矩阵, 矩阵的对角线和非对角线上的数据代表不同的转移类型, 需要分开进行可视化预处理. 以下将从对角线上的同类型转换和非对角线上的异类型转换两方面阐述.

定义转移矩阵E为M×M的矩阵:

$E = [{E_{ij}}](0 \le i < M,0 \le j < M)$

(1)

生态系统转移矩阵中, 通常转变的面积远小于未转变的面积, 所以一般情况下矩阵的对角线上的值远大于非对角线上的值. 在可视化的过程中, 如果将对角线上的元素考虑进去, 会严重影响可视化效果. 所以, 本文将矩阵进行了拆分, 非对角线上的元素构成的矩阵定义为异类型转换矩阵R, 即转换前后的生态系统种类不同; 对角线上的元素构成的矩阵定义为同类型转换矩阵I, 即转换前后的生态系统种类相同, 亦可理解为特定生态系统保留下来的面积构成的矩阵, 定义如下:

$\left\{\begin{split} & E = I + R \\ & I = [{E_{ij}}](i = j) \\ & R = [{E_{ij}}](i \ne j) \\ \end{split} \right.$

(2)

同类型转换矩阵体现的是生态系统中的不变面积, 这部分数值较大, 用来体现不同生态类型之间的不变面积数值的对比; 异类型转换矩阵体现的是不同种类生态类型相互之间的转换数量, 生态系统转移矩阵的研究价值主要体现在这一类数据上, 所以异类型转换矩阵是本文研究的重点. 同类型转换矩阵和异类型转换矩阵之间的关联体现在, 特定生态类型变化前的面积 $ {S}_{p} $ 等于该生态类型的不变面积 $ {S}_{t} $ 加上转移走的面积 ${S}_{\rm out}$ , 特定生态类型变化前的面积 $ {S}_{q} $ 等于该生态类型的不变面积 $ {S}_{t} $ 加上转移进的面积 ${S}_{\rm in}$ , 如下所示:

${S_p} = {S_t} + {S_{\rm out}},\;\;{S_q} = {S_t} + {S_{\rm in}}$

(3)

其中, 不变面积 $ {S}_{t} $ 在异类型矩阵中在变化前后保持不变, 且数值较大, 可以将这部分内容从异类型矩阵分析过程中统一排除, 所以将这两部分数据区分开来进行分析不影响转移矩阵的整体研究.

2.2.1 同类型转换

同类型转换矩阵I, 只有对角线上的值不为0, 可以将之转换为一个一维数组 $ {I}_{A} $ , 映射关系为:

${I_A} = \left[ {{\delta _i}} \right](0 \le i < M,{\delta _i} = {I_{ii}})$

(4)

该数组主要体现各个生态类型之间不变数值之间的对比关系和不同生态类型保留面积中的最大值、最小值和均值等特征信息. 所以针对同类型转换矩阵的可视化, 转化为对一维数组 $ {I}_{A} $ 的可视化, 本文采用柱状图的可视化方式, 并附加最大值、最小值和均值的展示即可将该类型矩阵的信息进行较好地可视化.

2.2.2 异类型转换

异类型转换矩阵R, 对角线上的数值都为0, 其他元素数值表示不同类型之间的转换值, 且为不对称矩阵. 以下用转入和转出描述同一生态类型, 从其他类型转移过来的面积和转移到其他类型的面积. 矩阵中包含了两类数据特征: 特定生态类型变化前后的面积总和变化情况和特定生态类型变化前后的源类型(转换前)和目的类型(转换后)的比例.

针对特定生态类型变化前后的面积总和变化情况, 定义R中元素为 $ {r}_{ij}\in R $ , 可以定义按行累加和按列累加的变量, 有如下关系式:

$C = \left\lfloor {\sum\nolimits_{i = 0}^{M - 1} {{r_{ij}}} } \right\rfloor (0 \le j < M,{r_{ij}} \in R)$

(5)

$D = \left\lfloor {\sum\nolimits_{j = 0}^{M - 1} {{r_{ij}}} } \right\rfloor (0 \le i < M,{r_{ij}} \in R)$

(6)

其中, C为按列累加得到的数组, 该数组中元素表示第j类生态类型的转出面积总和; D为按行累加得到的数组, 该数组中元素表示第i类生态类型的转入面积总和. 根据式(3)可知, 同一生态类型, 变化总面积 $ {S}_{c} $ 即等于转入的面积与转出的面积之差.

$ {S_C} = {S_q} - {S_p} = ({S_t} + {S_{\rm in}}) - ({S_t} + {S_{\rm out}}) = {S_{\rm in}} - {S_{\rm out}} $

(7)

设计分组柱状图的方式, 每一组有3个柱状图, 表示一种生态类型的转出面积 $ {S}_{p} $ 、变化总面积S_c和转入面积 $ {S}_{q} $ , 用冷色调表示转出面积, 暖色调表示转入面积, 暗色调表示变化总面积. 通过暗色调柱状图是否在0刻度以上, 可以判断该类型面积是否增加, 以及增加的量. 同时可以对比多个生态类型之间的转入面积和转出面积.

针对特定生态类型变化前后的源类型和目的类型的比例情况, 可以对 ${S}_{\rm out}$ 和 ${S}_{\rm in}$ 分别展开, 动态的分析两者的来源类型和目的类型, 根据式(5)和式(6)分别计算 ${S}_{\rm out}$ 和 ${S}_{\rm in}$ 的来源类型和目的类型的占比情况, 通过双向交互的方式为用户呈现实时的转换比例特征.

2.3 评价模型

通过转移矩阵的性质, 我们了解到, 其本身反映了各种类型之间的转移关系. 本文定义生态转移变化空间为K, 每一种类型的转移都表示变化空间的一个维度 $ {k}_{p} $ , 该维度表示一种覆被类型转移到另外一个覆被类型的概率. 基于此, 我们定义变化空间模型如下:

$\left\{\begin{split} & K = \left\{ {{k_p}} \right\},\;p = 1,2,3,4, \cdots, n \\ & {k_p} = {E_{ij}},i \ne j \end{split} \right.$

(8)

针对某一特定区域, 其评价模型计算步骤如算法1.

算法1. 评价模型计算算法

输入: 区域内两个时刻覆被点阵 $ \scriptstyle{A}_{1} $ 和 $\scriptstyle {A}_{2} $ .

输出: 评价模型K.

1) 读入点阵 $\scriptstyle {A}_{1} $ 和 $\scriptstyle {A}_{2} $ , 统计每一类转移发生的数量, 得到转移矩阵E.

2) 计算得到评价模型.

　　2.1) 转移矩阵中非对角线元素 $\scriptstyle {E}_{ij} $ 顺序赋值给第p个维度 $\scriptstyle {k}_{p} $ .

　　2.2) 对评价模型K进行单位化.

由于用户感兴趣区域的面积不同, 覆被类型也多种多样, 为了消除由于这两个因素导致的数值差异过大, 将评价模型进行单位化, 最终模型向量分布在一个n维空间球上, 其方向表示不同的生态转移类型. 数据评价模型空间的维度数量与覆被类型的数量的平方成正比, 所以该模型是一个高维数据, 需要进行降维处理.

3 可视化设计 3.1 评价模型可视化

通过分析, 我们得到本文定义的评价模型可以用一个双向图来表示, 通过对比分析多种可视化方式, 选定Sankey图(桑吉图)作为模型的可视化基础. 桑吉图有左右两侧两个堆叠柱状图和中间的双向弦图组成, 每个柱状图表示两个时刻不同类型的比例, 弦图的宽度表示该变化所占的比例, 宽度越宽, 比例越大. 该图的优势是既可以全局直观的呈现各个转移覆被类型的比例, 还可以通过鼠标悬浮的方式查看某一类型的出度和入度详情.

为了对类型转化实现更直观的可视化, 本文改进了桑吉图的显示方式, 在右侧柱状图上附加了一个柱状图, 用来表示转移后的覆被类型, 这样可以非常直观的对转移前后的覆被类型做对比可视化, 如图1 所示.

3.2 基于T-SNE的区域评价可视化

对特定区域进行生态变化评价, 需要对多个区域的评价模型进行对比可视化. 但是生态覆被类型繁多, 假设只按照7个大类计算, 维度也高达42个, 为了实现多区域对比可视化, 本文基于T-SNE算法设计了区域评价可视化方法. 表2为一种评价矩阵示例.

数据样本集为多个预设定区域的评价模型集合 $ P:\left\{{p}_{i}\right\} $ , 其中 $ {p}_{0} $ 为用户自定义区域的评价模型, 其他 $ {p}_{i} $ 为预定义评价模型, 是指用户根据自己的需要设定特定覆被类型之间的转移评价模型. 比如, 当用户需要对森林覆盖面积的变化情况进行判定分析时, 可以设计森林转移到城镇用地、森林转移到荒漠、森林转移到农田等维度较高的模型向量和其反向量, 即城镇用地转移到森林、荒漠转移到森林、农田转移到森林等维度较高的模型向量. 模拟生成多组符合用户意图的评价模型向量, 与用户自定义区域模型向量一同进行T-SNE降维运算, 得到二维平面上的分布. 由于覆被类型多样, 某种覆被类型转出的变化向量不一定都集中在一起, 用户自定义区域的向量投影所在的团簇判定采用二维空间距离的方法, 以自定义区域向量投影点为中心, 直线距离最近的预定义评价模型向量点即看做当前用户自定义区域向量的特征相似点. 根据该算法可以通过可视化的方式快速得出用户自定义区域的生态转移情况.

评价算法步骤如算法2.

算法2. 区域评价算法

输入: 自定义区域格点集 $\scriptstyle {A}_{1} $ 、 $\scriptstyle {A}_{2} $ 和评价矩阵M.

输出: 自定义区域生态转移评价.

1) 读入区域格点集 $\scriptstyle {A}_{1} $ , $\scriptstyle {A}_{2} $ 和评价矩阵M.

2) 计算自定义区域评价模型.

2.1) 根据格点集计算转移矩阵.

2.2) 利用转移矩阵计算评价模型并单位化.

3) 计算预定义评价模型.

3.1) 根据评价矩阵随机生成评价模型向量集.

3.2) 对向量集单位化, 并与自定义向量合并.

4) 对合并向量集做T-SNE降维可视化.

5) 根据欧氏距离判定当前点的评价.

通过以上评价算法, 可以实现对用户感兴趣区域的生态覆被类型变化情况的定性评价. 将以上算法流程绘制为流程图如图2所示.

4 用户交互设计 4.1 评价矩阵交互设计

设计评价矩阵的交互方式, 可以方便用户快速制定评价矩阵, 并根据评价矩阵计算预定义模型向量, 矩阵中选中项表示该项是否被作为评价指标, 每一项都被随机赋予一个0.5~1的值, 并对向量做单位化. 如图3所示为重点研究森林变化的评价矩阵, 预定的评价矩阵分类可以辅助队用户自定义区域的矩阵的评价.

4.2 界面交互设计 4.2.1 自定义区域交互

区域分析需要有区域选择方法, 基于GIS地图, 设计了矩阵和多边形两类区域选择方法, 通过鼠标交互在地图上获取区域范围, 然后将该区域范围的生态覆被数据按照数据模型计算得到数据样本点进行后续分析和可视化, 如图4所示.

4.2.2 分类图框选交互

基于T-SNE降维后的向量点, 通过散点图绘制出来, 增加框选的交互方式, 一方面可以快速选中一个团簇, 得到相关团簇的详细信息; 另一方面, 可以与地图进行联动, 选中的点对应在地图上的区域将被高亮显示, 如图5所示.

5 案例分析 5.1 转移矩阵分析

表1所示为2000年到2010年生态系统的转移矩阵数据, 从表中可见, 对角线上的数据是非对角线上的数据的近1000倍. 利用本文方法拆分后得到同类型转换矩阵, 并转化为数组

$\begin{split} {I_A} =& [{\rm{198\;739,189\;132,172\;365\;1,119\;405,}} \\ &{\rm{126\;419,105\;88,752\;210,386\;43}}] \\ \end{split} $

(9)

利用柱状图对该数组进行可视化, 如图6.

从图6可以得出, 2010年相比2000年, 我国的生态系统保留下来的面积中, 草地的面积最大, 城镇的面积最小. 草地和荒漠的保留面积在平均线以上, 其他类型面积在平均线以下.

同理可以得到:

$\left \{ \begin{split} {S_{\rm in}} =& [{\rm{524,1587,6112,5372,12\;818,}} \\ & {\rm{3244,3484,250}}] \\ {S_{\rm out}} =& [{\rm{619,2564,16612,3067,4671,}} \\ &{\rm{223,4811,824}}] \\ {S_C} =& [{\rm{ - 95, - 977, - 10500,2305,8147,}} \\ &{\rm{3021, - 1327, - 574}}] \\ \end{split} \right.$

(10)

利用分组柱状图对该数据进行可视化, 如图7.

从图7可以得到, 2010年相比2000年, 我国的生态系统中, 湿地、农田和城镇面积总体有所增长, 增长量最大的是农田, 城镇转入的面积远大于转出的面积, 说明城镇面积很少转移成其他类型面积, 却较大的从其他类型获取面积. 其他类型总体面积都呈减少趋势, 其中草地的减少量最多, 结合图6可知, 草地也是保留面积最多的类型. 森林的减少量最少, 其转入转出面积基本持平.

对R矩阵在横纵两个方向进行多对多的对应和比例计算, 可以得到不同类型转出面积中新生态类型的比例和转入面积中旧生态类型的比例情况. 通过图8进行交互式可视化.

如图8所示, 左侧为2000年转移走的各个生态类型的面积比例堆叠图, 右侧为2010年转移来的各个生态类型的面积比例堆叠图, 中间部分为左侧转出变换到右侧转入的比例连接线, 线的宽度表示转换的面积值. 其中不同的颜色表示不同的类型, 通过色彩让人直观的感受到不同的生态类型的特征, 比如用黑色表示荒漠, 用草绿色表示草地, 用黄绿色表示农田等. 右侧的堆叠图细分为左右两部分, 左侧部分对应了来源生态类型的颜色, 右侧部分表示了当前生态类型的颜色, 便于区分当前的生态类型种类.

从图8中可得到, 2000年流失的生态类型中, 草地占了一半, 城镇的流失面积最少. 2010年转移来的生态类型中, 农田占到了38%, 转移成其他(冰川或裸地)的较少, 通过两边的对比, 可以对发生变化的生态类型进行全局的对比分析. 鼠标悬浮在左右区域, 都会动态扩展当前鼠标选定的类型的详细比例信息. 以下对几类较典型的分析结果分别阐述.

从图9可见, 2000年发生转移的荒漠面积中, 转为森林的有0%, 有22%转为了农田, 其他两类较多的为湿地和草地, 而且转移的荒漠比例从2000年的14%下降到了2010年的10%, 说明对荒漠的治理初见成效.

5.2 森林变化分析

设定评价矩阵为森林相关系数, 并生成数据样本. 此时用户可以利用交互方式在地图上选择自己感兴趣的区域, 并做可视化分析.

分别对中国中部两个区域进行的选择, 如图9所示, 两个区域的森林变化特征呈现出截然相反的情况, 同时设定评价矩阵为森林模式, 如图3所示.

5.2.1 森林增加

根据图10左侧选定区域计算得到的可视化结果如图11所示. 左侧图是T-SNE降维数据结果, 蓝色点是用户选择点, 红色点是森林增加的点, 色点是森林减少的点, 从图中可以看到, 选中区域属于森林增加点的范围, 所以选中区域的森林变化出现了上涨, 对比右侧的图可以印证森林增加的结论, 同时可以对森林转移的比例进行可视化.

5.2.2 森林减少

如上节中所述同样的判定方法作用的图10右侧选中区域上, 得到了森林减少的判定结果, 对比桑基图验证, 可得到森林转移的比例情况(如图12).

5.3 城镇化分析

在生态分析中, 与森林变化同样重要的是城乡用地变化, 如图13所示, 为框选的上海和京津区域, 并对其进行城镇化生态判定.

图14是图13中左侧选中的上海区域的城镇化判定, 图15是图13右侧选中的京津区域的城镇化判定. 通过两个图左侧的降维可视化效果可以得到两个区域的城乡用地都出现了较明显的增长, 即2000年有大量的土地转换成了2010年的城镇用地.

5.4 模式聚类分析

生态覆被的转移区域是生态研究的重点, 但是没有发生转移的区域往往远远大于发生转移的区域, 通过一般的热度图对差异区域进行可视化, 得到的变化格点非常稀疏, 很难从空间上对生态变化进行可视化. 本文基于以上算法设计模式聚类分析方法. 该方法流程为将格点数据分为等距离的网格, 每一个网格作为一个数据采样, 计算转移矩阵, 并得到单位化的数据评价模型, 如图16.

设定好网格尺度后, 对其进行区域评价可视化(如图17), 右侧的图为降维后的可视化图, 对其中左上侧的团簇进行框选, 可以看到, 他们都集中在我国的东部地区, 而且集中在京津冀、山东和江浙一带, 说明这一带的生态覆被变化是类似的. 同理如图18 所示的团簇集中在我国的北部区域.

图19右侧框选的是中部的团簇, 得到的是东北的4个区域和西南的一个区域, 说明这5个区域在生态覆被变化上具有共性, 可以根据该可视化结果进行进一步的生态分析.

5.5 用户反馈

本文采用调查问卷的方式收集行业领域专家对该系统的使用评价. 评价主要集中在改进的桑基图方法和模式聚类分析方法两部分. “该方法中改进的桑基图非常高效的将生态覆被转移情况展示出来并进行交互, 极大的方便了对生态变化的快速分析”. “模式聚类分析方法提供了一种比较新颖的方式去分析生态覆被类型变化, 希望能对降维后团簇的特征做更深入的判别方法研究. ”

通过以上案例分析, 我们实现了对多个类型生态覆被转移变化进行整体和局部的对比分析; 在整体变化分布图中对某一生态类型变化进行判定; 对具有相同变化趋势的区域进行空间聚类, 研究不同空间的相似性.

6 结论与展望

本文针对当前生态系统数据及其研究中重要的指标性特征生态系统转移矩阵进行详细分析, 在此基础上设计数据模型, 利用降维算法处理后, 设计可视化视图, 实现生态数据的可交互的数据可视分析. 通过案例验证系统在时序分析和空间划分上的可视分析作用, 为评价生态系统变化提供了一种新的技术分析手段.