基于Bayes-Copula方法的用于植被脆弱性评估-栅格版

有梦想的Frank博士

已于 2025-03-10 11:26:38 修改

阅读量1.3k

点赞数 9

CC 4.0 BY-SA版权

分类专栏：数据分析数据绘图文章标签： copula 贝叶斯阈值敏感性

于 2025-02-25 12:08:12 首次发布

本文链接：https://blog.youkuaiyun.com/frankgis/article/details/145840893

数据分析同时被 2 个专栏收录

26 篇文章

订阅专栏

数据绘图

21 篇文章

订阅专栏

前言

随着气候变化的加剧，干旱等极端气候事件对生态系统，特别是植被生长的影响日益显著。干旱不仅直接影响水资源的可用性，还对农业生产、生态系统服务及人类社会的可持续发展构成了巨大的威胁。在全球范围内，植被脆弱性评估成为了研究气候变化影响的重要领域，尤其是在评估干旱对植被的影响方面，科学的评估方法对于制定应对策略至关重要。

传统的植被脆弱性评估方法多侧重于单一因素的分析，忽视了多种气候因子之间复杂的相互作用。近年来，随着统计学和机器学习方法的发展，基于联合概率分布的评估方法逐渐成为干旱影响评估的热门工具。Copula-Bayes方法作为一种结合了Copula函数与贝叶斯推理的统计方法，能够有效地处理不同气候因子之间的复杂关联，尤其是在涉及多变量联合分布的情境下表现出强大的灵活性与精准性。

Copula函数通过描述多变量之间的依赖关系，突破了传统相关分析的局限，使得干旱条件下植被损失的评估更加精细。贝叶斯推理则能够基于先验知识与观测数据，进行动态更新与优化，从而提高模型的预测能力和鲁棒性。在这一框架下，结合栅格数据的空间分布特征，可以在高分辨率的尺度上准确识别不同地区、不同条件下植被的脆弱性，进而为干旱的预警与应对提供科学依据。

本文基于Copula-Bayes方法，构建了一个用于植被脆弱性评估的框架，并进一步应用于栅格数据分析。通过精确评估不同干旱条件下植被的损失概率及其空间分布，本文不仅揭示了植被对气象干旱的脆弱性，还为全球气候变化背景下的生态风险评估与决策提供了有力支持。

一、原理

Copula函数是一种能够将两个变量之间的联合分布与各自的边际分布连接起来的数学工具。该方法已被广泛应用于评估植被生产力或农作物产量损失的概率。本研究基于Copula方法构建了干旱阈值评估框架，用以触发植被损失的发生。

我们分别界定了干旱的累积效应和滞后效应的时间尺度。进一步建立了与Rmax_cum和Rmax_lag时间尺度对应的SPEI与植被数据的联合分布模型。其计算公式如下：
在这里插入图片描述

其中，SPEIR表示在每个像素点上对应时间尺度的干旱累积值或滞后干旱值，FSPEI,VIs为SPEI与VIs之间的联合分布函数，而FSPEI和FVIs则分别表示SPEI和VIs的边际分布函数。在植被数据方面，我们选用了高斯分布、伽马分布和逻辑斯蒂分布作为候选的边际分布类型；对于SPEI数据，我们选择了高斯分布。C()表示在每个像素点上拟合的最佳Copula函数。为了找到最优的Copula函数，我们选择了高斯型、斯图登特t型、Clayton型、Gumbel型、Frank型和Joe型等六种候选函数，并通过最大似然估计法对参数进行了估算。使用Akaike信息准则（AIC）和Kolmogorov-Smirnov检验对各候选函数的分布类型进行了选择，以确定最佳模型。

接着，我们将选定的Copula函数与贝叶斯条件概率公式相结合，计算干旱导致植被损失的条件概率，并通过百分位数来量化植被损失的严重程度。具体地，我们定义VIs < 40百分位为轻度损失，VIs < 10百分位为重度损失。然后，我们计算在不同干旱条件下植被发生轻度或重度损失的概率。公式如下：
在这里插入图片描述
其中，X1和X2分别是给定SPEI区间的上限和下限。推导在某一程度干旱情景下，植被状态低于某一特定阈值（例如：低于长期生长季序列的第30、20、10百分位数）时的损失概率。这一推导将为后续的脆弱性评估提供可靠的数据支持。比如本研究的工作流程详见下图。SPEI值的通常范围为−3到3，当SPEI值低于−0.5时，表示干旱发生。因此，我们将X1的初始值设为−0.5，作为该范围的上限，并以0.1为间隔进行迭代，计算每个区间内植被损失的概率。如果某个区间内植被损失的概率大于0.5，则认为该区间的干旱阈值触发植被损失；如果所有区间的概率都小于0.5，则认为该像素点没有干旱阈值。在这一过程中，得出了因干旱的累积效应与滞后效应导致的植被损失的干旱阈值空间分布。在这里插入图片描述
用于计算导致植被损失的干旱阈值的工作流程图
考虑到植被与气候因子之间的紧密关联，植被的变化不仅受到当前气候状况的影响，也受到前期气候条件的影响，因此植被对气候的响应具有一定的滞后性。此外，不同地区的植被类型与气候条件各异，导致植被对气象因子的响应时段不同。一般而言，植被生长所需的水分主要来自降水和土壤水，尤其是在植被生长季节（5至9月）。因此，本研究重点探讨了气象干旱和农业干旱胁迫下，植被生长季节的损失。

二、代码

部分条件概率代码

clc; 
clear all;
%matlab是R2024b 版本 可以在https://mp.weixin.qq.com/s/5IutbJm-5o_SZLwpH1iA8w安装

% 修改为SPEI数据路径
imgdirX = 'D:\Codes\Geo_SPEI\'; 
addpath(genpath(imgdirX));

% 读取SPEI数据
filenamesX = dir([imgdirX '*.tif']);
for ii = 1:numel(filenamesX)
    A = imread(filenamesX(ii).name);  % 读取tif文件
    X(:,:,ii) = A;  % 将每张图像按顺序存储到三维矩阵X中
end

% 修改为NDVI数据路径
imgdirY = 'D:\Codes\GEO_NDVI\'; 
addpath(genpath(imgdirY));

% 读取NDVI数据
filenamesY = dir([imgdirY '*.tif']);
for ii = 1:numel(filenamesY)
    [da,R]= geotiffread(filenamesY(ii).name);
    A = imread(filenamesY(ii).name);  % 读取tif文件
    Y(:,:,ii) = A;  % 将每张图像按顺序存储到三维矩阵Y中
end

% 去除异常值
X(X == -9999) = nan;
Y(Y == -9999) = nan;
X(isinf(X)) = nan;
Y(isinf(Y)) = nan;

[rows, cols, depth] = size(X);

% 请确认数据的行列是否正确，若反了需要自己翻转
% 初始化一个空矩阵，用于存储结果数据
dryhot = nan(rows, cols);  % 初始化结果矩阵

% % 计算前5%的数据
row_end = round(rows);   
col_end = round(cols);  

% 循环遍历的像元
%这里是并行计算
parfor i = 1:row_end
    for j = 1:col_end

        p = X1(:);p=double(p);
        t = Y1(:);t=double(t);

        if all(isnan(p)) || all(isnan(t))||numel(p) == sum(p == p(1)) || numel(t) == sum(t == t(1))
            PP = nan;  % 如果全部为NaN，则结果为NaN
        else
            index = ~isnan(p) & ~isnan(t);  % 过滤掉NaN值
            X1 = p(index);
            Y1 = t(index);

            % 阈值计算,改这里即可
            % 表示中度干旱-1.5 < SPEI < -1
            phighth = -1;
            plowth = -1.5;
            tth = quantile(Y1, 0.4);  % 表示植被的40%分位数
  try
            %[D_U1, PD_U1] = marginfitdist(X1);  % 拟合X1的边缘分布

            
            % 拟合最佳边缘分布
            EP1 = cdf(PD_U1{1}, X1);  % 计算X1的CDF
            EP2 = cdf(PD_U2{1}, Y1);  % 计算Y1的CDF

            % 阈值的分布计算
            phighth = cdf(PD_U1{1}, phighth);
            plowth = cdf(PD_U1{1}, plowth);
            tth = cdf(PD_U2{1}, tth);

            EP12 = [EP1 EP2];
            EP12(EP12 <= 0) = 1e-4;  % 防止CDF为0，避免后续计算错误
            [Family1, thetahat1, loglik1] = doublecopulaselect(EP12, 'aic');  % 选择最佳的copula模型

            % 设置阈值
            EPpthigh = [phighth tth];
            EPptlow = [plowth tth];
            EPpthigh(EPpthigh <= 1e-4)=1e-4;
            EPptlow(EPptlow <= 1e-4)=1e-4;
            % 计算二维的条件概率
            %表示大于一个数和小于一个数的范围比如中度干旱（-1.5 < SPEI ≤ -1.0）
            ptor = (copulascdf(EPpthigh, Family1, thetahat1) - copulascdf(EPptlow, Family1, thetahat1)) / (phighth - plowth);
            %表示小于一个数范围，比如 SPEI ≤ -1.5
            %ptor =copulascdf(EPpthigh, Family1, thetahat1)/ phighth;
            %表示大于一个数范围，比如SPEI > 1.5
            %ptor = (EPpthigh-copulascdf(EPpthigh, Family1, thetahat1)) / (1 - plowth);
            
          % 如果计算失败，设置PP为NaN
                if isnan(ptor) || isinf(ptor)
                    PP = nan;  % 如果ptor为NaN或无穷大，跳过该像元
                else
                    PP = ptor;  % 计算结果
                end

            catch
                % 如果在try块中发生错误，设置PP为NaN并跳过该像元
                PP = nan;
            end
        end
        %disp(PP)  % 显示计算结果
        dryhot(i, j) = PP;  % 将计算结果存储到结果矩阵
    end
    
    % 更新进度信息
    fprintf('Processing row %d of %d...\n', i, row_end);  % 显示当前行处理进度
end

% 将结果矩阵dryhot保存为tif格式
outputDir = 'D:\desktop\result\';  % 修改为你希望保存结果的路径
outputFileName = 'Result.tif';  % 保存的文件名
outputPath = fullfile(outputDir, outputFileName);

% 确保输出文件夹存在
if ~exist(outputDir, 'dir')
    mkdir(outputDir);  % 如果文件夹不存在，创建它
end

% 由于dryhot可能包含NaN值，首先用nan填充，使用uint16格式保存数据
dryhot(isnan(dryhot)) = 0;  % 将NaN替换为0，避免保存为tif时出现问题
% 保存为.tif文件
geotiffwrite(outputPath,dryhot,R);
disp(['Result saved to: ', outputPath]);

2.部分阈值求解代码

clc; 
clear all;
%matlab是R2024b 版本 可以在https://mp.weixin.qq.com/s/5IutbJm-5o_SZLwpH1iA8w安装

% 修改为SPEI数据路径
imgdirX = 'D:\Codes\Geo_SPEI\'; 
addpath(genpath(imgdirX));

% 读取SPEI数据
filenamesX = dir([imgdirX '*.tif']);
for ii = 1:numel(filenamesX)
    A = imread(filenamesX(ii).name);  % 读取tif文件
    X(:,:,ii) = A;  % 将每张图像按顺序存储到三维矩阵X中
end

% 修改为NDVI数据路径
imgdirY = 'D:\Codes\GEO_NDVI\'; 
addpath(genpath(imgdirY));

% 读取NDVI数据
filenamesY = dir([imgdirY '*.tif']);

for ii = 1:numel(filenamesY)
    [da,R]= geotiffread(filenamesY(ii).name);
    A = imread(filenamesY(ii).name);  % 读取tif文件
    Y(:,:,ii) = A;  % 将每张图像按顺序存储到三维矩阵Y中
end

% 去除异常值
X(X == -9999) = nan;
Y(Y == -9999) = nan;
X(isinf(X)) = nan;
Y(isinf(Y)) = nan;

[rows, cols, depth] = size(X);

% 请确认数据的行列是否正确，若反了需要自己翻转
% 初始化一个空矩阵，用于存储结果数据
drythreshold = nan(rows, cols);  % 初始化结果矩阵
% 
row_end = round(rows);   
col_end = round(cols);  

parfor i = 1:row_end
    for j = 1:col_end
        p = X1(:); p = double(p);
        t = Y1(:); t = double(t);

        % 检查是否全部为NaN或者所有值相同
        if all(isnan(p)) || all(isnan(t)) || numel(p) == sum(p == p(1)) || numel(t) == sum(t == t(1))
            PP = nan;  % 如果全部为NaN，或者数据中所有值相同，则PP设为NaN
        else
            index = ~isnan(p) & ~isnan(t);  % 过滤掉NaN值
            X1 = p(index);
            Y1 = t(index);

            % 阈值计算
            % phighth = -0.2;
            % plowth = -5;
            mean_X1 = mean(X1, 'omitnan');
            std_X1 = std(X1, 'omitnan');
            mean_Y1 = mean(Y1, 'omitnan');
            std_Y1 = std(Y1, 'omitnan');

            % 设置阈值范围为数据的均值±标准差范围
            plowth = mean_X1 - 2*std_X1;
            phighth = mean_X1 + 2*std_X1;
            tth = quantile(Y1, 0.1);  % 植被的40%分位数
            try
                % 使用 fitdist 进行正态分布拟合
                pd_X1 = fitdist(X1, 'Normal');  % 拟合 X1 数据为正态分布
                D_U1 = struct('DistName', 'Normal', 'Params', pd_X1.Params);  % 保存拟合结果
                PD_U1 = {pd_X1};  % 将拟合后的正态分布对象保存到 PD_U1

                [D_U2, PD_U2] = marginfitdist(Y1);  % 拟合Y1的边缘分布

                % 拟合最佳边缘分布
                EP1 = cdf(PD_U1{1}, X1);  % 计算X1的CDF
                EP2 = cdf(PD_U2{1}, Y1);  % 计算Y1的CDF

                % 阈值的分布计算
                phighthcdf = cdf(PD_U1{1}, phighth);
                plowthcdf = cdf(PD_U1{1}, plowth);
                tth = cdf(PD_U2{1}, tth);

                EP12 = [EP1 EP2];
                EP12(EP12 <= 0) = 1e-4;  % 防止CDF为0，避免后续计算错误
      

                for spei = plowth:0.1:phighth+0.1
                    % 设置阈值
                    phighthcdf = cdf(PD_U1{1}, spei);
                    EPpthigh = [phighthcdf tth];
                    EPptlow = [plowthcdf tth];
                    EPpthigh(EPpthigh <= 1e-4) = 1e-4;
                    EPptlow(EPptlow <= 1e-4) = 1e-4;

                    % 计算二维的条件概率
                    ptor = (copulascdf(EPpthigh, Family1, thetahat1) - copulascdf(EPptlow, Family1, thetahat1)) / (phighthcdf - plowthcdf);

                    % 如果 ptor > 0.5，则退出循环并存储SPEI
                    if ptor >= 0.5
                        PP = spei;  % 找到符合条件的SPEI
                        found = true;
                        break;
                    else
                        % 找到最接近0.5的ptor
                        diff = abs(ptor - 0.5);
                        if diff < min_diff
                            min_diff = diff;
                            closest_spei = spei;  % 存储最接近的SPEI值
                        end
                    end
                end

                % 如果没有找到符合条件的值，则使用最接近0.5的SPEI值
                if ~found
                    PP = closest_spei;
                end

            catch
                % 如果还是没找到，设定PP为NaN
                PP = nan;
            end
        end

        disp(PP)  % 显示计算结果
        drythreshold(i, j) = PP;  % 将计算结果存储到结果矩阵
    end
    
    % 更新进度信息
    fprintf('Processing row %d of %d...\n', i, row_end);  % 显示当前行处理进度
end


% 将结果矩阵dryhot保存为tif格式
outputDir = 'D:\desktop\result\';  % 修改为你希望保存结果的路径
outputFileName = 'Resultthreshold.tif';  % 保存的文件名
outputPath = fullfile(outputDir, outputFileName);

% 确保输出文件夹存在
if ~exist(outputDir, 'dir')
    mkdir(outputDir);  % 如果文件夹不存在，创建它
end

% 由于dryhot可能包含NaN值，首先用nan填充，使用uint16格式保存数据
drythreshold(isnan(drythreshold)) = 0;  % 将NaN替换为0，避免保存为tif时出现问题
% 保存为.tif文件
geotiffwrite(outputPath,drythreshold,R);
disp(['Result saved to: ', outputPath]);

三.出图效果

以上代码只是片段，完整代码获取请关注公众号趣品科研回复BC
在这里插入图片描述
在中度气象或农业干旱胁迫下不同程度植被损失的条件概率

累积干旱（a，b）和时滞干旱（c，d）触发的植被损失阈值的空间分布

总结

本研究主要探讨了卫星衍生的太阳辐射驱动因子（SIF）和三种植被指数（VI）对气象干旱的累积效应与时间滞后效应的响应。研究发现，干旱对植被的影响呈现出显著的累积和滞后效应，其中，植被对干旱的累积效应和滞后效应的平均响应时间尺度主要集中在4至6个月之间。

值得注意的是，干旱的时间滞后效应比累积效应对植被的影响更为显著，这表明植被的响应不仅受到当前干旱状况的影响，还受到前期干旱的延续性作用。这种滞后效应在干旱的影响下，导致植被损失发生的时间和强度出现一定的滞后。

此外，研究还发现植被在不同水分梯度下对干旱的响应呈现出非线性特征。即随着水分状况的变化，植被对干旱的反应方式也发生了变化，这表明水分对植被生长的影响是多样而复杂的。

在干旱阈值方面，研究结果表明，干旱累积效应和滞后效应共同作用时，触发轻度植被损失的干旱阈值为SPEI值在−0.7至−1.0之间，而触发严重植被损失的干旱阈值则位于SPEI值的−2.0至−2.5区间。这些阈值为干旱的影响程度提供了量化的参考，能够帮助评估和预测不同强度干旱对植被的潜在影响。

综上所述，本研究揭示了干旱对植被的复杂影响机制，尤其是累积效应与滞后效应的不同作用，以及干旱阈值在触发不同程度植被损失中的重要作用。

参考文献
Shixian Xu, Yonghui Wang, Yuan Liu, Jiaxin Li, Kaixuan Qian, Xiuyun Yang, Xiaofei Ma,
Evaluating the cumulative and time-lag effects of vegetation response to drought in Central Asia under changing environments,
Journal of Hydrology,
Volume 627, Part B,
2023,
130455,
ISSN 0022-1694,
https://doi.org/10.1016/j.jhydrol.2023.130455.
韩知明.中国多类型干旱时空演变特征及其传播过程研究[D].西安理工大学,2022.DOI:10.27398/d.cnki.gxalu.2022.001707.