生态学报  2022, Vol. 42 Issue (20): 8331-8342

文章信息

罗爽, 刘会玉, 龚海波
LUO Shuang, LIU Huiyu, GONG Haibo
1982-2018年中国植被覆盖变化非线性趋势及其格局分析
Nonlinear trends and spatial pattern analysis of vegetation cover change in China from 1982 to 2018
生态学报. 2022, 42(20): 8331-8342
Acta Ecologica Sinica. 2022, 42(20): 8331-8342
http://dx.doi.org/10.5846/stxb202108252372

文章历史

收稿日期: 2021-08-25
网络出版日期: 2022-06-15
1982-2018年中国植被覆盖变化非线性趋势及其格局分析
罗爽1,2,3,4,5,6 , 刘会玉1,3,4,5,6 , 龚海波1,3,4,5,6     
1. 南京师范大学地理科学学院, 南京 210023;
2. 南京大学地理与海洋科学学院, 南京 210023;
3. 虚拟地理环境教育部重点实验室(南京师范大学), 南京 210023;
4. 江苏省地理环境演化国家重点实验室培育建设点, 南京 210023;
5. 江苏省地理信息资源开发与利用协同创新中心, 南京 210023;
6. 江苏省环境演变与生态建设重点实验室, 南京 210023
摘要: 探究植被覆盖变化是评估陆地生态系统环境变化的重要手段, 但现有研究多采用线性趋势来表达植被覆盖的变化情况而忽略了趋势的非线性。本文使用GLASS FVC数据, 利用BFAST方法和格局分析, 探讨了1982-2018年我国植被覆盖变化的非线性趋势及其分布格局。结果表明: (1)与线性趋势方法的对比发现, BFAST的检测结果揭示了四川盆地、黄土高原等地的植被覆盖显著增加趋势其实存在中断, 青海和东北等地植被覆盖经历了由退化到改善的过程而并非简单的线性增加, 而青藏高原中东部等地则由原先的改善趋势变为了退化趋势。(2)将非线性趋势结果进行分类, 其中单调型增加类型占比最多, 达到33.58%, 主要分布在内蒙古、陕西及河南等地; 单调型减少占比1.82%, 主要分布在东南沿海地区; 中断型增加占比22.91%, 主要分布在四川盆地东部和华北地区; 中断型减少占比2.68%, 主要分布在青藏高原东南部; 由增到减占比4.20%, 主要分布在青海等地; 由减到增占比14.62%, 主要分布在吉林等地。大范围的植被覆盖增加趋势充分反映了我国过去几十年植被的改善, 但同时存在的减少趋势表明潜在的植被退化风险仍不可忽视。(3)不同趋势类型发生改变的时间有所差异, 总体上1988-1999年间发生的改变较少, 而2000-2011年间发生的改变较多, 我国21世纪以来实施的大规模生态保护和恢复工程对植被的改善过程有重要影响。(4)分布格局上, 植被覆盖改善趋势类型(单调型增加, 中断型增加, 由减到增)呈现大聚集, 小分散的特点, 具有复杂的形状; 退化趋势类型(单调型减少, 中断型减少, 由增到减)的面积均较小, 分布也相对离散。全国尺度上趋势空间格局呈现一定规律但分布的异质性较大, 区域尺度上植被覆盖经受的干扰显著, 变化过程实际也是较为复杂的。本研究表明, 使用非线性趋势方法和格局分析, 可以更准确地评估植被覆盖的时空变化, 从而为生态环境相关工作的开展提供科学的参考。
关键词: BFAST    植被覆盖变化    非线性趋势    空间格局    
Nonlinear trends and spatial pattern analysis of vegetation cover change in China from 1982 to 2018
LUO Shuang1,2,3,4,5,6 , LIU Huiyu1,3,4,5,6 , GONG Haibo1,3,4,5,6     
1. College of Geography Science, Nanjing Normal University, Nanjing 210023, China;
2. School of Geography and Ocean Science, Nanjing University, Nanjing 210023, China;
3. Key Laboratory of Virtual Geographic Environment (Nanjing Normal University), Ministry of Education, Nanjing 210023, China;
4. State Key Laboratory Cultivation Base of Geographical Environment Evolution (Jiangsu Province), Nanjing Normal University, Nanjing 210023, China;
5. Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing Normal University, Nanjing 210023, China;
6. Jiangsu Key Laboratory of Environmental Change and Ecological Construction, Nanjing Normal University, Nanjing 210023, China
Abstract: The study of vegetation cover change is an important tool to evaluate the change of terrestrial ecosystem, but the previous studies mostly used linear trend to express the change of vegetation cover, ignoring the change of trend over time. Using GLASS FVC data, BFAST method and landscape pattern analysis, the nonlinear trend and distribution pattern of vegetation cover change in China were explored from 1982 to 2018. The results showed that: (1) compared with the results of linear trend method, although the vegetation cover in Sichuan Basin and the Loess Plateau presented a significantly linear increase trend, the detection results of BFAST showed that the increase trend was interrupted. Meanwhile, in Qinghai and Northeast China, the vegetation cover experienced a shift from degradation to improvement rather than simple linear increase, while it changed from the first improvement trend to later degradation in the middle east of Qinghai-Tibet Plateau. (2) Among the nonlinear trend types, the monotonic increase trend accounted for 33.58%, mainly distributed in the Inner Mongolia, Shaanxi, Henan and other places; the monotonic decrease trend accounted for 1.82%, mainly distributed in the southeast coastal area; the interrupted increase trend accounted for 22.91%, mainly distributed in the east of Sichuan Basin and the North China; the interrupted decrease trend accounted for 2.68%, mainly distributed in the southeast of Qinghai Tibet Plateau; the trend shifted from increase to decrease type accounted for 4.20%, mainly distributed in Qinghai; the trend shifted from the decrease to increase type accounted for 14.62%, mainly distributed in Jilin. The large-scale increasing trend of vegetation cover fully reflected the improvement of vegetation in China during the past few decades, but the monotonic decrease and the increase to decrease trends showed that the potential risk of vegetation degradation could not be ignored. (3) The time of change of different trend types were different, less occurred from 1988 to 1999, but mainly occurred from 2000 to 2011, which has been probably affected by the large-scale ecological protection and restoration projects since the 21st century. (4) The patches of the vegetation improvement trend type (monotonic increase, interrupted increase, decrease to increase) were characterized by large aggregation, small dispersion, and complex shape driven by strong interference, indicating that the process of vegetation cover change was actually very complex with high diversity on the regional scale; the patches of degradation trend type were small and dispersedly distributed. This study shows that the nonlinear trend and pattern analysis can more accurately evaluate vegetation cover change, so as to provide a scientific reference for the protection and recovery of ecological environment.
Key Words: BFAST    vegetation cover change    nonlinear trend    spatial pattern    

地表植被是陆地生态系统的重要组成部分, 它不仅是连接大气、土壤和水体等自然元素的“纽带”, 更是全球气候变化和生态系统变化的重要指示器[ 1 2]。近年来, 随着全球气候变暖和人类活动的加剧, 监测、理解和预测植被覆盖对全球变化的响应成为地球系统科学的核心活动[ 3]。了解植被覆盖的发展变化不仅有助于全面了解区域生态环境的演变历史, 也有利于评估区域生态环境系统质量从而制定可持续发展的生态环境保护策略, 因此对于植被覆盖长期动态变化的研究也已经成为生态学领域重点关注的热点之一[ 4]。植被覆盖度(Fractional Vegetation Cover, FVC)是指单位面积内植被地上部分(包括叶、茎、枝)在地面的垂直投影面积占统计区总面积的百分比[ 5], 植被覆盖的变化可以运用植被覆盖度的变化非常直观地指示, 而植被覆盖度的变化则能代表陆地生态系统一定程度上的波动或变化[ 6 7]。与此同时, 各种植被覆盖状态的空间排布与组合特征构成其分布格局[ 8], 研究植被覆盖分布格局能够充分反映植被空间分布及其在环境异质性和干扰状况综合控制下的动态变化特征[ 9]

目前, 针对我国植被覆盖动态变化的研究多数采用较为简单的线性趋势分析方法来提取其趋势特征, 包括一元线性回归分析方法、Theil-Sen斜率与Mann-Kendall检验结合的趋势分析法等[ 10], 如Mu等[ 11], 马梓策等[ 12], 金凯等[ 13]使用上述方法研究了中国植被覆盖长期变化趋势的时空特征, 发现了过去几十年来我国大部分地区的植被覆盖呈现线性增长趋势。但是普通线性分析方法认为植被覆盖变化趋势在整个时期是固定不变的, 有研究指出长时间总趋势可以由多个不同甚至相反的阶段性趋势组成[ 14 15]。针对线性趋势分析方法的不足, 我们更需要采取非线性趋势分析方法来检测植被覆盖的变化[ 16 17]。趋势分解法是一种有效的非线性趋势分析方法, 其中以季节性和趋势加断点检测法(Breaks For Additive Seasonal and Trend, BFAST)应用较普遍[ 10]。BFAST方法通常被用来检测长时间序列过程中趋势的突变, 从而可以将总趋势分为若干片段[ 18 20], 但目前也有研究使用该方法的变体BFAST01(即定义至多一个突变点), 来对植被NDVI长时间序列变化进行趋势分类[ 21 22]。为了强化植被覆盖长期变化趋势的提取从而全面认识植被覆盖动态变化特征, 对其进行非线性趋势研究并进行归类分析是很有必要的。

我国幅员辽阔, 自然条件复杂多样, 而随着近年来全球气候变化和人类活动的加剧, 尤其是在全球变暖、城市化迅速发展及20世纪末以来开展了大规模的生态恢复工程的大背景下, 我国的植被覆盖经历了复杂的变化, 但是针对植被覆盖长期动态变化趋势的研究仍存在一些不足。首先, 关于全国尺度植被覆盖变化的宏观格局认识仍较为缺乏, 尤其是长时间尺度的时空变化特征更是鲜有报道[ 23]。其次, 现有研究较少关注到植被覆盖长期变化中阶段性细节过程, 因此往往会忽略其中可能存在的突变和趋势转换, 最终会低估或者高估趋势变化带来的植被退化的风险。最后, 现有研究多针对分级的植被覆盖度进行分布格局动态变化分析[ 24], 少有研究探究不同植被覆盖变化趋势的空间格局。基于以上认识, 我们更加关注过去几十年间我国植被覆盖的非线性变化过程, 以识别其在何时何地发生了何种具体变化。本研究使用BFAST01方法结合GLASS FVC产品数据, 检测了1982—2018年期间全国尺度的植被覆盖非线性变化趋势, 并将一元线性回归分析法得到的结果与之对比, 辅以空间格局分析, 从而全面地评估了我国地表植被覆盖长期变化趋势特征, 可以为我国生态环境保护工作和生态文明建设的开展提供相应参考。

1 数据与方法 1.1 数据

本文数据来自于全球陆表特征参量(GLASS)系列数据中的植被覆盖度产品( http://www.glass.umd.edu/), 时间分辨率为8d, 空间分辨率为0.05°。GLASS FVC产品利用广义回归神经网络算法生产而成[ 25], 它具有覆盖全球的完整长时间系列数据并且适合植被覆盖趋势变化的研究[ 11], 其精度也被证明高于同类产品[ 26]。下载得到1982—2018年共37年GLASS FVC数据, 采用均值法得到年度FVC值以供后续方法的代入。 图 1为我国1982—2018年FVC多年均值分布示意图。

图 1 中国1982—2018年植被覆盖度多年均值空间分布 Fig. 1 Spatial distribution of annual mean FVC in China from 1982—2018 FVC: 植被覆盖度Fractional vegetation cover
图选项
1.2 研究方法 1.2.1 线性回归趋势分析法

基于最小二乘法的一元线性回归方法是常用的线性趋势分析法, 可以逐像元计算时序范围内的空间分异特性来反映植被的整体空间变化规律[ 4], 但它只能够分析存在于长时间序列中的植被总体变化趋势和变化量, 往往会忽略植被长期变化过程中的细节信息, 不能够用来检测植被的阶段性变化过程、季节性变化或者突变情况, 更不能够确定植被变化发生的时刻[ 10]。其计算公式如下:

$ \text { Slope }=\frac{n \times \sum\nolimits_i^n i \times \mathrm{FVC}_i-\sum\nolimits_i^n i \sum\nolimits_i^n \mathrm{FVC}_i}{n \times \sum\nolimits_i^n i^2-\left(\sum\nolimits_i^n i\right)^2} $ (1)

式中, Slope为植被覆盖变化趋势斜率;n为研究年份数目, 本文为37;FVCi为第i年的植被覆盖度值。Slope>0表示植被覆盖增加趋势, 即呈现改善状态;Slope < 0表示植被覆盖减少趋势, 即呈现退化状态。

1.2.2 BFAST方法

BFAST利用分段线性趋势和季节性模型将时间序列分解为趋势、季节性和残差组分, 并能检测趋势和季节性成分的突变[ 27]。假设一个加法的分解模型可以迭代出与趋势和季节性相匹配的分段线性模型[ 20], 该模型的算法形式为:

$ Y_t=T_t+S_t+e_t, \quad t=1, \ldots, n $ (2)

式中, Ytt时观测到的值, Tt为趋势组分, St为季节性组分, et为残差组分。

BFAST01则是改进过的BFAST方法, 它检测趋势中最有影响力的一个突变以将趋势分为前后两段而不是多个较小趋势[ 15]。在代入年度时序数据时, 模型不用考虑季节组分而仅有趋势组分和残差组分, 并且使用不同的结构变化试验以检测趋势和突变的显著性, 如果任何一个检验发现了结构性显著变化(P < 0.05)则检测为趋势的突变点[ 21]。BFAST及BFAST01程序需要在R环境中运行, 相关程序包下载于 https://mirrors.tuna.tsinghua.edu.cn/CRAN/。将1982—2018年共计37a的中国年度FVC栅格数据构成一个时间序列后代入BFAST01程序进行检测, 设置突变检验方式为4种, 并将最小趋势段长度设置为7a, 其他参数均为默认。在进行分析时, 考虑植被覆盖变化趋势的显著性, 将计算得到前后段趋势都不显著的栅格定义为非显著趋势并单独归类。

为了进一步呈现BFAST01对于植被覆盖变化的检测效果, 本研究重点关注趋势突变点和前后两段趋势的特征, 对像元计算结果进行分类分析。综合de Jong[ 15]和Higginbottom[ 21]的研究, 最终将BFAST01检测得到的植被覆盖变化趋势类型分为6类, 如 表 1 图 2所示。

表 1 BFAST01检测得到植被覆盖变化趋势类型结果 Table 1 The different types of vegetation cover change detected by BFAST01
趋势类型名称Type name 含义Meanings
单调型增加Monotonic increase 未检测出明显突变或检测出1个明显突变, 趋势整体表现为单调性增加
单调型减少Monotonic decrease 未检测出明显突变或检测出1个明显突变, 趋势整体表现为单调性减少
中断型增加Interrupted increase 检测出1个明显突变, 趋势表现为在增加过程中受到一次明显负面干扰
中断型减少Interrupted decrease 检测出1个明显突变, 趋势表现为在减少过程中受到一次明显正面干扰
由增到减Increase to decrease 检测出1个明显突变, 趋势表现从增加向减少的趋势转换
由减到增Decrease to increase 检测出1个明显突变, 趋势表现从减少向增加的趋势转换
表选项

图 2 BFAST01检测得到植被覆盖变化趋势类型示意图 Fig. 2 The figures of different types of vegetation cover change detected by BFAST01
图选项
1.2.3 空间分布格局分析

对于BFAST01所检测出来的趋势类型结果, 为了进一步分析其像元尺度分布的规律性, 可以借用景观生态学相关的景观指数来计算和分析格局。分布格局通常指的是不同形状、大小、属性的景观斑块的空间结构特征, 是景观异质性在空间上的综合表现, 也是不同类型的生态过程在不同尺度上作用的结果[ 28]。景观指数是高度量化的可以表达分布格局的信息, 因此景观指数可以充分反映景观斑块的变化情况, 能够在杂乱的斑块镶嵌中探索出一定的分布规律。参照有关研究[ 29], 在景观指数计算软件Fragstats4.2中选取相应指标, 所选指标及其取值范围、景观生态学含义如 表 2所示。本研究借鉴了景观格局指数从景观和趋势类型两个层次来分析植被覆被变化趋势的空间格局, 在类型层面选取的指数着重考虑不同植被覆盖趋势类型的斑块特质, 而在景观层面注重分析植被覆盖变化趋势整体的分布格局和多样性。

表 2 选取景观指数的取值范围和景观生态学含义 Table 2 The selected indexs′ value range and the meanning in landscape ecology
名称
Indexs
取值范围
Value range
景观生态学含义
Implications of landscape ecology
趋势类型层面
Class metrics
NP ≥1 某一类型斑块的数量(个)
AREA_MN ≥0 某一类型斑块的平均面积(km2)
LPI 0≤LPI≤100 某一类型斑块中最大斑块占景观面积的比例(%)
PD ≥0 某一类型斑块的密度(10-2个/km2)
LSI ≥0 某一类型斑块的形状复杂程度
AI 0≤AI≤100 某一类型斑块的聚合程度(%)
景观层面
Landscape metrics
NP ≥1 各类斑块总数量
SPLIT 1≤SPLIT≤NP2 低值表示景观破碎度小, 高值表示景观破碎度大
COHESION 0≤COHESION≤100 低值表示景观较分散, 高值表示景观较聚合
CONTAG 0≤CONTAG≤100 低值表示景观连通度小, 高值表示景观连通度大
SHDI ≥0 低值表示景观类型较少, 高值表示景观类型较多
SHEI 0≤SHEI ≤1 低值表示景观优势度高, 高值表示景观优势度低
NP: 斑块个数Number of patches; AREA_MN: 平均斑块面积Mean patch area; LPI: 最大斑块指数The largest patch index; PD: 斑块密度Patch density; LSI: 景观形状指数Landscape shape index; AI: 聚合度指数Aggregation index; SPLIT: 分离度指数Splitting index; COHESION: 凝聚度指数Patch cohesion index; CONTAG: 蔓延度指数Contagion; SHDI: Shannon′s多样指数Shannon′s diversity index; SHEI: Shannon′s均匀指数Shannon′s evenness index
表选项
2 结果与分析 2.1 我国长时间尺度植被覆盖变化趋势分析 2.1.1 线性趋势

采用一元线性回归趋势分析法对我国1982—2018年FVC数据进行处理。结果如 图 3所示, 69.82%的像元呈现显著改善趋势, 而仅有4.61%的像元呈现显著退化趋势(P < 0.05)。同时, 还有25.57%的像元呈现非显著趋势。植被覆盖变化显著改善趋势在我国的北方地区呈现大范围连续分布, 而在南方地区虽然分布较为分散但也占较大面积;显著退化趋势集中分布在东南沿海地带和中东部平原地区, 尤以京津冀、长三角、闽东南和珠三角这些城市化高度发展的地区为代表;非显著趋势主要集中在内蒙古高原中东部和大兴安岭南部地区, 此外在西藏的东南部、四川盆地中部及我国南方部分地区植被覆盖变化也不显著。

图 3 1982—2018年我国FVC线性变化趋势空间分布及占比 Fig. 3 Spatial distribution and proportion of FVC linear trends in China from 1982 to 2018
图选项
2.1.2 非线性趋势

通过使用BFAST01算法对我国1982—2018年植被覆盖度数据进行计算和检验, 得到了不同的植被覆盖非线性变化趋势类型空间分布示意图, 结果如 图 4所示。其中, 20.19%的地区趋势类型为非显著。显著趋势中, 区域面积占比达到33.58%的单调型增加最多, 其次为占比22.91%的中断型增加, 第三为占比14.62%的由减到增;单调型减少、中断型减少和由增到减占比相对较少, 依次为1.82%、2.68%和4.20%。由此可见, 非线性趋势检测结果分类明确且不同类型之间也存在较大差异, 若仅仅考虑线性趋势变化, 将不能正确评估趋势发生的变化。同时, 单调型增加、中断型增加和由减到增依次是占比最多的三种趋势类型, 充分表明了我国1982—2018年期间大部分地区植被覆盖的改善状况。但是, 也有部分地区存在植被覆盖减少的趋势, 潜在的植被退化风险不可忽视。

图 4 1982—2018年我国FVC非线性变化趋势空间分布及占比 Fig. 4 Spatial distribution and proportion of FVC nonlinear trends in China during 1982 to 2018
图选项

从空间分布来看, 和线性趋势相似的是, 内蒙古中东部地区趋势集中表现为非显著, 此外在南方地区非显著趋势也存在不同程度的镶嵌分布。单调型增加是最主要的趋势类型, 其重点分布于华中地区如湖北省、河南省和陕西省, 此外青藏高原的中西部、内蒙古高原南部以及天山山脉也有集中分布, 这些地区的植被覆盖未经受明显负面干扰甚至有过短期迅速绿化从而表现出长期改善的趋势。单调型减少趋势多集中分布于东南沿海和海南沿海地带, 而这些地区的植被覆盖未曾表现出改善趋势。中断型增加趋势, 即植被覆盖在长期逐渐改善的情况下经受了某些负面干扰因而出现短期的退化, 该趋势呈现大范围集中分布的特点, 在北方地区形成以黄土高原东部、山西省、河北省和东北平原的连续分布, 在南方地区的四川盆地和云贵高原地区形成另一连续分布。中断型减少趋势类型较少, 分散分布于其他趋势内, 在藏东南地区存在小规模集中分布, 说明该地区的植被覆盖曾因外界环境条件好转出现了短期的变绿。由增到减趋势多分布于青藏高原中东部地区, 在苏、皖地区也有部分集中。由减到增趋势在青海省东部和陕西省中部存在部分集中, 在东北地区形成大范围集中分布, 该种反转实则代表了植被生态环境条件由差变好。

2.1.3 非线性趋势与线性趋势对比

在空间分布上, 可以发现两种检测方式都指示了内蒙古中东部地区植被长期未发生显著改变, 并且东南沿海地区有着明显的植被退化进程。但是线性趋势结果所指示的显著改善区域对应的非线性趋势实际上多为单调型增加、中断型增加和由减到增, 如在四川盆地的东部、重庆、华中地区和黄土高原等地的显著改善趋势则经历了一次重要突变, 整体增长趋势遭遇短期中断。而在青海、吉林和黑龙江等地区的植被改善过程则可能是原先的退化趋势转变而来, 因此用单一的线性趋势来描述植被覆盖变化掩盖了其经历的具体细节性过程。

针对线性趋势分析法和BFAST01方法得到的结果, 为了进一步展现我国长时间尺度植被覆盖非线性变化趋势与线性变化趋势检测结果的差异, 现对两者不同趋势类型的转化情况进行分析, 如 表 3所示。

表 3 BFAST01与线性趋势分析法所得结果转化表/% Table 3 Transformation table of BFAST01 and linear trend analysis
非线性趋势Nonliear trend
单调型
增加
单调型
减少
中断型
增加
中断型
减少
由增到减 由减到增 非显著趋势
线性趋势
Linear trend
非显著趋势 0.9 0.1 15 8 5 13 58
显著退化 19 36 3 10 6 12 14
显著改善 44 5 25 0.4 3 14 8
表选项

表 3可以看出, 线性趋势分析方法和BFAST01非线性趋势分析方法检测结果存在明显的差异。线性趋势中的非显著趋势只有58%仍为BFAST01方法的非显著趋势, 其次则有15%和13%被检测为中断型增加和由减到增;而显著退化趋势只有36%被BFAST01方法定义为单调型减少, 这是线性减少趋势转化最多的非线性趋势类型;与此同时, 显著改善趋势转化非线性趋势类型当中, 尽管最多的是单调型增加达到了44%, 然而中断型增加和由减到增的面积也分别高达25%和14%。因此, BFAST01方法不仅检测出了总体趋势, 还包含了中断和转换等突变现象, 这是普通线性回归方法所不能做到的, 即线性趋势方法无法揭示植被覆盖长期变化中的细节。

2.2 植被覆盖趋势突变年份

在得到我国1982—2018年植被覆盖变化不同的非线性趋势类型的同时, 提取趋势突变(即趋势中断和转换)的年份( 图 5), 并进一步对各趋势类型发生突变的年份进行统计( 表 4)。本文所定义的单调型变化包括未检测出明显突变和检测出1个明显突变的两种子类型, 其中有突变的子类型在数据处理时跟其他类型一同提取了改变时间, 并绘制在 图 5中, 对照 图 4可确定单调型增加/减少地区其趋势发生改变的时间。由 图 5可知, 植被覆盖变化趋势发生突变的时间分布呈现一定的规律性。如在南方地区趋势改变多发生于2000年之后;山西省、陕西省及黄土高原一带趋势改变多发生在20世纪90年代末期, 东北地区情况与之类似但同时也存在着发生于2000—2005年期间的大规模趋势改变;四川盆地东部及云贵高原地区在2009年左右发生了一个大范围集中性植被覆盖变化趋势改变过程。

图 5 植被覆盖非线性变化趋势改变发生的年份 Fig. 5 Year of change of the vegetation cover nonlinear trend
图选项

表 4 植被覆盖非线性变化趋势发生改变的年份分布/% Table 4 The year distribution of nonlinear change trend of vegetation cover
趋势类型
Trend type
趋势发生改变年份(每6年)占比Proportion of 6-year trend change
1988—1993 1994—1999 2000—2005 2006—2011
单调型增加Monotonic increase 1.26 30.25 18.65 49.85
单调型减少Monotonic decrease 7.62 19.87 63.58 8.94
中断型增加Interrupted increase 18.99 17.08 35.85 28.07
中断型减少Interrupted decrease 11.83 20.18 13.57 54.42
由增到减Increase to decrease 3.97 17.75 17.23 61.05
由减到增Decrease to increase 27.96 37.87 19.36 14.81
每6年各趋势发生改变总和占比Total of proportion 17.50 24.99 26.22 31.29
表选项

表 4可知, 单调型增加发生趋势突变的时间多集中于1994—1999年和2006—2011年间, 而单调型减少发生趋势突变的时间多集中于21世纪初期;中断型增加发生趋势突变的时间多集中于2000年后, 而中断型减少趋势发生突变的时间多集中于2006—2011年间;由增到减趋势发生转换的时间明显集中于2005年后, 而由减到增趋势发生转换的时间集中于20世纪90年代后期。由每6年各趋势发生改变占比可以看出, 尽管2000—2005年的改变量(26.22%)同1994—1999年(24.99%)的相差不大, 但整体上2000年后发生的改变(57.51%)多于2000年前的(42.49%), 尤其是2006—2011年间的改变量(31.29%)接近总体的三分之一, 结合各趋势占比情况可知中断型增加、由减到增和含突变的单调型增加趋势对2000年后我国植被覆盖发生的趋势改变有主导作用。

2.3 植被覆盖非线性变化趋势类型分布格局

为了分析通过BFAST01所检测出来的我国植被覆盖非线性变化趋势类型的空间分布情况, 在Fragstats 4.2软件中分别从景观和趋势类型两个不同层面计算了不同的景观指数( 表 5)。

表 5 植被覆盖非线性变化趋势的格局分析计算结果 Table 5 Results of spatial pattern analysis in nonlinear trends of vegetation cover change
类型层面
Class metrics
值Value 景观层面
Landscape metrics

Value
1 2 3 4 5 6 7
NP 8388 2572 9812 4624 5381 9935 10317 NP 51029
AREA_MN/km2 351.2419 64.6775 164.8775 40.9199 55.0892 103.8756 158.5646 SPLIT 62.5236
LPI/% 11.5349 0.0328 2.0594 0.0197 0.0219 1.7543 2.5139 COHESION 95.2971
PD/(10-2/km2) 0.1092 0.0335 0.1278 0.0602 0.0701 0.1294 0.1343 CONTAG 21.5226
LSI 139.6711 56.3099 134.0418 73.2849 82.2911 122.7132 136.4628 SHDI 1.6111
AI/% 54.0879 20.0692 46.6770 14.5028 23.3772 38.8621 42.1535 SHEI 0.8280
1、2、3、4、5、6、7依次代表单调型增加、单调型减少、中断型增加、中断型减少、由增到减、由减到增和非显著趋势类型斑块
表选项

首先, 从景观层面来看, 总共有51029个斑块, 对应着我国面积约790万km2存在植被覆盖的土地;景观分离度指数(SPLIT)较低而凝聚度指数(COHESION)接近100说明整体植被趋势类型的分布格局较聚合, 但同时蔓延度指数(CONTAG)不高说明不同趋势类型之间的连通度低, 且趋势类型分布有一定破碎度, 格局异质性较高;Shannon′s多样指数(SHDI)不高仅为1.6111, 而Shannon′s均匀指数(SHEI)0.8280接近于1, 表明了趋势类型的多样性较小, 并且主导趋势类型的优势度不高。

其次, 从趋势类型层面来看, 非显著趋势的斑块个数最多而单调减少的最少;平均斑块面积(AREA_MN)最大的是单调型增加斑块, 是中断型减少斑块的8.75倍, 说明单调型增加分布比较集中斑块面积比较大, 中断型斑块的面积小, 分布离散;最大斑块指数(LPI)表明了某类斑块中面积最大的斑块个体占总景观的比例, 单调型增加斑块的该指数为11.5349远大于其他类型斑块, 说明存在一个大范围集中连片的单调型增加趋势;斑块密度(PD)最大的是由减到增, 同时它的平均斑块面积和最大斑块面积比较大, 说明了它存在大面积连片的斑块外, 还存在大量的小斑块分布;景观形状指数(LSI)当中, 单调型增加趋势大于单调型减少, 中断型增加大于中断型减少, 由减到增大于由增到减, 说明近几十年来植被增加趋势斑块形状趋于复杂化;聚合度指数(AI)最高的是单调型增加斑块, 可知其斑块较聚合连通性最好, 相反, 单调减小的聚合度指数最小其分布最不紧凑。结合 图 4, 还可以发现北方地区不同趋势的聚集性较强, 而南方地区分布则较为分散景观离散程度十分显著, 这应该与聚合度指数较高的单调型增加、中断型增加等类型斑块在北方地区集中分布有关。

从总体趋势变化类型来看, 又可将植被覆盖的趋势归成两种类型, 即改善趋势类型(单调型增加, 中断型增加, 由减到增)和退化趋势类型(单调型减少, 中断型减少, 由增到减)。一方面, 改善趋势类型不仅其斑块数量是退化趋势类型的2—4倍, 且其平均斑块面积也远大于后者, 说明了植被覆盖发生改善斑块呈现大面积集中的分布, 而发生退化斑块呈现出小范围零散分布。另一方面, 改善型斑块的斑块密度和聚集度约是退化型斑块的两倍, 表明除了一些聚集的大斑块外, 仍然存在大量镶嵌分布的小面积斑块, 使格局具有一定破碎分布特征。同时, 改善型斑块的形状指数都比退化型斑块的高, 这表明改善型斑块的形状较为复杂, 可能受外界的干扰比较强烈。

3 讨论

自20世纪80年代以来, 中国的植被覆盖呈现显著的绿化[ 30], 尤其有研究指出2000年后仅中国的森林面积增加量就贡献了世界植被覆盖变绿量的10%[ 31], 而本研究中使用线性趋势方法相应得到植被覆盖主要呈现显著改善(69.82%)的态势。但是, 不能据此认定我国的植被覆盖一直呈现出变绿的增长形势, 因为单一的线性方法仍然无法持续真实地体现植被覆盖的变化。全球尺度上已有利用不同非线性趋势方法的研究[ 15, 32]证明, 若简单假定长时间序列内植被覆盖的增长速率不变则会掩盖其阶段性退化的事实, 而采用BFAST01方法就是为了进一步揭示我国过去几十年间植被覆盖变化过程中存在着的中断以及趋势转换等现象。首先, 先前的研究存在一定局限性, 例如金凯等[ 13]对我国1982—2015年的植被生长季NDVI值进行了研究, 将植被覆盖变化的线性增减趋势按速率进行了分级并发现72.7%的地区呈现增加趋势, 本文研究则发现其结果中的华北和南方部分地区呈现的增加趋势实际上存在一次短期显著退化, 对应为本文的中断型增加趋势, 而在东北地区呈现的减少趋势在本研究结果中实际上检测为由减到增趋势。刘宪锋等[ 23]使用了分段线性回归趋势来研究1982—2012年中国区域的植被覆盖变化情况, 注意到了1997年前后我国植被覆盖总体变化趋势的差异并以此为界探讨了不同地区在两个时间段内的增减情况, 但是将全国尺度趋势的转变点代入像元尺度趋势并不一定符合区域的实际情况, 因此BFAST01方法逐像元提取趋势断点进行趋势分解可以证明是有效的。其次, 与前人对植被覆盖分段线性提取趋势所得结果[ 23, 32 33]相比, 本研究不仅检测了植被增长和退化趋势之间的转换, 还关注到了长期变化过程中的趋势中断性突变, 最后呈现出类别化的结果更有利于我们准确定位在何时何地发生何种变化。

从本研究得到的不同非线性变化趋势类别结果来看, 单调型增加趋势(33.58%)是主导的植被覆盖变化类型并且分布最广, 这主要因一方面全球变化背景下我国植被对气候条件的改善发生显著响应, 另一方面我国实施的生态环境工程措施对于2000年来大规模的植被恢复和改善有重要影响[ 31]。中断型增加趋势占比(22.91%)也超过20%, 形成黄土高原至东北平原、四川盆地和云贵高原地区的连续分布, 发生时间多为2006—2011年期间, 有研究显示这可能是由多年降水变化引起的干扰导致[ 4]。区域尺度气候条件的不利改变会导致植被覆盖出现退化, 也就意味着气候变化会造成植被覆盖持续趋势发生中断[ 15], 若仅使用线性趋势( 图 3)来看的话则更关注其前后接续增长而忽略掉这一次短期减少过程。中断型减少(2.68%)占比多于最少的单调型减少类型(1.82%), 空间分布和趋势改变时间都较为分散, 其变化过程中的突增现象可能是水热条件改善所致[ 14], 比如黄河流域河套地区降水的阶段性增多引起植被的增加[ 34]。仍然不可忽视的是, 植被覆盖的确存在由增到减的趋势反转(4.20%), 例如在青藏高原的三江源地区植被覆盖检测到发生在2009年左右的改变, 这可能是极端天气导致了该地区的植被覆盖波动下降[ 35]。与此同时, 近年来人类活动因素对于植被覆盖变化的影响越来越显著[ 23], 人类对地表进行改造使得植被覆盖发生剧烈变化。经济发展下的快速城市化和城市用地扩张进程导致了植被覆盖的破坏, 这体现为人类侵占耕地、毁林开荒以及环境破坏等造成了植被覆盖长期的持续退化和短时间的快速衰减, 同时也是东南沿海地区植被覆盖持续减少、华中及苏皖部分地区植被覆盖2000年来由增到减的主要人为原因。不过, 为了应对不断出现的环境问题我国开展了一系列大规模生态建设和修复工程, 这遏制了植被退化趋势并使植被覆盖出现阶段性、范围性的改善。20世纪90年代末及21世纪以来, 我国开展了“三北”防护林建设、退耕还林(草)、天然林保护和退牧还草等生态工程。结合 图 4 图 5可知, 我国西北和内蒙古东部地区的植被在2000年左右发生了明显改善, 相应地东北、黄土高原、青海等地区植被覆盖也发生了较为聚集的由减到增类型(14.62%)的趋势转变, 表明生态保护、恢复工程的开展对这些增长有重要贡献[ 31, 36]

整体格局上来看, 由于面积占比较大的植被覆盖改善趋势型的斑块数量、平均斑块面积、形状指数以及聚合度等指数都大于退化型, 表明了改善型斑块呈现大范围聚集分布、形状复杂且受外界条件影响较大的特点, 这证实气候变化和人类活动的作用对我国植被覆盖影响是显著的。像元尺度上来看, 各种趋势类型破碎分布, 致使主导趋势类型的优势度下降。已有研究[ 37 38]指出外界驱动力不仅使得过去几十年来我国植被有大范围改善的情况, 区域性事件如夏季降水减少和土地利用方式的改变对植被覆盖造成干扰, 这可能促使各趋势类型形成分布不均的现象。所以, 在改善型趋势呈现大面积集中分布的特征下, 众多小面积且零散的趋势使整体分布格局形成一定的空间异质性, 显示了区域尺度上植被覆盖经历的变化实际更加多样, 相关因素的作用机制也更为复杂。相比于前人针对我国长时间植被覆盖变化趋势的研究[ 13, 23], 通过景观生态学中的方法定量分析了趋势的格局, 更客观的评估了其空间上的分布规律。

综上所述本文认为, 一方面, 在全球变暖背景下, 气温的升高使植被生长周期变长从而有利于植被绿化, 因此植被覆盖会出现持续性的增加过程, 这也是单调型增加、中断型增加和由减到增趋势类型占比大的主要原因。此外, 气候条件的波动还会引起植被覆盖变化出现突变, 导致趋势的持续性中断。另一方面, 近年来人类活动因素对于植被覆盖变化的影响越来越显著, 人类对地表进行改造从而使得植被覆盖发生变化, 城市化过程加剧了部分地区植被的持续性或阶段性退化, 而我国开展的植树造林、退耕还林等大规模生态环境工程则对于其改善有着重要影响。

4 结论

本文从植被覆盖非线性变化趋势角度出发, 通过BFAST方法结合GLASS FVC数据针对我国1982—2018年的植被覆盖情况进行了研究分析, 重点关注其长期变化过程中存在的非线性趋势及突变现象, 经过前文结果分析和讨论, 得到的主要结论如下:

(1) BFAST方法的检测结果揭示了我国植被覆盖存在的非线性变化趋势。如黄土高原等地经历了一次短期退化、青藏高原中东部等地从增加到减少趋势以及青海和东北等地区的从退化到改善过程。将所得结果进行类型划分, 其中:单调型增加类型占比最多, 达到33.58%, 主要分布在内蒙古、陕西、河南等地;单调型减少类型占比1.82%, 主要分布在东南沿海地区;中断型增加类型占比22.91%, 主要分布在四川盆地东部和华北地区;中断型减少类型占比2.68%, 主要分布在青藏高原东南部;由增到减类型占比4.20%, 主要分布在青海、陕西等地;由减到增类型占比14.62%, 主要分布在青海、吉林等地。大范围的植被覆盖增加趋势充分反映了我国过去几十年植被的改善, 但同时存在的植被覆盖减少趋势表明潜在的植被退化风险仍不可忽视。

(2) 不同趋势类型发生突变和趋势反转的时间有所差异。总体上1988—1999年间发生的改变较少, 而2000—2011年间发生的改变较多。单调型增加发生趋势突变的时间多集中于1994—1999年和2006—2011年间, 而单调型减少发生趋势突变的时间多集中于21世纪初期;中断型增加发生趋势突变的时间多集中于2000年后, 而中断型减少趋势发生突变的时间多集中于2006—2011年间;由增到减趋势发生转换的时间明显集中于2005年后, 而由减到增趋势发生转换的时间集中于20世纪90年代后期。

(3) 趋势分布格局在全国尺度上呈现一定规律性, 但整体又有一定破碎度且趋势分布异质性较大。相比于植被覆盖退化趋势类型斑块(单调型减少, 中断型减少, 由增到减), 改善趋势类型斑块(单调型增加, 中断型增加, 由减到增)呈现大聚集, 小分散的特点, 且具有复杂的形状。因此, 区域尺度上植被覆盖经历的变化实际是复杂多样的, 经受的干扰较为显著, 反映了区域性事件会深刻影响其变化过程。

本文将我国1982—2018年植被覆盖的变化情况从非线性趋势角度进行了探讨, 但仍然存在一些不确定性。有研究指出在植被覆盖呈现出长期较强的退化状态时植被指数才能检测出显著的减少趋势[ 39], 因此本研究使用BFAST01方法进行长时间植被覆盖变化趋势检测时, 可能存在将不显著的减少趋势被定义为持续性增加趋势或是一段增加趋势的情况, 有待进一步揭示植被覆盖变化的原因。

参考文献
[1]
Cui L L, Shi J. Temporal and spatial response of vegetation NDVI to temperature and precipitation in Eastern China. Journal of Geographical Sciences, 2010, 20(2): 163-176. DOI:10.1007/s11442-010-0163-4
[2]
Bohovic R, Dobrovolny P, Klein D. The spatial and temporal dynamics of remotely-sensed vegetation phenology in central Asia in the 1982-2011 period. European Journal of Remote Sensing, 2016, 49(1): 279-299. DOI:10.5721/EuJRS20164916
[3]
Piao S L, Yin G D, Tan J G, Cheng L, Huang M T, Li Y, Liu R G, Mao J F, Myneni R B, Peng S S, Poulter B, Shi X Y, Xiao Z Q, Zeng N, Zeng Z Z, Wang Y P. Detection and attribution of vegetation greening trend in China over the last 30 years. Global Change Biology, 2015, 21(4): 1601-1609. DOI:10.1111/gcb.12795
[4]
李美丽, 尹礼昌, 张园, 苏旭坤, 刘国华, 王晓峰, 奥勇, 伍星. 基于MODIS-EVI的西南地区植被覆盖时空变化及驱动因素研究. 生态学报, 2021, 41(3): 1138-1147.
[5]
穆少杰, 李建龙, 陈奕兆, 刚成诚, 周伟, 居为民. 2001—2010年内蒙古植被覆盖度时空变化特征. 地理学报, 2012, 67(9): 1255-1268.
[6]
Gitelson A A, Kaufman Y J, Stark R, Rundquist D. Novel algorithms for remote estimation of vegetation fraction. Remote Sensing of Environment, 2002, 80(1): 76-87. DOI:10.1016/S0034-4257(01)00289-9
[7]
Parmesan C, Yohe G. A globally coherent fingerprint of climate change impacts across natural systems. Nature, 2003, 421(6918): 37-42. DOI:10.1038/nature01286
[8]
王新闯, 刘文锴, 杨会军, 王世东, 张合兵. 河南省植被覆盖度及其景观格局时空变化. 水土保持通报, 2015, 35(6): 241-247, 254.
[9]
陈赛赛, 孙艳玲, 杨艳丽, 王中良. 三北防护林工程区植被景观格局变化分析. 干旱区资源与环境, 2015, 29(12): 85-90.
[10]
王巨. 基于时序NDVI植被变化检测与驱动因素量化方法研究——以河西地区为例[D]. 兰州: 兰州大学, 2020.
[11]
Mu B H, Zhao X, Wu D H, Wang X Y, Zhao J C, Wang H Y, Zhou Q, Du X Z, Liu N J. Vegetation cover change and its attribution in China from 2001 to 2018. Remote Sensing, 2021, 13(3): 496. DOI:10.3390/rs13030496
[12]
马梓策, 于红博, 曹聪明, 张巧凤, 侯丽丽, 刘月璇. 中国植被覆盖度时空特征及其影响因素分析. 长江流域资源与环境, 2020, 29(6): 1310-1321.
[13]
金凯, 王飞, 韩剑桥, 史尚渝, 丁文斌. 1982—2015年中国气候变化和人类活动对植被NDVI变化的影响. 地理学报, 2020, 75(5): 961-974.
[14]
de Jong R, Verbesselt J, Schaepman M E, Bruin S. Trend changes in global greening and browning: contribution of short-term trends to longer-term change. Global Change Biology, 2012, 18(2): 642-655. DOI:10.1111/j.1365-2486.2011.02578.x
[15]
de Jong R, Verbesselt J, Zeileis A, Schaepman M. Shifts in global vegetation activity trends. Remote Sensing, 2013, 5(3): 1117-1133. DOI:10.3390/rs5031117
[16]
Phillips P C B, Moon H R. Linear regression limit theory for nonstationary panel data. Econometrica, 1999, 67(5): 1057-1111. DOI:10.1111/1468-0262.00070
[17]
Wan L, Liu H Y, Gong H B, Ren Y J. Effects of climate and land use changes on vegetation dynamics in the Yangtze River Delta, China based on abrupt change analysis. Sustainability, 2020, 12(5): 1955. DOI:10.3390/su12051955
[18]
Geng L Y, Che T, Wang X F, Wang H B. Detecting spatiotemporal changes in vegetation with the BFAST model in the Qilian Mountain region during 2000-2017. Remote Sensing, 2019, 11(2): 103. DOI:10.3390/rs11020103
[19]
Ma J N, Zhang C, Guo H, Chen W L, Yun W J, Gao L L, Wang H. Analyzing ecological vulnerability and vegetation phenology response using NDVI time series data and the BFAST algorithm. Remote Sensing, 2020, 12(20): 3371. DOI:10.3390/rs12203371
[20]
刘宝柱, 方秀琴, 何祺胜, 荣祁远. 基于MODIS数据和BFAST方法的植被变化监测. 国土资源遥感, 2016, 28(3): 146-153.
[21]
Higginbottom T P, Symeonakis E. Identifying ecosystem function shifts in Africa using breakpoint analysis of long-term NDVI and RUE data. Remote Sensing, 2020, 12(11): 1894.
[22]
Bernardino P N, de Keersmaecker W, Fensholt R, Verbesselt J, Somers B, Horion S. Global-scale characterization of turning points in arid and semi-arid ecosystem functioning. Global Ecology and Biogeography, 2020, 29(7): 1230-1245.
[23]
刘宪锋, 朱秀芳, 潘耀忠, 李宜展, 赵安周. 1982—2012年中国植被覆盖时空变化特征. 生态学报, 2015, 35(16): 5331-5342.
[24]
苏艳琴, 赖日文, 闫琦, 余莉莉, 李红霖, 赖炽敏. 植被覆盖度提取及景观格局分析. 森林与环境学报, 2018, 38(2): 164-170.
[25]
梁顺林, 程洁, 贾坤, 江波, 刘强, 刘素红, 肖志强, 谢先红, 姚云军, 袁文平, 张晓通, 赵祥. 陆表定量遥感反演方法的发展新动态. 遥感学报, 2016, 20(5): 875-898.
[26]
Jia K, Yang L Q, Liang S L, Xiao Z Q, Zhao X, Yao Y J, Zhang X T, Jiang B, Liu D Y. Long-term global land surface satellite (GLASS) fractional vegetation cover product derived from MODIS and AVHRR data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(2): 508-518.
[27]
Verbesselt J, Hyndman R, Newnham G, Culvenor D. Detecting trend and seasonal changes in satellite image time series. Remote Sensing of Environment, 2010, 114(1): 106-115.
[28]
傅伯杰, 陈利顶, 马克明. 景观生态学原理及应用. 北京: 科学出版社, 2001.
[29]
陈文波, 肖笃宁, 李秀珍. 景观指数分类、应用及构建研究. 应用生态学报, 2002, 13(1): 121-125.
[30]
Tucker C J, Slayback D A, Pinzon J E, Los S O, Myneni R B, Taylor M G. Higher northern latitude normalized difference vegetation index and growing season trends from 1982 to 1999. International Journal of Biometeorology, 2001, 45(4): 184-190.
[31]
Chen C, Park T, Wang X H, Piao S L, Xu B D, Chaturvedi R K, Fuchs R, Brovkin V, Ciais P, Fensholt R, Tømmervik H, Bala G, Zhu Z C, Nemani R R, Myneni R B. China and India lead in greening of the world through land-use management. Nature Sustainability, 2019, 2(2): 122-129.
[32]
Pan N Q, Feng X M, Fu B J, Wang S, Ji F, Pan S F. Increasing global vegetation browning hidden in overall vegetation greening: insights from time-varying trends. Remote Sensing of Environment, 2018, 214: 59-72.
[33]
Chu H S, Venevsky S, Wu C, Wang M H. NDVI-based vegetation dynamics and its response to climate changes at Amur-Heilongjiang River Basin from 1982 to 2015. Science of the Total Environment, 2019, 650: 2051-2062.
[34]
黄建平, 张国龙, 于海鹏, 王闪闪, 管晓丹, 任钰. 黄河流域近40年气候变化的时空特征. 水利学报, 2020, 51(9): 1048-1058.
[35]
徐嘉昕, 房世波, 张廷斌, 朱永超, 吴东, 易桂花. 2000—2016年三江源区植被生长季NDVI变化及其对气候因子的响应. 国土资源遥感, 2020, 32(1): 237-246.
[36]
唐见, 曹慧群, 陈进. 生态保护工程和气候变化对长江源区植被变化的影响量化. 地理学报, 2019, 74(1): 76-86.
[37]
Li Y, Piao S L, Li L Z X, Chen A P, Wang X H, Ciais P, Huang L, Lian X, Peng S S, Zeng Z Z, Wang K, Zhou L M. Divergent hydrological response to large-scale afforestation and vegetation greening in China. Science Advances, 2018, 4(5): eaar4182.
[38]
Piao S L, Wang X H, Ciais P, Zhu B, Wang T, Liu J. Changes in satellite-derived vegetation growth trend in temperate and boreal Eurasia from 1982 to 2006. Global Change Biology, 2011, 17(10): 3228-3239.
[39]
Wessels K J, van den Bergh F, Scholes R J. Limits to detectability of land degradation by trend analysis of vegetation index data. Remote Sensing of Environment, 2012, 125: 10-22.

PHP网站源码荷坳优化同乐网页制作深圳网站搭建福永网络广告推广同乐关键词排名包年推广民治网站关键词优化西乡关键词按天扣费吉祥seo优化沙井网站设计广州阿里店铺托管横岗关键词按天收费龙华网站排名优化惠州百度网站优化吉祥关键词按天收费广州高端网站设计双龙百搜标王罗湖网页设计吉祥外贸网站建设盐田百姓网标王推广坪山百度网站优化同乐模板网站建设宝安网站制作设计丹竹头百搜标王宝安阿里店铺托管深圳外贸网站制作福田网站优化永湖如何制作网站布吉百搜词包大鹏网站优化按天计费福永网站设计歼20紧急升空逼退外机英媒称团队夜以继日筹划王妃复出草木蔓发 春山在望成都发生巨响 当地回应60岁老人炒菠菜未焯水致肾病恶化男子涉嫌走私被判11年却一天牢没坐劳斯莱斯右转逼停直行车网传落水者说“没让你救”系谣言广东通报13岁男孩性侵女童不予立案贵州小伙回应在美国卖三蹦子火了淀粉肠小王子日销售额涨超10倍有个姐真把千机伞做出来了近3万元金手镯仅含足金十克呼北高速交通事故已致14人死亡杨洋拄拐现身医院国产伟哥去年销售近13亿男子给前妻转账 现任妻子起诉要回新基金只募集到26元还是员工自购男孩疑遭霸凌 家长讨说法被踢出群充个话费竟沦为间接洗钱工具新的一天从800个哈欠开始单亲妈妈陷入热恋 14岁儿子报警#春分立蛋大挑战#中国投资客涌入日本东京买房两大学生合买彩票中奖一人不认账新加坡主帅:唯一目标击败中国队月嫂回应掌掴婴儿是在赶虫子19岁小伙救下5人后溺亡 多方发声清明节放假3天调休1天张家界的山上“长”满了韩国人?开封王婆为何火了主播靠辱骂母亲走红被批捕封号代拍被何赛飞拿着魔杖追着打阿根廷将发行1万与2万面值的纸币库克现身上海为江西彩礼“减负”的“试婚人”因自嘲式简历走红的教授更新简介殡仪馆花卉高于市场价3倍还重复用网友称在豆瓣酱里吃出老鼠头315晚会后胖东来又人满为患了网友建议重庆地铁不准乘客携带菜筐特朗普谈“凯特王妃P图照”罗斯否认插足凯特王妃婚姻青海通报栏杆断裂小学生跌落住进ICU恒大被罚41.75亿到底怎么缴湖南一县政协主席疑涉刑案被控制茶百道就改标签日期致歉王树国3次鞠躬告别西交大师生张立群任西安交通大学校长杨倩无缘巴黎奥运

PHP网站源码 XML地图 TXT地图 虚拟主机 SEO 网站制作 网站优化