生物科学2001—2022年黄土高原不同植被对气象干旱的响应

杨梅焕 ,  蒋永涛 ,  王涛 ,  张靖铷 ,  王博 ,  高波 ,  聂磊 ,  陶冰玉 ,  罗状 ,  杨昊伟

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (18) :7601-7611

PDF (19034KB)
科学技术与工程 ›› 2026, Vol. 26 ›› Issue (18) :7601-7611 DOI: 10.12404/j.issn.1671-1815.2505934
生物科学

生物科学2001—2022年黄土高原不同植被对气象干旱的响应

作者信息 +

Responses of Different Vegetation Types to Meteorological Drought in the Loess Plateau from 2001 to 2022

Author information +
文章历史 +

摘要

频繁的干旱事件已对植被正常生长与生态系统稳定构成威胁。为解析干旱对植被影响的复杂性,以2001—2022年黄土高原逐月植被总初级生产力(gross primary productivity, GPP)和标准化降水蒸散发指数(standardized precipitation and evapotranspiration index, SPEI)为基础,采用Copula函数利用去趋势化后的植被GPPdt和SPEI构建了植被生长季气象干旱胁迫下黄土高原植被总初级生产力损失的条件概率模型,并确定了导致植被总初级生产力损失的干旱触发阈值。结果表明:2001—2022年黄土高原植被生长季GPPdt与SPEI以正相关为主(95.05%),植被对干旱的响应时间集中在2~6个月,响应时间由短到长依次为稀疏植被(3.89个月)、草本植被(5.78个月)、落叶阔叶林(5.80个月)和常绿针叶林(6.38个月);对相同植被GPPdt损失发生阈值,干旱等级越高,植被GPPdt损失发生的条件概率越大。相同干旱等级下,植被GPPdt损失发生阈值越高,植被GPPdt损失发生的条件概率越大;黄土高原西北部和中部干旱触发阈值较高,对干旱响应更为敏感;以2001—2022年植被生长季GPPdt40th作为植被GPPdt损失发生阈值,植被GPPdt损失的平均概率由大到小依次为稀疏植被(52.21%)、草本植被(48.49%)、落叶阔叶林(47.77%)和常绿针叶林(42.78%),干旱触发阈值依次为稀疏植被(SPEI=-1.486)、草本植被(SPEI=-1.573)、落叶阔叶林(SPEI=-1.601)和常绿针叶林(SPEI=-1.676)。研究结果为黄土高原的植被生态保护和干旱防御策略提供科学支撑。

Abstract

Frequent drought events are recognized as threats to the normal growth of vegetation and the stability of ecosystem in the Loess Plateau. To elucidate the intricate effects of drought on vegetation, monthly GPP (gross primary productivity) and SPEI (standardized precipitation evapotranspiration index) from 2001 to 2022 were used to explore vegetation responses to meteorological drought. Using GPPdt (detrended GPP) and SPEI, copula functions were applied to construct a conditional probability model for vegetation GPP loss under growing-season meteorological drought stress, and drought-trigger thresholds that lead to GPPdt loss were determined. The results show as follows. During 2001—2022, GPPdt exhibits a predominantly positive correlation with SPEI (95.05%) during the growing season. Vegetation responses to drought are concentrated within 2~6 months, with the response time lengthening from sparse vegetation (3.89 months) to herbaceous vegetation (5.78 months), deciduous broadleaf forest (5.80 months), and evergreen needleleaf forest (6.38 months). Given the same GPPdt-loss threshold, more severe drought conditions lead to higher conditional probabilities of GPPdt loss. Under identical drought level, higher GPPdt-loss thresholds also correspond to increased conditional probabilities. Elevated drought-trigger thresholds are observed in the northwest and central region of the Loess Plateau, which indicate heightened drought sensitivity. Employing the 40th percentile of growing-season GPPdt from 2001 to 2022 as the loss threshold, the average probability of GPPdt loss declines from sparse vegetation (52.21%) to herbaceous vegetation (48.49%), deciduous broadleaf forest (47.77%), and evergreen needleleaf forest (42.78%). The drought-trigger thresholds for these vegetation types are SPEI = -1.486, -1.573, -1.601, and -1.676 respectively. The research findings provide scientific support for vegetation ecological protection and drought prevention strategies on the Loess Plateau.

Graphical abstract

关键词

条件概率 / 脆弱性评价 / 总初级生产力 / 标准化降水蒸散指数 / 黄土高原

Key words

conditional probability / vulnerability assessment / gross primary productivity / standardized precipitation and evapotranspiration index / Loess Plateau

引用本文

引用格式 ▾
杨梅焕, 蒋永涛, 王涛, 张靖铷, 王博, 高波, 聂磊, 陶冰玉, 罗状, 杨昊伟. 生物科学2001—2022年黄土高原不同植被对气象干旱的响应[J]. 科学技术与工程, 2026, 26(18): 7601-7611 DOI:10.12404/j.issn.1671-1815.2505934
YANG Mei-huan, JIANG Yong-tao, WANG Tao, ZHANG Jing-ru, WANG Bo, GAO Bo, NIE Lei, TAO Bing-yu, LUO Zhuang, YANG Hao-wei. Responses of Different Vegetation Types to Meteorological Drought in the Loess Plateau from 2001 to 2022[J]. Science Technology and Engineering, 2026, 26(18): 7601-7611 DOI:10.12404/j.issn.1671-1815.2505934

登录浏览全文

4963

注册一个新账户 忘记密码

全球气候变暖及其导致的极端气候事件频发问题是中外学者关注的热点问题[1-2]。随着大气升温,大气饱和水汽压增大,可容纳水量增加,可能导致气象干旱强度和频率进一步上升,从而对植被生长产生影响[3]。干旱对植被的影响包括直接影响和间接影响,其中直接影响主要表现为干旱通过减少水分供应,抑制植物光合作用、生长和生物量积累[4-6],严重时导致植物损伤甚至死亡[7-8];间接影响表现为干旱通过降低光合作用效率和减少地上及地下生物量,削弱生态系统碳汇能力[9]。同时,植被生产力下降可能影响全球物质循环和能量流动过程,从而降低生物群落稳定性和生态系统恢复力[10]。因此,量化干旱胁迫下植被脆弱性概率(植被遭受干旱破坏的可能性)[11]对于植被生态系统保护和碳循环具有重要意义。
随着遥感技术快速发展,多种基于遥感数据的植被指数和干旱指数广泛应用于植被脆弱性研究[12]。植被脆弱性研究中,传统研究主要采用确定性方法[13],即通过相关性将植被指数与干旱指数联系起来,评估长期气候条件干旱对植被的影响[14-15]。如王小红等[16]利用植被总初级生产力(gross primary productivity, GPP)和标准化降水蒸散指数(standardized precipitation and evapotranspiration index, SPEI),评估了秦巴山区GPP的干旱脆弱性和干旱风险。Ding等[12]利用归一化差异植被指数(normalized difference vegetation index, NDVI)和帕默尔干旱指数(Palmer drought severity index, PDSI)的最大相关系数,评估了中国不同植被覆盖对干旱的响应。李阳阳等[17]利用NDVI和SPEI数据,结合线性回归模型,分析了干旱对黄河流域植被的影响。尽管已有研究为植被与干旱之间的关系提供了新的见解,但仅通过相关性、线性统计等方法,难以刻画变量间的非线性依赖关系。因此,将确定性和概率方法相结合,定量评估干旱胁迫下植被脆弱性概率成为当前研究的主要趋势[18-19]。在概率论与统计学中,Copula函数是用来描述多元随机变量依赖结构的一种数学工具。它将多元随机变量的边缘分布“连接”成联合概率分布,既能描述变量之间的整体相关结构,还能刻画不同变量之间的尾部依赖性[20-22],从而达到量化不同等级极端事件下植被脆弱性的目的。当前,Copula函数已广泛应用于植被对气候干旱响应阈值的研究[23-25]。如Won等[18]利用Copula函数构建了标准化降水指数(standardized precipitation index, SPI)和植被健康指数(vegetation health index, VHI)的双变量联合概率分布,量化不同气象干旱情景下植被的响应,并评估干旱胁迫下植被脆弱性。Gao等[26]基于NDVI和SPEI,采用Copula函数和条件概率函数,定量评估了蒙古高原不同NDVI损失情景下气象干旱触发阈值及其动态变化。Li等[27]运用Copula方法,基于NDVI和SPEI,量化了伊犁河流域植被干旱脆弱性。
黄土高原区是我国重要的农业区之一[28],也是同纬度干旱灾害最频发的地区之一[29]。干旱频发对植被生态可持续性和粮食安全构成极大威胁[30-31]。黄土高原大部分地区植被变化与干旱密切相关,且不同类型植被对干旱响应具有一定的时空分异性[32]。已有黄土高原地区植被对干旱响应研究多聚焦于二者的相关性方面,基于Copula函数构建黄土高原植被损失的条件概率模型,并从概率角度评估从气象干旱到植被损失的干旱触发阈值的研究报道较少。
因此,现以2001—2022年黄土高原逐月GPP和SPEI为基础,采用Copula函数利用去趋势化后的植被GPPdt和SPEI构建植被生长季(4—9月)[33]气象干旱胁迫下黄土高原植被总初级生产力损失的条件概率模型,并确定导致植被总初级生产力损失的干旱触发阈值,以期为黄土高原的植被生态保护和干旱防御策略提供科学支撑。

1 数据与方法

1.1 研究概况

黄土高原位于黄河中游(107°28'E~111°15'E, 35°21'N~39°34'N),涉及青海、甘肃、宁夏、内蒙古、山西、陕西、河南7个省区,总面积为62.46×104 km2,海拔介于85~8 210 m,如图1(a)所示。黄土高原由东南向西北气候类型依次为温带大陆性半湿润季风气候、温带大陆性半干旱气候、高原气候等,夏季炎热多雨,冬季寒冷干燥。多年平均年气温5.8~15.2 ℃,多年平均年降水量为108~919 mm,年降水量由东南向至西北递减。黄土高原植被类型主要以稀疏植被、草本植被(草原、草丛、草甸)、落叶阔叶林、常绿针叶林和栽培植被为主,共占总面积的87.95%,如图1(b)所示。其中稀疏植被、草本植被、落叶阔叶林、常绿针叶林和栽培植被分别占总面积的7.02%、46.21%、11.22%、2.70%和20.80%。由于栽培植被受人类活动影响较大,在一定程度上能缓冲干旱对植被的胁迫,本研究主要分析稀疏植被、草本植被、落叶阔叶林和常绿针叶林植被对干旱的响应。

1.2 数据来源

(1)气象数据:2001—2022年分辨率为0.01°的中国逐月降水量和潜在蒸散发量栅格数据,下载自国家冰川冻土沙漠科学数据中心(https://s.ncdc.ac.cn/portal/metadata)。该数据是基于气象站点观测数据、再分析资料和卫星遥感数据制作而成[34]。

(2)遥感数据:2001—2022年500 m分辨率8 d中国GPP栅格数据,来源于GEE平台的PML_V2(CAS/IGSNRR/PML/V2_v018)数据产品。PML_V2通过冠层电导理论耦合GPP,其8 d平均GPP与实测值具有良好的一致性(R2=0.77,RMSE=1.99gC/m2·d)[35],在精度与稳定性方面表现优异,被广泛应用于陆地碳循环监测[36]。

(3)DEM数据:覆盖黄土高原的90 m分辨率SRTM DEM数据,下载自地理空间数据云(https://www.gscloud.cn/home)。

(4)植被类型数据:2000年、2005年、2010年、2015年、2020年和2022年共6期30 m分辨率中国土地覆盖分类栅格数据产品(GLC_FCS30)[37],下载自可持续发展大数据国际研究中心(https://data.casearth.cn)。依据黄土高原边界将研究区植被类型划分为稀疏植被、草本植被(草原、草丛、草甸)、落叶阔叶林、常绿针叶林、栽培植被和其他植被共6类。

(5)黄土高原矢量边界数据:黄土高原矢量边界数据下载自全球变化科学研究数据出版系统(https://www.geodoi.ac.cn)[38]。

利用黄土高原矢量边界,对以上数据进行裁切并重采样为0.01°,分别得到2001—2022年黄土高原逐月平均降水量、潜在蒸散发量和GPP数据,及2000年、2005年、2010年、2015年、2020年和2022年植被类型数据。

1.3 研究方法

1.3.1 标准化降水蒸散发指数(SPEI)

SPEI是通过降水量与潜在蒸散发量(potential evapotranspiration, PET)之间的平衡关系反映干旱情况,基于2001—2022年黄土高原逐月降水量和潜在蒸散发量,计算得到了2001—2022年黄土高原1~12个月SPEI数值,即SPEI1~SPEI12,具体计算公式参考文献[39]。根据国家标准化管理委员会发布的气象干旱等级划分标准[40],将SPEI划分为轻度干旱(-1<SPEI≤-0.5)、中度干旱(-1.5<SPEI≤-1)、重度干旱(-2<SPEI≤-1.5)和极端干旱(SPEI≤-2)。

1.3.2 植被GPPdt对干旱的响应时间确定

根据已有研究[41],GPP数据可能受长期趋势影响,如CO2浓度变化、土地利用变化或植被演替引起的植被生产力变化。为了更好地观察和分析短期或季节性变化,需对GPP数据进行去趋势化处理,记为GPPdt。并以2001—2022年黄土高原植被生长季GPPdt与不同时间尺度SPEI(SPEI-1~SPEI-12)的相关系数绝对值最大值所对应的SPEI时间尺度k,作为植被对干旱的响应时间[42-43]。计算公式为

$\begin{array}{l}{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{m}}=\frac{\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}({\mathit{x}}_{\mathit{i}}-\stackrel{-}{\mathit{x}})({\mathit{y}}_{\mathit{i}}^{\mathit{m}}-{\stackrel{-}{\mathit{y}}}^{\mathit{m}})}{\sqrt{\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}({\mathit{x}}_{\mathit{i}}{-\stackrel{-}{\mathit{x}})}^{2}}\sqrt{\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}({\mathit{y}}_{\mathit{i}}^{\mathit{m}}-{\stackrel{-}{\mathit{y}}}^{\mathit{m}}{)}^{2}}},\\ 1\le \mathit{m}\le 12\end{array}$
$\left|{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k}}\right|$=max($\left|{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-1}\right|$,$\left|{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-2}\right|$,…,$\left|{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-12}\right|$)

式中:${\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{m}}$为2001—2022年黄土高原植被生长季GPPdt月值与时间尺度为m的SPEI(SPEI-m)之间的相关系数;xi为2001—2022年黄土高原植被生长季第i个月GPPdt;$\stackrel{-}{\mathit{x}}$为2001—2022年黄土高原植被生长季GPPdt序列的均值;${\mathit{y}}_{\mathit{i}}^{\mathit{m}}$为2001—2022年黄土高原植被生长季第i个月SPEI-m值;${\stackrel{-}{\mathit{y}}}^{\mathit{m}}$为2001—2022年黄土高原植被生长季SPEI-m序列的均值;n为2001—2022年黄土高原植被生长季总月数,n=132;m为SPEI的时间尺度,即1个月、2个月、…、12个月;$\left|{\mathit{r}}_{\mathit{x}\mathit{y}}^{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k}}\right|$为2001—2022年黄土高原植被生长季GPPdt与不同时间尺度SPEI序列的最大相关系数的绝对值;SPEI-k为植被对干旱的响应时间k所对应的SPEI序列。

1.3.3 联合分布函数

Copula函数是由Sklar[44]于1959年提出的多元联合分布函数,它能够将任意两个或多个随机变量的边缘函数进行连接并组成联合分布。利用2001—2022年黄土高原植被生长季GPPdt和SPEI-k构建二维联合分布模型。

$\begin{array}{l}{\mathit{F}}_{\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{d}\mathrm{t}},\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k}}(\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{d}\mathrm{t}},\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k})=\mathit{C}\left[{\mathit{F}}_{\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{d}\mathrm{t}}}\right(\mathrm{G}\mathrm{P}{\mathrm{P}}_{\mathrm{d}\mathrm{t}}),\\ \mathit{ }\mathit{ }{\mathit{F}}_{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k}}(\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}-\mathit{k})]\end{array}$

式(3)中:FGPPdt,SPEI-k为GPPdt和SPEI-k的二维联合分布;FGPPdt和FSPEI-k分别为GPPdt和 SPEI-k序列的边缘分布;C(·)为Copula函数。

选用正态分布(Normal)、伽马分布(Gamma)、对数正态分布(Log-normal)、Logistic分布(Logistic)和威布尔分布(Weibull)分别拟合GPPdt序列。由于SPEI是采用标准正态分布变换计算得到的指数,故用正态分布(Normal)对其进行拟合。利用柯尔莫哥洛夫斯米尔诺夫检验(Kolmogorov-Smirnov, K-S)对分布函数拟合优度进行检验,选择通过K-S检验且贝叶斯信息准则(Bayesian information criteria, BIC准则)最小的分布作为最佳边缘分布。选用Gaussian、Student's t、Clayton、Frank、Gumbel、Joe分别构建GPPdt与SPEI-k的联合分布函数。采用极大似然估计法确定Copula函数的参数,并通过赤池信息准则(Akaike information criterion, AIC)选择最优Copula函数,AIC值越小表明联合分布越好。

1.3.4 条件概率函数

基于Copula函数的联合分布函数和条件概率函数,以2001—2022年黄土高原植被生长季GPPdt第40、第30、第20、第10百分位数作为植被GPPdt损失发生阈值,结合不同干旱等级,计算不同干旱等级和不同植被GPPdt损失发生阈值下的条件概率,计算公式[45]为

=$\begin{array}{l} P\left(\mathrm{GPP}_{\mathrm{dt}} \leqslant q_{p} \mid u_{1} \leqslant \mathrm{SPEI}-k \leqslant-u_{2}\right) \\ =\frac{P\left(\mathrm{GPP}_{\mathrm{dt}} \leqslant q_{p}, u_{1} \leqslant \mathrm{SPEI}-k \leqslant u_{2}\right)}{P\left(u_{1} \leqslant \mathrm{SPEI}-k \leqslant u_{2}\right)} \\ =\frac{F_{\mathrm{GPP}_{\mathrm{dt}}, \mathrm{SPEI}_{-k}}\left(q_{p}, u_{2}\right)-F_{\mathrm{GPP}_{\mathrm{dt}}, \mathrm{SPEI}-k}\left(q_{p}, u_{1}\right)}{F_{\mathrm{SPEI}-k}\left(u_{2}\right)-F_{\mathrm{SPEI}-k}\left(u_{1}\right)} \\ =\frac{C\left[F_{\mathrm{GPP}_{\mathrm{dt}}}\left(q_{p}\right), F_{\mathrm{SPEI}-k}\left(u_{2}\right)\right]-C\left[F_{\mathrm{GPP}_{\mathrm{dt}}}\left(q_{p}\right), F_{\mathrm{SPEI}-k}\left(u_{1}\right)\right]}{F_{\mathrm{SPEI}-k}\left(u_{2}\right)-F_{\mathrm{SPEI}-k}\left(u_{1}\right)} \end{array}$

式(4)中:u1和u2分别为1.3.1中不同干旱等级SPEI区间的上限和下限;qp为2001—2022年黄土高原植被生长季GPPdt的第p百分位数,p=40、30、20、10。

1.3.5 植被GPPdt损失的干旱触发阈值确定

选择2001—2022年黄土高原植被生长季GPPdt第40、30、20、10百分位数作为植被GPPdt损失发生阈值,以-3.00和-0.50作为导致植被GPPdt损失的干旱触发阈值区间的下限和上限[46],令SPEI从-0.50开始,以0.1为步长进行迭代,计算植被GPPdt损失发生时(GPPdt低于2001—2022年植被GPPdt由小到大排序的第10、20、30和40百分位数,即GPPdt≤GPPdt10th、GPPdt≤GPPdt20th、GPPdt≤GPPdt30th、GPPdt≤GPPdt40th),不同SPEI间隔内发生的条件概率,并将条件概率≥50%时,对应的SPEI间隔的下限作为植被GPPdt损失的干旱触发阈值。计算公式[36]为

$\begin{array}{l} P\left(\mathrm{GPP}_{\mathrm{dt}} \leqslant q_{p} \mid x_{i-0.1} \leqslant \operatorname{SPEI}-k \leqslant x_{i}\right) \\ =\frac{F_{\mathrm{GPP}_{\mathrm{dt}}, \mathrm{SPEI}_{-} k}\left(q_{p}, x_{i}\right)-F_{\mathrm{GPP}_{\mathrm{dt}}, \mathrm{SPEI}-k}\left(q_{p}, x_{i-0.1}\right)}{F_{\mathrm{SPEI}-k}\left(x_{i}\right)-F_{\mathrm{SPEI}-k}\left(x_{i-0.1}\right)} \\ =\frac{C\left[F_{\mathrm{GPP}_{\mathrm{dt}}}\left(q_{p}\right), F_{\mathrm{SPEI}-k}\left(x_{i}\right)\right]-C\left[F_{\mathrm{GPP}}\left(q_{\mathrm{dt}}\right), F_{\mathrm{SPEI}-k}\left(x_{i-0.1}\right)\right]}{F_{\mathrm{SPEI}-k}\left(x_{i}\right)-F_{\mathrm{SPEI}-k}\left(x_{i-0.1}\right)} \end{array}$

式(5)中:xi-0.1和xi分别为SPEI间隔的下限和上限。

2 结果分析

2.1 植被GPPdt对干旱的响应时间

2001—2022年黄土高原植被生长季GPPdt与SPEI的最大相关系数以正相关关系为主,占像元总数95.05%(不含水体、裸地等无数据像元)。GPPdt与SPEI最大相关系数集中分布在0.1~0.5,共占像元总数的85.98%,其中0.1~0.2占像元总数的29.04%,0.2~0.3占像元总数的29.67%,0.3~0.4占像元总数的17.13%,0.4~0.5占像元总数的10.15%,如图2(a)所示。空间分布上,2001—2022年黄土高原植被生长季GPPdt与SPEI的最大相关系数高值区集中在甘肃北部、东部,内蒙古西部边缘地带和宁夏中部;低值区(负相关的区域)主要分布在黄土高原西部和南部边缘,如陕西南部、青海西部和甘肃南部。2001—2022年黄土高原植被对干旱的响应时间集中在2~6个月,共占像元总数的61.85%,其中响应时间为2个月、3个月、4个月、5个月和6个月的像元分别占像元总数的9.14%、14.60%、13.45%、8.42%和16.23%,如图2(b)所示。空间分上,植被对干旱的响应时间≤3个月的区域集中在黄土高原北部的内蒙古北部、宁夏北部,反映出该区域植被对短期水分变化响应迅速;植被对干旱的响应时间≥10个月的区域主要分布在黄土高原西部的甘肃北部、宁夏西部,反映出这些地区植被对干旱的相应相对滞后,可能与不同植被类型的生理适应性差异有关。

选择2000—2022年共6期植被类型数据中,植被类型(稀疏植被、草本植被(草原、草丛、草甸)、落叶阔叶林和常绿针叶林)未发生改变的像元,分析2001—2022年不同植被类型GPPdt与SPEI的最大相关系数绝对值及对干旱的响应时间,如表1所示。结果表明,2001—2022年黄土高原稀疏植被GPPdt与SPEI最大相关系数绝对值的像元平均值最高,为0.343;其次为草本植被,为0.248;落叶阔叶林和常绿针叶林最低,分别为0.110和0.104。植被GPPdt对干旱的响应时间由短到长依次为稀疏植被(3.89个月)、草本植被(5.78个月)、落叶阔叶林(5.80个月)和常绿针叶林(6.38个月)。由于稀疏植被和草本植被多以一年生草本植物为主,根系较浅,易受短期干旱影响,落叶阔叶林和常绿针叶林等乔木林根系较深,可通过吸收深层土壤水维持生长,对干旱具有更强的适应能力从而响应迟缓。

2.2 不同干旱等级下植被GPPdt损失发生的条件概率

随着干旱等级的增加,相同植被GPPdt损失发生阈值下,植被GPPdt损失发生的条件概率逐渐增大。具体表现为:GPPdt≤GPPdt10th(以植被GPPdt第10百分位数作为植被GPPdt损失发生阈值)情景下,不同干旱等级(轻度、中度、重度和极度干旱)下植被GPPdt损失发生的条件概率平均值分别为12.05%、14.28%、16.77%和21.35%,如图3(a)~图3(d)所示。GPPdt≤GPPdt20th(以植被GPPdt第20百分位数作为植被GPPdt损失发生阈值时)情景下,不同干旱等级下植被GPPdt损失发生的条件概率平均值分别为22.07%、25.05%、28.20%和33.41%,如图3(e)~图3(h)所示。GPPdt≤GPPdt30th(以植被GPPdt第30百分位数作为植被GPPdt损失发生阈值)情景下,不同干旱等级下植被GPPdt损失发生的条件概率平均值分别为31.55%、35.00%、38.45%和43.78%,如图3(i)~图3(l)所示。GPPdt≤GPPdt40th(以植被GPPdt第40百分位数作为植被GPPdt损失发生阈值)情景下,不同干旱等级下植被GPPdt损失发生的条件概率平均值分别为42.54%、46.19%、49.70%和54.81%,如图3(m)~图3(p)所示。

而随着植被GPPdt损失发生阈值的增加(从第10百分位数增加至第20、30、40百分位数),相同干旱等级下,植被GPPdt损失发生的条件概率逐渐增大。轻度干旱等级下,GPPdt≤GPPdt10th、GPPdt≤GPPdt20th、GPPdt≤GPPdt30th、GPPdt≤GPPdt40th情景下植被GPPdt损失发生的条件概率平均值分别为12.05%、22.07%、31.55%和42.54%,如图3(a)、图3(e)、图3(i)和图3(m)所示。中度干旱条件下,GPPdt≤GPPdt10th、GPPdt≤GPPdt20th、GPPdt≤GPPdt30th、GPPdt≤GPPdt40th情景下植被GPPdt损失发生的条件概率平均值分别为14.28%、25.05%、35.00%和46.19%,如图3(b)、图3(f)、图3(j)和图3(n)所示。重度干旱条件下,GPPdt≤GPPdt10th、GPPdt≤GPPdt20th、GPPdt≤GPPdt30th、GPPdt≤GPPdt40th情景下,植被GPPdt损失发生的条件概率平均值分别为16.77%、28.20%、38.45%和49.70%,如图3(c)、图3(g)、图3(k)和图3(o)所示。极端干旱条件下,GPPdt≤GPPdt10th、GPPdt≤GPPdt20th、GPPdt≤GPPdt30th、GPPdt≤GPPdt40th情景下,植被GPPdt损失发生的条件概率平均值分别为21.35%、33.41%、43.78%和54.81%,如图3(d)、图3(h)、图3(l)和图3(p)所示。

总体上,相较于干旱等级,植被GPPdt损失发生阈值的增加对植被GPPdt损失发生的条件概率影响更大。

2.3 从气象干旱到植被GPPdt损失的干旱触发阈值

当干旱触发阈值较小时(较小的SPEI),表明只有发生更严重的干旱才会导致植被GPPdt损失,反映出植被具有较强的抗旱能力。相反,当干旱触发阈值较高时(较大的SPEI)表明即使较轻的干旱也可能导致植被GPPdt损失,即植被对干旱事件较为敏感。参考气象干旱等级划分标准,本研究将干旱触发阈值划分为高干旱触发阈值(-1<SPEI≤-0.6)、中干旱触发阈值(-2<SPEI≤-1)和低干旱触发阈值(-3≤SPEI≤-2)。结果表明,随着植被GPPdt损失阈值的增加(从第10百分位数增加至第20、30、40百分位数),干旱触发阈值的像元数量逐渐变多,如图4所示,且随着植被GPPdt损失发生阈值的增加干旱触发阈值四分位距(interquartile range, IQR)向高SPEI值方向移动,如图5所示。

GPPdt≤GPPdt10th(以植被GPPdt第10百分位数作为植被GPPdt损失发生阈值)情景下,干旱触发阈值平均值为-2.67,IQR介于-3.0~-2.5,以低干旱触发阈值为主,占像元总数的99.67%。空间上主要分布在甘肃北部、宁夏中部、内蒙古西部、陕西南部和河南东部,如图4(a)和图5所示。GPPdt≤GPPdt20th(以植被GPPdt第20百分位数作为植被GPPdt损失发生阈值)情景下,干旱触发阈值平均值为-2.51,IQR介于-2.8~-2.3,中、低干旱触发阈值像元分别占像元总数的5.51%和94.49%。中干旱触发阈值像元主要分布在甘肃北部和东部、内蒙古西部、陕西南部和河南东部。低干旱触发阈值像元主要分布中干旱触发阈值周围,如图4(b)和图5所示。GPPdt≤GPPdt30th(以植被GPPdt第30百分位数作为植被GPPdt损失发生阈值)情景下,干旱触发阈值平均值为-2.22,IQR介于-2.6~-1.9。高、中和低干旱触发阈值像元分别占像元总数的0.04%、31.83%和68.13%,中干旱触发阈值主要分布在甘肃北部和东部、宁夏中部、内蒙古西部和陕西北部。低干旱触发阈值主要分布在陕西北部和内蒙古东部,如图4(c)和图5所示。GPPdt≤GPPdt40th(以植被GPPdt第40百分位数作为植被GPPdt损失发生阈值)情景下,干旱触发阈值平均值为-2.51,IQR介于-2.0~-1.2,高、中和低干旱触发阈值像元分别占像元总数的11.62%、65.10%和23.28%。高干旱触发阈值主要分布在甘肃北部和东部、陕西北部和内蒙古西部,低干旱触发阈值主要分布在山西东部和内蒙古北部,如图4(d)和图5所示。

总体上,2001—2022年黄土高原高干旱触发阈值主要分布在区域西北部和中部,即甘肃北部和东部、宁夏中部、内蒙古西部和陕西北部,这些区域较其他地区对干旱响应更为敏感。

2.4 不同植被类型GPPdt损失概率

以2001—2022年黄土高原植被生长季GPPdt第40百分位所对应数值GPPdt40th作为植被GPPdt损失发生阈值(GPPdt≤GPPdt40th),统计不同植被类型(稀疏植被、草本植被(草原、草丛、草甸)、落叶阔叶林和常绿针叶林)在不同干旱等级下GPPdt损失概率及导致GPPdt损失的干旱触发阈值。

通过表2看出,2001—2022年黄土高原稀疏植被在不同干旱条件下GPPdt损失概率平均值为52.21%,其中在轻度、中度、重度和极端干旱下植被GPPdt损失概率分别为42.74%、48.60%、54.39%、63.10%。随着干旱等级增加,稀疏植被GPPdt损失概率呈快速上升趋势,极端干旱较轻度干旱情景下,植被GPPdt损失概率增加20.36%。草本植被在不同干旱条件下GPPdt损失概率平均值为48.49%。其中在轻度、中度、重度和极端干旱下植被GPPdt损失概率分别为42.51%、46.31%、49.93%、55.22%。落叶阔叶林在不同干旱条件下植被GPPdt损失概率平均值为47.77%。其中在轻度、中度、重度和极端干旱下植被GPPdt损失概率分别为46.64%、47.47%、48.11%、48.84%。随着干旱等级的增加,落叶阔叶林GPPdt损失概率变化幅度较小(相较于轻度干旱,极端干旱植被GPPdt损失概率仅增加2.20%),呈相对稳定的特征。常绿针叶林在不同干旱条件下GPPdt损失概率平均值为42.78%,其中在轻度、中度、重度和极端干旱下植被GPPdt损失概率分别为41.60%、42.44%、43.11%、43.98%,在4种植被类型中GPPdt损失概率最小,且极端干旱较轻度干旱的植被GPPdt损失概率仅增加2.38%。

总体上,2001—2022年黄土高原植被GPPdt损失的干旱触发阈值由小到大排序依次为稀疏植被(SPEI=-1.486)、草本植被(SPEI=-1.573)、落叶阔叶林(SPEI=-1.601)和常绿针叶林(SPEI=-1.676),反映出稀疏植被对干旱事件最为敏感,其次为草本植被和落叶阔叶林,常绿针叶林干旱触发阈值较低,表明对干旱具有较强的抵抗力。

3 结论

(1)2001—2022年黄土高原植被生长季GPPdt与SPEI以正相关关系为主,占像元总数95.05%。植被对干旱的响应时间主要集中在2~6个月,共占像元总数的61.85%。不同植被类型对干旱的响应时间由短到长依次为稀疏植被(3.89个月)、草本植被(5.69个月)、落叶阔叶林(5.83个月)和常绿针叶林(6.38个月)。

(2)2001—2022年黄土高原相同植被GPPdt损失发生阈值下,植被GPPdt损失发生的条件概率随着干旱等级的增加而增加。相同干旱等级下,植被GPPdt损失发生的条件概率也随植被GPPdt损失发生阈值的增加而增加。

(3)2001—2022年黄土高原高干旱触发阈值主要分布在区域西北部和中部,即甘肃北部和东部、宁夏中部、内蒙古西部和陕西北部,这些区域较其他地区对干旱响应更为敏感。

(4)以2001—2022年黄土高原植被生长季GPPdt第40百分位所对应数值GPPdt40th作为植被GPPdt损失发生阈值(GPPdt≤GPPdt40th),不同植被类型GPPdt损失概率由高到低分别为稀疏植被(52.21%)、草本植被(48.49%)、落叶阔叶林(47.77%)和常绿针叶林(42.78%)。导致植被GPPdt损失的干旱触发阈值由小到大排序依次为稀疏植被(SPEI=-1.486)、草本植被(SPEI=-1.573)、落叶阔叶林(SPEI=-1.601)和常绿针叶林(SPEI=-1.676)。

参考文献

[1]

Wang S P, Zhang Q, Liu Y Z, et al. Effects of global warming on drought onset in China[J]. Journal of Hydrology, 2024, 632: 130964.

[2]

曲芷程, 黄绍普, 刘司博, 等. 1963—2020年蒙古锡林河流域极端气候及其对水文干旱的影响[J]. 生态学报, 2025, 45(11): 5386-5397.

[3]

Qu Zhicheng, Huang Shaopu, Liu Sibo, et al. Extreme climate and its impact on hydrological drought in the Xilin River Basin of Inner Mongolia over the past 60 years[J]. Acta Ecologica Sinica, 2025, 45(11): 5386-5397.

[4]

Zeng Z, Wu W, Li Y, et al. Increasing meteorological drought under climate change reduces terrestrial ecosystem productivity and carbon storage[J]. One Earth, 2023, 6(10): 1326-1339.

[5]

张璐, 畅建霞, 杨贵斌, 等. 干旱胁迫下植被生产力与冠层结构易损性的非一致性评估[J]. 农业工程学报, 2025, 41(6): 145-153.

[6]

Zhang Lu, Chang Jianxia, Yang Guibin, et al. Assessing the non-uniform vulnerability of vegetation productivity and canopy structure under drought stress[J]. Transactions of the Chinese Society of Agricultural Engineering, 2025, 41(6): 145-153.

[7]

范虎, 党星海, 赵健赟, 等. 基于温度植被干旱指数的西北干旱区干旱时空变化[J]. 科学技术与工程, 2024, 24(27): 11537-11546.

[8]

Fan Hu, Dang Xinghai, Zhao Jianyun, et al. Temporal and spatial variation of drought in northwest arid region based on temperature vegetation drought index[J]. Science Technology and Engineering, 2024, 24(27): 11537-11546.

[9]

Buttlar J V, Zscheischler J, Rammig A, et al. Impacts of droughts and extreme-temperature events on gross primary production and ecosystem respiration: a systematic assessment across ecosystems and climate zones[J]. Biogeosciences, 2018, 15(5): 1-39.

[10]

杨肖丽, 崔周宇, 任立良, 等. 1966—2015年长江流域水文干旱时空演变归因[J]. 水科学进展, 2023, 34(3): 349-359.

[11]

Yang Xiaoli, Cui Zhouyu, Ren Liliang, et al. Patterns and attributions of hydrological drought in the Yangtze River basin from 1966 to 2015[J]. Advances in Water Science, 2023, 34(3): 349-359.

[12]

Yang G B, Chang J X, Wang Y M, et al. Understanding drought propagation through coupling spatiotemporal features using vine copulas: a compound drought perspective[J]. Science of the Total Environment, 2024, 921: 171080.

[13]

Li L, Zheng Z, Biederman J A, et al. Drought and heat wave impacts on grassland carbon cycling across hierarchical levels[J]. Plant, Cell & Environment, 2021, 44(7): 2402-2413.

[14]

Li Y, Zhuang Q, Zhao H, et al. Evaluation of the resistance and resilience of terrestrial ecosystems to drought in southwest China[J]. Journal of Hydrology, 2025, 646: 132318.

[15]

Fang W, Huang S, Huang Q, et al. Probabilistic assessment of remote sensing-based terrestrial vegetation vulnerability to drought stress of the Loess Plateau in China[J]. Remote Sensing of Environment, 2019, 232: 111290.

[16]

Ding Y, Xu J, Wang X, et al. Spatial and temporal effects of drought on Chinese vegetation under different coverage levels[J]. Science of the Total Environment, 2020, 716: 137166.

[17]

Yao T, Wu C, Yeh P J F, et al. Assessing the response lag and vulnerability of terrestrial vegetation to various compound climate events in Mainland China under different vegetation types[J]. Stochastic Environmental Research and Risk Assessment, 2025, 39(5): 2231-2250.

[18]

Gouveia C M, Trigo R M, Beguería S, et al. Drought impacts on vegetation activity in the Mediterranean region: an assessment using remote sensing data and multi-scale drought indicators[J]. Global and Planetary Change, 2017, 151: 15-27.

[19]

史峥, 崔婷茹, 孙小诺, 等. 基于1990—2015年归一化植被指数的河北植被时间变化特征及其对干旱的响应[J]. 科学技术与工程, 2020, 20(25): 10347-10355.

[20]

Shi Zheng, Cui Tingru, Sun Xiaonuo, et al. Temporal variations of vegetation and its response to drought based on normalized difference vegetation index data in Hebei Province in 1990—2015[J]. Science Technology and Engineering, 2020, 20(25): 10347-10355.

[21]

王小红, 刘宪锋, 孙高鹏, 等. 2001—2020年秦巴山区植被生产力对干旱的响应[J]. 应用生态学报, 2022, 33(8): 2105-2112.

[22]

Wang Xiaohong, Liu Xianfeng, Sun Gaopeng, et al. Response of vegetation productivity to drought in the Qinling-Daba Mountains, China from 2001 to 2020[J]. Chinese Journal of Applied Ecology, 2022, 33(8): 2105-2112.

[23]

李阳阳, 董国涛, 薛华柱, 等. 1982—2020年黄河流域多时间尺度气象干旱对植被影响[J]. 水土保持学报, 2024, 38(1): 187-196.

[24]

Li Yangyang, Dong Guotao, Xue Huazhu, et al. Effects of multi-timescale meteorological drought on vegetation in the Yellow River Basin from 1982 to 2020[J]. Journal of Soil and Water Conservation, 2024, 38(1): 187-196.

[25]

Won J, Seo J, Lee J, et al. Vegetation drought vulnerability mapping using a copula model of vegetation index and meteorological drought index[J]. Remote sensing, 2021, 13(24): 5103.

[26]

Yuan Y, Bao A, Jiang P, et al. Probabilistic assessment of vegetation vulnerability to drought stress in Central Asia[J]. Journal of Environmental Management, 2022, 310: 114504.

[27]

周辉, 张新松, 郭晓丽, 等. 基于Copula理论的风电出力相关性建模及在输电网规划中的应用[J]. 科学技术与工程, 2018, 18(4): 273-278.

[28]

Zhou Hui, Zhang Xinsong, Guo Xiaoli, et al. Modeling of wind power correlation based on Copula theory and its application in transmission network expansion planning[J]. Science Technology and Engineering, 2018, 18(4): 273-278.

[29]

任青阳, 施俭, 王彦丁, 等. 不同拟合指标下模拟节理岩体结构面产状的Copula函数方法[J]. 科学技术与工程, 2025, 25(10): 4274-4283.

[30]

Ren Qingyang, Shi Jian, Wang Yanding, et al. Copula function method for simulating the occurrence of jointed rock mass structural plane under different fitting indices[J]. Science Technology and Engineering, 2025, 25(10): 4274-4283.

[31]

杨蕊, 吴时强, 高学平, 等. 基于Vine Copula函数的河湖连通水环境多因子联合风险识别研究[J]. 水利学报, 2020, 51(5): 606-616.

[32]

Yang Rui, Wu Shiqiang, Gao Xueping, et al. A vine copula-based study on identification of multivariate water environmental risk under different connectivity of rivers and lakes[J]. Journal of Hydraulic Engineering, 2020, 51(5): 606-616.

[33]

包青格乐, 吴志祥, 刘文杰, 等. 复合干热事件对橡胶林总初级生产力的影响研究[J]. 热带作物学报, 2025, 46(9): 2190-2203.

[34]

Bao Qinggele, Wu Zhixiang, Liu Wenjie, et al. The impact of compound drought-hot events on the gross primary productivity of rubber plantations[J]. Chinese Journal of Tropical Crops, 2025, 46(9): 2190-2203.

[35]

Poonia V, Jha S, Goyal M K. Copula based analysis of meteorological, hydrological and agricultural drought characteristics across Indian river basins[J]. International Journal of Climatology, 2021, 41(9): 4637-4652.

[36]

Jha S, Das J, Sharma A, et al. Probabilistic evaluation of vegetation drought likelihood and itsimplications to resilience across India[J]. Global and Planetary Change, 2019, 176: 23-35.

[37]

Gao S, Huang S, Singh V P, et al. Dynamic response of vegetation to meteorological drought and driving mechanisms in Mongo-lian Plateau[J]. Journal of Hydrology, 2025, 650: 132541.

[38]

Li Y, Yang J, Wu J, et al. Spatiotemporal responses and vulnerability of vegetation to drought in the Ili river transboundary basin: a comprehensive analysis based on copula theory, SPEI, and NDVI[J]. Remote Sensing, 2025, 17(5): 801.

[39]

Gao X, Zhao Q, Zhao X, et al. Temporal and spatial evolution of the standardized precipitation evapotranspiration index (SPEI) in the Loess Plateau under climate change from 2001 to 2050[J]. Science of the Total Environment, 2017, 595: 191-200.

[40]

刘哲琼, 佘敦先, 夏军, 等. 黄土高原地区气象干旱动态格局演变及其对植被的影响[J]. 地理科学, 2023, 43(9): 1659-1671.

[41]

Liu Zheqiong, She Dunxian, Xia Jun, et al. Impacts of meteorological drought on vegetation and their dynamic pattern evolutions on the Loess Plateau[J]. Scientia Geographica Sinica, 2023, 43(9): 1659-1671.

[42]

苏文壮, 杨永崇, 王涛, 等. 近20年黄土高原干旱变化特征及其影响因素[J]. 科学技术与工程, 2023, 23(11): 4551-4560.

[43]

Su Wenzhuang, Yang Yongchong, Wang Tao, et al. Variation characteristics and influencing factors of drought in the Loess Plateau in recent 20 years[J]. Science Technologyand Engineering, 2023, 23(11): 4551-4560.

[44]

She D, Xia J. Copulas-based drought characteristics analysis and risk assessment across the Loess Plateau of China[J]. Water Resources Management, 2018, 32(2): 547-564.

[45]

李家誉, 佘敦先, 张利平, 等. 黄土高原植被变化对气象干旱多尺度响应特征与机制[J]. 水土保持学报, 2022, 36(6): 280-289.

[46]

Li Jiayu, She Dunxian, Zhang Liping, et al. Multi-scale response characteristics and mechanism of vegetation to meteorological drought on the Loess Plateau[J]. Journal of Soil and Water Conservation, 2022, 36(6): 280-289.

[47]

张含玉, 方怒放, 史志华. 黄土高原植被覆盖时空变化及其对气候因子的响应[J]. 生态学报, 2016, 36(13): 3960-3968.

[48]

Zhang Hanyu, Fang Nufang, Shi Zhihua. Spatio-temporal patterns for the NDVI and its responses to climatic factors in the Loess Plateau, China[J]. Acta Ecologica Sinica, 2016, 36(13): 3960-3968.

[49]

Zhang L, Li X, Zheng D H, et al. Merging multiple satellite-based precipitation products and gauge observations using a novel double machine learning approach[J]. Journal of Hydrology, 2021, 594:125969.

[50]

Zhang Y, Kong D, Gan R, et al. Coupledestimation of 500 m and 8-day resolution global evapotranspiration and gross primary production in 2002—2017[J]. Remote sensing of environment, 2019, 222: 165-182.

[51]

Wang Y, Gao G, Huang Y, et al. Ecosystem water use efficiency and carbon use efficiency respond oppositely to vegetation greening in China's Loess Plateau[J]. Science of the Total Environment, 2025, 964: 178575.

[52]

Zhang X, Zhao T, Xu H, et al. GLC_FCS30D: The first global 30 m land-cover dynamic monitoring product with a fine classification systemfrom 1985 to 2022 using dense time-series Landsat imagery and continuous change-detection method[J]. Earth System Science Data Discussions, 2023, 1-32.

[53]

王正兴. 黄土高原地区[J]. 全球变化数据学报(中英文), 2017, 1(1): 113,236.

[54]

Wang Zhengxing. Boundary data of the Loess Plateau region[J]. Journal of Global Change Data & Discovery, 2017, 1(1): 113,236.

[55]

Vicente-Serrano S M, Beguería S, López-Moreno J I. A multiscalar drought index sensitive to global warming: the standardized precipitation evapotranspiration index[J]. Journal of Climate, 2010, 23(7): 1696-1718.

[56]

中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会. 气象干旱等级: GB/T 20481—2017[S]. 北京: 中国标准出版社, 2017.

[57]

General Administration of Quality Supervision,Inspection and Qua-rantine of the People's Republic of China, National Standardization Administration. Classification of meteorological drought: GB/T 20481—2017[S]. Beijing: Standards Press of China, 2017.

[58]

张更喜, 王慧敏, 粟晓玲, 等. 复合干热胁迫下黄土高原夏季植被脆弱性评估[J]. 农业工程学报, 2024, 40(6): 339-346.

[59]

Zhang Gengxi, Wang Huimin, Su Xiaoling, et al. Assessing the vegetation vulnerability of Loess Plateau under compound dry and hot climates[J]. Transactions of the Chinese Society of Agricultu-ral Engineering, 2024, 40(6): 339-346.

[60]

Zhao A, Yu Q, Feng L, et al. Evaluating the cumulative and time-lag effects of drought on grassland vegetation: a case study in the Chinese Loess Plateau[J]. Journal of Environmental Management, 2020, 261: 110214.

[61]

Wei X, He W, Zhou Y, et al. Global assessment of lagged and cumulative effects of drought on grassland gross primary production[J]. Ecological Indicators, 2022, 136: 108646.

[62]

Sklar A. Fonctions de répartition à n dimensions et leurs marges[M]. Paris: Publication Institute Statistics University, 1959.

[63]

Guo W, Huang S, Huang Q, et al. Drought trigger thresholds for different levels of vegetation loss in China and their dynamics[J]. Agricultural and Forest Meteorology, 2023, 331: 109349.

[64]

Yao N, Li Y, Dong Q, et al. Influence of the accuracy of refe-rence crop evapotranspiration on drought monitoring usingstandar-dized precipitation evapotranspiration index in Mainland China[J]. Land Degradation & Development, 2020, 31(2): 266-282.

基金资助

国家自然科学基金(42271309)

陕西省哲学社会科学研究专项(2025YB0239)

陕西省自然科学基础研究项目(2024JC-YBQN-0306)

中国地质调查局地质调查项目(DD20220882)

PDF (19034KB)

48

访问

0

被引

详细

导航
相关文章

/

〈 〉