基于 Python 的栅格数据 Theil-Sen 与 Mann-Kendall 趋势分析及显著性检验

在气候变化、生态环境等研究中，如何准确分析长时间序列数据的趋势至关重要。Theil-Sen Median 趋势分析（Sen 分析）和 Mann-Kendall 显著性检验（MK 检验）是两种常用的非参数方法，特别适合处理含有异常值或缺失值的遥感数据。本文将梳理这两种方法的原理，并提供完整的 Python 代码实现。

核心方法解析

1. Theil-Sen Median 趋势分析（Sen 分析）

这是一种稳健的非参数统计方法。不同于传统的最小二乘法，Sen 分析通过计算数据的中位数斜率来评估时间序列的趋势，对异常值和极端值不敏感。

具体计算以蒸散发（ET）为例：

其中，ETi 和 ETj 分别为时间序列中不同时间点的数据值，Sen 为中位数斜率。若 Sen > 0，表示上升趋势；若 Sen < 0，表示下降趋势。

2. Mann-Kendall 显著性检验（MK 检验）

MK 检验用于判断时间序列的趋势是否显著，不要求数据服从特定分布。它通过对每对数据点进行比较，计算趋势的符号，并利用统计量评估显著性。

基本步骤如下：

1. 比较所有数据点对，得到符号（升序为 1，降序为 -1，相等为 0）；
2. 计算统计量 S，反映总体趋势；
3. 根据 S 及其方差计算 Z 值，判断显著性。

Z 值越大，趋势越显著。符号函数累加和得到统计量 S：

3. MK 显著性结果划分

在 0.05 置信水平下，依据 Z 值判断：

• 显著变化：上升趋势（Z > 1.96）或下降趋势（Z < -1.96）
• 不显著变化：-1.96 ≤ Z ≤ 1.96

4. 方法优势

• 稳健性强：有效处理含异常值或缺失值的数据。
• 显著性检验：确保趋势结果具有统计学意义。
• 适用范围广：适用于无法满足正态分布假设的时间序列数据。

代码实现

以下代码实现了逐像元栅格进度可视化，只需修改标注位置即可运行。代码逻辑上优化了空值处理和内存占用，适合长时间序列分析。

 os
 rasterio
 numpy  np
 tqdm  tqdm
 matplotlib.pyplot  plt
 pathlib  Path


base_path = os.getcwd()


data_path = os.path.join(base_path, )
result_path = os.path.join(base_path, )


os.makedirs(result_path, exist_ok=)


start_year = 
end_year = 
cd = end_year - start_year +   

:
    
    first_file = os.path.join(data_path, )
     rasterio.(first_file)  src:
        a = src.read()  
        transform = src.transform  
        metadata = src.meta.copy()  
        
        
        metadata.update({
            : ,
            : np.nan  
        })
        
        
        m, n = a.shape
        ()
        
        
        
        all_data = np.full((m, n, cd), np.nan)
        
        
        ()
         i, year  (tqdm((start_year, end_year + ), desc=)):
            filename = os.path.join(data_path, )
             rasterio.(filename)  src:
                all_data[:, :, i] = src.read()
        
        
        sen_result = np.full((m, n), np.nan)
        
        
        
        ()
        valid_pixels = 
        
        valid_mask = np.(all_data > , axis=)
        total_valid = np.(valid_mask)
        
         i  tqdm((m), desc=):
             j  (n):
                 valid_mask[i, j]:
                    data = all_data[i, j, :]
                     np.(~np.isnan(data))  np.(data) > :  
                        slopes = []
                         k1  (, cd):
                             k2  (k1):
                                
                                slope = (data[k1] - data[k2]) / (k1 - k2)
                                slopes.append(slope)
                        
                         slopes:
                            sen_result[i, j] = np.median(slopes)
                            valid_pixels += 
        
        ()
        
        
        sen_output_path = os.path.join(result_path, )
        
        
        sen_result = sen_result.astype()
        
        
         rasterio.(sen_output_path, , **metadata)  dst:
            dst.write(sen_result, )
        (, sen_output_path)
        
        
        ()
        mk_result = np.full((m, n), np.nan)
        
         i  tqdm((m), desc=):
             j  (n):
                 valid_mask[i, j]:
                    data = all_data[i, j, :]
                     np.(~np.isnan(data))  np.(data) > :
                        sgnsum = []
                         k1  (, cd):
                             k2  (k1):
                                
                                sgn = np.sign(data[k1] - data[k2])
                                sgnsum.append(sgn)
                        
                        mk_result[i, j] = np.(sgnsum)
        
        
        ()
        vars_mk = cd * (cd - ) * ( * cd + ) / 
        z_scores = np.full((m, n), np.nan)
        
        
        z_scores[~np.isnan(mk_result) & (mk_result == )] = 
        z_scores[~np.isnan(mk_result) & (mk_result > )] = (mk_result[~np.isnan(mk_result) & (mk_result > )] - ) / np.sqrt(vars_mk)
        z_scores[~np.isnan(mk_result) & (mk_result < )] = (mk_result[~np.isnan(mk_result) & (mk_result < )] + ) / np.sqrt(vars_mk)
        
        
        mk_output_path = os.path.join(result_path, )
        
        z_scores = z_scores.astype()
         rasterio.(mk_output_path, , **metadata)  dst:
            dst.write(z_scores, )
        (, mk_output_path)
        
        
        ()
        
        S2 = np.full((m, n), np.nan)
        M2 = np.full((m, n), np.nan)
        
        
        S2[np.isnan(sen_result)] = -
        S2[~np.isnan(sen_result) & (sen_result <= -)] = -
        S2[~np.isnan(sen_result) & (sen_result >= )] = 
        S2[~np.isnan(sen_result) & (sen_result > -) & (sen_result < )] = 
        
        
        M2[z_scores > ] = 
        M2[~np.isnan(z_scores) & (z_scores <= )] = 
        
        
        reclassify = (S2 * M2).astype(np.int16)
        
        
        reclass_output_path = os.path.join(result_path, )
        
        
        metadata.update({
            : ,
            : -
        })
        
        
         rasterio.(reclass_output_path, , **metadata)  dst:
            dst.write(reclassify, )
        (, reclass_output_path)
        
 Exception  e:
    ()