APP下载

冻结季沱沱河源多年冻土区活动层土壤水分含量分析

2022-06-19李智斌刘广岳邹德富汪凌霄杜二计胡国杰周华云幸赞品赵建婷殷秀峰

冰川冻土 2022年1期
关键词:土壤水分含水量植被

李智斌, 赵 林,2, 刘广岳, 邹德富, 汪凌霄, 杨 斌, 杜二计,胡国杰, 周华云, 王 翀, 幸赞品, 赵建婷, 殷秀峰,

迟鸿飞5, 谭昌海4, 陈 文4

(1.南京信息工程大学地理科学学院,江苏南京 210044; 2.中国科学院大学,北京 100049; 3.中国科学院西北生态环境资源研究院,甘肃兰州 730000; 4.中国地质调查局自然资源综合调查指挥中心,北京 100055; 5.中国科学院青藏高原研究所,北京 100101)

0 引言

多年冻土是寒冷气候的产物,其中存储了大量的地下冰,地下冰的变化能够在一定程度上调节地区的水循环过程[1]。青藏高原现有多年冻土的面积约为106×104km2,约占整个高原面积的40%[2]。活动层是指位于多年冻土之上夏季融化、冬季冻结的土层[3],是多年冻土与大气之间进行水热能量交换的过渡层,其水热变化能在一定程度上反映多年冻土区地气能量交换状况[4-5]。土壤水分的变化增加了陆面和大气相互作用的复杂性,对气候变化的研究与预测起着非常重要的作用[6-9]。

多年冻土区活动层内部较高的土壤水分条件不仅是维系高寒生态系统稳定的基本条件,也通过其水分的变化和冻结融化,形成多年冻土区独特的地气能水交换过程和径流形成过程。准确地描述活动层内部水分分布情况无论是对于高原上的水文研究、生态研究还是陆面过程研究,都是不可忽略的重要一环[10-12]。长期以来,多年冻土区水文及气象数据的获取都是青藏高原研究的难点,近年来被广泛应用于陆面过程模式和遥感研究的几个主要土壤水分数据分别来自于玛曲、那曲和青藏公路沿线地区,其中前两个地区属于季节冻土区,而后者的监测多分布于青藏公路沿线,对于整个青藏高原来说,已有站点数量较少且分布不均,基于站点尺度的青藏高原多年冻土区土壤水分遥感反演和模式模拟研究不能满足精细化的空间数据研究和验证需求。

实测数据的缺乏限制了青藏高原多年冻土区活动层土壤水分的研究,通过卫星遥感手段反演的土壤水分数据及陆面过程模式同化得到的再分析土壤水分数据也由于缺乏实测数据而存在一定的误差[13-15]。全球陆面过程模式同化数据GLDAS 与ERA5 再分析数据近些年在高原地区应用广泛,Cheng 等[16]在玛曲、那曲及阿里等地区将ERA5 再分析土壤水分产品与多源卫星土壤水分产品进行对比评估,发现目前ERA5 的土壤水分精度要低于卫星遥感产品,但是ERA5 再分析资料具有良好的时空连续性;王婉昭等[17]利用GLDAS 数据分析了青藏高原及其周边地区气温和降水的时空分布和变化特征,检验了其在高原地区的适用性。邓明珊等[18]选取青藏高原中部那曲地区试验点的土壤水分数据与GLDAS 中陆面过程模型模拟的土壤水分产品进行对比分析,发现NOAH 陆面模式资料在青藏高原地区的适用性较好。总的来说,虽然两种土壤水分产品都能一定程度上表达区域或全球尺度土壤水分分布状况,但是对于较小区域尺度来说,其低空间分辨率的数据难以描述土壤水分分布的空间分异性,也无法表示不同下垫面对于土壤水分分布的影响。再分析数据的不足,与实测数据的缺乏有较大的关系,因此,扩充多年冻土区的土壤水分实测数据对于更精确的青藏高原区域陆面生态水文过程模拟具有重要意义。

本研究基于“第二次青藏高原综合科学考察研究——多年冻土对亚洲水塔的影响”专题,于2020年10 月12 日至11 月22 日开展了长江源多年冻土区野外科学考察工作。科考区域主要位于长江源区的沱沱河源头、唐古拉山脉格拉丹冬北坡。气候变化背景下,长江源区气温增温显著,蒸发量、径流量总体呈增加趋势[19],多年冻土也正在退化,活动层内部水分状态和动态过程也随之发生了变化,准确获取活动层的水分情况是评估和预估其未来变化的前提。本次科考通过钻探与坑探两种方式采集了大量土壤样品来获取土壤含水量数据,本文基于对不同植被类型、不同坡位、不同坡向以及不同冻结类型的活动层总含水量、剖面含水量分布及其空间异质性进行了分析和讨论。该研究可以为长江源地区的水文研究、生态研究以及陆面过程研究提供参考与验证。此外,本研究选用GLDAS-Noah 同化数据和ERA5-Land 再分析资料与实测资料进行对比分析,评估两种资料的土壤水分数据在此区域的适用性,为青藏高原土壤水分的后续研究提供参考。

1 数据和方法

1.1 研究区概况

青藏高原唐古拉山脉格拉丹冬雪山北坡的沱沱河源区以及通天河源区属于长江源区。沱沱河源区位于青海省的西南部、唐古拉山主峰格拉丹冬雪山的北侧,流域面积约为1.58×104km2[20](图1)。沱沱河流域内山体多,地势高,海拔高度在4 489~6 468 m之间[19],气候干冷,土壤冻结期长,属于连续多年冻土区[21]。流域内地表植被类型主要以高寒草甸和高寒草原为主,部分区域高寒草甸退化明显,高寒沼泽草甸分布较少。

1.2 数据

1.2.1 实测数据

野外采样过程主要参考《多年冻土调查手册》,利用钻探与坑探两种方法进行[22]。本次野外考察共布设钻孔32 个(钻孔编号为TTH01~TTH32),探坑20个。钻探采样深度范围为0~20 cm、20~50 cm、50~80 cm、80~100 cm、100~150 cm、150~200 cm、200~300 cm、300~400 cm、400~500 cm;坑探采样深度范围为0~10 cm、10~20 cm、20~30 cm、30~50 cm、50~70 cm、70~100 cm、100~150 cm、150~200 cm;本研究选取数据以钻探数据为主,坑探为辅。科考过程中采用物探与钻探相结合的方式来确定多年冻土基本特征,结果显示,26 个钻孔所在位置存在多年冻土,其他6 个钻孔位于融区内,属于季节冻土(图1)。

图1 研究区及采样点Fig.1 Study area and sampling sites

样品的采集是环刀法进行的,为避免土壤水分损失,在采样现场进行称重,获取湿土重量及湿容重。随后带回实验室在烘箱内105 ℃烘干24 h,最后称重计算质量含水量及干容重。如式(1):

式中:W为质量含水量;Mw为湿土质量;Md为干土质量。

土壤干容重计算如式(2):

式中:Bd为土壤干容重;Vs为土体积。

土壤体积含水量V的计算可以由质量含水量与干容重计算而来,如式(3):

由于此次野外考察区域属于高山高寒地区,表层土壤碎石含量较高,样品采集较为困难,环刀采样法在部分采样点无法使用,因此只在12个采样点采集了容重数据。在分析活动层含水量及剖面土壤水分变化时,均采用质量含水量,而与再分析数据对比分析时,则使用体积含水量。

1.2.2 SRTMTPI 90 m坡位数据

坡位指数(TPI),由Weiss[23]提出,其目的是用来描述地形部位的一个地形参数。坡位指数基本思路是:用某点高程与其周围一定范围内平均高程的差,结合该点的坡度,来表示其在坡面上所处的部位,从而确定研究目标点与其周围地形的位置关系。TPI 与坡位等级对应情况如表1 所示,SD 代表邻域高程标准差。

表1 坡位类型判断标准Table 1 The Criteria for judging TPI

本研究采用的坡位数据来源于中国科学院计算机网络信息中心地理空间数据云平台,是由SRTM V4.1 90 m 分辨率数字高程数据加工而成,在此数据集中,坡位类型已采用分类代码(1~6)进行标识。

1.2.3 SRTM

SRTM 由美国太空总署和国防部国家测绘局联合发布,通过对获取的雷达影像数据进行处理,制成了数字地形高程模型,即现在的SRTM 地形产品数据。本研究采用SRTM 30 m 数据集,分辨率为30 m,用于对研究区进行阴坡与阳坡的分类。

1.2.4 ERA5-Land再分析数据

ERA5 是欧洲中期天气预报中心发布的第五代再分析资料。相较于第四代ERA-Interim 数据,该数据将更多的历史观测数据利用到先进的数据同化模式系统中,用以估计更为准确的大气状况,并采用了4DVar 数据同化方法,其数据更新速度、时间分辨率与空间分辨率均有了较大提升[24-25]。

ERA5-Land 是在ERA5 的基础上对陆地区域进行重新模拟而成的再分析数据集,相较于ERA5 具有更高的空间分辨率。本研究根据土壤样品采集时间提取研究区2020 年10 月12 日至11 月22 日的ERA5-Land 四层再分析土壤水分数据:Layer1(0~7 cm)、Layer2(7~28 cm)、Layer3(28~100 cm)以及Layer4(100~289 cm),空间分辨率为0.1°×0.1°。

1.2.5 GLDAS-Noah同化数据

全球陆面数据同化系统GLDAS 的数据产品目前在陆面水文过程研究方面应用广泛,该数据集更新速度较快,可提供1979 年至今的同化数据资料。GLDAS 调用了4 个陆面过程模型:Noah(0.25°×0.25°)、CLM、VIC 和Mosaic(1°×1°)。本研究筛选使用Noah 模式驱动GLDAS 的4 层土壤水分资料:Layer1(0~10 cm)、Layer2(10~40 cm)、Layer3(40~100 cm)以及Layer4(100~200 cm),时间段为2020年10月12日至11月22日,时间分辨率为3 h。

1.3 方法

本研究以野外考察实测数据为基础,结合SRT⁃MTPI 90 m 坡位因子数据、SRTM 数据及采样点环境因子记录,分析区内不同植被类型、不同坡位、不同坡向以及不同冻土类型的活动层土壤含水量情况,同时与GLDAS-Noah 同化产品、ERA5-Land 再分析资料进行对比,评估两种土壤含水量数据在此区域的适用性,方法流程如图2所示。

图2 研究方法流程Fig.2 The flowchart of research methodology

本研究在分析不同植被类型、坡位类型、坡向类型以及冻土类型的活动层土壤含水量特征时采用了箱线图展示。箱线图主要包括6 个部分:下边缘、下四分位数、中位数、上四分位数、上边缘以及超出上下边缘的异常值。箱线图可用于反映原始数据分布的特征,还可以进行多组数据分布特征的比较。由于分类原因,部分类别数据量较少,故而箱线图的几个部分会呈现出重合现象。

1.3.1 植被类型分类

植被状况是地表土层中的水分、热量和养分动态平衡过程的一种外在表现,不同植被类型下的土壤持水能力和土壤水热动态运输过程也不同[26]。该研究考虑了植被和土壤含水量的协同作用关系,对比分析了不同植被类型活动层土壤含水量的差异。此次野外考察以现场识别的方式判别了研究区内各个采样点所属植被类型,其中,高寒沼泽草甸采样点2个、高寒草甸11个、高寒草原12个,不同植被类型的孔位分布如图3所示。

1.3.2 坡位与坡向分类

坡位通过影响土壤水分下渗迁移路径从而改变土壤含水量[27-28]。该研究区域内山地众多、地形复杂度高,坡位因子对活动层土壤含水量的影响也尤为重要。本研究利用SRTMTPI 90 m 数据集,结合钻探与坑探编目,将采样点的坡位分为四类:上坡位、中坡位、平坡位和下坡位。其中,属于上坡位类型的采样点2个,中坡位3个,平坡位18个,下坡位2个。不同坡位类型的采样点分布如图3所示。

图3 采样点的植被类型和坡位类型Fig.3 Slope position type and vegetation type of sampling sites

不同坡向接收的太阳辐射量和风力作用强度不同,导致降水量、蒸散发量乃至土壤物质组成均会存在一定差异,进而对活动层水分特征形成影响。本文对研究区采样点的坡向进行计算,把坡向为0°~90°和270°~359°所有偏北的坡定义为阴坡,90°~270°之间所有偏南的坡定义为阳坡。采样点的归类结果为,阳坡10个、阴坡15个。

2 结果与分析

2.1 不同植被类型活动层水分差异分析

三种植被类型活动层总含水量及活动层剖面含水量分布特征显示如图4。

图4 不同植被类型的活动层土壤含水量Fig.4 Soil moisture content of active layer under different vegetation types:the total soil moisture content of active layer under different vegetation types(a);the soil moisture content of active layer profile in alpine marsh meadow(b);the soil moisture content of active layer profile in alpine meadow(c);the soil moisture content of active layer profile in alpine steppe(d)

对采样点探坑编目记录和质量含水量测量结果的统计分析显示,该区域采样点的活动层平均厚度约为2.72 m,活动层土壤含水量平均约为14.0%,不同植被类型活动层含水量随深度的变化存在一定差异。

在空间分布上,不同植被类型活动层总含水量差异较大[图4(a)],高寒沼泽草甸活动层土壤含水量最高,平均为17.7%,在10 cm 处土壤含水量达到43.0%,其次为高寒草甸植被,约为16.8%,高寒草原植被的土壤最干,约为10.4%。

在垂直剖面上,高寒沼泽草甸区的活动层表层土壤含水量高于深层;高寒草甸的活动层含水量呈现表层高、中间低和底部高的特征,不同观测点10 cm、175 cm 及250 cm 深度处土壤含水量的空间差异性较大;高寒草原植被类型活动层土壤水分在垂直剖面上的分布与高寒草甸相似,但土壤含水量在垂直剖面上的变化幅度较小,这可能与该类区域的活动层平均厚度较大有关(约为3 m)。

2.2 多年冻土区与季节冻土区水分差异分析

多年冻土区与季节冻土区的土壤含水量也存在着较大差异[29]。通过对比研究区内高寒草原类型下多年冻土区与融区(季节冻土)浅表层0~350 cm 深度范围的质量含水量发现,多年冻土区含水量相对较高[图5(a)],约为10.5%,季节冻土区则较低,约为6.1%。在垂直剖面上,随着深度的增加,两者含水量均呈现出先增大后减小再增大的趋势,但多年冻土区的这种变化趋势更为显著[图5(b),5(c)]。

图5 不同冻土类型土壤含水量Fig.5 Soil moisture content of different frozen ground types:0~350 cm total soil moisture content of different types of frozen soil(a);the soil moisture content of soil profile in permafrost region(b);the soil moisture content of soil profile in seasonally frozen ground(c)

2.3 不同坡位活动层水分差异分析

对比不同坡位类型活动层质量含水量[图6(a)]表明,上坡位活动层含水量最高,约21.9%;中坡位和下坡位的含水量相近,分别为14.2%和14.4%,平坡位最低,为13.1%。

在垂直剖面上,上坡位[图6(b)]和中坡位[图6(c)]表层含水量的空间差异性较大,而平坡位[图6(d)]类型区整个剖面的含水量都存在一定的空间差异性,下坡位[图6(e)]的空间差异性最小。上坡位与下坡位类型的活动层含水量变化趋势相似,均呈现出由表层向深部逐渐增加的趋势;中坡位10 cm 处含水量最高,约达22.9%,随深度的增加呈明显下降趋势;平坡位类型的多年冻土活动层含水量随深度变化不大。

图6 不同坡位类型的活动层土壤含水量Fig.6 Soil moisture content of active layer with different slope types:the total soil moisture content of active layer with different TPI(a);the soil moisture content of active layer profile in upper slope(b);the soil moisture content of active layer profile in middle slope(c);the soil moisture content of active layer profile in flat slope(d);the soil moisture content of active layer profile in lower slope(e)

2.4 不同坡向活动层水分差异分析

对比不同坡向活动层质量含水量[图7(a)]表明,位于阳坡的活动层整体含水量为15.6%,高于阴坡的活动层含水量,为13.1%。阳坡位置的多年冻土活动层厚度平均约为2.73 m,阴坡约为2.71 m,对比阴坡与阳坡的活动层剖面含水量分布情况[图7(b),7(c)],发现阴坡与阳坡活动层剖面含水量变化情况相似,均呈现出先减小后增大再减小再增大的变化趋势。

图7 不同坡向类型的活动层土壤含水量Fig.7 Soil moisture content of the active layer in different slope directions:the total soil moisture content of active layer on shady slope and sunny slope(a);the profile soil moisture content of active layer on sunny slope(b);the profile soil moisture content of active layer on shady slope(c)

2.5 GLDAS-Noah 同 化 产 品 与ERA5-Land 再 分析资料的适用性评估

本研究将两种数据产品(GLDAS-Noah、ERA5-Land)与实测资料不同深度土壤含水量进行精度对比,由于两种数据产品的空间分辨率相较于实测数据来说均较高,故在此选用12个采样点不同深度体积含水量的均值来与两种数据产品进行对比,根据深度进行加权平均计算得到四个深度范围(0~50 cm,50~100 cm,100~200 cm,0~200 cm)的土壤含水量。发现GLDAS-Noah 土壤水分产品在此区域的适用性更好,无论是不同剖面深度的土壤含水量,还是0~200 cm 土壤总含水量,均优于ERA5-Land 再分析资料。就0~200 cm 土壤总含水量而言(表2),实测数据的含水量约为23.2%,GLDAS-No⁃ah 约为24.2%,ERA5-Land 约为50.6%,由此可见,与实测数据相比,GLDAS-Noah 同化水分数据的误差较小,与实测数据对比的误差在10%以内,而ERA5-Land再分析资料存在明显高估。

表2 不同类型数据土壤体积含水量Table 2 The comparison of soil volumetric moisture content for different types of data

对于0~200 cm 土壤剖面含水量变化而言(图8),ERA5-Land 与GLDAS-Noah 的土壤含水量均呈现出随深度增加而增大的趋势,而实测数据则先增大后减小。

图8 实测资料与再分析资料的土壤含水量对比Fig.8 The comparison of soil moisture content between actual measured data and reanalysis data:the comparison of total soil moisture content of measured data and reanalysis data(a);the soil moisture content of shallow soil profile of measured data(b);the soil moisture content of shallow soil profile of GLDAS-Noah(c);the soil moisture content of shallow soil profile of ERA5-Land(d)

图9 展示了两种土壤水分产品数据在0~30 cm及100~200 cm 深度范围的空间分布特征。ERA5-Land 再分析土壤水分数据与实测水分数据空间分布存在较大区别,大部分地区都存在显著的高估。GLDAS-Noah 水分产品虽然在0~30 cm 处也存在一定的高估与低估,但误差均较小,且在100~200 cm处具有较高的一致性。

图9 实测资料与再分析资料土壤含水量空间分布Fig.9 The spatial distribution of soil moisture content of measured data and reanalysis data:0~30 cm measured data and ERA5-Land(a);100~200 cm measured data and ERA5-Land(b);0~30 cm measured data and GLDAS-Noah(c);100~200 cm measured data and GLDAS-Noah(d)

图10 统计了两种土壤水分产品数据在高寒草甸和高寒草原两种植被类型区与实测数据相比情况。相较于实测数据,ERA5-Land 数据均存在明显的高估,GLDAS-Noah 数据在这两种植被类型区也存在一定的高估,但是误差范围较小,且在高寒草甸地区的准确性要高于高寒草原地区。

图10 两种植被类型下的实测数据与再分析水分产品对比Fig.10 The comparison of measured data and reanalysis of moisture products under two vegetation types

3 讨论

沱沱河源区土壤含水量与植被发育情况存在明显的正相关关系,这主要由于植被发育越好,土壤持水性也越好,从而使得表层土壤含水量较高。冻结季活动层内土壤含水量的垂直剖面分布特征主要受植被类型、土壤质地和水分迁移的影响[30]。研究区内高寒草甸和高寒草原植被类型的活动层内部土壤含水量呈现表层和底部高与中间低的特征,表明在区域尺度,由于气候条件相似,水分迁移是影响活动层内部土壤水分分布的主要原因;但植被发育状况与活动层内部的土壤水分是相互影响的,高寒沼泽草甸植被类型由于其持水性较好,表层土壤含水量大,影响冻结季水分迁移的量级和时间,因此高寒沼泽草甸植被类型的活动层内部土壤水分分布特征主要是植被类型影响[31-32]。整体上,研究区内多年冻土活动层平均厚度约为2.72 m,高于青藏高原2.30 m 的平均活动层厚度[33],结合采样点的环境编目记录,表明沱沱河源区多年冻土正处于退化状态,一定程度上影响了活动层的土壤含水量[34]。

对比不同坡位的活动层水分分布特征,发现该研究区的坡位水分特征与目前大部分研究得出的结论有所差异,此区域的坡位与土壤含水量并没有明显的相关性[35]。研究区内地形以山地为主,受冻结季水分迁移的影响,土壤水分以自下而上先向冻结锋面迁移为主,山体坡面活动层土壤水分侧向流动有所减弱[36],因此上坡位土壤含水量高于下坡位[27]。此外,土壤质地可能也是影响不同坡位活动层含水量的因素,由于重力作用,上坡位粗颗粒土壤及砾石会向中坡位及平坡位滚动迁移,而中坡位区域的大部分砾石也会向平坡位迁移,所以平坡位土壤粒径大于中坡位,其持水性弱于中坡位,导致其含水量也低于中坡位[37];而下坡位地区坡度较小,地表起伏度也较小,上坡位及中坡位区域的砾石及粗颗粒土壤无法达到该区,但坡上细颗粒土壤部分会随着流水作用到达该区域,导致该区域土壤质地较为紧密,持水性较强,因此,下坡位的土壤含水量也相对较高。

GLDAS-Noah 同化产品和ERA5-Land 再分析数据与实测数据在剖面水分变化方面相比均存在误差,其原因可能是冻结季活动层内部的水热迁移以及不同坡位类型土壤质地存在差异[38],从而引起同化系统模式的冻融参数化方案存在一定的问题,造成其对土壤的冻融过程模拟得不准确,难以捕捉到土壤含水量随深度的变化。同时,研究区内各地理要素具有较强的空间异质性,对两种土壤水分产品的精度评估也会造成一定的影响。

4 结论

本研究基于“青藏高原第二次综合科学考察研究——多年冻土对亚洲水塔的影响”专题,利用野外综合调查所采集的多年冻土区土壤水分数据,从不同植被类型、不同坡向类型、不同坡位类型以及不同冻土类型四个方面分析了该区域冻结季多年冻土活动层及浅表层0~350 cm 土壤剖面水分的分布特征与整体含水量。

(1)研究表明,此区域多年冻土存在退化现象,活动层厚度平均约2.72 m,活动层土壤质量含水量约为14.0%。活动层土壤含水量与植被类型存在明显的相关性,高寒沼泽草甸地区土壤含水量最高,其次为高寒草甸地区,高寒草原地区最低。同时,除高寒沼泽草甸外,高寒草甸和高寒草原类型下的多年冻土区活动层剖面含水量分布特征均呈现出先增大后减小的趋势。

(2)研究区内多年冻土区与季节冻土区的土壤含水量存在明显区别,就高寒草原地区浅表层土壤含水量而言,多年冻土区高于季节冻土区。多年冻土区浅表层土壤剖面含水量呈现随着深度的增加先增大后减小再增大的趋势,季节冻土区也存在此趋势但没有多年冻土区显著。

(3)坡位和坡向对于活动层土壤含水量的影响也较为明显。坡位水分的空间分布一定程度上受到冻结季水分迁移与土壤质地的影响,就活动层整体含水量而言,上坡位>下坡位>中坡位>平坡位。对比不同坡位类型的活动层土壤含水量变化情况,上坡位与下坡位土壤含水量均呈现出由浅入深的增加趋势,中坡位主要呈下降趋势,平坡位类型下的多年冻土活动层土壤含水量整体变化不大。阳坡活动层含水量高于阴坡,且两者的活动层剖面水分变化趋势一致。(4)评估GLDAS-Noah 同化水分数据与ERA5-Land 再分析水分数据在该区的适用性情况,发现前者的准确度更高,误差在10%以内,且表层土壤含水量的模拟精度最高,模拟精度随深度的增加而降低,但两种数据的剖面含水量变化情况描述均与实测不符,造成误差的原因可能与同化系统的模式冻融参数化方案不准确以及冻结期水分迁移等因素有关。

猜你喜欢

土壤水分含水量植被
基于高分遥感影像的路域植被生物量计算
喀斯特坡耕地块石出露对土壤水分入渗的影响
基于根系加权土壤水分有效性的冬小麦水分生产函数
磷素添加对土壤水分一维垂直入渗特性的影响
北京土石山区坡面土壤水分动态及其对微地形的响应
基于Landsat-8影像数据的北京植被覆盖度时空特征与地形因子的关系
眼干的人别选高含水量隐形眼镜
追踪盗猎者
第一节 主要植被与自然环境 教学设计
数字说