基于Landsat时序数据的1995~2024年深圳湾红树林空间分布与NDVI时空变化特征研究

Characterization of Spatial Distribution and Spatial-Temporal Variation of NDVI of Mangrove Forests in Shenzhen Bay from 1995 to 2024 Based on Landsat Time Series Data

摘要: 本文使用Mann-Kendall + Theil-Sen联合分析了深圳湾红树林从1995至2024年的NDVI时空变化，并使用Pearson分析了红树林NDVI的变化与海平面和海表面温度之间的关系。使用Hurst指数联合Theil-Sen分析对红树林NDVI的未来趋势进行了预测。研究表明1995~2024年中国红树林NDVI空间平均值上升了0.21，线性变化斜率为0.0073。NDVI与海平面高度的Pearson相关系数为0.201，呈现正相关，关系不显著。NDVI与海表面温度Pearson相关系数为0.146。预计在未来深圳湾红树林NDVI值呈上升趋势的面积为658.8 ha，占比89.39%，呈下降趋势的面积为42.93 ha，占比5.82%。

Abstract: In this paper, spatial and temporal changes in NDVI of mangrove forests in Shenzhen Bay from 1995 to 2024 were analyzed using Mann-Kendall + Theil-Sen joint analysis, and the relationship between NDVI changes and sea-level and sea-surface temperatures was analyzed using Pearson analysis. The future trend of mangrove NDVI was predicted using the Hurst index combined with Theil-Sen analysis. The study showed that the spatial mean of mangrove NDVI in China increased by 0.21 from 1995 to 2024, with a linear trend slope of 0.0073. The Pearson correlation coefficient between NDVI and sea level height was 0.201, and 0.146 between NDVI and sea surface temperature—both positive but insignificant. It is projected that NDVI values in the mangrove forests of Shenzhen Bay will continue to increase, with 658.8 ha (89.39%) of the area showing an upward trend and only 42.93 ha (5.82%) showing a downward trend.

1. 引言

红树林是一种分布于热带与亚热带滨海陆地和海洋环境之间动态界面上的特殊植物群落，在海岸保护、养分循环以及维持生物多样性等方面发挥着重要作用。此外红树林生态系统因其卓越的碳封存能力而得到认可，其单位面积的碳封存量是热带雨林的四倍，是全球生态过程和气候变化调节计划中的关键部分。但近几十年来由于经济快速发展、城市化、全球气候变化和海平面上升等因素，红树林生长活动及其生境经历着剧烈变化。据统计，全球范围内有35%的红树林已经消失，目前还在以每年1%~2%的速率减少，67%的红树林遭受了不可逆转的破坏。监测红树林生长活动的空间模式和演化特征，对海岸带生态保护、海洋资源管理和可持续发展具有重要意义。

遥感技术因其广域空间尺度、信息获取速度快、重返周期短的特点，极大提升了监测红树林生长活性的能力。随着卫星遥感技术发展，多个开源遥感数据集如Landsat系列，MODIS，SPOT和Sentinel-2等，为捕捉全球或区域的红树林在年、季、月尺度上动态变化提供了可能。相比于其他卫星遥感数据，Landsat系列(1972~至今)具有更长时间序列连续性观测优势兼具中低空间分辨率优势(30 m, 15 m)，而更适于捕捉生境复杂地区红树林长期动态变化。

2. 材料与方法

2.1. 研究区域

深圳湾(中国内地以外称后海湾)，地理位置为北纬22˚28.6'，东经113˚57.5'，为中华人民共和国香港特别行政区和广东省深圳市之间的一个海湾，位于珠江口伶仃洋东侧中部，湾口朝向西南[1]。深圳湾的东北部为“福田国家级自然保护区”，东南部为“香港米埔自然保护区”。气候属于亚热带季风气候，雨热同期。共有“深圳河”、“元郎河”以及“大沙河”三条河流流入。本文的研究区域来源于贾明明等创建的10 m分辨率全球红树林数据集(HGMF_2020)，筛选出深圳湾红树林区域，如图1所示。

Figure 1. Study area

图1. 研究区域

2.2. 数据来源与处理

2.2.1. Landsat系列卫星数据

本文通过谷歌地球引擎(Google Earth Engine)获取美国航空航天局(NASA)和美国地质调查局(USGS)共同管理的Landsat系列卫星数据，并在GEE上完成各种指数运算。由于Landsat不同系列的卫星有不同的寿命，以及每一代卫星传感器有所提升，因此本研究采用Landsat5、Landsat7以及Landsat8这三个系列卫星数据，其空间分辨率为30m。其中1995年至2011年使用的是Landsat5 TM，2012、2013使用的是Landsat7 ETM+，2014年至2024年使用的是Landsat8 OLI。

2.2.2. 平均海平面高度数据

海平面高度数据来源于国家科技资源共享服务平台——国家海洋科学数据中心(https://mds.nmdis.org.cn/)的“水位综合数据集”。该数据涵盖了国内外多个大型合作项目和计划，包括全球海平面观测系统(GLOSS)、IOC/IODE (政府间海洋学委员会/国际海洋资料与信息交换委员会)、东北亚海洋观测系统(NEAR-GOOS)以及香港天文台等发布全球海域范围内的水位资料。海平面高度的单位为1厘米(cm)，时间分辨率为1 h，数据格式为xls表格。使用均值合成的方法获取平均海平面高度。

2.2.3. 平均海表温度数据

海表面高度数据来源于IAPv4海洋温度和海洋热含量网格化数据集，数据单位为1摄氏度(℃)。使用ArcGis Pro中的Extract by Mask工具，选中距离本文研究区最近的海表温度数据，使用均值合成的方法获取平均海表面温度数据。

2.3. 研究方法

2.3.1. Theil Sen Median Slope + Mann Kendall趋势分析

Theil Sen中值斜率是经过计算时间序列所有斜率的中值，它可以很好地减少噪声的干扰，但其本身不能实现序列趋势显著性判断。Mann-Kendall检验不需要时间序列服从标准正态分布，且对异常值不敏感。本文研究区属于亚热带季风气候，获取的1995~2024年Landsat时序数据会不可避免地受到不稳定大气的影响，因此引入该方法可以完成对时间序列趋势的显著性判断[2]。

Sen中值斜率的计算方法为：

$\text{Se<msup> n <mo>′</mo> </msup> Slope}=\text{Median}\left(\frac{y_j-y_i}{j-i}\right)\text{for}j>i$ (1)

式中，Median是一个统计函数，表示取序列中的中值，$y_j$ 和$y_i$ 分别表示第j年和第i年的NDVI值。当Sen’Slope > 0时序列呈上升趋势，当Sen’Slope < 0时序列成下降趋势。

Mann-Kendall检验的计算方法为：

$S={\displaystyle {\sum }_i^{n-1}{\displaystyle {\sum }_j^n\mathrm{sgn}\left(y_j-y_i\right)}}$ (2)

$\mathrm{sgn}\left(y_j-y_i\right)=\{\begin{array}{ll}1,\hfill & y_j-y_i>0\hfill \\ 0,\hfill & y_j-y_i=0\hfill \\ -1,\hfill & y_j-y_i<0\hfill \end{array}$ (3)

$Z=\{\begin{array}{ll}\frac{S-1}{\sqrt{\frac{n\left(n-1\right)\left(2n+5\right)}{18}}},\hfill & S>0\hfill \\ 0,\hfill & S=0\hfill \\ \frac{S+1}{\sqrt{\frac{n\left(n-1\right)\left(2n+5\right)}{18}}},\hfill & S<0\hfill \end{array}$ (4)

在本研究中选取显著水平$\alpha =0.05$ ，临界值$Z_{1-\alpha /2}$ 对应标准正态分布表的数值为1.96，当$\left|Z\right|\ge Z_{1-\alpha /2}$ 时认为趋势显著。

将Sen与Mann-Kendall结合，分级依据如表1所示。可获得深圳湾红树林NDVI在1995~2024的30年变化趋势。

MK突变检验：对30年NDVI时间序列$y_1,y_2,y_3,\cdots ,y_n$ ，构造一秩序列：

$S_k={\displaystyle {\sum }_{i=1}^k{R}_i},k=2,3,4,\cdots ,n$ (5)

式中，R_i表示$y_i>y_j$ 的累计数，定义统计量：

$U{F}_k=\frac{S_k-E\left(S_k\right)}{\sqrt{Var\left(S_k\right)}}$ (6)

式中，当k = 1时，UF₁ = 0，$E\left(S_k\right)$ ，$Var\left(S_k\right)$ 表示S_k的均值和方差，将时间序列倒序，重复以上步骤可以得到UB_k。分析UB_k和UF_k的曲线，两条曲线的交点即为突变点。

Table 1. Significance grading of Sen + MK trends

表1. Sen + MK趋势显著性分级

显著性水平α	Sen’Slope	$\left|Z\right|$	变化趋势
α = 0.05	Sen’Slope > 0.0005	$\left|Z\right|>1.96$	显著上升
	Sen’Slope > 0.0005	$\left|Z\right|\le 1.96$	不显著上升
	−0.0005 < Sen’Slope < 0.0005		无变化
	Sen’Slope < −0.0005	$\left|Z\right|\le 1.96$	不显著下降
	Sen’Slope < −0.0005	$\left|Z\right|>1.96$	显著下降

2.3.2. Hurst指数分析

Hurst指数由英国水文专家Hurst提出，是用以衡量未来时间序列数据变化趋势为持续性或反转性的有效方法，在水文、地质等方面具有广泛应用[3]。当0 < H < 0.5时，未来的趋势与过去的趋势呈负相关，意味着如果时间序列在过去呈现上升趋势，那么在未来下降的可能性较大；当H = 0.5时，时间序列是随机的，即未来的趋势与过去的趋势不存在长期相关性；当0.5 < H < 1时，则未来的趋势与过去的趋势保持一致，且H值越大，正持续性越强[4]。本研究采用重标极差分析法(Rescaled Range Analysis, R/S)来计算Hurst指数。

Hurst指数计算公式：

将总长度为N的时间序列划分为A个长度为n的连续子区间$I_a\left(a=1,\cdots ,A\right)$ ，$I_a$ 中的每一点记作$R_{k,a}$ ；

针对不同的子区间长度n，计算A个子区间的平均重标极差：

${\left(R/S\right)}_n=\frac{1}{A}{\displaystyle {\sum }_{k=1}^A{\left(R/S\right)}_k}$ (7)

对每个子区间上的重标极差进行以下计算：

对每个长度为n的子区间，在对相应时间序列作零均质化处理后，计算其累计离差：

$X_{N,a}={\displaystyle {\sum }_{k=1}^N\left(R_{k,a}-\frac{1}{N}{\displaystyle {\sum }_{k=1}^n{R}_{k,a}}\right)}$ (8)

定义单个子区间上的极差：

$R_a=\max \left(X_{k,a}\right)-\min \left(X_{k,a}\right),k=1,2,\cdots ,n$ (9)

计算各个子区间上的重标极差值

${\left(R/S\right)}_a=R_a/S_a$ (10)

其中，

$S_i=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{\left(R_{k,i}-\frac{1}{N}{\displaystyle {\sum }_{k=1}^n{R}_{k,a}}\right)}^2}},i=1,2,\cdots ,A$ (11)

由于样本的平均重标极差值与样本长度之间存在标度关系：

${\left(R/S\right)}_n=c\times n^H$ (12)

因此，对不同的时间尺度重复以上步骤，并将所得的平均重标极差值${\left(R/S\right)}_n$ 对n进行双对数回归：

$\log {\left(R/S\right)}_n=\log \left(C\right)+H\cdot \log \left(n\right)$ (13)

本文进一步将Hurst指数结果与Sen趋势分析结果进行联合分析，分析结果如表2。

Table 2. Division of mangrove growth trends

表2. 红树林生长变化趋势划分

Hurst指数	Sen趋势	变化趋势	描述
0 ≤ H < 0.5	Sen < 0	反持续性降低	未来变化趋势与过去持续降低的趋势相反
	Sen > 0	反持续性升高	未来变化趋势与过去持续升高的趋势相反
0.5 < H ≤ 1	Sen < 0	持续性降低	未来呈持续降低趋势
	Sen > 0	持续性升高	未来呈持续升高趋势
H = 0.5		无持续特征	无法预测未来趋势

2.3.3. Pearson相关性分析

本文采用Pearson相关关系，计算1995~2024年红树林NDVI与平均海平面高度和平均海表温度之间的关系，计算公式如下：

$r=\frac{{\displaystyle {\sum }_{i=1}^n\left(x_i-x^{\prime }\right)\left(y_i-y^{\prime }\right)}}{\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(x_i-x^{\prime }\right)}^2}}\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(y_i-y^{\prime }\right)}^2}}}$ (14)

式中，n为NDVI的时间序列长度，$x_i$ 为第i年NDVI的值，$y_i$ 为第i年的平均海平面高度或海表温度，x′表示多年平均NDVI，y′表示多年平均海平面高度或海表温度。相关系数绝对值越大，变量之间线性相关程度越高。

2.3.4. Random Forest分类

本研究在GEE上采用随机森林(Random Forest) [5]模型对深圳湾1995~2024年的红树林范围进行提取，随机森林的优点使得它能够处理复杂的遥感图像分类[6]。模型输入特征有：SR_B1、SR_B2、SR_B3、SR_B4、SR_B5、NDVI、NDWI、EVI、IBI，且使用QA_PIXEL波段去除含云量大于3%的影像，使用平均值函数对影像进行合成。训练样本数据来自Landsat5的真彩色合成影像，通过目视解译的方式进行标注，并且采用二分类的方式，仅生成红树林与非红树林的标注，结合红树林数据集(HGMF_2020)，随机森林分类的总体精度与Kappa系数均能达到1的效果。

3. 结果

3.1. NDVI时空变化趋势特征

3.1.1. NDVI总体变化趋势

由图2可知，在过去30年中，深圳湾红树林的NDVI值整体呈波动上升的趋势，且增速明显，总体上升值为0.21，线性趋势斜率为0.0073，R² = 0.993。UF与UB两条线的交点出现在临界线外(±1.96)，表明1995~2024年间深圳湾红树林NDVI变化不存在突变，属于稳定变化。

将NDVI划分为低(0~0.25)、中(0.25~0.5)、高(0.5~0.75)三个等级，由图3发现深圳湾红树林低等级NDVI在经历一年的快速上升后，开始逐渐下降。中等级NDVI呈现先上升后下降的变化趋势。高等级NDVI则是先经历波动上升后，变化为稳定上升趋势。

Figure 2. Annual mean NDVI values and MK mutation test

图2. 年均NDVI值与MK突变检验

Figure 3. Changes in NDVI area for the three classes

图3. 三个等级的NDVI面积变化

3.1.2. NDVI全局变化趋势

由表3与图4可知，在1995~2024年间深圳湾红树林NDVI总体呈显著上升趋势，面积为757.35 ha，占比95%，不显著上升面积为22.86 ha，占比2.86，无明显变化面积为0.54 ha，占比0.06%，不显著下降面积为12.06 ha，占比1.5%，显著下降面积为4.32 ha，占比0.5。合计NDVI呈增加趋势的面积为780.21 ha，NDVI呈下降趋势的面积为16.38 ha，分别占97.86%与2%。

3.1.3. NDVI变化的持续性特征

由表4与图5可知，在1995~2024年间深圳湾红树林NDVI值持续变化面积为665.88 ha，占比90.76%，反持续面积为32.85 ha，占比4.45%，无持续性面积为35.19 ha，占比4.77%。

Table 3. Percentage of NDVI classification trend

表3. NDVI分级趋势占比

分级类型	面积/ha	面积占比/%
显著上升	757.35	95
不显著上升	22.86	2.86
无明显变化	0.54	0.06
不显著下降	12.06	1.5
显著下降	4.32	0.5

Figure 4. NDVI trend test of mangrove forests in Shenzhen Bay

图4. 深圳湾红树林NDVI趋势检验

Table 4. NDVI Hurst persistence classification

表4. NDVI Hurst持续性分级

分级类型	面积/ha	面积占比/%
持续性升高	658.17	89.31
反持续性升高	32.22	4.37
持续性降低	10.71	1.45
反持续性降低	0.63	0.08
无持续性特征	35.19	4.77

通过Hurst指数结合Sen中值斜率分析，可知在未来，深圳湾红树林NDVI值呈上升趋势的面积为658.8 ha，占比89.39%，呈下降趋势的面积为42.93 ha，占比5.82%。

3.2. NDVI相关性分析

将深圳湾红树林多年平均NDVI值，与深圳湾多年平均海平面和多年平均海表温度分别进行Pearson相关性分析。结果显示NDVI与平均海平面高度的Pearson相关系数r = 0.201，NDVI与平均海表温度的Pearson相关系数r = 0.146。参照Pearson相关性表格，由图6可知，深圳湾红树林多年平均NDVI值与平均海平面高度呈弱相关，与平均海表温度呈极弱相关或不相关。

3.3. 深圳湾红树林面积变化

根据随机森林提取红树林的结果，由图7与图8可知，深圳湾红树林面积在1995~2024这30年间增加了478.89 ha，主要增加区域为香港米埔自然保护区的深圳河入湾口处。

4. 讨论

4.1. 深圳湾红树林的30年变化趋势

虽然深圳湾处于深圳与香港两个超一线城市之间，但深圳福田自然保护区与香港米埔自然保护区分别于1986和1976年建立，在两大保护区的共同努力下，1995~2024年深圳湾红树林面积与NDVI值均呈上升趋势，并且高等级NDVI面积变化量达177%，低等级NDVI面积变化量达−280%。NDVI与净初级生产力(NPP)、冠层覆盖面积、生物量等具有明显相关性，通过NDVI的变化可以间接反映出，深圳湾红树林的生态环境质量正在逐渐提高。

4.2. 深圳湾红树林时空变化的影响因子

在本研究中，深圳湾红树林NDVI值与平均海平面高度呈0.201的低相关性，与平均海表面温度呈0.146的极低相关性。在褚艳玲等人的研究中，发现在海平面升高0.55 m的情景下，深圳湾红树林将面积增加22%，而在海平面升高0.87 m的情景下，深圳湾红树林面积将会减少22.4%。这表明红树林与海平面高度之间并非简单的线性关系[7]。

Figure 5. Analysis of NDVI Hurst index of mangrove forests in Shenzhen Bay

图5. 深圳湾红树林NDVI Hurst指数分析

Figure 6. Multi-year mean sea level and sea surface temperature changes in Shenzhen Bay

图6. 深圳湾多年平均海平面与海表温度变化

Figure 7. Changes in the area of mangrove forests in Shenzhen Bay between 1995 and 2024

图7. 深圳湾红树林1995~2024年间面积变化

Figure 8. Area range of mangrove forests in Shenzhen Bay in 1995, 2009 and 2024

图8. 深圳湾红树林1995、2009、2024年面积范围

5. 结论

本文研究结果表明：

1995~2024年，深圳湾红树林年均NDVI值在整体上呈波动上升趋势，且增速明显，总提升值为0.21，高等级NDVI面积增加了206.19 ha，低等级NDVI面积减少了129.69ha，NDVI显著增长面积占比最高，且绝大部分地区的红树林NDVI增长迅速，仅有少部分的红树林NDVI值缓慢下降。

在未来，深圳湾红树林NDVI值呈上升趋势的面积为658.8 ha，占比89.39%，呈下降趋势的面积为42.93 ha，占比5.82%。

海平面高度和海表面温度等因子的变化对红树林的生长具有微弱影响，在海平面高度和海表面温度越高的地区NDVI的值越高，其中海平面高度对NDVI的影响稍微更加显著。

1995~2024年，30年间深圳湾红树林面积增加了478.89 ha，主要增加区域为香港米埔自然保护区的深圳河入湾口处。

NOTES

^*通讯作者。