PotatoTom Blog

The distribution pattern and driving factors of vegetation resilience differences inside and outside protected areas at different levels in China

思路

中国不同层次保护区内外植被韧性差异分布格局及其驱动因子：实验完整流程

0. 研究目标与总体思路

你的文章可以围绕一个主问题展开：

中国不同层次保护区内部、外部及其周边反事实区域的植被韧性是否存在系统差异？这种差异是否表现为保护区外部韧性溢出或泄漏？其空间格局和驱动因子是什么？

本文建议借鉴 Protected area management has significant spillover effects on vegetation 的核心设计：将空间单元划分为保护区内部 A、潜在溢出区 B、独立反事实区 C，并通过 A vs B、B vs C、A vs C 三类对比判断保护区外部效应。该论文强调，判断溢出不能只比较保护区外部和远处区域，还需要同时判断保护区内部与外部近邻区是否相似；其中 B vs C 是溢出效应大小的核心对比。

1. 实验总体框架

1.1 空间框架

每个保护区 (PA_i) 周围构建三个空间单元：

A：保护区内部稳定自然植被像元
B：保护区外部缓冲区稳定自然植被像元
C：远离保护区、但环境条件相似的非保护区稳定自然植被像元

其中外部缓冲区 B 分为：

B1：0–1 km
B2：1–5 km
B3：5–10 km
B4：10–20 km

1.2 核心指标

每个像元计算植被韧性：

[ R=-AR(1) ]

其中：

  • (AR(1)) 越大，表示植被异常状态持续性越强，恢复越慢；
  • (R=-AR(1)) 越大，表示植被韧性越高，恢复越快。

1.3 三类关键差异

保护区内外差异

[ \Delta R_{in-out}=R_{in}-R_{out} ]

外部韧性溢出或泄漏

[ Spillover_R=R_{near\ outside}-R_{matched\ control} ]

保护区内部相对于反事实背景的差异

[ Protection_R=R_{in}-R_{matched\ control} ]

2. 数据准备

2.1 必备数据

数据用途空间分辨率建议时间范围
kNDVI 或由 NDVI 计算的 kNDVI计算植被韧性250 m2000–2020
保护区边界构建 A、B、C 空间单元矢量最新可用
土地覆盖数据剔除农田、城市、水体、非植被和土地覆盖变化像元250 m 或重采样至 250 m2000–2020
植被类型图分组统计和匹配250 m 或更高静态或多年稳定
生态区/气候区分组统计和匹配矢量或栅格静态
DEM海拔、坡度30–250 m静态
气候数据年均温、年降水、干旱指数等1 km 或更高2000–2020
人类压力数据夜间灯光、人口、道路距离、居民点距离250 m–1 km2000–2020 或多年平均

3. 数据预处理

3.1 建立统一空间基准

所有数据必须统一到同一空间参考、分辨率和范围。

操作步骤

  1. 选择中国 Albers 等面积投影。
  2. 以 250 m kNDVI 栅格作为基准模板。
  3. 所有栅格重采样到 250 m。
  4. 所有矢量数据转为同一投影。
  5. 保护区边界栅格化到 250 m。

输出文件

template_250m.tif
pa_id_250m.tif
landcover_annual_250m/
climate_250m/
terrain_250m/
human_pressure_250m/

3.2 kNDVI 构建

如果你已有 250 m kNDVI 产品，可以直接进入质量控制。

如果你从 NDVI 计算 kNDVI，可使用：

[ kNDVI=\tanh(NDVI^2) ]

其中：

[ NDVI=\frac{NIR-Red}{NIR+Red} ]

操作步骤

  1. 读取 16 天或月尺度 NDVI。
  2. 对 NDVI 进行质量控制。
  3. 将 NDVI 转换为 kNDVI。
  4. 合成为月尺度 kNDVI。

3.3 质量控制

对每个时间步的 kNDVI，剔除以下像元：

云
云阴影
雪
低质量观测
异常反射率
明显离群值

建议规则

对每个像元 (p)：

如果 2000–2020 年有效观测比例 < 80%，剔除；
如果连续缺测 > 3 个月，标记为低质量；
如果某个月所有年份有效值少于 15 个，剔除该像元；
如果 kNDVI < 0 或 kNDVI > 1，剔除该值。

3.4 月尺度合成

如果原始数据为 16 天尺度，建议合成为月尺度。

推荐公式

对像元 (p)、年份 (y)、月份 (m)：

[ kNDVI_{p,y,m}=median(kNDVI_{p,y,m,valid}) ]

即使用当月所有有效观测的中位数。

说明

  • 中位数比最大值更稳健；
  • 如果某月只有一个有效观测，可以保留；
  • 如果某月没有有效观测，可用时间序列平滑补齐，但要记录缺测比例。

3.5 时间序列平滑

建议使用 Savitzky–Golay 滤波。

操作步骤

  1. 对每个像元的月尺度 kNDVI 序列按时间排序。
  2. 对短缺测进行线性插值。
  3. 对插值后的序列进行 SG 平滑。
  4. 输出平滑后的月尺度 kNDVI。

推荐参数

时间步长：月
窗口长度：5 或 7
多项式阶数：2

敏感性分析

建议至少比较：

未平滑
SG 窗口 = 5
SG 窗口 = 7

如果主要结论一致，说明平滑方法不会主导结果。

4. 稳定自然植被像元筛选

这是你的研究中非常关键的一步。因为人为土地利用突变会显著影响 AR(1)，导致韧性指标失真。

4.1 剔除非自然植被像元

剔除以下类型：

农田
城市建设用地
水体
永久冰雪
裸地
人工表面
湿地中长期水体主导区域

保留：

森林
灌丛
草地
草甸
荒漠草地
稀疏自然植被
其他自然植被

4.2 剔除土地覆盖变化像元

对每个像元统计 2000–2020 年土地覆盖类型。

规则 1：严格稳定规则

如果 2000–2020 年期间土地覆盖类型始终为同一种自然植被类型，则保留。

2000 = 森林
2001 = 森林
...
2020 = 森林
→ 保留

规则 2：宽松稳定规则

如果 21 年中至少 80% 年份属于同一自然植被大类，且从未转为农田、城市、水体，则保留。

[ StableRatio=\frac{N_{dominant}}{N_{total}} ]

若：

[ StableRatio \geq 0.8 ]

且没有转为农田、城市、水体，则保留。

推荐主分析

使用严格稳定规则。

推荐敏感性分析

使用宽松稳定规则。

4.3 输出稳定自然植被掩膜

最终生成：

stable_natural_vegetation_250m.tif

像元值建议为：

1 = 稳定自然植被
0 = 剔除

同时保留主导植被类型：

stable_veg_type_250m.tif

5. 植被韧性计算

5.1 构建每个像元的月尺度序列

对每个稳定自然植被像元，得到：

[ kNDVI_t ]

其中：

t = 1, 2, 3, ..., 252

2000–2020 年共 21 年，月尺度为：

[ 21 \times 12 = 252 ]

5.2 去季节性

对每个像元分别计算每个月的多年平均值。

[ Season_m=mean(kNDVI_{m}) ]

其中 (m=1,2,…,12)。

然后得到去季节序列：

[ A_t=kNDVI_t-Season_m ]

其中 (A_t) 是 kNDVI 异常值。

5.3 去趋势

对去季节序列拟合线性趋势：

[ A_t=a+bt+\varepsilon_t ]

其中：

  • (a)：截距；
  • (b)：长期趋势；
  • (\varepsilon_t)：去季节、去趋势后的残差。

最终用于计算韧性的序列是：

[ \varepsilon_t ]

5.4 平稳性检验

对每个像元的 (\varepsilon_t) 进行 ADF 检验。

保留规则

如果 ADF 检验 p < 0.05，保留；
如果 p ≥ 0.05，剔除或进入敏感性分析。

输出

adf_pvalue_250m.tif
stationary_mask_250m.tif

5.5 计算 AR(1)

对每个像元，用残差序列计算一阶自相关。

简单公式

[ AR(1)=cor(\varepsilon_t,\varepsilon_{t-1}) ]

其中：

  • (\varepsilon_t)：当前月残差；
  • (\varepsilon_{t-1})：前一个月残差。

5.6 计算植被韧性

[ R=-AR(1) ]

解释：

(AR(1))(R)含义
高低扰动后恢复慢，韧性低
低高扰动后恢复快，韧性高
接近 1接近 -1强持续性，可能接近临界减速
接近 0接近 0异常不持续，恢复较快
小于 0大于 0异常快速反转，韧性较高

输出文件

AR1_2000_2020_250m.tif
R_resilience_2000_2020_250m.tif

6. 保护区数据处理

6.1 保护区属性整理

每个保护区至少需要以下字段：

字段含义
PA_ID保护区唯一编号
PA_name保护区名称
PA_level保护层级
PA_type保护区类型
Province所属省份
Area_km2面积
Year建立年份
Main_veg主导植被类型
Ecozone生态区
Climate_zone气候区

6.2 处理保护区重叠

中国保护区之间可能存在重叠，必须处理。

推荐主方案：层级优先

优先级示例：

国家公园 > 国家级自然保护区 > 省级自然保护区 > 地方级自然保护区 > 自然公园

重叠区域归入最高层级保护区。

推荐敏感性方案：保护区复合体

如果多个保护区边界相连或相距小于 1 km，可合并为一个保护区复合体。

6.3 剔除不适合分析的保护区

剔除或单独标记：

面积过小的保护区
内部稳定自然植被像元过少的保护区
外部 20 km 内几乎没有稳定自然植被的保护区
边界严重重叠、无法构建独立外部缓冲区的保护区
主要保护对象为水域或地质景观而非植被生态系统的保护区

推荐最小样本量

单元最小像元数
保护区内部 A≥ 30 个像元
每个外部缓冲带 B≥ 30 个像元
候选反事实区 C≥ 300 个像元
匹配后 C至少与 B 数量相当

如果某个保护区不满足要求：

保留在全国地图展示中；
不进入严格统计模型；
在补充材料中报告剔除原因。

7. 构建保护区内外空间单元

7.1 构建保护区内部 A

对每个保护区 (PA_i)，定义：

[ A_i={p:p \in PA_i,\ p为稳定自然植被像元} ]

建议进一步划分：

A_all：保护区内全部稳定自然植被像元
A_core：距离保护区边界 > 1 km 的内部核心像元
A_edge：距离保护区边界 0–1 km 的内部边缘像元

主分析建议

使用：

A_core

作为保护区内部韧性代表。

原因是 A_core 可减少边界混合像元、道路、居民点和缓冲区误差影响。

7.2 构建保护区外部 B

对每个保护区构建外部缓冲区：

B1：0–1 km
B2：1–5 km
B3：5–10 km
B4：10–20 km

定义：

[ B_{i,d}={p:p \in buffer_{i,d},\ p \notin PA,\ p为稳定自然植被像元} ]

其中 (d) 表示距离带。

外部像元筛选规则

每个 B 区像元必须满足：

不是任何保护区内部像元；
是稳定自然植被；
不是农田、城市、水体、裸地；
不是土地覆盖明显变化像元；
最近的保护区是当前 PA_i；
不位于多个保护区缓冲区重叠区，或按最近保护区归属。

推荐主方案

如果外部像元同时位于多个保护区缓冲区，剔除。

推荐敏感性方案

如果外部像元同时位于多个保护区缓冲区，归属于最近保护区。

7.3 构建反事实区 C

C 区是最重要的对照区。它不是简单的“远处区域”，而是：

远离保护区、但环境条件与保护区外部缓冲区相似的非保护区稳定自然植被像元。

定义：

[ C_i={p:p \notin PA,\ Distance(p,nearestPA)>D,\ p为稳定自然植被像元} ]

其中 (D) 是远离保护区的距离阈值。

推荐主阈值

[ D=30km ]

敏感性阈值

D = 20 km
D = 50 km

8. 反事实匹配

8.1 为什么必须做匹配

保护区不是随机分布的，往往位于高海拔、低人口密度、低开发强度或生态敏感区。如果直接比较保护区附近和远处区域，可能把环境背景差异误认为保护区效应。

Cumming 的论文也使用匹配方法控制局地生物物理条件和保护区选址偏差，并用反事实区域估计溢出效应。

8.2 匹配对象

对每个保护区 (PA_i) 和每个距离带 (d)，单独匹配。

处理组

[ T=1: p \in B_{i,d} ]

即保护区外部某一距离带的稳定自然植被像元。

对照组

[ T=0: p \in C_i ]

即远离保护区的候选反事实像元。

8.3 匹配变量

建议匹配变量分为 5 类。

类型变量
气候年均温、年降水、干旱指数、降水季节性
地形海拔、坡度、坡向、地形起伏度
土壤土壤有机碳、土壤质地、土壤水分、土壤深度
初始植被状态2000–2002 年平均 kNDVI、kNDVI 季节振幅
人类压力距道路距离、距居民点距离、人口密度、夜间灯光、距耕地距离

不要用于匹配的变量

不要把最终结果变量用于匹配：

不要用 R
不要用 AR(1)
不要用 2000–2020 全期平均 kNDVI 趋势

否则可能把需要解释的差异提前控制掉。

8.4 匹配方法

推荐使用两步匹配。

第一步：精确分层

先要求 B 和 C 属于相同或相近类别：

同一生态区
同一主导植被类型
同一气候区

如果样本不足，可以放宽为：

同一生态区 + 同一植被大类

第二步：距离匹配

对连续变量做标准化：

[ Z=\frac{X-\bar{X}}{SD} ]

然后计算 B 像元和 C 像元之间的环境距离：

[ Dist(p,q)=\sqrt{\sum (Z_p-Z_q)^2} ]

其中：

  • (p)：B 区像元；
  • (q)：候选 C 区像元；
  • (Z)：标准化后的匹配变量。

为每个 B 像元选择最相似的 (K) 个 C 像元。

推荐参数

K = 3 或 5
最大环境距离 caliper = 0.5 或 1.0

如果没有满足 caliper 的 C 像元，该 B 像元不参与匹配分析。

8.5 匹配质量检验

对每个匹配变量计算标准化均差：

[ SMD=\frac{\bar{X}_B-\bar{X}_C}{SD} ]

判断标准：

| (|SMD|) | 质量 | |---:|---| | < 0.1 | 可接受 | | < 0.05 | 较好 | | < 0.01 | 非常好 |

保留规则

如果所有主要变量 |SMD| < 0.1，保留该保护区—距离带组合；
如果部分变量 |SMD| ≥ 0.1，重新匹配或扩大候选 C 区；
如果仍不合格，剔除该组合。

9. 计算保护区内外韧性差异

9.1 计算保护区内部韧性

对每个保护区：

[ R_{in,i}=median(R_p), \quad p \in A_i ]

建议同时计算：

mean
median
standard deviation
25% 分位数
75% 分位数
像元数

主文使用中位数，补充材料报告均值。

9.2 计算外部缓冲区韧性

对每个保护区和每个距离带：

[ R_{out,i,d}=median(R_p), \quad p \in B_{i,d} ]

其中：

d = 0–1 km
d = 1–5 km
d = 5–10 km
d = 10–20 km

9.3 计算内外差异

[ \Delta R_{in-out,i,d}=R_{in,i}-R_{out,i,d} ]

解释：

(\Delta R_{in-out})含义
> 0保护区内部韧性高于外部
< 0外部韧性高于保护区内部
≈ 0内外韧性相近

9.4 输出保护区级结果表

建议生成：

PA_resilience_inner_outer.csv

字段包括：

字段含义
PA_ID保护区编号
PA_name保护区名称
PA_level保护层级
PA_type保护区类型
Province省份
Ecozone生态区
Veg_type主导植被类型
Climate_zone气候区
Area_km2面积
R_in内部韧性
R_out_0_10–1 km 外部韧性
R_out_1_51–5 km 外部韧性
R_out_5_105–10 km 外部韧性
R_out_10_2010–20 km 外部韧性
DeltaR_0_1内部 - 0–1 km
DeltaR_1_5内部 - 1–5 km
DeltaR_5_10内部 - 5–10 km
DeltaR_10_20内部 - 10–20 km
N_in内部像元数
N_out_0_10–1 km 像元数
N_out_1_51–5 km 像元数
N_out_5_105–10 km 像元数
N_out_10_2010–20 km 像元数

10. 计算外部韧性溢出或泄漏

10.1 计算匹配反事实区韧性

对每个保护区和距离带，计算匹配后的 C 区韧性：

[ R_{control,i,d}=mean(R_q), \quad q \in C_{matched,i,d} ]

如果匹配方法产生权重，则使用加权平均：

[ R_{control,i,d}=\sum w_qR_q ]

其中：

  • (w_q)：匹配权重；
  • (R_q)：匹配对照像元韧性。

10.2 计算韧性溢出效应

[ Spillover_{R,i,d}=R_{out,i,d}-R_{control,i,d} ]

解释：

(Spillover_R)含义
> 0保护区外部近邻区韧性高于匹配反事实区，可能存在正向韧性溢出
< 0保护区外部近邻区韧性低于匹配反事实区，可能存在韧性泄漏或负向边缘效应
≈ 0没有明显外部韧性效应

10.3 建立 A–B–C 三对比指标

借鉴 Cumming 论文的三空间单元对比逻辑，你的韧性研究建议计算以下三个指标。该论文的 Fig. 2 使用 A、B、C 三组样本区分正向溢出、负向溢出和无溢出，并指出 B vs C 是溢出幅度的核心对比，但严格确认还要看 A vs B。

Contrast 1：内部 vs 外部

[ C1_{i,d}=R_{in,i}-R_{out,i,d} ]

表示保护区内部与外部近邻区是否相似。

Contrast 2：外部 vs 反事实

[ C2_{i,d}=R_{out,i,d}-R_{control,i,d} ]

这就是主溢出指标：

[ Spillover_R=C2 ]

Contrast 3：内部 vs 反事实

[ C3_{i,d}=R_{in,i}-R_{control,i,d} ]

表示保护区内部是否高于相似环境背景下的非保护区自然植被。

11. 溢出、泄漏和无效应分类

11.1 设置阈值

因为 (R=-AR(1))，其理论范围大致在 -1 到 1 之间。建议设置一个简单阈值：

[ \delta=0.05 ]

即：

差异绝对值小于 0.05，认为生态意义较弱；
差异大于 0.05，认为有较明确差异。

敏感性分析

同时测试：

δ = 0.02
δ = 0.05
δ = 0.10

11.2 正向韧性溢出

如果：

[ C2_{i,d}>\delta ]

则认为外部近邻区相对于反事实区韧性更高。

更严格的正向溢出要求：

[ C2_{i,d}>\delta ]

且：

[ C3_{i,d}>\delta ]

且：

[ |C1_{i,d}|<\delta ]

解释：

保护区内部韧性高于反事实区；
保护区外部近邻区也高于反事实区；
保护区外部近邻区与保护区内部相似；
因此可判断为较严格的正向韧性溢出。

11.3 韧性泄漏或负向边缘效应

如果：

[ C2_{i,d}<-\delta ]

则认为外部近邻区韧性低于匹配反事实区。

如果同时：

[ C3_{i,d}>\delta ]

且：

[ C1_{i,d}>\delta ]

说明：

保护区内部韧性较高；
但保护区外部近邻区韧性较低；
可能存在边界外人为压力集中、资源利用转移、道路干扰、旅游开发或管理断裂。

11.4 无明显外部效应

如果：

[ |C2_{i,d}| \leq \delta ]

则认为没有明显韧性溢出或泄漏。

11.5 保护区内部有效但无外部溢出

如果：

[ C3_{i,d}>\delta ]

且：

[ |C2_{i,d}| \leq \delta ]

说明：

保护区内部韧性高于相似非保护区背景；
但这种优势没有扩展到外部缓冲区。

12. 距离衰减分析

12.1 构建距离序列

每个保护区都有四个距离带的溢出效应：

Spillover_R_0_1
Spillover_R_1_5
Spillover_R_5_10
Spillover_R_10_20

为每个距离带设置中点：

距离带中点
0–1 km0.5 km
1–5 km3 km
5–10 km7.5 km
10–20 km15 km

12.2 计算简单距离衰减斜率

对每个保护区拟合：

[ Spillover_R=a+bD ]

其中：

  • (D)：距离带中点；
  • (b)：距离衰减斜率。

解释：

(b)含义
< 0离保护区越远，正向溢出越弱
> 0离保护区越远，韧性越高，近边界可能存在负效应
≈ 0溢出效应随距离变化不明显

12.3 全国距离衰减曲线

按保护层级、生态区、植被类型分别计算平均溢出：

[ \overline{Spillover}{d}=mean(Spillover{R,i,d}) ]

绘制：

横轴：距离带
纵轴：Spillover_R
分组：保护层级 / 生态区 / 植被类型 / 气候区

13. 全国格局分析

13.1 按保护层级统计

分组计算：

[ \overline{R}_{in} ]

[ \overline{R}_{out,d} ]

[ \overline{\Delta R}_{in-out,d} ]

[ \overline{Spillover}_{R,d} ]

分组包括：

国家公园
国家级自然保护区
省级自然保护区
地方级自然保护区
自然公园
其他保护地

13.2 按生态区统计

每个生态区统计：

保护区数量
平均 R_in
平均 R_out
平均 DeltaR
平均 Spillover_R
正向溢出比例
负向泄漏比例
无效应比例

正向溢出比例：

[ P_{positive}=\frac{N(Spillover_R>\delta)}{N} ]

负向泄漏比例：

[ P_{negative}=\frac{N(Spillover_R<-\delta)}{N} ]

无效应比例：

[ P_{none}=\frac{N(|Spillover_R|\leq \delta)}{N} ]

13.3 按植被类型统计

分组包括：

森林
灌丛
草地
高寒草甸
荒漠草地
稀疏植被

重点比较：

哪类植被内部韧性最高？
哪类植被最容易出现正向溢出？
哪类植被最容易出现边界外韧性泄漏？
草地和森林的距离衰减是否不同？

13.4 按气候区统计

分组包括：

湿润区
半湿润区
半干旱区
干旱区
高寒区

重点判断：

干旱区保护区是否更容易出现韧性溢出？
湿润区森林保护区是否主要表现为内部优势？
高寒区保护区内外是否差异较小？

14. 驱动因子分析

14.1 构建保护区级驱动因子表

每个保护区—距离带一行。

输出文件：

PA_spillover_driver_table.csv

字段包括：

字段含义
PA_ID保护区编号
Distance_band距离带
Spillover_R韧性溢出效应
DeltaR内外差异
R_in内部韧性
R_out外部韧性
PA_level保护层级
PA_type保护区类型
Area_km2面积
Year建立年份
Veg_type植被类型
Ecozone生态区
Climate_zone气候区
MAT年均温
MAP年降水
AI干旱指数
Elevation海拔
Slope坡度
Soil土壤变量
Road_dist距道路距离
Settlement_dist距居民点距离
Nightlight夜间灯光
Population人口密度
Cropland_ratio_20km周边农田比例
Built_ratio_20km周边建设用地比例
NDVI_initial初始 kNDVI
kNDVI_trendkNDVI 长期趋势

14.2 主模型：RF / XGBoost + SHAP

响应变量

主响应变量：

[ Y=Spillover_R ]

也可以分别建模：

Y = DeltaR
Y = R_in
Y = R_out

预测变量

包括：

保护层级
保护区面积
建立年份
植被类型
生态区
气候区
年均温
年降水
干旱指数
海拔
坡度
土壤变量
夜间灯光
人口密度
道路距离
居民点距离
周边农田比例
周边建设用地比例

操作步骤

  1. 读取 PA_spillover_driver_table.csv。

  2. 删除缺失严重变量。

  3. 类别变量做 one-hot 编码。

  4. 连续变量标准化。

  5. 按保护区 ID 划分训练集和测试集，避免同一保护区不同距离带同时进入训练和测试。

  6. 训练 RF 或 XGBoost。

  7. 使用交叉验证选择参数。

  8. 计算模型性能：

    • (R^2)
    • RMSE
    • MAE

  9. 计算 SHAP 值。

  10. 输出变量重要性和 SHAP 依赖图。

模型评价公式

[ RMSE=\sqrt{mean((Y_{obs}-Y_{pred})^2)} ]

[ MAE=mean(|Y_{obs}-Y_{pred}|) ]

14.3 解释模型：GAMM

GAMM 用于解释非线性关系，并控制保护区重复观测。

推荐模型

[ Spillover_R=s(MAP)+s(MAT)+s(AI)+s(Elevation)+s(HumanPressure)+PA_level+VegType+DistanceBand+u_{PA} ]

其中：

  • (s())：平滑函数；
  • (u_{PA})：保护区随机效应。

解释重点

降水是否存在阈值？
干旱指数是否调节溢出效应？
人类压力升高是否导致 Spillover_R 由正转负？
不同保护层级是否仍然显著？

14.4 机制模型：SEM

SEM 用于表达假设机制。

推荐机制路径

气候背景 → 植被韧性
地形条件 → 人类压力 → 外部韧性
保护层级 → 人类压力下降 → 外部韧性上升
保护区面积 → 边界干扰降低 → 正向溢出增强
周边农田/建设用地比例 → 外部韧性降低

简化路径公式

[ HumanPressure = Climate + Terrain + PA_level ]

[ R_{out}=HumanPressure + Climate + Terrain + PA_level ]

[ Spillover_R=R_{out}-R_{control} ]

实际 SEM 可设置为：

PA_level, Area → HumanPressure
Climate, Terrain → R_in, R_out
HumanPressure → R_out
R_in, R_out, R_control → Spillover_R

14.5 空间异质性：MGWR 或滑动窗口

MGWR 响应变量

[ Y=Spillover_R ]

解释变量

人类压力
降水
干旱指数
海拔
保护区面积
保护层级
周边农田比例
周边建设用地比例

目标

判断不同地区驱动因子的作用方向和强度是否不同。

例如：

在东部人口密集区，人类压力可能主导负向泄漏；
在西北干旱区，降水和干旱指数可能主导韧性；
在青藏高原，海拔和温度限制可能更重要。

15. 统计检验

15.1 内外差异是否显著

以保护区为单位，使用配对检验：

[ R_{in,i} \quad vs. \quad R_{out,i,d} ]

如果数据近似正态：

配对 t 检验

如果不正态：

Wilcoxon 符号秩检验

15.2 不同保护层级是否不同

响应变量：

DeltaR
Spillover_R

检验方法：

Kruskal-Wallis 检验
Dunn 事后检验

或使用混合效应模型：

[ Spillover_R=PA_level+DistanceBand+VegType+ClimateZone+u_{PA} ]

15.3 溢出比例是否不同

例如比较不同保护层级的正向溢出比例：

[ P_{positive}=\frac{N(Spillover_R>\delta)}{N} ]

可使用：

卡方检验
二项 logistic 回归

logistic 回归形式：

[ Positive=PA_level+VegType+ClimateZone+HumanPressure ]

其中：

Positive = 1，如果 Spillover_R > δ
Positive = 0，否则

16. 不确定性估计

16.1 像元重采样

对每个保护区和每个距离带：

从 B 区像元中有放回抽样；
从匹配 C 区像元中有放回抽样；
重复 1000 次；
每次计算 Spillover_R；
取 2.5% 和 97.5% 分位数作为 95% CI。

16.2 空间块 Bootstrap

因为像元存在空间自相关，建议使用空间块而非单像元重采样。

操作步骤

  1. 将全国划分为 5 km × 5 km 或 10 km × 10 km 网格。
  2. 每个像元归入一个空间块。
  3. Bootstrap 时抽样空间块。
  4. 重新计算 (R_{in})、(R_{out})、(Spillover_R)。

16.3 置信区间判定

若：

[ CI_{95%}(Spillover_R)>0 ]

则认为正向溢出较稳健。

若：

[ CI_{95%}(Spillover_R)<0 ]

则认为负向泄漏较稳健。

若置信区间跨 0，则认为不确定。

17. 敏感性分析

至少做 8 类敏感性分析。

17.1 韧性指标敏感性

比较：

R = -AR(1)
R_alt = 1 - AR(1)
AR(1) 原值

主文使用：

[ R=-AR(1) ]

补充材料报告其他指标。

17.2 平滑方法敏感性

比较：

未平滑
SG 窗口 5
SG 窗口 7

17.3 稳定植被掩膜敏感性

比较：

严格稳定自然植被
80% 年份稳定自然植被
90% 年份稳定自然植被

17.4 反事实距离阈值敏感性

比较 C 区定义：

远离所有保护区 >20 km
远离所有保护区 >30 km
远离所有保护区 >50 km

17.5 外部缓冲区敏感性

比较：

0–1 km
0.5–1 km
1–5 km
5–10 km
10–20 km

尤其检查 0–1 km 是否受边界误差影响。

17.6 匹配方法敏感性

比较：

最近邻匹配
1:3 匹配
1:5 匹配
全匹配
倾向得分匹配
环境距离匹配

17.7 阈值敏感性

比较：

[ \delta=0.02 ]

[ \delta=0.05 ]

[ \delta=0.10 ]

17.8 保护区重叠处理敏感性

比较：

重叠缓冲区剔除
最近保护区归属
保护区复合体合并

18. 结果图件设计

图 1：技术路线图

展示：

kNDVI 预处理 → 韧性计算 → A/B/C 空间单元 → 匹配 → 内外差异 → 溢出效应 → 驱动因子

图 2：中国植被韧性空间分布图

显示全国稳定自然植被像元的：

[ R=-AR(1) ]

图 3：保护区内部韧性分布图

每个保护区用：

[ R_{in} ]

着色。

图 4：保护区内外韧性差异图

每个保护区用：

[ \Delta R_{in-out} ]

着色。

可分别绘制：

0–1 km
1–5 km
5–10 km
10–20 km

图 5：韧性溢出/泄漏空间图

每个保护区或保护区边界外缓冲区用：

[ Spillover_R ]

着色。

颜色建议：

正值：正向韧性溢出
负值：韧性泄漏
接近 0：无明显效应

图 6：不同保护层级箱线图

展示：

R_in
DeltaR
Spillover_R

分组：

国家公园
国家级自然保护区
省级自然保护区
地方级自然保护区
自然公园

图 7：距离衰减曲线

横轴：

0.5 km
3 km
7.5 km
15 km

纵轴：

[ Spillover_R ]

分组：

保护层级
植被类型
气候区
生态区

图 8：驱动因子 SHAP 图

包括：

变量重要性条形图
SHAP beeswarm 图
关键变量依赖图

图 9：GAMM 非线性响应图

展示：

Spillover_R 对降水的响应
Spillover_R 对干旱指数的响应
Spillover_R 对人类压力的响应
Spillover_R 对保护区面积的响应

图 10：MGWR 空间异质性图

显示不同驱动因子的局部回归系数：

人类压力系数空间分布
降水系数空间分布
保护区面积系数空间分布
道路距离系数空间分布

19. 结果表格设计

表 1：数据来源与预处理

数据来源分辨率时间用途

表 2：保护区样本统计

保护层级原始数量保留数量剔除数量平均面积稳定自然植被比例

表 3：不同保护层级韧性差异

保护层级(R_{in})(R_{out,0-1})(R_{out,1-5})(\Delta R)(Spillover_R)

表 4：不同植被类型韧性溢出

植被类型正向溢出比例负向泄漏比例无效应比例平均 (Spillover_R)

表 5：模型性能

模型响应变量(R^2)RMSEMAE

表 6：敏感性分析结果

敏感性方案正向溢出比例负向泄漏比例平均 (Spillover_R)与主结果是否一致

20. 推荐文件结构

project/
│
├── data_raw/
│   ├── kNDVI/
│   ├── protected_areas/
│   ├── landcover/
│   ├── climate/
│   ├── terrain/
│   ├── human_pressure/
│
├── data_processed/
│   ├── kNDVI_monthly_250m/
│   ├── kNDVI_sg_250m/
│   ├── stable_vegetation_mask/
│   ├── resilience/
│   ├── PA_raster/
│   ├── buffers/
│   ├── matched_controls/
│
├── scripts/
│   ├── 01_preprocess_kNDVI.py
│   ├── 02_stable_vegetation_mask.py
│   ├── 03_calculate_resilience.py
│   ├── 04_prepare_PA_buffers.py
│   ├── 05_matching_controls.py
│   ├── 06_calculate_spillover.py
│   ├── 07_statistical_analysis.R
│   ├── 08_driver_RF_XGBoost.py
│   ├── 09_GAMM_SEM_MGWR.R
│
├── results/
│   ├── tables/
│   ├── figures/
│   ├── maps/
│   ├── model_outputs/
│
└── manuscript/

21. 可以直接执行的实验步骤清单

阶段一：数据预处理

1. 下载或准备 2000–2020 年 250 m kNDVI 数据。
2. 统一投影为中国 Albers。
3. 对 kNDVI 进行质量控制。
4. 合成为月尺度 kNDVI。
5. 对月尺度 kNDVI 进行 SG 平滑。
6. 准备 2000–2020 年土地覆盖数据。
7. 剔除农田、城市、水体、裸地等非自然植被。
8. 剔除 2000–2020 年土地覆盖发生明显变化的像元。
9. 得到稳定自然植被掩膜。

阶段二：韧性计算

10. 对每个稳定自然植被像元提取 2000–2020 年月尺度 kNDVI。
11. 计算每个月多年平均值。
12. 用原始序列减去月平均值，得到去季节异常。
13. 对异常序列拟合线性趋势。
14. 提取残差序列。
15. 对残差序列进行 ADF 平稳性检验。
16. 对通过 ADF 检验的像元计算 AR(1)。
17. 计算 R = -AR(1)。
18. 输出全国植被韧性栅格。

阶段三：保护区空间单元构建

19. 清理保护区边界。
20. 处理保护区重叠。
21. 给每个保护区赋予 PA_ID 和保护层级。
22. 栅格化保护区边界。
23. 构建保护区内部 A。
24. 构建内部核心区 A_core。
25. 构建外部缓冲区 B1、B2、B3、B4。
26. 剔除 B 中非稳定自然植被像元。
27. 剔除 B 中位于其他保护区内部的像元。
28. 处理多个保护区缓冲区重叠像元。

阶段四：反事实区构建与匹配

29. 构建远离所有保护区 >30 km 的候选 C 区。
30. 只保留 C 中稳定自然植被像元。
31. 给 B 和 C 像元提取气候、地形、土壤、人类压力和初始 kNDVI 变量。
32. 对每个保护区和每个距离带单独匹配。
33. 先按生态区、植被类型、气候区分层。
34. 在分层内用环境距离或倾向得分匹配。
35. 检查匹配质量 SMD。
36. 剔除匹配质量不合格的保护区—距离带组合。

阶段五：内外差异和溢出效应计算

37. 计算每个保护区 R_in。
38. 计算每个保护区每个距离带 R_out。
39. 计算 DeltaR = R_in - R_out。
40. 计算匹配反事实区 R_control。
41. 计算 Spillover_R = R_out - R_control。
42. 计算 C1 = R_in - R_out。
43. 计算 C2 = R_out - R_control。
44. 计算 C3 = R_in - R_control。
45. 按 δ = 0.05 分类为正向溢出、负向泄漏、无效应等类型。
46. 用 bootstrap 估计 Spillover_R 的置信区间。

阶段六：空间格局分析

47. 绘制全国 R 分布图。
48. 绘制保护区 R_in 分布图。
49. 绘制 DeltaR 分布图。
50. 绘制 Spillover_R 分布图。
51. 按保护层级统计 R_in、R_out、DeltaR、Spillover_R。
52. 按生态区统计正向溢出比例和负向泄漏比例。
53. 按植被类型统计溢出距离衰减。
54. 按气候区比较韧性差异。

阶段七：驱动因子分析

55. 构建保护区—距离带驱动因子表。
56. 使用 RF 或 XGBoost 建模 Spillover_R。
57. 使用 SHAP 解释变量重要性。
58. 使用 GAMM 分析非线性响应。
59. 使用 SEM 检验机制路径。
60. 使用 MGWR 或滑动窗口分析空间异质性。

阶段八：稳健性检验

61. 更换平滑方法。
62. 更换稳定植被定义。
63. 更换 C 区距离阈值。
64. 更换匹配方法。
65. 更换 δ 阈值。
66. 更换保护区重叠处理规则。
67. 比较主结果是否稳定。

22. 论文方法部分可直接使用的表述

22.1 韧性计算表述

本研究基于 2000–2020 年月尺度 250 m kNDVI 时间序列计算植被韧性。首先对 kNDVI 进行质量控制、月尺度合成和 SG 平滑；随后剔除非自然植被、农田、城市、水体以及 2000–2020 年土地覆盖发生明显变化的像元，仅保留长期稳定自然植被。对每个稳定自然植被像元，先去除月尺度季节循环，再去除长期线性趋势，并对残差序列进行 ADF 平稳性检验。通过检验的像元计算残差序列的一阶自相关系数 (AR(1))，并定义 (R=-AR(1)) 作为植被韧性指标。该指标数值越大，表示植被异常状态持续性越弱、恢复速度越快、韧性越高。

22.2 A–B–C 空间设计表述

为识别保护区内外植被韧性差异及外部韧性溢出效应，本研究借鉴保护区溢出效应研究中的 A–B–C 反事实框架，将每个保护区划分为保护区内部稳定自然植被区 A、保护区外部不同距离缓冲区 B，以及远离保护区但环境条件相似的非保护区反事实区 C。外部缓冲区包括 0–1、1–5、5–10 和 10–20 km。反事实区 C 由远离所有保护区至少 30 km 的稳定自然植被像元构成，并通过生态区、植被类型、气候区分层以及气候、地形、土壤、人类压力和初始 kNDVI 等变量匹配获得。

22.3 溢出效应表述

对每个保护区及其每个外部缓冲区，本研究计算三类差异：(C1=R_{in}-R_{out})、(C2=R_{out}-R_{control})、(C3=R_{in}-R_{control})。其中 (C2) 被定义为外部韧性溢出效应，即 (Spillover_R=R_{out}-R_{control})。当 (Spillover_R>0) 时，表示保护区外部近邻区韧性高于匹配反事实区，可能存在正向韧性溢出；当 (Spillover_R<0) 时，表示保护区外部近邻区韧性低于匹配反事实区，可能存在韧性泄漏或负向边缘效应。为避免将微弱差异过度解释，本研究以 (\delta=0.05) 作为主要生态意义阈值，并使用 (\delta=0.02) 和 (\delta=0.10) 进行敏感性分析。

23. 最终核心公式汇总

kNDVI

[ kNDVI=\tanh(NDVI^2) ]

去季节异常

[ A_t=kNDVI_t-Season_m ]

去趋势残差

[ A_t=a+bt+\varepsilon_t ]

一阶自相关

[ AR(1)=cor(\varepsilon_t,\varepsilon_{t-1}) ]

植被韧性

[ R=-AR(1) ]

保护区内部韧性

[ R_{in}=median(R_A) ]

保护区外部韧性

[ R_{out}=median(R_B) ]

内外差异

[ \Delta R_{in-out}=R_{in}-R_{out} ]

反事实区韧性

[ R_{control}=mean(R_C) ]

韧性溢出效应

[ Spillover_R=R_{out}-R_{control} ]

A–B–C 三对比

[ C1=R_{in}-R_{out} ]

[ C2=R_{out}-R_{control} ]

[ C3=R_{in}-R_{control} ]

距离衰减

[ Spillover_R=a+bD ]

24. 最推荐的主实验设定

项目主设定
时间范围2000–2020
时间尺度月尺度
空间分辨率250 m
植被指数kNDVI
平滑方法SG 滤波，窗口 5 或 7
韧性指标(R=-AR(1))
稳定植被定义2000–2020 年长期稳定自然植被
内部区保护区内距离边界 >1 km 的稳定自然植被
外部区0–1、1–5、5–10、10–20 km
反事实区远离所有保护区 >30 km
匹配变量气候、地形、土壤、人类压力、初始 kNDVI
匹配质量(SMD<0.1)
溢出指标(Spillover_R=R_{out}-R_{control})
主阈值(\delta=0.05)
主驱动模型RF/XGBoost + SHAP
解释模型GAMM
机制模型SEM
空间异质性MGWR 或滑动窗口

25. 最终论文逻辑

你的文章可以形成如下逻辑链：

长期稳定自然植被 kNDVI 时间序列
        ↓
去季节、去趋势、平稳性检验
        ↓
计算 AR(1) 和 R = -AR(1)
        ↓
构建保护区内部 A、外部缓冲区 B、匹配反事实区 C
        ↓
计算 R_in、R_out、R_control
        ↓
分析 ΔR_in-out：保护区内外韧性差异
        ↓
分析 Spillover_R：保护区外部韧性溢出或泄漏
        ↓
按保护层级、生态区、植被类型、气候区统计空间格局
        ↓
用 RF/XGBoost + SHAP、GAMM、SEM、MGWR 识别驱动机制
        ↓
提出保护区外围治理、生态廊道建设和分区管理建议