不同地表条件下青藏公路对多年冻土的热影响差异研究

周 保， 魏 刚， 张永艳， 魏赛拉加， 蒋观利

（1.青海省地质环境监测总站，青海西宁810008；

2.青海省环境地质勘查局，青海西宁810008；

3.中国科学院西北生态环境资源研究院冻土工程国家重点实验室，甘肃兰州730000）

青藏公路东起青海省西宁市，西止于西藏自治区拉萨市，穿越约550 公里青藏高原连续多年冻土区，是世界上海拔最高、线路最长的柏油公路。青藏公路始建于1950 年，于1954 年通车，并于1974—1986 年间开始第一次全面提升改建，加铺沥青路面。由于对多年冻土缺乏足够的了解，在青藏公路最初的设计和建设中并没有将多年冻土与公路的相互作用充分纳入考虑，因此在公路开通运营之后由于多年冻土变化而导致公路病害频发。近年来，青藏高原的持续发展带来更多的社会经济活动，同时对道路等基础设施的安全稳定标准提出了更高的要求，这就使得道路工程与多年冻土间的相互作用问题变得极为重要［1-3］。

目前，在气候变化背景下，青藏高原多年冻土显著退化，其地温升高、活动层厚度增加［4-5］，引起广泛关注［6-9］。与此同时，随着经济社会的发展，青藏高原多年冻土区各种人类活动越加频繁，对于多年冻土区工程设施的稳定运营要求更高，这就使得各种相应措施被应用在工程建设中，以减少对多年冻土的扰动，保障工程设施的稳定。在多年冻土区修建道路时，通过考虑路基结构和路面形式，可以有效减少道路和下伏多年冻土之间的相互作用，例如采用块碎石路基［10-11］、通风管路基［12-13］、隔热路基［14］和水泥混凝土路面［15-16］等。此外，在道路修建过程中，确保合理的选线和精心的施工，也可以有效减少工程病害，保证工程质量［17-18］。在青藏高原有诸多关于工程活动对多年冻土产生影响的研究，这些研究侧重于工程和多年冻土的稳定性，探讨二者间的相互作用［19-21］。研究结果表明，工程活动对多年冻土环境具有重大影响，导致其退化。然而，随着地表条件的不同，工程活动对多年冻土环境的影响也各不相同，这主要是因为地表条件，尤其是植被条件，与多年冻土的热状态之间存在着密切的关系［22-27］。

不同的地表植被条件，通过影响下伏多年冻土与大气间的水热交换过程，进一步使得近地表范围内的能量平衡产生差异，这将显著影响多年冻土热状态及活动层水热过程特征［28-29］。不同的高寒生态系统下，多年冻土的地温、活动层厚度和冻融指数都不相同［30］；
地表植被盖度越小，多年冻土的冻结和融化起始时间越早、相变量越大、冻结深度的积分越小，且负值等温线和融化期间的正值等温线的最大侵入深度明显增加［31-33］。此外，在不同的地表植被条件下，近地表层的土壤含水率和融冻持续时长等也都不相同［34-35］。关注植被与多年冻土热状态相互作用的研究较多，且研究较为充分；
但在不同的地表植被条件下，有关工程设施，尤其是道路工程设施对多年冻土的热影响存在何种差异的问题，则相对研究得较少。

在多年冻土区，道路工程设施会显著改变地表和大气之间的热交换，从而使多年冻土的热状态受到影响［36-38］。随着多年冻土区社会经济活动的日益频繁，青藏高原工程走廊内各线性工程早已密布，且不断有新的工程开工建设，因此对道路工程设施与多年冻土之间的相互关系开展研究，无论是对道路工程设施的稳定运营，还是对保持多年冻土环境的相对稳定，都是必要且具有现实意义的。在这一研究领域，加强对不同地表条件差异的认识，有助于加深对多年冻土环境的了解，进而可以更好地指导工程实践，例如，对多年冻土区各工程之间的合理间距设定，以及即将建设的青藏高速公路双向路基间的合理距离设计都可以提供帮助。目前，关于青藏高原不同地表条件下多年冻土的相关研究，主要关注气候变化如何改变它们之间的相互作用［39-42］。关于道路工程对不同地表条件下多年冻土的热影响的大多数研究都基于室内模型或数值模拟［43-45］，通过现场监测数据进行的研究相对较少［46］。尽管模型或模拟方法具有其自身优势，但考虑到现实情况的复杂性，我们决定通过现场监测数据来开展研究。因此，在青藏高原腹地多年冻土区，我们沿青藏公路选择了两个不同地表条件的地点，布设监测场地，以研究不同地表条件下多年冻土对青藏公路热影响的响应差异。

在青藏高原腹地，我们沿青藏公路选择了两个具有不同地表条件的地点：斜水河和北麓河，作为监测场地。这两个监测场地的位置示意图如图1所示。两个监测场地地势都比较平坦，几乎没有坡度；
周边都建有距离较近的气象站，气候条件接近。斜水河监测场地的公路走向约为西南30°，北麓河监测场地的公路走向约为西南60°。如图1所示，青藏公路较南向的一侧称为左侧（L），较北向的一侧称为右侧（R）。

图1 监测场地位置Fig.1 Locations of monitoring sites

斜水河监测场地和北麓河监测场地均位于高寒干旱气候区，海拔4 500～4 700 m，年平均气温为-5～-3 ℃，年平均降水量约300 mm。两个监测场地所在区域多年冻土的主要类型是多冰和富冰冻土，多年冻土厚度约为50～80 m［47-49］。根据目视调查，斜水河监测场地的植被覆盖率低于20%。在北麓河监测场地，地表植被类型为高寒草甸，根据目视调查，其植被覆盖率约为70%～90%（图2），优势物种为矮嵩草和藏嵩草。通常情况下，北麓河监测场地的植被高度小于15 cm。根据钻孔数据，两个监测场地的地层和土壤含水量如图3 所示，土壤干密度如表1所示。

表1 监测场地土壤干密度Table 1 Soil bulk densities of monitoring sites

图2 监测场地地表特征Fig.2 Surface conditions of monitoring sites

图3 监测场地地层情况与土壤含水率Fig.3 Stratigraphy and soil water contents of monitoring sites

两个监测场地的传感器布置方式相同。钻孔以青藏公路为中心，在公路两侧对称分布，每侧有6个深度为15 m 的钻孔，钻孔位置及间距如图4 所示。土壤测温线缆布设于各钻孔中，每条土壤测温线缆长度为15 m，线缆编号如图4所示；
每条土壤测温线缆上的土壤温度探头的编号及深度如表2 所示。斜水河监测场地青藏公路的路基高度为3 m，北麓河监测场地青藏公路的路基高度为2 m。

表2 土壤温度探头布设深度Table 2 Depths for soil temperatures sensors

在每个监测场地，青藏公路左右两侧各有一个气温探头和两组土壤热通量板。气温探头编号为AL3 和AR3，架设高度距离地面2 m，分别位于L3和R3 土壤测温线缆上方。每组土壤热通量板由两块热通量板组成，安装在地表以下5 cm 和10 cm 深度。根据热通量板的监测数据，可以通过计算获得通过热通量板的总热量。各组热通量板的编号为HL1、HL5、HR1和HR5，分别与L1、L5、R1和R5的土壤测温线缆相邻。为便于叙述，将土壤测温线缆与左侧公路坡脚之间的距离定义为正值，如2 m、7 m，… ，67 m，右侧的相应距离定义为负值，如-2 m、-7 m，…，-67 m，如图4所示。考虑到左右坡脚处的距离都为0 m，为对其进行区分，将左侧坡脚处的距离标记为L0 m，左侧坡脚处的距离标记为R0 m。

图4 监测场地探头布设示意图Fig.4 Layout for sensors in monitoring sites

土壤温度探头由冻土工程国家重点实验室设计开发，其精度为±0.05 °C。气温探头为Campbell Scientific 公司（USA）的HMP45C-L11。热通量板为Hukseflux 公司（Netherlands）的HFP01sc。监测数据从2014 年9 月开始收集，并由数采仪（Campbell Scientific 公司，USA，CR3000）每30 分钟自动采集一次。日平均值是当天所有数据值的算术平均值。

3.1 气温和地温变化

表3显示了两个监测场地和五道梁国家气象站2015—2021 年的年平均气温（监测场地的气温值为左右两侧气温数据的平均值）。五道梁国家气象站位于两个监测场地中间，是距离最近的国家气象站。数据显示，北麓河监测场地的年平均气温接近五道梁国家气象站的年平均气温，高于斜水河监测场地的年平均气温。

表3 年平均气温Table 3 Mean annual air temperatures for monitoring sites and Wudaoliang National Weather Station

两个监测场地5 m 和15 m 深度处的月平均地温数据表明，在2015 年至2021 年期间，两个监测场地左侧距离公路最近的位置（L1）和最远的位置（L6）处的地温都在持续升高（图5）。这种地温升高的趋势，使得北麓河监测场地L1位置处和斜水河监测场地L6 位置处的5 m 深度，分别在2019 年（北麓河L1 5 m）（图5～6）和2016 年（斜水河L6 5 m）由持续冻结的多年冻土层转变为周期性融冻的活动层。图6显示了两个监测场地中左侧距离公路最近的位置（L1）处的活动层厚度的变化情况。数据表明，北麓河监测场地的活动层厚度小于斜水河监测场地，且从2015 年到2021 年，这两个监测场地的多年冻土上限变得越来越深。

图5 两个监测场地L1处2015—2021年5 m和15 m深度月平均地温变化情况Fig.5 Mean monthly soil temperature variations at 5 m and 10 m in L1 for both sites from 2015 to 2021

图6 两个监测场地L1处2015—2021年活动层厚度变化情况Fig.6 Active layer thickness variations in L1 for both sites from 2015 to 2021

3.2 公路的热影响范围

根据地温监测数据，两个监测场地的多年冻土活动层均大约在10 月中旬接近最大融化深度，且整个活动层基本处于融化状态；
在3 月中旬达到最大冻结状态，整个活动层基本完全冻结。基于这一趋势，选择2015 年和2021 年的10 月15 日和3 月15日来比较监测场地的土壤热状态，以检视公路的热影响范围。

图7显示了两个监测场地在2015 年和2021 年10月15日，多年冻土活动层融化深度接近最大融化深度时的地温状态，根据图中土壤温度等值线的分布，可以查看青藏公路的热影响范围。目前，在北麓河监测场地，青藏公路热影响的水平范围在左右两侧都约为15 m；
在斜水河监测场地，青藏公路两侧的热影响水平范围约为37 m。在这两个监测场地，随着时间的推移，青藏公路的水平热影响范围在逐渐扩大。图8 显示了2015 年和2021 年3 月15日两个监测场地地温状态，此时整个活动层基本完全冻结。从图中可以看到，靠近路基的位置地温最高，并朝远离路基的方向，地温向两侧依序降低。基于图中土壤温度的分布情况来看，在北麓河监测场地，青藏公路热影响的水平范围在左右两侧约为15 m；
在斜水河监测场地，青藏公路的水平热影响水平范围在左右两侧约为37 m，与图7 所反映的情况基本一致。

图7 2015年和2021年10月15日两个监测场地地温状态Fig.7 Soil temperature on October 15th of 2015 and 2021 for both sites

图8 2015年和2021年3月15日两个监测场地地温状态Fig.8 Soil temperature on March 15th of 2015 and 2021 for both sites

3.3 年平均地温和活动层厚度比较

年平均地温和活动层厚度是评价多年冻土状态的重要参数。在图9 和图10 中，对两个监测场地的上述两个参数进行了比较。图9 显示了2015年至2021 年两个监测场地公路两侧不同距离的年平均地温情况（年平均地温为15 m 深度处的地温）。数据表明，在北麓河监测场地中，距离公路越近，年平均地温越高；
同等位置左侧（较南向侧）的年平均地温高于右侧（较北向侧）。在斜水河监测场地，在公路左侧（较南向侧），距离公路越近，年平均地温越高，且这种相关性较为显著；
但在右侧（较北向侧），距离和年平均地温之间的相关性不如另一侧显著。数据同时表明，斜水河监测场地的年平均地温显著高于北麓河监测场地。此外，图9 同时也反映两个监测场地不同距离处的年平均地温都在持续升高。从2015 年至2021 年，北麓河监测场地两侧各距离处的年平均地温大约升高了0.10～0.14 ℃；
斜水河监测场地左侧各距离处的年平均地温大约升高了0.04～0.06 ℃，右侧各距离处的年平均地温大约升高了0.03～0.04 ℃。

图9 2015—2021年多年冻土年平均地温对比Fig.9 Comparison of mean annual ground temperatures from 2015 to 2021

在图10 中，对比了2015 年至2021 年两个监测场地多年冻土活动层的年度最大厚度。数据表明，在北麓河监测场地中，在两侧距离路基坡脚最近的两个位置（分别为L0 m、R0 m 和±2 m），活动层的年度最大厚度相对较大，最大值为5.95 m（2021 年，L1）。对于距离更远的其他位置，活动层厚度相对较小，最大值仅为3.01 m（2021 年，L6）。在斜水河监测场地左侧，距离路基坡脚最近的四个位置（分别为L0 m、2 m、7 m 和17 m），活动层厚度较大，最大厚度为8.32 m（2021 年，L1）。其余两个位置（分别为37 m 和67 m）的活动层厚度相对较小，最大厚度为5.84 m（2021 年，L6）。相比之下，在斜水河监测场地右侧，只有两个距离路基坡脚最近的位置（分别为R0 m 和-2 m），有较大的活动层厚度，最大值为6.63 m（2021 年，R1）。在其余的四个位置，活动层厚度都不超过6 m。当与青藏公路的距离相同时，斜水河监测场地的多年冻土活动层厚度明显大于北麓河监测场地。图10 同样也反映了活动层持续增厚的趋势。在北麓河监测场地，从2015 年至2021 年，左侧距离坡脚最近的L1 处，活动层厚度从4.00 m 增加到了5.95 m，为该场地在此期间的最大增幅。但在斜水河监测场地，活动层厚度的最大增幅并不出现在最靠近路基的位置，而是出现在右侧的R5 处；
该处活动厚度从2015 年的3.05 m 增加至2021年的4.69 m。

图10 2015—2021年各年活动层最大厚度对比Fig.10 Comparison of annual maximum active layer thicknesses from 2015 to 2021

图11显示了两个监测场地2015—2021 年多年冻土年平均地温总增幅及活动层厚度年平均增幅情况。数据表明，从2015 年至2021 年，北麓河监测场地年平均地温总增幅显著高于斜水河监测场地，且在两个监测场地，年平均地温的总增幅与距离青藏公路的远近程度都不相关。但活动层厚度年平均增幅却呈现不同的情况。在北麓河监测场地，随着与青藏公路距离的增加，活动层厚度的年平均增幅大致呈现递减的情况。然而在斜水河监测场地，活动层厚度的年平均增幅与距离青藏公路的远近程度并没有呈现出一个清晰的相关性关系。

图11 2015—2021年多年冻土年平均地温总增幅及活动层厚度年平均增幅对比Fig.11 Comparisons of the increase of the mean annual ground temperature and the active layer thicknesses from 2015 to 2021

3.4 年度总热量比较

多年冻土的形成是大气与地层之间热交换的结果。在北麓河与斜水河监测场地周边各有一个气象站。在北麓河监测场地，气象站位于其西面约600 m 处；
在斜水河监测场地，气象站位于其北面2 km 处。两处气象站都安装有辐射仪，可通过其测量的下行短波辐射近似估计周围地区的太阳辐射量。图12 显示了这两处气象站的年度下行短波辐射总量值。数据表明，北麓河与斜水河监测场地周边的年度下行短波辐射总量接近，这表明两个监测场地周边的太阳辐射量基本接近，其中北麓河监测场地略低。

图12 监测场地各年下行短波辐射总量对比Fig.12 Comparison of annual total downward solar radiations

在两个监测场地，热通量板每30分钟采集一次热通量数据。每平方米每30 分钟的总热量可通过公式（1）计算，其中q是热通量测量值（W·m-2），Q是每平方米每30 分钟的总热量（J）。通过热通量板的日总热量和年度总热量分别是时限内所有热量值的总和。

图13显示了这两个监测场地2015 年至2021 年各年通过所布设的热通量板的年度总热量数据。在北麓河监测场地，10 cm 深度处的年度总热量几乎都低于5 cm 深度处的年度总热量，仅有三处例外，但差值都很小。在斜水河监测场地，总的来说2015 年至2021 年的年度总热量在数值上小于北麓河监测场地。在该监测场地，靠近青藏公路的两组热通量板处（HL1 和HR1），10 cm 深度处的年度总热量常常高于5 cm 深度处的年度总热量；
甚至有时当5 cm 深度处的热传递方向是从地层到大气时，10 cm深度处的热传递仍然是从大气到地层方向。

图13 监测场地通过各热通量板的年度总热量对比（当年度总热量值为负时，代表通过热通量板的总热量方向为从土壤到大气）Fig.13 Comparison of total heat passing through each heat flux board（the direction of the heat transfer is from the soil to the atmosphere when the value of the heat is negative）

多年冻土的热状态必然是外部影响因素（如气温、降水、太阳辐射等）与其自身特性（如地表条件、土壤性质、含水量等）相互作用的结果。虽然这两处监测场地的太阳辐射量相当，但这两处场地的多年冻土状态因其他因素差异而有所不同。斜水河监测场地的气温低于北麓河监测场地、青藏公路的路基也高于北麓河监测场地。然而，即使有这两个有利条件，斜水河监测场地的多年冻土退化状态却比北麓河监测场地更显著，这可能是因为两处监测场地的自身特性存在差异。由于路基正下方没有地温探头，但从监测断面地温等值线的形状推测，青藏公路对路基正下方的多年冻土的热影响最大。在图7 和图8 中，我们可以看到地温等值线方向的变化：靠近青藏公路的位置比远离青藏公路的位置，其地温等值线的倾斜程度要大得多；
但青藏公路两侧地温等值线的倾斜程度并不对称，这应该是由路基的阴阳坡效应所造成的。在斜水河监测场地，这种不对称性要比北麓河监测场地明显，这无疑是因为两处监测场地不同的地表条件在发挥作用。

植被对土壤具有显著的影响，其通过遮蔽和蒸散作用影响土壤热状态［50-51］。在植被生长季节，一半以上的进入地层的能量会通过蒸散作用以潜热的形式消耗［52］，当水分充足的时候，这一比例可能超过80%［53］。在北麓河监测场地，地表的植被覆盖率和土壤含水量都高于斜水河监测场地，这可能导致更大比例的进入地层的能量以蒸散发作用由地层向大气传递。通过植被的蒸散和水分的蒸发过程，这种能量传递以潜热的形式发生。热通量板只能用于测量显热，但显热只是地气之间能量传递的一部分。这可能解释了即使北麓河监测场地的年总热量相对较高，但该监测场地多年冻土的年平均地温却显著低于斜水河监测场地，且多年冻土活动层厚度也相对较小。不仅如此，地表条件的差异，除了植被覆盖度情况不同之外，还有与之相关联的其他不同之处。在青藏高原，相关研究表明，植被、浅层土壤特征和浅层土壤含水量密切相关［54-56］。在植被覆盖度较高的北麓河监测场地，其浅层土壤为亚黏土，且含水量较高；
而在植被覆盖度较低的斜水河监测场地，其浅层土壤为砂砾石土，含水率相对较低（图3）。在这种差异背景下，较高的植被覆盖度、较细的浅层土壤颗粒和较高的浅层土壤含水量这些因素，将使得北麓河监测场地植被郁蔽作用、地面土壤蒸发作用和植被蒸腾作用都要比斜水河监测场地强烈，消耗更多的浅层土壤热量，有效降低地面温度，同时也减少向地层深部的热量传递，进一步使多年冻土地温得以保持在相对较低的状态。

图11中的数据表明，在两个监测场地，年平均地温的总增幅与距离青藏公路的远近程度都不相关，这或许是因为在较短时间尺度内，如仅仅几年内，青藏公路对于年平均地温（即15 m 深度处的地温）的直接影响并不显著，而年平均地温的影响因素仅仅取决于土壤性质等客观原因，由此导致在相同的监测场地内，年平均地温的总增幅都比较接近，但并不与距离青藏公路的远近程度相关。但也需注意到，在较长的时间尺度上，青藏公路对于年平均地温应该仍然是有影响的。如图9 所示，靠近公路的位置，其年平均地温要高于远离公路的位置。这或许是因为青藏公路开通运营以来，对多年冻土年平均地温数十年的影响结果。对于2015—2021 年间活动层厚度年平均增幅所呈现出来的情况（图11），或许可以通过气候变化和青藏公路共同的影响来进行解释。在北麓河监测场地，活动层厚度相对较小，除L1 位置的活动层厚度接近6 m 外，其他位置的活动层厚度都小于4 m（图10）。因此，对于相对较小的活动层厚度来说，气候变化和青藏公路都会对其产生影响，使得在北麓河监测场地，活动层厚度的年平均增幅随着与青藏公路距离的增加，而大致呈现递减的情况。但在斜水河监测场地，活动层厚度相对较厚，尤其是在L1、L2、L3 和L4，以及R1 和R2 位置，活动层厚度都超过6 m，某些位置甚至超过8 m。对于较大的活动层厚度而言，在较短时间尺度内，即使其距离青藏公路较近，公路所能发挥的热影响也相对有限，这或许是该场地活动层厚度的年均增幅与距离青藏公路的远近并无明确相关性的原因。但气候变化的影响依然对斜水河监测场地的活动层厚度变化产生了影响，尤其是对活动层厚度相对较小的位置，其影响相对更大，因此也使得在该监测场地，活动层厚度相对更小的位置，其活动层厚度年平均增幅更大。但同时也需注意，在较长的时间尺度上，青藏公路对于活动层厚度仍然是有影响的。如图10所示，在更靠近公路的位置，其活动层厚度要比远离公路的位置更大。这或许也是青藏公路数十年来对活动层厚度产生影响的结果。

此外，通常情况下在浅层土壤中，当热量从大气进入地层后，随着土壤对热量的吸收，热量的传递过程一般存在自上而下递减的趋势。当出现图11 中所示的相反的情况时，或许是由于热通量板数据测定的不确定性导致，而这一不确定性来源于两方面因素：一方面，斜水河监测场地地面仅有极低的植被覆盖度，几乎是裸地，这就导致5 cm 深度处的热通量板的测量数据极易受到外界因素的干扰，如风速等，可能会带走浅层土壤的热量，增加5 cm深度处的热量散失；
另一方面，由于青藏公路主要以大型货运车辆为主的巨大车流量，对于紧邻公路的HL1 和HR1 这两处热通量板安装位置而言，车辆频繁经过所导致的气流扰动，可能也会带走浅层土壤的热量，增大这两处热通量板、尤其是较浅的5 cm 深度处测量数据的不确定性。在上述两方面因素共同作用之下，使得5 cm 深度处土壤的年度吸收总热量与10 cm 深度处土壤的年度吸收总热量差距缩小，并出现5 cm 深度处土壤的年度总热量小于10 cm、甚至为负值的情况。而在北麓河监测场地，尽管也存在同样的风速和车流的扰动，但由于此处地表具有更高的植被覆盖度，从而减少了上述因素的扰动程度，保持了相对稳定的热量传递过程。

地表条件对多年冻土的稳定性以及多年冻土地区道路基础设施的安全运行至关重要。根据目前的研究结果，得出以下主要结论：

（1）青藏公路对其两侧的多年冻土有明显的热影响。距离公路越近，多年冻土退化越显著。与植被覆盖率更高的北麓河监测场地相比，植被覆盖率较差的斜水河监测场地多年冻土年平均地温较高，多年冻土活动层更厚。

（2）青藏公路对植被覆盖度高的多年冻土的热影响范围相对较小。在植被覆盖率较高的北麓河监测场地，青藏公路热影响的水平范围在左右两侧都约为15 m。在斜水河监测场地，地表植被覆盖率较差，青藏公路两侧的热影响水平范围约为37 m。

（3）多年冻土的热状态受地表条件的显著影响。随着地表条件，如植被覆盖度的不同，大气和地层间的热交换模式也存在差异。热交换过程中，地表植被必然起到了非常重要的影响，且在地表植被背后的土壤性质和土壤含水量等关联因素也发挥着重要的作用。

猜你喜欢 水河植被厚度 基于高分遥感影像的路域植被生物量计算中国交通信息化(2022年7期)2022-10-27大厚度填土场地勘察方法探讨建材发展导向(2022年6期)2022-04-18鉴水河之歌歌海(2021年2期)2021-06-22追踪盗猎者课外生活·趣知识(2021年2期)2021-05-24第一节 主要植被与自然环境 教学设计学校教育研究(2021年24期)2021-03-28与众不同的河水和湖泊作文周刊·小学二年级版(2020年36期)2020-10-15杭锦旗植被遥感分析农业与技术(2020年12期)2020-07-04脸皮究竟有多厚读者·校园版(2020年4期)2020-06-01诗要有温度，有厚度中华诗词(2019年8期)2020-01-06家乡的赧水河作文评点报·小学三、四年级(2019年48期)2019-12-30