东北地区归一化植被指数与环境因子间的尺度依存关系1)

李继红 焦裕欣

(东北林业大学，哈尔滨，150040)

自然体系是多维生态要素在不同时空尺度上耦合形成的非线性耗散系统，其不仅具有内生演替周期，还受外部生态梯度效应影响，进而产生复杂、多样的时序特征与空间格局[1]。归一化植被指数(NDVI)富含植被功能与结构等生态、物理信息，是景观分析、水土保持、环境调查、气候演变等全球变化研究的重要监测因子[2]，国内外学者已将归一化植被指数数据成功应用于资源环境监测、生态建模、遥感地学反演等领域中，并关注尺度转换过程中归一化植被指数与时空位置的依赖关系[3]。

关于对归一化植被指数尺度性研究，学者们进行了时空双维度分析。从时间维度上，采用了回归、信息熵、小波分析等工具，研究归一化植被指数的时序变化、物候节律、循环周期等特性[4-6]；从空间维度上，应用了地统计学、相关分析、谱分析、直方图等方法，分析了其与地形、气候、土壤、土地利用等自然人文因素空间关联、尺度耦合规律[7-9]。其中，小波分析具有多分辨率分析的优点，被广泛应用于归一化植被指数动态特征挖掘、时空多尺度格局的研究中。小波分析提供的多尺度分解手段，能够准确识别生态序列空间自相关、生态过程-时空相关性，为研究归一化植被指数时空分布结构性规律与环境相关性提供了良好途径。但是，已有相关研究多集中于省域[2]、小流域[10]、黄土高原[11]等中小尺度单元，缺乏大尺度地理空间上归一化植被指数与生态因子序列空间多尺度耦合的实证，且很少展开对归一化植被指数与地形、气候因子的定量分析。为此，本研究以东北地区为研究区，构建了西南—东北方向、南—北方向的两条生态样带，以归一化植被指数为目标变量，以海拔、坡向、坡度等地形因子和气温、降水等气候因子为解释变量，应用小波分析识别归一化植被指数与环境因子间尺度特征相关性，并进一步研究其空间相关性，运用冗余分析(RDA)厘定各因子对归一化植被指数空间分布的影响，旨在为区域生态资源保护、景观规划和林业管理提供参考。

1 研究区概况

中国东北地区地处亚欧大陆东端中高纬度带，地理坐标115°33′35″～134°58′28″E，38°43′15″～53°33′35″N，面积147万km2，包括黑龙江省、吉林省、辽宁省、内蒙古自治区中东部。研究区域为亚热带季风性气候，季候分明、雨热同期，自南向北分属中温带、寒温带，年平均气温在-1.3～10.5 ℃之间，降水量为250～1 000 m。地形分异明显，区域西部为大兴安岭、内蒙古高原，东北、东南侧为小兴安岭、长白山，松嫩、辽河平原分布于中南部，高程介于0～2 638 m。东北地区是我国最大森林覆盖区，拥有完整的寒温带森林生态系统，极具物种多样性，典型植被为温带落叶阔叶林、针阔混交林。

2 研究方法

2.1 样带设置

为揭示归一化植被指数空间尺度性，以典型样带为直观描述对象。样带设置不仅考虑区域空间主轴，还应反映自然环境差异性[3，5]。为此，构建了南—北向、西南—东北向的样带(见图1)。样带A——西起111°59′56″E、42°38′13″N，东至134°9′36″E、47°8′9″N，纵贯1 750 km，跨越内蒙古高原(浑善达克沙地)、大兴安岭、科尔沁沙地、松嫩平原、小兴安岭、三江平原地区；样带B——北起123°14′38″E、53°31′40″N，南终123°26′42″E、39°47′45″N，延伸1 475 km，纵穿大兴安岭、松嫩平原、辽河平原、辽东丘陵区。东北地区跨度较大，地形复杂，环境多样化，选择样带不仅要考虑区域空间主轴，还应反映自然环境差异性。样带A内蒙古高原-松嫩平原-三江平原，主要考虑海拔因素，样带自西向东涵盖整个东北地区且最能反映海拔变化信息，该样带能显著反映高程变化给归一化植被指数带来的影响。样带B主要取决于纬度因素，纬度决定了气候特征，气候从南到北逐渐变化，其植被长势也逐渐变化，该样带能显著反映气候的环境因素变化给归一化植被指数带来的影响。

2.2 数据来源与处理

本研究中归一化植被指数数据来源于国家航空航天局提供的2000—2017年中分辨率成像光谱仪(MODIS)标准合成产品(见图2)，产品序列号为MOD13Q1，数据版本为V005，时间序列是2017年的植被生长季节(6—9月份)。该产品地面分辨率是250 m的16 d合成图像，呈正弦投影。MOD13Q1版本V005产品提供了每个像素的植被指数值。有两个主要植被层：第一植被层是归一化植被指数，它被称为现有国家海洋与大气管理局的连续性指数；第二植被层是增强型植被指数(EVI)衍生的归一化植被指数，在高生物量区域具有更高的敏感性。另外，还包括来自每日图像的最高质量像素，像元的质量控制文件QA，对于质量控制信息的MODLAND QA字段，我们转换为0(代表最高质量)和3(最低质量)之间的整数等级。这些图像是使用USGS MODIS重投影工具Web界面(MRTWeb)从美国地质调查局网站(http://mrtweb.cr.usgs.gov/)免费下载的。在这个平台上，图像是以子集的形式存储，并使用WGS84参考椭球投影到地理坐标。遥感影像预处理过程见图3，使用软件Envi 5.3、MRT对图像进行投影转换和重采样，得到不同分辨率的归一化植被指数图像；再用与植被指数投影方式一样的东北地区矢量图进行研究区裁剪；依据Matlab2018b使用最大值合成法合成归一化植被指数；在Envi5.3下，使用时间序列滤波软件对归一化植被指数进行时序滤波。最大值复合方法可以消除大部分云、大气、太阳高度角的部分影响和干扰：Imi=maxIij。式中：Imi为16 d周期内第一个i的归一化植被指数最大组合值；Iij为第i天的16 d周期内第j天的归一化植被指数值。此外，还对遥感图像处理软件(ENVI)中提取的图像进行了无缝拼接、投影变换，最终转换为WGS84坐标系统，并制作感兴趣区域(ROI)数据。

尽管对植被指数数据已经使用最大值符合方法进行了预处理，但仍然存在影响数据质量的干扰，例如地理位置误差、角度变化、残留云、大气干扰。为了最大程度地减少这些影响，使用时间序列滤波软件包的卷积平滑算法对数据进行了平滑处理，并为每个粒度采样点生成了5 a归一化植被指数时间序列曲线。中分辨率成像光谱仪(MODIS)测得的像素可用于对时间序列中的每个点进行加权：值为0(良好数据)的权重为最大(1.0)，值为1～2(边际数据，雪/冰)的权重为一半(0.5)，值为3(多云)的权重最小(0.1)。时间序列滤波中使用的函数拟合参数为卷积平滑滤波程序，2个拟合步骤上的4点窗口，适应强度为2.0，无尖峰或幅值截止，季节截止为0，季节开始和结束阈值为20%。时间序列滤波的输出是一组文件，这些文件包含年度季节性参数，长度(季节的长度)、基值(左右最小值的平均值)、峰值时间(季节中拟合函数的最大数据值对应的时间)、峰值(季节中拟合函数的最大数据值)、振幅(季节性振幅)等；包含原始数据的平滑表示的拟合函数文件。生成平均归一化植被指数值的时间序列后，计算了归一化植被指数的统计平均值、标准偏差、最小值、最大值、范围，以该年该序列产品进行最大值合成，以合成之后的归一化植被指数数据为计算标准(见图2)。

图2 东北地区不同位置地形与归一化植被指数

样带空间跨度大、地理与生态单元复杂，为揭示归一化植被指数空间分布与环境因子关系，以地形、气候等2个维度环境要素作为影响因子。其中：地形因素包含高程、坡向、坡度等因子；气候因素主要是气温、降水。高程数据来源于美国宇航局和日本经济产业省推出的航天飞机雷达地形测绘所得高程数据(http://dat-amirror.csdb.cn/dem/)，其空间分辨率约为90 m，坡度、坡向则以高程数据为基础进行提取。其中将坡向值按照下式处理：λ=cos(|180°-θ|)。式中，θ为实际坡度；λ为坡度计算值，其值越大，表明坡度越向阳。多年平均气温和降水数据，由地理国情监测云平台提供(http://www.dsac.cn/DataProduct/Detail/)，其中气温和降水数据时间序列为2000—2017年，其空间分辨率均为1 km。

前述各类栅格数据的空间信息与投影系统存在差异，为便于后续小波相干分析，均转换为兰伯特(Lambert)投影系统，将像元大小重采样为250 m，进而以像元为单位提取归一化植被指数及地形、气候因子属性值。样带A、B分别有10 654、4 162个数据点，包含归一化植被指数与环境因子信息，应用Matlab R2017b计算归一化植被指数与环境因子的多尺度关系，并绘制相干性图谱；应用Rstudio的vegan程序包定量分析环境因子对归一化植被指数空间分布的影响[12]。

图3 归一化植被指数影像预处理

2.3 小波分析原理

小波变换是J. Morlet[13]提出的对具有时空或频率域变量的分析方法，其通过一维伸缩、平移、提取等方法进行不同尺度的细化分析，从中获取相应尺度下结构信息。对于一维变量，小波函数描述为[5-6]：

式中：f(x)为分析小波函数；ψ(x)为小波基函数，ψ(x)的收缩或扩张因子即为尺度参数，用a表示；b则表征小波中心位置。将小波基与待分析函数通过上式运算，便可获取多重尺度性小波系数。

小波方差为变量序列尺度效应的度量，用小波系数模利差平方和表示。一般认为，在一定的尺度窗口下，小波方差越大，则其承载的高频信息越丰富、越能体现变量在此处的细节特征。小波方差定义为：

式中：W(a,b)为信号f(n)在尺度为a、位置为b处的小波变换系数；V(a)为小波方差；n为样本总数[14]。

小波功率谱：表征不同尺度对应的能量密度，其值越大，表明该尺度越强烈。小波功率谱公式如下[16]：

式中：N为级数X的长度。通常也希望提取单个小波尺度的结果，特别是如果小波功率位于有限数量的标度中。

3 结果与分析

3.1 归一化植被指数多尺度空间分布

3.1.1 样带A归一化植被指数多尺度空间分布

小波功率谱能够提取归一化植被指数空间分布的隐含显著尺度，进而显现其空间异质格局(见图4)。由图4可见：样带A存在显著的尺度结构，其中小尺度结构2～4 km在800～1 100 km距离上显著，主要由于此处为松嫩平原地区，该区土地利用结构复杂、地表植被分别连续性不强，因而归一化植被指数在小尺度上表现出较强的异质性。在300～1 500 km距离上(大兴安岭、科尔沁沙地、松嫩平原、小兴安岭、三江平原)，归一化植被指数的尺度结构>32 km，原因是该方向上除了西北内蒙古局部沙地区植被覆盖稀疏、归一化植被指数空间破碎化程度高外，大兴安岭以东地区降水较多、林草生长较好，归一化植被指数空间连续性较好，在大尺度上变异性较小。在500～1 000 km距离上(松嫩平原)、1 100～1 300 km距离上(松嫩平原向小兴安岭延伸地带)存在显著尺度结构，依次为32、40 km，其体现了区域耕地用地类型上归一化植被指数的空间分布一致性。在0～400 km距离上(内蒙古高原)存在64、35 km的尺度结构。综合看，东北地区西南—东北走向上，植被生长在不同区域空间分布特性不一致，以松嫩平原区归一化植被指数尺度较小，说明其空间变异大、分布复杂；该方向上存在4个明显的结构尺度，分别为4～8、8～16、16～50、50～70 km，其在50～70 km尺度上功率谱值达到峰值(为4.3)，表明该方向上归一化植被指数以该尺度结构分布。

3.1.2 样带B归一化植被指数多尺度空间分布

样带B的归一化植被指数在540 km距离(松嫩平原)的尺度为4、8 km，并达到显著水平，表明该局部植被分布连续性差、归一化植被指数空间破碎化。在500～900、600～1 200 km距离上(辽河平原)，归一化植被指数存在32、32～64 km的显著尺度结构，主要反映了平坦地形区农用地类地块归一化植被指数整体连续性好、差异不大。在0～300 km距离上(大兴安岭地区)，归一化植被指数尺度结构较大，为64 km，表明植被在该地段分布良好、空间异质性小。样带B的归一化植被指数功率谱在不同尺度存在多个峰值，具体相关尺度为16～32、32～70、70～160 km，其中尺度约为128 km时达到峰值(为9.6)，表明其是东北地区北—南向归一化植被指数空间分布主要尺度。

对比可知(见图4、图5)，样带A内归一化植被指数的尺度较多而复杂，并存在较小的结构尺度；样带B则小尺度较少，主要特征尺度较大，表明东北地区北—南向植被长势与空间分布差异性小，而西南—东北向异质性大。

3.2 归一化植被指数与环境因子多尺度相干性

3.2.1 样带A归一化植被指数与环境因子相干性

图6为样带A归一化植被指数与环境因子的小波方差和相干性图谱，图谱中线圈内暖色调为通过们Monte Carlo检验的区域。样带不同位置上归一化植被指数随尺度变化表现出不同相干性，为便于统一分析其相干性显著水平，结合小波功率谱显示的尺度周期(见图4)，将样带划分为3个尺度区域：小尺度4～10 km、中尺度15～65 km、大尺度区>200 km，并计算各尺度区域内通过95%水平信度显著检验的面积率。

结果表明：在小尺度上，归一化植被指数与海拔、坡度、坡向、气温、降水量的相干性系数，依次为0.461、0.408、0.378、0.392、0.481，相干显著性通过率为5.87%、4.72%、5.23%、3.95%、7.64%，大部分空间地带上未通过显著性检验，表明在小尺度上归一化植被指数与环境因子关系不显著。在中尺度上，归一化植被指数与海拔、坡度、坡向、气温、降水量的相干性系数，依次为0.625、0.547、0.518、0.438、0.674，其相干显著性通过率逐渐增多，分别达到16.37%、12.51%、9.27%、8.19%。在大尺度上，归一化植被指数与海拔、坡度、坡向、气温、降水量的相干性系数，依次为0.618、0.466、0.627、0.419、0.683，相干性通过率为4.69%、5.24%、21.65%、2.17%、14.68%。由此可见，环境梯度与归一化植被指数的尺度性，在随空间位置、空间尺度大小而异，二者之间关系复杂、多样。

归一化植被指数与海拔在样带300 km以西(即样带内蒙古高原、大兴安岭地区)、1 300 km以东(样带小兴安岭、三江平原地区)与海拔的相干性最突出，与海拔呈正相关关系。这是因为西部内蒙古高原地区沙漠化严重，林草植被主要集中于大兴安岭南部山地；样带东部三江平原地区植被受人为影响较大，导致该地区的归一化植被指数降低，小兴安岭地区植被覆盖度则较好而归一化植被指数增大。坡度是海拔的派生地形因子，其与归一化植被指数的相干性，与海拔相一致。在样带600～1 200 km位置上，归一化植被指数与坡向的相干性最显著，松嫩平原整体地形由北向南倾向，坡向算子值主要介于0～1之间，其与归一化植被指数呈负相关关系；主要由于本区阳坡为迎风坡，水热资源丰富，植被长势相对较好。降水量与归一化植被指数在150 km位置(样带A西端内蒙古高原)、550～800 km位置(大兴安岭向松嫩平原延伸区)相干性显著，样带A西端距离海洋远，降水量是植被生长的限制因子；而大兴安岭东侧是迎风坡，地形抬升汇聚大量水汽，便利了成云致雨，利于植被生长。气温在较大尺度上与归一化植被指数相干性显著，其在宏观尺度上对植被长势起作用。

3.2.2 样带B归一化植被指数与环境因子相干性

根据归一化植被指数自相关尺度(见图5)及其与环境因子相干性图谱位置分布，将尺度划分为5个尺度区域：0～15、15～30、30～80、80～150、>150 km。

在0～15 km尺度上，归一化植被指数与坡度、坡向的相干系数为0.431、0.415，高于其他环境因子；表明在小尺度上，坡度、坡向对植被空间分异的作用大于其他环境因子。

在15～30 km尺度上，海拔等地形因子与归一化植被指数相干性显著，通过率相近，海拔、坡度、坡向通过率分别为12.85%、14.77%、12.95%，降水量对归一化植被指数的影响最小(5.16%)，温度则最大(16.54%)；表明在该尺度上，降水量对归一化植被指数空间分布异质性的作用最大。

在30～80 km尺度上，1 300 km以南(辽东丘陵)地区，归一化植被指数与海拔相干性显著，随着海拔升高，植被覆盖度增加；归一化植被指数与坡度、坡向的相干性，与海拔相似，并在1 000～1 100 km距离上相干性显著；归一化植被指数与降水量，在<200 km(样带B大兴安岭地区)、>1 350 km(辽东丘陵区)位置上相干性较强。该尺度上，归一化植被指数与海拔、坡度、坡向、降水、气温相干性显著，通过率依次为13.36%、17.68%、14.23%、16.63%、12.57%。

在80～150 km尺度上，归一化植被指数与地形因子的相关显著性，主要集中于样带南端辽中南和辽东地区，该区地形由平原向山地丘陵过渡，水热资源受地形梯度分配，引起植被覆盖空间分异，而降水量与归一化植被指数的相干性特征与地形因子相似。该尺度上，海拔、坡度、坡向、降水、气温对归一化植被指数相干性，显著性图谱面积率分别为5.93%、8.72%、8.71%、17.88%、27.86%。

在>150 km尺度上，地形因子、气温对归一化植被指数的相干性，在<600 km的位置上(松嫩平原北部—大兴安岭北部地区)相干性显著，且显著性通过率相差不大，介于8.72%～10.23%之间；而降水量与归一化植被指数相干性，显著通过率达36.54%，表明在该尺度上降水是植被长势的控制因子。

图7 样带B各因子随尺度变化的小波方差和因子间相干系数(主图右侧的渐变图为相关关系)

3.3 归一化植被指数空间分布的影响因素

应用RDA排序分析得到各环境因子对归一化植被指数空间分布的影响(见图8、表1)。对于样带A中归一化植被指数地带性分布，地形、气候因素分别独立贡献了11%、17%；地形-气候组合效应共同贡献了24%；各因素组综合贡献能力分别为气候因素(41%)、地形因素(35%)。从单因子看(见表1)，降水量、海拔的贡献程度最高，分别为33.52%、29.27%，并在0.001水平上达到显著性，表明其是东北地区西南—东北向上植被生长空间分异的主导因素。东北地区该方向上降水量自东向西呈递减分布，样带A依次横穿湿润区、半湿润区、半干旱区，相应地距海洋距离由近及远，其地表覆被由森林向草原、荒漠呈纬向演替。另外，该样带自西向东分属高原、山地、平原等不同地形区，域内地形差异引起地表覆被类型、人为活动强度差异明显。

对于样带B，地形、气候因素，分别独立贡献了13%、21%的归一化植被指数分布信息，二者共同贡献了42%；另有24%的信息未能贡献，是因为与其他因素有关。各因素组综合贡献能力分别为地形因素(55%)、气候因素(63%)。其中，气温、海拔是东北地区南—北方向上植被长势空间分异的重要因素，其贡献了33.67%、30.17%，在0.001水平上显著。样带B区，气温、降水量由南向北减少分布，这与归一化植被指数空间格局相逆(归一化植被指数与气温相关系数为-0.55)，而与高程分布相一致(相关系数为0.68)。海陆位置差异引起水热资源分配不均衡，进而对植被长势产生影响，东北地区纵跨中温带、寒温带，较低纬度地区气候温暖、人为活动强烈，用地类型以建设用地、耕地为主，其植被长势相对于北部高纬地区弱。高海拔地区气温较低，便利了动植物残体的矿化积累，因而土壤肥厚、促进植被生长，加之其地形复杂，人为活动较少，总体看，北部地区归一化植被指数高于南部，并在大小兴安岭北端发育形成中国最大的亚寒带针叶林区。该方向上降水量，从南至北由700 mm逐渐递减为400 mm，归一化植被指数与降水量呈负相关关系(相关系数为-0.36)，贡献了22.5%的分布信息，是归一化植被指数的次要影响因素。东北地区南北地段降水量整体属于半湿润，能够满足植被生长需求，受气温影响，区域蒸发量呈南高北低分布，大气湿度整体差异不大，因而降水量对归一化植被指数空间分异的影响较小。

对比分析可知，样带A跨越的地理单元较多，域内地形、气候、人文环境差异更大，归一化植被指数空间分布的影响因素较多，故5种环境因子对其贡献能力有限，另有48%的贡献率由其他因素引起；与样带B相比，其地理单元相对单调、环境梯度差异较小，因而该环境因子对其贡献能力较好，残差仅为0.24。同时，这也表明，东部地区西南—东北方向上，植被长势相对于南—北方向上变异性大。

图8 不同环境因子对归一化植被指数空间变异影响的方差分解

表1 环境因子重要性排序和显著性检验

4 结论

植被覆盖与长势在空间上呈非线性分布，并具有对空间位置的依赖性，根据小波能量谱可以清晰识别归一化植被指数空间尺度结构。在西南—东北样带上，其特征尺度为4～8、8～16、16～50、50～70 km，其在50～70 km尺度上功率谱值最大，是其主要尺度。在南—北方向上，结构尺度为16～32、32～70、70～160 km，其中70～160 km是其主要尺度。在特征尺度上，植被覆盖与环境变量的结构信息丰富、空间异质性高，西南—东北方向上，归一化植被指数空间特征尺度较多且复杂，其变异性较南—北方向强。

东北地区植被覆盖和归一化植被指数空间分布格局，是多维环境因素综合作用的结果。方差分析表明，地形、气候因子，对西南—东北方向样带归一化植被指数空间分布的综合贡献率依次为35%、41%，对南—北方向样带上归一化植被指数的贡献率为55%、63%；其中，西南—东北方向环境梯度差异大，环境因子总体贡献率为52%，南—北方向上地理单元相对单调，环境因子对其贡献能力达76%。对于西南—东北方向归一化植被指数，各单一环境因子中，降水量、海拔的贡献能力达33.54%、29.27%，其中降水是其空间分异的首要影响因子；对于南—北方向上，气温、海拔、降水的贡献能力较好，依次为33.67%、30.17%、22.25%，其中气温是其空间分布的主导因子。

本研究采用小波相干分析，能够直观揭示归一化植被指数与环境因子随空间尺度、空间位置变化的依赖关系。在小尺度上，地形因子与归一化植被指数相关系数的相干性图谱显著通过率，大于气候因子；上述现象也充分表明，在小尺度上，地形对植被空间分布起着主导作用。大尺度上，气候因子与归一化植被指数的相干系数与显著性通过率高，表明大尺度上气候资源分布是植被生长和归一化植被指数的主控因素。

采用小波分析东北地区归一化植被指数与环境因子的多尺度空间的相关性，能够体现出小波变换在多尺度分解生态时间序列上的方法优势；通过小波分析，对于识别主导的生态因素、预测生态环境变化的发展趋势，具有一定的可行性。但是，在实际的研究中需要进行空间重采样，将研究中的连续空间数据转变为离散的空间序列数据。小波分析技术在实际运用中，核函数的不同也会对最终的分析结果产生影响，因此，在下一步的研究中，将会针对不同的核函数进行选择，优化小波分析的最终结果。

小波相干能够成功分析归一化植被指数与环境因子的尺度依赖性，是空间分析强有用的工具，但是，针对一维样带结果会受到样带的地理位置、分布方向、长度以及采样密度等因素的影响，所以，将小波相干分析延伸应用于二维数据或者影像的整体分析是很好的突破口。同时，对归一化植被指数与环境因子相互作用的多尺度检验，有助于确定环境因子对生态系统影响的有效规模，这对预测生态分布很有帮助。