1. 引言
2021年10月,中共中央、国务院印发《黄河流域生态保护和高质量发展规划纲要》,为有序开展黄河流域生态保护和高质量发展等工作指明了方向。“十四五”规划和2035远景目标纲要中明确提出:“加强长江、黄河等大江大河和重要湖泊湿地生态保护治理,加强重要生态廊道建设和保护” [1] 。上述文件表明,流域生态保护工作已成为当前我国生态环境建设的工作重点。作为流域的重要组成部分,河口三角洲区土质肥沃、人口稠密,常面临着由于过度开发所导致的生态环境恶化、人地矛盾突出等困境,因此妥善协调好河口三角洲区的人地关系,对开展流域生态保护工作具有充分的引领和示范作用 [2] 。
目前,已有对黄河三角洲地区黄河三角洲NDVI变化特征及其影响因素的分析研究,如岳喜元等人对其调查研究,通过对黄河三角洲NDVI值的演化计算,得出目前人类活动对黄河三角洲的自然保护区具有正向影响,由此对其生态环境格局的影响也可以通过NDVI值进行计算 [3] 。刘新旭等人通过对黄河三角洲NDVI进行反演、归因分析,对黄河三角洲湿地蒸发量进行估算以及时空演变规律分析,研究总结了黄河三角洲蒸发量对植被覆盖度的影响 [4] 。
采用人工量测法进行森林三维结构、监测林木长势需要较高的人力物力成本,难以大面积开展。传统的光学遥感影像可实现大尺度无缝监测,但具有易饱和的缺点,只能得到林冠表层的绿度信息,无法穿透树冠,获得林木的三维垂直结构。GEDI激光雷达高度探测器,可测量植被的垂直结构,并提供更精确的植被高度数据;Sentinel-2高分辨率遥感卫星则可以提供多光谱遥感数据,包括红、绿、蓝和近红外波段。融合两者的数据,能够更全面地描述植被生长状态,有助于更全面地了解和管理地球表面的植被分布和生长状态。研究土地利用类型对于保障土地资源的合理利用、加强土地资源保护、推进经济社会发展等具有重要意义,土地作为人类进行生产、生活活动的空间载体,其利用变化映射出人类进行各种活动的过程和结果,研究区域土地利用转型有助于摸清土地利用规律、促进区域合理高效利用土地资源、加快实现区域可持续发展 [5] 。
基于以上研究,本研究采用黄河三角洲2000~2020年植被覆盖度NDVI值,空间分辨率为30 m进行研究黄河三角洲的植物生长趋势。同时结合土地覆盖类型和植被高度分析研究黄河三角洲总体生态质量。
2. 材料与方法
2.1. 研究区
黄河三角洲的定义边界一直都在变化。古代时期,这个地区是由黄河流经而形成的三角洲。1855年开始,黄河改道至大清河入海口后,至1934年黄河分流点从宁海下移至渔洼之前,形成了近代的黄河三角洲。而自1934年以来,黄河三角洲以渔洼为顶点,呈现出了现代的形态。黄河三角洲集中分布在山东省东营市和滨州市境内,由于黄河三角洲93%的面积主要集中在东营市境内,因此本文主要研究东营市内黄河三角洲,为了数据获取的完整性,本文选取东营市的4个行政区——河口区、垦利区、东营区和利津县作为研究整体,即近现代黄河三角洲 [6] ,见图1。
2.2. 数据源与处理方法
本课题所用的植被归一化数据来源于国家科技资源共享服务平台,采取了其2000~2020年中国30米年最大NDVI数据集(http://nesdc.org.cn/sdo/detail?id=60f68d757e28174f0e7d8d49),共21年的NDVI数据 [7] [8] ,本数据集基于Google Earth Engine云计算平台,利用全年所有Landsat 5/7/8遥感数据,通过系列数据预处理和数据平滑等方法,得到2000~2020年每个像元一年中的NDVI最大值。该数据集空间分辨率为30 m,时间分辨率为每年。
地表覆盖分布是研究地球表面不同类型的地形和植被覆盖情况的重要指标,对气候变化、生态环境评估和国土资源和地理国情监测具有重要作用。本文土地栅格利用数据来源全球30米地表覆盖精细分类产品(地球大数据科学工程数据共享服务系统(https://data.casearth.cn/sdo/detail/5fbc7904819aec1ea2dd7061))得到2000~2020年每五年的数据产品 [9] ,通过图谱分析后再通过人工目视解译获取所需相关信息。
植被高度数据来源于https://mp.weixin.qq.com/s/1zXB99wy_5QdnDspe4PcJg,ETH Global Sentinel-2 10米冠层高度(2020年) [10] ,通过融合星载激光雷达GEDI (Global Ecosystem Dynamics Investigation)和Sentinel-2影像数据获取。管理陆地生态系统、缓解气候变化和防止生物多样性的丧失,需要地理空间上的明确信息,而且最好是高度解析的信息。本研究使用了一种全球首个连续逐栅格、地面采样距离为10米的树冠高度图,并取得了突破性的进展。与传统的卫星图像估算方法相比,该方法能够有效消除饱和效应,准确估算高大树冠的高度。这为高碳密度区域的精确测量提供了重要支持。
2.3. 研究方法
利用2000年到2020年最大NDVI数据集,进行一元线性回归和Sen + MK分析,得到黄河三角洲21年间植被覆盖区域NDVI的空间分布特征、时间变化特征和变化趋势特征,为黄河三角洲生态质量时空变化研究提供基础依据,加上20年间每五年的土地栅格数据制成土地利用转移图谱,得到土地利用类型变化规律,结合两者空间重复率定性分析植被绿都在土地利用的基础上的生态质量。再结合2020年黄河三角洲植被高度数据,分析该地区生态系统稳定性。以此综合分析土地覆盖和生态质量的时空变化。
2.3.1. 一元线性回归变化斜率法
回归分析是研究多个变量之间统计联系的一种重要方法,是研究植被长时序变化趋势的重要方法 [11] 。对一组时间自变量x与NDVI因变量y数据,可以用如下的数学模型来描述:
(1)
式中,a、k是未知常数;ε是随机误差。利用观测值(xi, yi) (
)可以求出未知参数k,而k也就是slope斜率。
(2)
对于NDVI长时间连续数据来说,采用一元线性最小二乘法拟合后得到的相应的线性方程式。式中的斜率k值也就是所求的slope (坡度)值,可以说明像元NDVI值的多年变化度的变化趋势 [12] ,主要关联见表1。
![](Images/Table_Tmp.jpg)
Table 1. Relationship between rate of change and vegetation greenness
表1. 变率与植被绿度关系
2.3.2. Theil-Sen Median趋势分析和Mann-Kendall检验方法(Sen + MK趋势分析)
在计算趋势值时,通常会使用Sen斜率估计法,并结合MK非参数检验来判断趋势是否显著。具体而言,先使用Sen斜率估计法计算趋势值,再使用MK方法进行显著性检验。
Theil-Sen Median方法又被称为Sen斜率估计,是一种稳健的非参数统计的趋势计算方法。Sen趋势度是计算序列的中值,虽然本身不能实现趋势显著性判断,但可以很少地避免其他各种因素的干扰,对于测量误差和离群数据不敏感。该方法计算效率高,常被用于长时间序列数据的趋势分析中。
Sen趋势度的计算公式(3)为:
(3)
式中:xi和xj为时间序列数据。Β大于0表示时间序列呈现上升趋势;β小于0表示时间序列呈现下降趋势。
Mann-Kendall检验是一种非参数检验方法,与其他参数检验方法相比,具有更强的适应性。它不需要样本数据遵从特定分布,适用于非线性、非正态分布以及存在异常值的数据。相比于参数检验,Mann-Kendall方法对数据序列分布不要求,且对异常值的影响较小,更适合处理顺序变量。在水文、气象趋势变化相关研究中,Mann-Kendall检验已经得到了广泛应用。例如,用于评估径流、降水、气候变化趋势的显著性 [13] 。
其检验统计量S由公式(4)计算为:
(4)
(5)
当n ≥ 10时,统计量S近似服从标准正态分布,使用定义统计量Z进行趋势检验,Z值由公式(6)计算:
,其中
(6)
式中,n表示时间序列的长度;sign是符号函数;而统计量Z的取值范围一般为(−∞, +∞)。在给定显著性水平p-value下,当|Z| > 1.96时,表示研究序列在p-value水平上存在显著的变化。一般取p-value = 0.05,本文判断在0.05置信水平上黄河三角洲NDVI值时间序列变化趋势的显著性。
本研究使用R语言Raster包进行栅格计算,利用trend包sen.slope函数进行Sen + MK的计算,输出结果包含三个波段,Z值、slope和p值。
2.3.3. 土地利用转移解读
土地利用变化是自然环境因子和其他人类生产活动环境因子之间综合相互作用的必然结果,是影响全球土地气候变化的重要形成动因,已发展成为当前的科学研究发展热点 [14] 。本研究采用了土地利用转移矩阵方法对黄河三角洲的土地利用动态进行监测和分析。这种方法不仅可以精确地描述不同土地利用类型之间的转化情况,还可以量化不同土地综合利用转移类型之间的土地转移速率 [15] 。通过这种方法,我们可以更全面地了解和分析不同土地利用类型之间的变化以及转移速率,为土地管理和规划提供科学依据和数据支持。通过这种方法,我们可以评估黄河三角洲的人工刺槐林生长状况,以实现有效的监测效果。本章从黄河三角洲土地利用的显性形态入手,运用ArcGIS 10.8的Spatial Analyst Tools的重分类工具进行栅格数据唯一值分类,将2000、2005、2010、2015、2020年5个时期的土地利用遥感数据,从数量结构、用途转换空间分布来分析黄河三角洲土地利用覆盖类型的时空特征变化规律,分类后的土地覆盖类型见表2。
随后运用ArcGIS 10.8的Spatial Analyst Tools进行地图代数运算,得到黄河三角洲2000~2020年土地利用转型信息图谱,以实现土地利用格局与演变过程的有效集成 [16] 。所建的土地利用转型图谱是基于已有研究并运用地图代数叠加运算得出的:通过比较两个时期的土地利用类型图谱单元编码,本研究成功地生成了新的图谱编码。具体来说,当土地利用类型数量小于或等于10时,使用了以下的合成公式:
(7)
式中,N表示新生成的图谱编码,A和B分别代表研究初期和末期的土地利用图谱单元编码。这种方法不仅为本研究提供了一种适用于未来研究的方法,而且还可以提高土地利用分析数据的准确性。由此得到2000~2005、2005~2010、2010~2015、2015~2020年四段黄河三角洲的土地利用转型图谱。从宏观角度准确把握黄河三角洲土地覆盖的演变趋势。
![](Images/Table_Tmp.jpg)
Table 2. Classification of major land cover types
表2. 主要土地覆盖类型分类
2.3.4. 植被高度数据解读
GEDI空间仪器提供了高度数据,尽管数据稀疏,但其覆盖范围是前所未有的。相对地,Sentinel-2光学卫星图像观测密集,但不能直接测量垂直结构。这两种技术的优点互补,共同为地表高度和垂直结构提供了更全面的观测和解决方案。通过融合GEDI和Sentinel-2,本数据源开发了一个概率深度学习模型,从地球上任何地方的Sentinel-2图像中检索树冠高度,并对这些估计的不确定性进行量化。这种方法可以为正在进行的森林保护工作服务,并有可能促进气候、碳和生物多样性模型的进步。而本文利用此模型融合的数据来观测分析在黄河三角洲NDVI值变化趋势和土地利用覆盖类型在时间和空间分布上的生态质量环境的变化趋势。
利用ArcGIS 10.8和研究区边框将黄河三角洲的植被高度数据裁剪出,利用con函数将0 m及0 m以下的地区转化为无值区,以便观测。若观测区植被高度较低会影响生态系统的许多方面。一方面,它可能会导致土地侵蚀,因为没有足够的根系来绑定土壤,使其易受风蚀和水蚀的影响。此外,植被高度低还会影响地球的水循环。由于没有足够的蒸发量,水分很少从土壤中蒸发到大气中,这可能会导致干燥地区的一系列问题,如干旱和低水位等。另一方面,植被高度低也会影响生物多样性和食物链。许多动物需要一定高度的植被来隐藏和寻找食物,如果植被高度太低,则这些生物的生存将受到严重威胁。因此,维持适当的植被高度对于生态系统的健康至关重要。
3. 结果与分析
3.1. 植被绿度趋势分析
本文采用了Theil-Sen Median方法和Mann-Kendall (MK)方法两种技术相互协作,来分析计算黄河三角洲21年间的植被归一化值趋势变化特征。我们利用R语言出图结果中的Z、slope和p-value趋势来研究植被长势,并进一步分析生态质量的时空变化。
3.1.1. Slope趋势分析
根据黄河三角洲21年数据综合分析可以得出图2黄河三角洲的坡度连续值,为方便人工目译观测分析,故将栅格图谱分为5段。见图2,可知2000~2020年黄河三角洲大部分地区slope值高于0,植被长势良好,北部和东部靠海地区有较大部分slope较低,植被有退化迹象。结合本文随后的可信度分析,黄河三角洲地区植被长势良好地区多为耕地,也就是人工农田生态区;滩地水体以及沼泽的植被长势趋于平缓或者退化;植被长势良好的非人工生态地区集中分布于利津区东北部、河口区西南部以及垦利的中部区域。而上述地区与|Z| > 1.96地区的重合度较高,也就是数据可信度较高。
![](//html.hanspub.org/file/15-2420562x17_hanspub.png?20230705100255132)
Figure 2. Analysis of Interannual variation trend of NDVI
图2. NDVI年际变化趋势分析
3.1.2. Z值图译分析
在Sen + MK趋势分析中,我们可以通过观察Z值的绝对值大小来进行显著性检验。当|Z|分别大于1.65、1.96和2.58时,通常表示趋势通过了90%、95%和99%的显著性检验。这种检验方法可以帮助我们判断分析结果的可靠性,并用于评估数据趋势的显著性水平。本文取Z值1.96为分界线,Z的绝对值大于1.96的部分为高可信度区域,以此来进行分类。从Z值结果图中可以看到,在垦利区中部,东营区北部,利津区与河口区交接部分,趋势分析的可信度较高。此外,在黄河口生态旅游区和广南水库等地区所计算出来的可信度极高。由此,可以推断此部分区域的植被趋势分析的准确度较高。同时,研究区内可信度高的区域一般与农田、草本覆盖、水塘等地区重叠度高,而人工建筑等建设用地的地区与可信度低的地区重叠度较高,但灌溉农田附近地域可信度也较低。具体分布见图3。
![](//html.hanspub.org/file/15-2420562x18_hanspub.png?20230705100255132)
Figure 3. Significance distribution based on Z-value
图3. 基于Z值的显著性分布
3.1.3. 显著性分析
见图4,从显著性结果分析中可以看出,p值小于0.05的区域与上图中可信度高的区域有多个地区重合,说明此区域植被变化度并不显著,也就是植被长势生态质量改善并不高。植被变化度大的区域在黄河三角洲靠近入海口处地区集中分布,内地零星分布,分布地与农田和植被长势较缓等地有多地重合。可以初步推测近20年黄河三角洲生态环境质量的时空变化,其环境质量并没有得到很好的保障并且有退化的趋势,因此科学治理黄河三角洲生态环境已是一个迫在眉睫的问题。
![](//html.hanspub.org/file/15-2420562x19_hanspub.png?20230705100255132)
Figure 4. Significance distribution based on p-value
图4. 基于p值的显著性分布
3.2. 土地利用转变
黄河三角洲土地利用时空变化分析
根据中国科学院资源环境与科学数据中心提供的山东省土地利用遥感解译数据,基于ArcGIS软件叠加相关行政界线,对所得土地利用栅格数据进行重分类,并依据数据手册对其11项分类进行合并,利用唯一值分类后结果见表3。将“旱地农业”和“灌溉农业”合并得到“耕地”,“草本覆盖”和“草地”合并得到“草地”,“落叶阔叶林”与“常绿针叶林”合并得到“林地”,“海洋”、“沼泽湿地”与“水体”合并为“水域”,“建设用地”保留,“裸地”与“稀疏植被”和“未利用地”合并为“未利用地”六项,制成黄河三角洲2020年土地利用现状图。
为更细致地从总体观察土地覆盖时空变化,保留了2020年栅格数据唯一值分类出来的所有土地利用类型,结果见图5。各土地利用类型中,草本覆盖面积最大,广泛分布在垦利、利津以及其它非沿海地区,但因为NDVI草地面积不到其1/10,在后续分析故将其合并分析;耕地以灌溉农田为主,主要分布在河道附近,其次为旱地农业,远离河道分布,本文为方便观测,将旱地农业与灌溉农业合并为耕地,其总面积为:水域面积包括了主要水体和沼泽湿地部分,主要为水库坑塘和滩涂用地,集中分布在沿海地区;城乡建设用地面积居第五位,零星分布于内陆地区,主要集中在东营与垦利交界地带;未利用地主要为裸地和稀疏植被,面积较小;林地面积较小,只有为落叶阔叶林的5,并且相较于往年,常绿针叶阔叶林完全消减至零。而黄河三角洲的森林资源以灌丛为主,其中包括河滩防护林、海滩林、淡水湿地灌丛等类型。地区内的湿地区域则主要以杨树、楸树等为主,而平原地区则以榆树等阔叶树种为主。乔木林带则主要种植了松树、柏树等树种。因为数据产品只有落叶阔叶林和常绿针叶林,所以大部分灌丛森林被分类到了草本覆盖地区,在本文后会利用植被高度数据进行观测。
![](Images/Table_Tmp.jpg)
Table 3. Land use and cover types in the Yellow River Delta in 2020
表3. 黄河三角洲2020年土地利用覆盖类型表
![](//html.hanspub.org/file/15-2420562x20_hanspub.png?20230705100255132)
Figure 5. All land use types in the Yellow River Delta in 2020
图5. 黄河三角洲2020年所有土地利用类型
黄河三角洲地区2000~2020年土地利用时空分布格局见图6,土地利用面积构成矩阵转移见图7所示。2000~2020年间,黄河三角洲地区的草地逐年呈下降趋势,从2000年的3809 km2到2020年仅余3108 km2,林地也是仅有的第二个逐年减少的土地类型。耕地面积在原有的基础上逐年增多,主要集中在黄河口、黄河故道等河道附近;垦利区东部靠海处也有增加。
![](//html.hanspub.org/file/15-2420562x21_hanspub.png?20230705100255132)
Figure 6. Land cover types from 2000 to 2020
图6. 2000年~2020年土地覆盖类型
总体上,黄河三角洲地区土地分布有一定规律性,主要类型还是以草本植物覆盖为主,水体多分布在河道和靠海地区,其余用地零星分布其中。耕地大多集中分布在河道离水源近的地区或者内陆平原地区。建筑用地集中分布在东营区和垦利区交界处。因为数据来源30 m土地栅格数据中不包含大部分当地林种的数据,所以本研究中所涉及到的林地不是指研究区图1的人工刺槐林一类,而是黄河三角洲南部零星接触到的阔叶林和针叶林。数据源中的产品并未包括乔木等高大树种的生态群落,所以才划分的草地分类中包含了灌丛木一类的树木林地。
在图7图例编号中,12代表耕地转变为草地,18代表耕地转变为水体,19则代表耕地转变为建设用地;而21表示草地转变为耕地,81则表示水体转变为耕地,其余编号同理。通过对比四个时间段的黄河三角洲土地利用类型转变情况,本研究发现,在20年的时间里,黄河三角洲草地转变为水体的速率最快。同时,水体的扩张速度也相当快,表明黄河三角洲的水流体系在不断改变。而耕地的面积也在不断扩大,侵占了草地和林地,高大乔木的混交林从本就不多的状态降至更少,当然这也与本地适合灌丛木等低矮树木生长有关。这些转变都反映了该地区的生态环境质量需要进一步改善。随着草地的减少和水域面积的扩大,我们需要采取更加有效的措施来保护该地区的环境健康,例如退耕还林还草等一系列措施,以及加大力度建设和拓展人工刺槐林面积,不仅可以提高生态坏境整体质量,还可以有效减缓当地的土地盐碱化趋势,保持入海口地区的土地环境质量。
![](//html.hanspub.org/file/15-2420562x22_hanspub.png?20230705100255132)
Figure 7. Land use change from 2000 to 2020
图7. 2000年~2020年土地利用转移
3.3. 植被高度分析
该数据其空间分辨率为10米*10米,使得观察该区域更加便捷且相对于其他数据的精准度更高,具体分布见图8。该区域大部分地区都有植被覆盖,但是植被高度较低,土地主要由草本和农田覆盖类型组成,以灌丛为主导,其次是阔叶树林和针阔混交林,只有少部分地区存在高植被分布。观察该区域的图像可以发现,黄河三角洲沿河道部分以及靠近内陆地区多为高大的树木,而靠近海岸和滩地沼泽附近地区则树种较为矮小,大部分地区的植被高度都较低。观察该区域的图像可以发现,高植被覆盖地程零星块状分布于行政区边界附近,与上文植被绿度变化且可信度高的地区重合度较低,不同数据之间有差值;植被高于均值的地区总数不过半,结合土地利用类型可以推测,大多数为城市道路旁植树林,且大面积人工林数量占比少。而植被高度较低代表了垂直生态系统成分较为简单,易于受到外界干扰。因此,从分析来看,黄河三角洲地区的生态系统稳定性较低,容易受到外界因素干扰和破坏,因此迫切需要加强人工对生态系统的干预和保护。由于该数据没有多年连续观测值,无法获取高度变化信息,未来研究中可通过多年高度反演值,分析高度变化情况。
![](//html.hanspub.org/file/15-2420562x23_hanspub.png?20230705100255132)
Figure 8. Spatial distribution map of vegetation height
图8. 植被高度空间分布图
4. 结论与讨论
4.1. 结论
本研究基于空间分析功能和可视化技术合成地学信息图谱,研究了2000~2020年黄河三角洲地区土地用地时空格局动态变化特征,分析植被NDVI数据Sen + MK斜率变化趋势,综合评价黄河三角洲地区生态环境质量变化;利用土地利用转移矩阵得到不同的类型转变,进而得到黄河三角洲地区土地覆盖的时空变化特征,以此来总结土地覆盖与生态质量时空变化趋势的影响因素。通过对以上数据的综合分析,发现土地覆盖、生态环境和人类活动之间的复杂关系。
研究发现,黄河三角洲地区的植被覆盖度总体上呈增长趋势,其中农业用地和城市用地对植被覆盖的影响较大,而森林和草地则有助于保护和促进植被的生长。此外,由于城市化的快速发展,建设用地近年来呈增加趋势,导致植被覆盖度下降。另外,与土地利用相比,植被高度的观测数据能够更加直观地反映出植被的生长情况,并且与土地利用数据相结合能够进一步揭示出土地利用和植被生长之间的关系。然而,测量植被高度是一项复杂而昂贵的任务,需要大量的人力和物力支持。
在综合分析以上数据的基础上,笔者认为,加强生态保护和推进生态文明建设是黄河三角洲地区可持续发展的重要措施。例如,合理规划和利用土地资源、加强对植被生长的监测和保护、改善土地生态环境等措施可以有助于促进该地区的可持续发展。
4.2. 讨论
笔者发现,本研究存在以下不足:1) 未全面考虑气候变化、人类活动等因素对植被覆盖的影响;2) 使用的数据可能存在误差;3) 未进行相关系数分析;4) 研究仅局限于部分地区,未考虑整体区域特征。未来的研究需要加强上述方面的完善,应该探寻经济与生态的协同发展模式,制定出符合当地特点的科学合理的土地利用政策、环保政策,并积极倡导绿色发展理念,以实现和谐共生的经济生态发展,为保护黄河三角洲地区的环境提供更可靠的科学依据。
基金项目
山东省高等学校省级大学生创新创业训练计划(S202110446141)。