企业官网建设流程全解析

做NPP数据的趋势分析，几乎是生态遥感里最常遇到的需求之一。Theil-Sen Median斜率估计搭配Mann-Kendall趋势分析，这套组合我已经用了很多年，处理过的多年NPP数据从省域到全国都有。这篇文章会把方法原理、数据预处理、Python完整实现、结果分类，以及我实测踩过的坑一次讲清楚，手头正在处理NPP、NDVI、LAI这类栅格时序的同学可以直接照着抄。

1. 为什么给多年NPP数据做趋势分析，我首选Theil-Sen+Mann-Kendall

1.1 NPP数据的“脾气”决定了常规回归容易翻车

NPP的中文全称是净初级生产力，指绿色植物在单位时间和单位面积内通过光合作用固定的有机碳，扣除自身呼吸消耗之后的部分。说人话就是：一片地一年到头真正积累了多少“干货碳”。它直接反映植被生产力水平，是生态遥感里衡量生态系统健康的核心指标之一，也是碳循环研究里最常用的输入量。

但真正的NPP栅格数据远没有课本里画的那么规矩。拿MODIS的MOD17A3HGF产品举例，它是500米分辨率、逐年的全球NPP数据。处理多了你就会发现几个鲜明特点。

第一个特点是非正态分布严重。干旱半干旱区的大片像元NPP值集中在低值区间，而森林、农田区域又拖出长长的右尾。你没法指望用均值、方差那一套经典统计假设去描述它，很多经典参数检验方法在这里并不适用。

第二个特点是异常值多。云残留、冰雪覆盖、传感器退化、气溶胶污染，都会让某一年某个像元的NPP出现离谱跳变。有些值比正常范围高好几倍，有些直接掉到0附近。这些异常值不是偶发，而是系统性的、年年都有。

第三个特点是时间序列短且端点敏感。常见产品也就二十几年，如果首尾年份出现一个异常值，普通最小二乘回归的斜率会被明显拽动。这种情况下，回归线几乎是被“端点绑架”的。

这才是问题的根源。最小二乘线性回归本质上是“均值回归”，对异常值没有免疫力，一个坏点就能改变整条趋势线的走向。我第一次用普通线性回归跑某省2000-2020年NPP趋势时，草原区有一大片像元显示出“显著退化”。后来一查原因，居然是2010年那期数据受严重云污染影响，当年NPP被系统性低估。这种伪趋势在生态结论里是非常致命的，也是从那之后我再也不敢拿普通回归直接交差了。

1.2 Theil-Sen + Mann-Kendall到底是一套什么样的组合

Theil-Sen Median斜率估计和Mann-Kendall趋势检验，这两个方法经常成对出现。一句话概括：前者算趋势的“大小”，后者检验趋势的“可信度”。

Theil-Sen斜率估计的思路非常朴素。对时序数据中所有点对做两两配对，计算每一对之间的斜率，最后取这些斜率的中位数作为整体趋势斜率。因为是取中位数而不是取均值，天然对异常值不敏感，抗差能力极强。统计上它的崩溃点（breakdown point）可以到大约29.3%，意思是在极端情况下，即使有近三成的数据是异常值，它给出的斜率依然不会彻底失灵。

Mann-Kendall检验则是一个非参数的趋势显著性检验方法。它不关心数据服从什么分布，也不要求方差齐性，只比较每个数据点之间的大小关系（上升、下降、持平），统计出一个S统计量，再用正态近似算出标准化Z值和p值，判断趋势是否显著。它和Kendall's tau系数在数学上是同一套逻辑，可以看作tau趋势检验的一种具体应用。

把两者搭配起来，逻辑非常清楚：Theil-Sen给你“趋势方向和速率”，Mann-Kendall给你“这个趋势到底可不可信”。两者结合，就能得到生态遥感里最常见的趋势分类图。这也是国内外文献里分析NPP、NDVI、LAI、GPP长时间序列变化的主流做法，比单纯回归要稳得多。

2. 原理拆解：Theil-Sen斜率估计与Mann-Kendall检验的底层逻辑

2.1 Theil-Sen斜率到底是怎么算出来的

假设你有n年的NPP时间序列，年份记为x，NPP值记为y。把所有满足 $i < j$ 的点对 $(x_i, y_i)$ 和 $(x_j, y_j)$ 都拿出来，计算它们之间的斜率：

$$ \beta_{ij} = \frac{y_j - y_i}{x_j - x_i} $$

把所有 $\beta_{ij}$ 从小到大排序，取中位数，就是Theil-Sen斜率：

$$ \beta = \operatorname{median}{\beta_{ij} \mid 1 \le i < j \le n} $$

举个例子，假设有6年的NPP数据（单位g C/m²/yr）：2000年500、2001年520、2002年490、2003年530、2004年560、2005年545。年份间隔统一为1年，点对数量是 $C_6^2 = 15$ 个。把这15个点对的斜率都算出来，按从小到大排列，取第8个（中位数），得到的就是这套时序的Theil-Sen斜率。如果你有20年数据，点对数是190个，计算量并不大。

需要注意一个细节：Theil-Sen斜率的单位是“NPP单位除以时间单位”。NPP常用g C/m²/yr时，斜率就是“每年增长或减少多少g C/m²”，比如β=3.2表示NPP平均每年增加3.2 g C/m²。这个数字可以直接用来表达NPP变化速率，也可以在像元尺度上累加得到区域总生产力变化量。

这里补充一个很容易犯的错误：如果年份不是逐年连续的，比如中间有断年，或者年份间隔不一致，分母必须除以实际年数差，不能想当然用序号1、2、3代替。我后面在代码部分会专门强调这一点。

2.2 为什么中位数能抗干扰

先做个生活化类比。五个人报身高：170、172、171、169、300。平均值是196.4，一看就不对劲；但取中位数是171，基本反映了真实情况。Theil-Sen用同样的哲学处理斜率：即使某一年因云污染出现一个离谱的NPP值，它只会影响跟这个点有关的那些点对斜率，而对整体中位数的影响非常有限。

这也正是它跟最小二乘回归的本质区别。最小二乘在目标函数里对每个点的误差做平方惩罚，异常值误差巨大、权重巨大，会“拉着”回归线往自己方向偏；Theil-Sen则把每个点对斜率都当作一次“投票”，异常值只有少数几个“投票权”，中位数又对极端投票不敏感。所以在处理现实NPP数据时，Theil-Sen的稳定性要明显优于普通回归。

很多人会问：那直接用中位数回归（quantile regression的0.5分位）行不行？理论上可以，但计算复杂度和稳定性不一定比Theil-Sen好。Theil-Sen的实现极其简单、结果可解释性强，这也是它在遥感领域长盛不衰的原因。

2.3 Mann-Kendall检验的S统计量与Z值

Mann-Kendall检验的核心是一个符号统计量S。对时序数据里所有的点对 $(i < j)$，比较 $y_j$ 和 $y_i$ 的大小：

• 如果 $y_j > y_i$，记+1；
• 如果 $y_j < y_i$，记-1；
• 如果 $y_j = y_i$，记0。

把所有点对的符号加起来：

$$ S = \sum_{i=1}^{n-1}\sum_{j=i+1}^{n} \operatorname{sign}(y_j - y_i) $$

S为正说明整体有上升趋势，S为负说明整体有下降趋势。S的绝对值越大，趋势倾向越强。

但光有S还不够，得判断它在统计上是否显著。在零假设（无趋势）下，S近似服从均值为0的正态分布，方差为：

$$ \operatorname{Var}(S) = \frac{n(n-1)(2n+5) - \sum_{p} t_p(t_p-1)(2t_p+5)}{18} $$

这里的 $t_p$ 是第p组相等值（ties）的个数。NPP数据经常会出现大量相同值，尤其是整数型产品，如果不做这个平局校正，方差会被系统性低估，导致显著性检验虚高，这一点非常关键。

然后标准化得到Z统计量：

$$ Z = \begin{cases} \frac{S-1}{\sqrt{\operatorname{Var}(S)}} & S > 0 \ 0 & S = 0 \ \frac{S+1}{\sqrt{\operatorname{Var}(S)}} & S < 0 \end{cases} $$

Z值近似服从标准正态分布。双侧检验下，$|Z| > 1.96$ 对应p < 0.05（95%置信水平），$|Z| > 2.58$ 对应p < 0.01（99%置信水平）。你也可以直接算出p值：p = 2 * (1 - scipy.stats.norm.cdf(abs(Z)))。

这里再提一个进阶注意事项。如果时间序列本身存在显著的自相关，Mann-Kendall检验会倾向于高估趋势的显著性，也就是把随机波动误判成趋势。NPP年度数据通常自相关不算强，但严谨起见可以在分析前画一下自相关图，或者使用趋势预白化的做法（TFPW-MK），即先估计并去除序列中的趋势成分，对残差做预白化后再重新检验。如果处理的是月尺度数据，还要考虑季节周期，那就更适合用Seasonal Mann-Kendall方法，不能直接用原始序列跑。

2.4 斜率和显著性是怎么配合使用的

有了Theil-Sen斜率β和Mann-Kendall检验的Z值（或p值），大多数研究会把两者联合起来给每个像元贴标签。

举个例子，如果β > 0且p < 0.05，说明NPP显著增加；如果β < 0且p < 0.05，说明NPP显著减少；如果p >= 0.05，无论β正负都只能算“不显著变化”。因为Z本身就是统计量，很多文章会直接用“显著改善/不显著变化/显著退化”这三级划分，也有文章进一步细分出“极显著变化”（p < 0.01）。

这里有一个非常容易犯的错误：只看斜率不看显著性。NPP的时序里如果只有单调但微弱的上升，斜率可能为正，但如果不显著，这种趋势在统计上没有说服力，审稿人一眼就能挑出问题。反过来，只看显著性不看斜率也没意义，因为一个“显著”的变化如果速率极小，生态学含义也有限。正确做法一定是两者结合。

3. 数据准备：多年NPP数据去哪里拿、如何预处理

3.1 常见NPP数据产品怎么选

做多年NPP趋势分析，第一步是选对数据。现在主流的产品有这么几类，我列个表给你参考：

我个人最常用的是MOD17A3HGF，理由有三：分辨率够细（500m），年份从2000年一直到当前年，更新稳定；单位是kg C/m²/yr，方便换算；已经有大量文献用同一产品做分析，结果可对标。

GIMMS NDVI反演的NPP产品虽然分辨率粗，但适合做1980年代以来的长时序分析。如果你要做“近40年”NPP趋势，选它更合适。GLASS NPP则是近几年的热门选择，时序长、质量控制做得好，但下载和数据格式处理相对繁琐。

3.2 数据预处理的四个关键环节

拿到NPP数据后，别急着跑趋势分析，先把下面四件事做好。

第一，检查投影和坐标系。全球产品通常用经纬度（WGS84），但在区域研究中，需要统一到目标区域的投影坐标系，否则面积计算和像元对应都会出错。我习惯把分析区域重投影到Albers等积投影或UTM，并且保证所有年份的像元范围、行列数完全一致。

第二，处理单位和无数据值。MOD17A3HGF的原始数据是整型（HDF格式里乘以了10000），需要除以10000转成kg C/m²/yr，再根据需要乘以1000转成g C/m²/yr。它的填充值通常是65533、65534、65535这类，不处理的话会把趋势计算彻底搞乱。更麻烦的是有些产品用0表示水域、用负数表示特殊状态，需要仔细看数据说明书。

第三，做时序上的筛选与掩膜。城市、水体、裸岩、永久冰雪这些非植被像元，NPP通常没有意义，最好用土地覆盖数据（如MCD12Q1）或NPP本身的多年均值阈值把它掩膜掉。否则建筑物像元NPP常年为0或极小，在趋势分类里会被误判为“显著退化”，非常影响区域统计结果。

第四，检查数据质量标志。MOD17A3HGF有对应的质量控制图层，建议把质量差的像元剔除或标记。如果某个像元在多个年份被标记为低质量，直接排除比强行插值更稳妥。因为插值会引入人为趋势，这是趋势分析里最忌讳的。

4. 代码实操：Python实现Theil-Sen斜率估计与Mann-Kendall检验

4.1 单像元时间序列的完整实现

先从一个像元讲起。假设你已经把某个像元2000-2021年共22年的NPP值读到了一维数组里，年份也是对应的整数数组。

我这里写一个既有Theil-Sen斜率、又有Mann-Kendall检验的完整函数。为了便于理解，先用numpy实现，代码里有注释说明每一步在干什么：

import numpy as np from scipy import stats def ts_slope_mk(years, values): """ 单像元Theil-Sen斜率估计 + Mann-Kendall显著性检验 返回: slope, z, p slope: Theil-Sen斜率，单位 = values单位/年 z : Mann-Kendall标准化统计量 p : 双侧p值 """ years = np.asarray(years, dtype=float) values = np.asarray(values, dtype=float) # 剔除无效值像元 valid = np.isfinite(values) if valid.sum() < 3: return np.nan, np.nan, np.nan years = years[valid] values = values[valid] n = len(values) # ---------- Theil-Sen斜率 ---------- # 上三角索引: i < j i_idx, j_idx = np.triu_indices(n, k=1) dy = values[j_idx] - values[i_idx] # y_j - y_i dx = years[j_idx] - years[i_idx] # x_j - x_i valid_slope = dx != 0 slopes = dy[valid_slope] / dx[valid_slope] if len(slopes) == 0: return np.nan, np.nan, np.nan slope = np.median(slopes) # ---------- Mann-Kendall ---------- # S统计量，diffs[i,j] = y_j - y_i diffs = values[None, :] - values[:, None] sign_matrix = np.sign(diffs) iu = np.triu_indices(n, k=1) s = np.sum(sign_matrix[iu]) # 方差，带平局校正 unique, counts = np.unique(values, return_counts=True) ties = counts[counts > 1] var_s = n * (n - 1) * (2 * n + 5) if len(ties) > 0: var_s -= np.sum(ties * (ties - 1) * (2 * ties + 5)) var_s /= 18.0 if var_s <= 0: return slope, 0.0, 1.0 # 标准化Z统计量 if s > 0: z = (s - 1) / np.sqrt(var_s) elif s < 0: z = (s + 1) / np.sqrt(var_s) else: z = 0.0 # 双侧p值 p = 2 * (1 - stats.norm.cdf(abs(z))) return slope, z, p

这个函数里有两个细节值得说。第一个是有效值筛选：如果某一年NPP是NaN，最简单的做法是直接剔除该年后再算，而不是用0填充。但要注意，如果缺失年份太多（比如22年里缺了8年），结果可靠性会大打折扣，建议这种情况下把该像元标记为“数据不足”。第二个是平局校正：很多NPP产品是整数型的，大量像元在多年间数值完全相同，如果不减去ties那一项，方差算小了，p值就会假性偏小，容易得出“假显著”的结论。

调用方式很简单：

years = np.arange(2000, 2022) # 2000-2021 values = np.array([500, 520, 490, 530, 560, 545, ...]) # 该像元实际NPP值 slope, z, p = ts_slope_mk(years, values) print(f"Theil-Sen斜率: {slope:.2f} g C/m2/yr") print(f"Mann-Kendall Z: {z:.3f}, p = {p:.4f}")

如果只是快速验证一两个像元，也可以直接借助现成库pymannkendall，一行代码就能拿到结果：

import pymannkendall as mk res = mk.original_test(values) print(res.slope, res.z, res.p, res.Tau)

不过要注意，pymannkendall的循环在栅格尺度上非常慢，几百万像元根本跑不动。它适合做单点验证或小样本分析，全栅格运算还是用下面的向量化方案更靠谱。

4.2 全栅格逐像元计算：从循环到向量化

真实项目里不可能只算一个像元。一片500米分辨率、覆盖一个省级区域的NPP数据，可能有几百万个有效像元。如果每个像元都调用上述Python纯循环函数，计算量会大到怀疑人生。

我的建议是采用“矩阵化”写法：一次性把多年栅格读入成一个三维数组（年份×行×列），重排成二维数组（年份×像元），然后把所有像元批量计算。核心思路是利用numpy的广播和索引机制，把点对运算向量化。

下面给出一个直接可用的实现，假设你已经把每年的NPP GeoTIFF准备好，文件名类似NPP_2000.tif、NPP_2001.tif：

import numpy as np import rasterio from scipy import stats def load_npp_stack(year_list, path_template): """读取多年NPP的GeoTIFF，返回三维数组和元数据""" stack = [] with rasterio.open(path_template.format(year=year_list[0])) as src: meta = src.meta.copy() stack.append(src.read(1).astype(np.float64)) for year in year_list[1:]: with rasterio.open(path_template.format(year=year)) as src: stack.append(src.read(1).astype(np.float64)) return np.stack(stack, axis=0), meta def trend_analysis_stack(stack, years, nodata=None): """ 全栅格Theil-Sen斜率 + Mann-Kendall检验（向量化版本） stack: (n_years, rows, cols) years: (n_years,) 年份数组 返回 slope, z, p，形状与stack单层一致 """ n_years, rows, cols = stack.shape flat = stack.reshape(n_years, -1) # 有效像元：所有年份都有限且不等于nodata if nodata is not None: valid = np.all(np.isfinite(flat) & (flat != nodata), axis=0) else: valid = np.all(np.isfinite(flat), axis=0) data = flat[:, valid] # (n_years, n_valid_pixels) # 构造上三角索引 i_idx, j_idx = np.triu_indices(n_years, k=1) dy = data[j_idx, :] - data[i_idx, :] # (n_pairs, n_pixels) dx = years[j_idx] - years[i_idx] # (n_pairs,) # Theil-Sen斜率：逐对斜率取中位数 pair_slopes = dy / dx[:, None] slope_valid = np.median(pair_slopes, axis=0) # Mann-Kendall S统计量 sign_pair = np.sign(dy) s_valid = np.sum(sign_pair, axis=0) # 方差（平局校正） n = n_years var_base = n * (n - 1) * (2 * n + 5) / 18.0 sorted_data = np.sort(data, axis=0) tie_correction = np.zeros(data.shape[1]) for col in range(data.shape[1]): count = 1 for row in range(1, n): if sorted_data[row, col] == sorted_data[row - 1, col]: count += 1 else: if count > 1: t = count tie_correction[col] += t * (t - 1) * (2 * t + 5) count = 1 if count > 1: t = count tie_correction[col] += t * (t - 1) * (2 * t + 5) var_valid = var_base - tie_correction / 18.0 # 标准化Z统计量 z_valid = np.zeros(data.shape[1]) pos = (s_valid > 0) & (var_valid > 0) neg = (s_valid < 0) & (var_valid > 0) z_valid[pos] = (s_valid[pos] - 1) / np.sqrt(var_valid[pos]) z_valid[neg] = (s_valid[neg] + 1) / np.sqrt(var_valid[neg]) # 双侧p值 p_valid = 2 * (1 - stats.norm.cdf(np.abs(z_valid))) p_valid[~(pos | neg)] = 1.0 # S=0或方差为0时无显著趋势 # 填回原位置 slope = np.full(rows * cols, np.nan) z = np.full(rows * cols, np.nan) p = np.full(rows * cols, np.nan) slope[valid] = slope_valid z[valid] = z_valid p[valid] = p_valid return slope.reshape(rows, cols), z.reshape(rows, cols), p.reshape(rows, cols)

这个向量化版本有个地方需要说明：平局校正我用了逐像元的小循环。如果像元数量上千万，这个小循环会成为瓶颈。但n通常只有20-40年，内部循环规模很小，实测几百万像元几分钟内也能跑完，具体看机器性能。如果你处理的栅格特别大，比如全国范围的500m数据，建议用rasterio.windows分块读取，再配合xarray加dask做延迟计算。把趋势分析函数包装成逐块处理，内存占用可以从几十G降到2-3G，这是处理大区域数据的标准姿势。

4.3 更快的方式：用Numba加速循环

如果你想进一步提速，Numba是个很好的选择。用@njit装饰一个逐像元计算的循环，编译器会把Python代码编译成机器码，速度能提升几十倍。对于标准NPP趋势分析，n较小，计算瓶颈主要在IO上，Numba版本可以有效处理大栅格。

from numba import njit from math import erf, sqrt @njit def ts_slope_mk_numba(years, values): n = len(values) # 有效值处理 valid_count = 0 for k in range(n): if not np.isnan(values[k]): valid_count += 1 if valid_count < 3: return np.nan, np.nan, np.nan y = np.empty(valid_count) x = np.empty(valid_count) cnt = 0 for k in range(n): if not np.isnan(values[k]): y[cnt] = values[k] x[cnt] = years[k] cnt += 1 # Theil-Sen斜率 m = valid_count n_pair = m * (m - 1) // 2 slopes = np.empty(n_pair) idx = 0 for i in range(m): for j in range(i + 1, m): if x[j] != x[i]: slopes[idx] = (y[j] - y[i]) / (x[j] - x[i]) idx += 1 if idx == 0: return np.nan, np.nan, np.nan slope = np.median(slopes[:idx]) # Mann-Kendall S统计量 s = 0 for i in range(m - 1): for j in range(i + 1, m): if y[j] > y[i]: s += 1 elif y[j] < y[i]: s -= 1 # 平局校正 y_sorted = np.sort(y) var_ties = 0.0 cnt_run = 1 for k in range(1, m): if y_sorted[k] == y_sorted[k - 1]: cnt_run += 1 else: if cnt_run > 1: var_ties += cnt_run * (cnt_run - 1) * (2 * cnt_run + 5) cnt_run = 1 if cnt_run > 1: var_ties += cnt_run * (cnt_run - 1) * (2 * cnt_run + 5) var_s = m * (m - 1) * (2 * m + 5) / 18.0 - var_ties / 18.0 if var_s <= 0: return slope, 0.0, 1.0 # 标准化Z统计量 if s > 0: z = (s - 1) / sqrt(var_s) elif s < 0: z = (s + 1) / sqrt(var_s) else: z = 0.0 # 标准正态分布近似 cdf = 0.5 * (1 + erf(abs(z) / sqrt(2.0))) p = 2 * (1 - cdf) return slope, z, p @njit(parallel=True) def trend_analysis_numba(stack, years): n_years, rows, cols = stack.shape slope = np.full((rows, cols), np.nan) z = np.full((rows, cols), np.nan) p = np.full((rows, cols), np.nan) for r in range(rows): for c in range(cols): vals = stack[:, r, c] if np.all(np.isnan(vals)): continue slope[r, c], z[r, c], p[r, c] = ts_slope_mk_numba(years, vals) return slope, z, p

需要注意，Numba的parallel=True对嵌套循环并行化有版本要求。实测中如果单像元计算量不大，并行收益未必明显，反而可能因为调度开销变慢。对大栅格，我建议优先用numba.prange对行循环做并行，再配合分块读取数据，这样的组合效果最好。

4.4 结果分类与制图

计算完slope、z、p三个栅格后，把结果组合成分类图。常用的五级分类方案如下：

用numpy就可以一句话完成分类：

def classify_trend(slope, p): out = np.zeros(slope.shape, dtype=np.int16) sig = p < 0.05 e_sig = p < 0.01 up = slope > 0 down = slope < 0 out[(up & e_sig)] = 2 out[(up & sig & ~e_sig)] = 1 out[(down & sig & ~e_sig)] = -1 out[(down & e_sig)] = -2 out[np.isnan(slope) | np.isnan(p)] = -32768 return out

分类完成后，用matplotlib出图。我习惯用离散色带，红色系代表退化、绿色系代表改善、浅色代表不显著变化。出图时保留研究区边界和经纬度坐标，比例尺、指北针、图例一个不能少，这是科研制图的基本要求。

import matplotlib.pyplot as plt import matplotlib.colors as mcolors from matplotlib.patches import Patch import numpy.ma as ma cat_colors = [ (0.6, 0.0, 0.0), # -2 极显著退化 (0.9, 0.6, 0.2), # -1 显著退化 (0.9, 0.9, 0.8), # 0 不显著变化 (0.5, 0.8, 0.2), # 1 显著改善 (0.0, 0.4, 0.0), # 2 极显著改善 ] labels = ['极显著退化', '显著退化', '不显著变化', '显著改善', '极显著改善'] cmap = mcolors.ListedColormap(cat_colors) norm = mcolors.BoundaryNorm([-2.5, -1.5, -0.5, 0.5, 1.5, 2.5], cmap.N) fig, ax = plt.subplots(figsize=(10, 8)) # 用掩膜数组，避免NaN被当作0参与配色 result_masked = ma.masked_invalid(result) im = ax.imshow(result_masked, cmap=cmap, norm=norm) legend_handles = [Patch(color=cat_colors[i], label=labels[i]) for i in range(5)] ax.legend(handles=legend_handles, loc='lower right', frameon=False) ax.set_title('2000-2021年NPP变化趋势分类') plt.savefig('npp_trend_class.png', dpi=300, bbox_inches='tight')

这里有个制图细节很容易被忽略：如果直接用imshow显示带缺失值的数组，缺失值会被当成0参与映射，显示成“不显著变化”的颜色，这会完全误导读图人。务必用numpy.ma.masked_invalid把无效值掩膜掉，或者填特殊值后用cmap.set_bad()单独设置颜色。

5. 结果解读：趋势分类、统计分析与科研表达

5.1 趋势分级怎么统计和描述

栅格趋势图出来后，不能只说一句“有显著变化”。审稿人和导师更关心的是：显著改善的面积占总研究区的百分之多少？空间上集中在哪些区域？不同土地覆盖类型上趋势有没有差异？

统计面积比例的标准做法是：按分类结果逐类统计有效像元数量，再乘以单个像元面积。如果用的是0.05°栅格，单个像元面积在每个纬度上是变化的，不能简单用固定面积乘，要先生成纬度面积权重栅格，或者把栅格转成等积投影后再统计。我在第一次做全国分析时直接用经纬度栅格数像元，得出的面积比例偏差了好几倍，后来重投影成Albers等积投影才纠正过来。

统计表一般长这样：

这种表格几乎是每篇NPP趋势分析论文的标配。除了面积统计，还可以进一步做分区统计（按省份、流域、生态区），或者按土地覆盖类型提取趋势值，观察哪一类生态系统的生产力在提升或退化，这一步对生态政策评估特别有用。

5.2 NPP趋势结果背后的生态学解释

得到趋势结果后，最难的一步其实是解释。NPP上升不一定全是好事，NPP下降也不一定全是坏事，必须结合研究区的气候背景和人类活动来讨论。

比如在退耕还林还草工程区，NPP显著增加的像元往往集中在坡耕地退耕区域，这反映植被恢复成效。在干旱区，如果某年降水异常偏多，NPP也会出现一次高位脉冲，但这不代表生态系统发生了根本性好转。同样的地形和气候条件下，灌溉农田的NPP趋势可能非常平稳，而天然草地则表现出强烈的年际波动，MK检验对“波动大的序列”尤其容易判为不显著，这是非参数检验的特性，解释时要特别小心。

还有一点很关键：NPP趋势分析是单一指标，不要过度解读成“生态系统健康状况”。NPP高可能只是意味着生物量大，不一定是生物多样性高或生态系统服务强。在论文里表述时，我一般写成“植被生产力呈上升趋势”，而不是“生态系统明显改善”。

5.3 论文里的常见表达参考

在科研论文里，方法描述部分我习惯这样写，你可以直接参考：

本研究运用Theil-Sen Median斜率估计方法计算每个像元的NPP变化速率。该方法通过对时间序列所有点对斜率取中位数，能有效抑制异常值对趋势估计的干扰。同时采用Mann-Kendall非参数检验评估趋势的统计显著性，其统计量S基于序列内所有数据对的大小比较构建，在零假设下近似服从正态分布，可据此计算标准化统计量Z和显著性水平p。当p < 0.05时认为趋势达到显著水平。

方法段落写清楚“用了什么、为什么用、显著性标准是什么”三部分，就足够规范了。很多期刊对方法部分的要求就是“可复现”，把参数写明白，别人才能照着做。

6. 常见问题与排错技巧：我实测踩过的那些坑

6.1 常见问题速查表

6.2 我踩过的三个坑

第一个坑：忘记平局校正。有一回我用某个省的NPP数据跑趋势分析，算了十个像元的p值，发现所有结果都异常地显著，p值普遍小于0.001。一查原因，是数据里大量年份NPP相同（产品做了整型量化），但我的方差计算公式里没有减平局校正项。补上校正后，p值立刻回归合理范围。这是教科书里很容易被忽略、但实操里最要命的细节。

第二个坑：把“年份”当成序号直接用。Theil-Sen斜率的计算公式里，分母是年份差，不是序号差。如果数据是2001、2003、2005这种间隔为2年的序列，直接用1、2、3当分母，斜率会被放大2倍。我一直强调，写代码时年份一定要作为真实数值传入，不要用range(len(years))替代。

第三个坑：掩膜顺序。我早期习惯先做趋势分析，再对结果做掩膜。这在大部分情况下没大问题，但在NPP常年为0的像元上会浪费大量计算，而且会产生“0值像元趋势为0、判定为不显著”的假象。正确做法是在数据读取阶段就把无效区域设成NaN，让趋势分析直接跳过，这样又省时间又避免污染结果。

6.3 最后再分享一个实用小技巧

做完整套分析后，我建议顺手输出一个“像元数-斜率”直方图和一个“显著像元占比柱状图”。这两个图放在结果第一页，能帮你快速检查趋势结果是否合理。比如斜率分布如果出现明显的双峰，可能说明研究区里有两种完全不同变化模式的地类；显著像元占比如果超过50%，就要怀疑是不是没有做平局校正或者掩膜没到位。

另外，如果你要发表论文，记得把所有中间结果（每年的NPP平均图、趋势斜率图、显著性图、分类图）都存成带坐标系的GeoTIFF，命名规则统一。过三个月再返工的时候，你会发现当初这点“举手之劳”能救你一条命。我自己就是因为当年随手存了带地理信息的中间文件，后来补分析、改配色、换分类阈值时省了整整两天时间。