地表反照率定义为地表反射太阳辐射与入射太阳辐射之比,是影响地表辐射能量收支和地-气相互作用的一个重要的地表参数(Dickinson,1995)。地表反照率在很大程度上决定了太阳辐射能量在地-气系统中的分配,影响大气环流和生物物理过程,对区域以至全球气候有着重要影响(Hu and Boos,2017;Planque et al.,2017),因而是不同时空尺度的气候模式(Haywood et al.,2016;Li et al.,2016;Levine and Boos,2017)、数值天气预报模式(Vahmani and Ban-Weiss,2016)以及陆面过程模式(谢志鹏等,2017)中重要的地表参数。此外,地表反照率还可以用于监测植被动态(Luke$\check{s }$et al.,2016;Wang et al.,2017)、沙漠面积变化(Liu et al.,2017),以及研究干旱对气候的影响(Evans et al.,2017)、火灾对生态系统的影响(Lyons et al.,2008;Saha et al.,2017),反演大气气溶胶状况及区域能见度变化(林楚勇等,2015;包云轩等,2018),城市化影响(储鹏等,2016)。

青藏高原平均海拔约4500 m,被称为“世界屋脊”、“第三极”,是世界海拔最高的高原。青藏高原大范围热动力作用不仅影响全球大气环流,而且对中国及亚洲地区的灾害性天气乃至全球气候变化有着重要影响(Ye,1981;Ye and Wu,1998)。研究并获得青藏高原地表反照率,对于促进青藏高原辐射能量收支以及地-气相互作用研究具有重要的科学意义。

青藏高原地形复杂、气候恶劣,卫星遥感是获得青藏高原高精度、高时空分辨率地表反照率的唯一有效途径(方宗义等,1996;徐兴奎和林朝晖,2002;陈爱军等,2009,2018)。目前,已有多种卫星遥感反演的地表反照率产品可以提供青藏高原的地表反照率,例如:美国MODIS(MODerate-resolution Imaging Spectroradiometer)地表反照率产品MCD43(Schaaf et al.,2010),我国北京师范大学制作的全球地表卫星遥感产品GLASS(Global LAnd Surface Satellites)地表反照率(Liang et al.,2013;Liu et al.,2013),以及欧洲航天局发布的全球地表反照率GlobAlbedo(Lewis et al.,2012)。

地表反照率作为重要的气候变量,全球气候观测系统GCOS(Global Climate Observing System)要求其绝对精度不低于0.05(https://public.wmo.int/en/programmes/global-climate-observing-system)。青藏高原卫星遥感反演地表反照率的精度一直备受关注,不少研究分析和评估了青藏高原MODIS地表反照率产品的精度(Wang,2004;余予等,2010;陈爱军等,2012a),也有不少研究分析和评估了青藏高原GLASS地表反照率产品的精度(Liu et al.,2009;Wang et al.,2010;齐文栋等,2014;王立钊等,2014;陈爱军等,2015)。由于观测站点和资料有限,只能分析和评估有限的卫星遥感反演结果的精度,少有研究深入分析影响卫星遥感反演结果精度的因素。反演质量是分析和评估大范围卫星遥感反演结果精度的一条重要途径(Jin,2003;陈爱军等,2016a,2016b)。此外,地表反照率卫星遥感反演结果的空间分布连续性,也是气候和陆面过程模拟与应用研究必须考虑的一个重要因素。

MODIS地表反照率是目前唯一能够近实时业务化提供青藏高原地表反照率的卫星遥感反演产品(Schaaf et al.,2010),GLASS地表反照率则是我国自主研发、时间序列最长的地表反照率卫星遥感反演产品(Liang et al.,2013;Liu et al.,2013),两者对于青藏高原地表辐射能量收支和地-气相互作用研究均有重要的应用价值。本文应用MODIS地表反照率和GLASS地表反照率反演质量数据,对比分析了两者在青藏高原大范围的空间分布连续性和高质量反演结果的比例,应用青藏高原CAMP/Tibet国际合作科学观测实验期间获得的高精度地表反照率观测结果评估它们的精度,通过人工目视解译MODIS地表反射率图像并结合MODIS积雪产品,分析和讨论了影响卫星遥感反演结果精度的因素。

1 研究数据

1.1 MODIS地表反照率数据

MODIS是搭载在美国地球观测系统EOS(Earth Observing System)上午星Terra和下午星Aqua的一个重要的探测器(https://modis.gsfc.nasa.gov/)。美国航空航天局NASA(National Aeronautics and Space Administration) 利用其观测数据开发了大量的大气、陆地、海洋和冰雪产品,MODIS地表反照率产品就是其中的陆地产品之一。

MODIS地表反照率采用半经验的线性RossThick-LiSparse-R核驱动的二向反射分布函数BRDF(Bidirectional Reflectance Distribution Function)模型,利用经过大气校正后的Terra和Aqua双星MODIS数据反演地表BRDF参数。反演时,根据反演周期内累积的高质量多角度晴空观测总数(记为N),采用全反演(the full inversion)或者当量反演(the magnitude inversion)获得反演结果,或者直接赋予填充值(无反演结果)(Lucht et al.,2000)。这里的高质量多角度晴空观测,是指卫星对地观测时,卫星-像元-太阳之间形成的观测几何角度尽可能地遍及整个半球空间。如果N≥7,首先尝试全反演;如果3≤N<7,首先尝试当量反演;如果N<3,直接赋予填充值。全反演失败,则尝试当量反演,当量反演失败则直接赋予填充值(Lucht et al.,2000)。所谓的当量反演,是用高质量多角度晴空观测数据修正先验的(a priori)BRDF参数,使它们更加符合地表的实际状态。然后,对BRDF积分获得地表黑空反照率BSA(Black-Sky Albedo)和白空反照率WSA(White-Sky Albedo)。

Aqua卫星发射前,MODIS地表反照率仅由Terra MODIS数据生成,每隔16 d反演一次,产品编号MOD43(Schaaf et al.,2002)。之后,则由16 d期间的Terra和Aqua双星MODIS数据每隔8 d反演一次,以便提高反演质量和效率,产品编号改为MCD43(Schaaf et al.,2010)。

MODIS地表反照率产品MCD43以500 m、1 km和0.05°气候网格等三种空间分辨率(分别编号A、B、C)提供地表BRDF模型参数、地表BRDF参数反演质量、地表反照率和BRDF校正后的天底反射率NBAR(Nadir BRDF Adjusted Reflectance)等四种产品(分别编号1、2、3和4)。空间分辨率为500 m和1 km的MODIS地表反照率产品采用正弦投影,空间分辨率为0.05°气候网格的产品则采用等经纬度投影。空间分辨率相同的MODIS地表反照率产品,具有完全一致的时空关系。

MODIS陆地正弦投影网格(http://landdb1.nascom.nasa.gov/developers/is_tiles/is_grid.html)编号为h24v05、h25v05、h26v05、h25v06和h26v06的5个区域正好覆盖青藏高原。本文所用数据为这5个区域2002—2004年空间分辨率为1 km的MCD43正弦投影数据,包括MCD43B3和MCD43B2。MCD43B3主要提供MODIS 1—7通道和可见光(0.3~0.7 μm)、近红外(0.7~5.0 μm)及短波(0.3~5.0 μm)等3个宽波段的地表黑空反照率BSA和白空反照率WSA(Schaaf et al.,2002)。MCD43B2主要提供地表BRDF参数反演质量,其中:质量标志0、1和255分别表示地表BRDF参数为“全反演结果”、“当量反演结果”和“填充值”(https://lpdaac.usgs.gov/products/modis_products_table/mcd43b2)。研究还采用了部分仅由Terra MODIS数据制作的MOD43B3和MOD43B2数据。

研究数据由美国地质调查局USGS(United States Geological Survey)地球资源观测与科学EROS(Earth Resource Observation and Science)数据中心提供,数据质量能够满足科学研究的需要(https://www.umb.edu/spectralmass/terra_aqua_modis)。

1.2 GLASS地表反照率

GLASS地表反照率是北京师范大学开发的5个全球地表卫星遥感产品之一(Liang et al.,2013;梁顺林等,2014)。该产品首先采用AB(Angular Bin)算法生成初级产品(逐日的地表反照率),然后采用基于统计学的时域滤波STF(Statistics-based Temporal Filter)算法对初级产品进行数据融合处理,包括数据的时间序列平滑,以及缺失数据的填补,获得多天合成的终级产品(梁顺林等,2014)。AB算法实际是一种直接反演算法:利用卫星观测的多波段方向性地表反射率,采用线性统计反演关系获得可见光、近红外和短波波段等三个宽波段的反照率(刘强等,2012)。具体而言,可以采用AB1算法或AB2算法,前者的多波段方向性地表反射率经过大气校正,后者则未经大气校正。

目前,对外公开发布的GLASS地表反照率产品的时间范围为1981—2013年。1981—1999年的GLASS地表反照率产品由AVHRR(Advanced Very High Resolution Radiometer)生成,仅采用AB1算法生成初级产品,终级产品GLASS02B05每隔8 d由33 d期间(产品标记日及其前后16 d)的初级产品采用STF算法融合而成。GLASS02B05采用等经纬度投影,空间分辨率为0.05°。2000—2013年的GLASS地表反照率产品由MODIS数据生成,每天由Terra和Aqua MODIS数据分别采用AB1算法和AB2算法获得4种初级产品(GLASS02A2X,X=1,2,3,4)。然后,每隔8 d采用STF算法融合17 d期间(产品标记日及其前后8 d)的4种初级产品,生成终级产品GLASS02A06。GLASS02A06采用与MCD43相同的正弦投影方式,空间分辨率为1 km(Liang et al.,2013;Liu et al.,2013;梁顺林等,2014)。

本文研究所用数据为2002—2004年的GLASS02A06(第1版),数据区域与研究所用MCD43数据对应。GLASS02A06主要提供短波波段(0.3~3.0 μm)的黑空反照率BSA和白空反照率WSA,同时提供反照率的总体质量及其他数据信息(Liu et al.,2013;王立钊等,2014)。GLASS02A06采用质量标记0、1、2和3分别表示产品的总体质量“好”、“可接受”、“较差”和“不可用”。研究所用数据由北京师范大学全球陆表特征参量产品免费提供(http://glassproduct.bnu.edu.cn)。

1.3 地面观测数据

本文所用地面观测数据为BJ、NPAM、TGL和XDT等四个站点在青藏高原CAMP/Tibet国际合作科学观测实验期间获得的高精度地表下行和上行太阳短波总辐射观测资料。这四个站点位于MODIS陆地正弦投影网格编号为h25v05的区域,处于青藏高原中部(表1)。观测站点所处位置都比较开阔,周围地势也比较平坦,下垫面能够代表藏北高原大范围的地表特征,具有较好的代表性。所有仪器在试验前都经过严格的标定和校准,试验期间也得到较好的维护和保养,因而观测资料的质量较好,具有较好的准确性和可比性。研究所用数据由中国科学院青藏高原研究所提供,其中:BJ站和NPAM站的观测数据为2002年资料,时间分辨率为1 h;XDT站和TGL站的观测数据为2006年,时间分辨率为30 min。

1.4 其他数据

研究还采用了MODIS陆地正弦投影网格编号为h25v05区域的另外两类数据:一类是空间分辨率为500 m经过大气校正的MODIS地表反射率产品,包括Terra MODIS数据生成的MOD09GA和Aqua MODIS数据生成的MYD09GA(https://modis.gsfc.nasa.gov/data/dataprod/mod09.php);另一类是空间分辨率为1 km的MODIS积雪产品,包括Terra MODIS数据制作的MOD10A1和Aqua MODIS数据制作的MYD10A1(https://modis.gsfc.nasa.gov/data/dataprod/mod10.php)。对于大气校正后的MODIS地表反射率数据,研究主要利用它们提供的MODIS 1—7通道经过大气校正的地表反射率数据及其质量信息。对于MODIS积雪产品,研究主要利用它们提供的MODIS云/雪判识结果及其质量信息。这些数据均由美国地质调查局USGS地球资源观测与科学EROS数据中心提供,数据质量能够满足科学研究的需要。

2 研究方法

2.1 地面辐射观测数据处理

首先,参考陈爱军等(2015)的处理方法,以下行太阳短波总辐射不低于600 W/m²为标准筛选出“晴空”辐射观测资料,以上行总辐射与下行总辐射之比计算时间分辨率为1 h或30 min的“晴空”地表反照率,然后采用算术平均计算日平均地表反照率。

然后,采用与研究所用GLASS地表反照率产品GLASS02A06完全一致的时域窗口,应用算术平均计算17 d平均地表反照率。由于积雪覆盖对地表反照率的影响较大,研究参考陈爱军等(2015)的处理方法,以日平均地表反照率0.4为标准将17 d分为“积雪日”(不低于0.4)和“无雪日”(低于0.4),根据17 d期间“积雪日”和“无雪日”的多寡,选择对应条件的日平均地表反照率计算17 d平均地表反照率。

2.2 卫星遥感反演产品处理

无论是MODIS地表反照率,还是GLASS地表反照率产品,它们提供的只是黑空反照率和白空反照率(Schaaf et al.,2002;刘强等,2012),也就是太阳辐射完全直射和完全漫射条件下的地表反照率,并不是真实的地表反照率。真实的地表反照率需要根据实际天空漫射光比例对两者进行线性加权平均(Pinty et al.,2005),即:

其中 _WSA和α_BSA分别是白空反照率和黑空反照率;r为实际天空漫射光比例,可以根据如下经验公式获得(Stokes and Schwartz,1994):

其中θ是正午时的太阳天顶角。

2.3 数据比较方法

采用“点对点”的比较方式对比分析地表反照率地面观测结果与卫星遥感反演结果之间的差异,即:将站点对应像元的卫星遥感反演结果与站点的地面观测进行比较。由于GLASS地表反照率GLASS02A06与MODIS地表反照率MCD43B3和MOD43B3的时域合成窗口存在差异,本文将时间标记为第9天的GLASS02A06与时间标记为第1天的MCD43B3或MOD43B3进行比较(图1),并依次类推让两种卫星遥感反演产品在时域上尽可能对应一致,同时忽略两者时域合成周期相差1 d引起的差异。另外,本文忽略地面观测结果、GLASS地表反照率短波波段、MODIS地表反照率短波波段的谱宽差异。

3 反演结果对比

3.1 反演效率

反演效率,即有效反演结果的比例,可以反映地表反照率反演结果的空间分布连续性。统计分析结果表明:2002—2004年青藏高原GLASS地表反照率的反演效率几乎达到100%,最低也超过99.9%(图2a);MODIS地表反照率的反演效率则相对偏低,最高约98%,一般在80%以上,最低时甚至不足60%(图2b)。这表明GLASS地表反照率的空间分布连续性更好,更加有利于气候和陆面过程模拟与应用研究,而MODIS地表反照率的空间分布连续性仍有较大的改进空间。分析其中的原因,与两种卫星遥感反演产品的算法有着密切关系。根据GLASS地表反照率算法,只要有晴空观测数据即可由AB算法获得日产品,而且终级产品由17 d期间的日产品采用STF算法融合而成,因而反演效率较高。然而,MODIS地表反照率反演算法对高质量多角度晴空观测总数有着较高的要求,即使当量反演

也要求至少3个,全反演则要求至少7个(Schaaf et al.,2010)。

3.2 高质量反演结果的比例

反演质量是评价大范围卫星遥感反演结果精度的一个重要指标(Jin,2003)。地表反照率高质量反演结果的精度更高,而低质量反演结果的精度具有较大的不确定性(陈爱军等,2012a,2012b,2015)。分析结果表明:2002—2004年,青藏高原MODIS地表反照率高质量反演结果的比例一般介于40%和80%,最高约90%(2003年第281天开始的8 d期间),最低不足1%(2002年第73天开始的8 d期间);GLASS地表反照率高质量反演结果的比例大多超过80%,最高约98%(2002年第313天开始的8 d期间),最低约56%(2002年第1天开始的8 d期间)(图3)。另外,不同反演周期的GLASS地表反照率和MODIS地表反照率高质量反演结果的比例存在较大的差异。总体而言,青藏高原GLASS地表反照率高质量反演结果的比例要高于MODIS

地表反照率高质量反演结果的比例,且相对稳定。其中的原因,主要与各反演周期的晴空观测总数有关:即使晴空观测总数不足,GLASS地表反照率反演算法也可以通过融合多天的反演结果最终获得质量较高的反演结果,而MODIS地表反照率反演算法不仅对晴空观测总数有要求,而且只有高质量多角度晴空观测总数不少于7个时才有可能通过全反演获得高质量反演结果,必然影响其高质量反演结果的比例。

3.3 精度分析

对比分析TGL、XDT、NPAM和BJ站的地面观测结果MODIS地表反照率和GLASS地表反照率,绝大多数时段内它们都有较好的一致性(图4)。尽管某些时段存在一定的差异,但是仍然能够保持较为一致的变化趋势,较好地反映地表因为经历降雪、融雪过程而引起的地表反照率异常变化过程,例如:在TGL站第65、73、81、89、97天开始的时段,以及第289、297、305、313、321、329、337天开始的时段(图4a)。

MODIS地表反照率或GLASS地表反照率在某些时段也与地面观测结果存在较大的差异,大致可以分为以下几种情形:第一种是MODIS地表反照率与地面观测结果基本一致或相差较小,而GLASS地表反照率明显高于地面观测结果,例如:TGL站第41天开始的8 d期间(图4a)和XDT站第41天开始的8 d期间(图4b);第二种是GLASS地表反照率与地面观测结果基本一致或者相差较小,而MODIS地表反照率明显低于地面观测结果,例如:XDT站2006年第297天开始的8 d期间(图4b);第三种是MODIS地表反照率和GLASS地表反照率均显著高于地面观测结果,例如:NPAM站2002年第25天开始的8 d期间(图4c);第四种是MODIS地表反照率和GLASS地表反照率均显著低于地面观测结果,例如:BJ站2002年第17、25、33、57、65天开始的8 d期间(图4d)。深入分析MODIS地表反照率或GLASS地表反照率与地面观测结果存在较大差异的原因,结果表明:云覆盖是导致第一种情形的主要原因。以TGL站2006年第41天开始的8 d期间为例,第33—49天期间,人工目视解译结果表明:TGL站及其周边区域除8次无雪之外,其余9次均为云覆盖。因此,导致GLASS地表反照率明显高于地面观测结果的唯一可能原因是将云覆盖误判为积雪,因而反演结果高达0.47。MODIS地表反照率虽然比GLASS地表反照率的输入数据少1 d,但是通过MODIS积雪产品提供的云雪判识结果对输入数据的质量进行了严格控制,只能以地表无雪状态的观测数据反演地表BRDF参数和反照率,反演结果与地面观测结果非常接近。这说明云雪判识是卫星遥感反演地表反照率必不可少的环节,有助于提高反演结果的精度,尤其是GLASS地表反照率。

第二种情形的原因,与地表处于积雪状态和MODIS地表反照率当量反演结果的不确定性有关。根据MODIS积雪产品提供的云雪判识结果,Terra和Aqua卫星于2006年第289—304天期间过境观测时,XDT站5次无雪、11次积雪、1次云覆盖。根据MODIS地表反照率反演算法,应该采用晴空积雪状态下的观测数据反演。理想情况下,完全可能通过全反演获得反演结果,然而实际反演结果不仅是当量反演结果,而且低至0.23(图4b),与地面观

测结果的绝对偏差高达0.36。这既说明MODIS地表反照率反演算法在积雪状态下存在一定的局限性,也说明MODIS地表反照率当量反演结果存在较大的不确定性。相反,GLASS地表反照率反演结果为0.56,与地面观测结果的绝对偏差只有0.03,显示出该反演算法的优势,尤其地表处于积雪状态。

对于第三种情形,较为合理的解释是地面观测站点无雪,而站点对应的像元存在积雪。根据MODIS积雪产品提供的云雪判识结果,2002年第17—33天期间NPAM站有9次积雪、3次无雪状态,另有5次云覆盖状态。人工目视解译结果也显示该站有10次积雪、3次无雪状态,其余4次云覆盖状态。在此期间,GLASS地表反照率和MODIS地表反照率分别高达0.54和0.73(图4c),都是较为合理的结果。然而,这17天期间的地表反照率地面观测结果表明该站并无积雪,只有第21天的地表反照率略高,为0.34,其余大多约0.22。因此,只有NAPM站对应的像元因为积雪覆盖引起卫星观测的地表反射率偏高,导致地表反照率反演结果偏高。比较而言,这期间GLASS地表反照率应该比MODIS当量反演结果的可靠性更高。另外,这期间MODIS地表反照率仍然为当量反演结果,也再次说明MODIS地表反照率反演算法在积雪状态下的局限性,以及当量反演结果的不确定性。

第四种情形则与站点对应像元存在局地积雪有着密切关系。以BJ站2002年第57天开始的8 d期间为例,人工目视解译结果表明:2002年第49—66天期间,BJ站及周边区域出现6次积雪,而且BJ站多次处于薄雪区或积雪边缘区,另有7次云覆盖和4次无雪。BJ站地面观测的日平均地表反照率普遍较高,最高约0.80,最低也有0.41,只有2 d的日平均地表反照率约为0.30。显然,BJ站出现了积雪和消融过程,因而17 d平均的地表反照率超过0.60(图4d)。然而,由于站点对应的像元经常处于薄雪或积雪边缘区,积雪覆盖比例较小,对卫星观测的地表反射率的影响较小,导致地表反照率卫星遥感反演结果偏低。其他几个类似的情形,也与2002年第57天开始的8 d期间的情况基本类似。

导致第三、四种情形的实质都是局地积雪的影响,根本原因在于卫星遥感反演结果与地面观测结果比较时存在“点”对“面”的差异。采用对应时次的高空间分辨率、高精度卫星遥感反演结果,无疑是消除这种差异、实现“面”对“面”比较的一条有效途径。然而,目前仍然是定量评估卫星遥感反演结果精度的难点之一。

综合TGL、XDT、NPAM和BJ等4个站点的地面观测结果,统计分析它们与卫星遥感反演结果的均方根误差,结果表明:MODIS地表反照率和GLASS地表反照率与地面观测结果的均方根误差分别高达0.0861和0.0771;即使只考虑高质量反演结果,也仍然高达0.0787和0.0775。尽管MODIS地表反照率的均方根误差略小,仍然无法达到全球气候观测系统GCOS对地表反照率的绝对精度不低于0.05的要求。不过,结合人工分析剔除错误的反演结果,MODIS地表反照率和GLASS地表反照率的均方根误差均减小至0.05以内,分别为0.0280和0.0329;如果只考虑高质量反演结果,则进一步减小至0.0206和0.0330。

4 结论

本文采用MODIS和GLASS地表反照率反演质量数据,分析了两种卫星遥感反演结果在青藏高原大范围的空间分布连续性和高质量反演结果的比例,应用青藏高原CAMP/Tibet试验期间BJ、NPAM、TGL和XDT等4个站点的高精度地面观测结果评估了两种遥感反演结果的精度,通过人工目视解译MODIS地表反射率图像并结合MODIS积雪产品提供的云雪判识结果,分析和讨论了影响卫星遥感反演结果精度的因素,结论如下:

1)GLASS地表反照率和MODIS地表反照率的反演效率和高质量反演结果的比例在不同的反演周期存在较大的差异:GLASS地表反照率的反演效率超过99.9%,高质量反演结果的比例大多超过80%,最高约98%;MODIS地表反照率的反演效率一般超过80%,但最低不足60%,高质量反演结果的比例一般介于40%和80%,最低不足1%。相对而言,GLASS地表反照率具有比MODIS地表反照率更好的空间分布连续性和更高的反演质量。

2)MODIS地表反照率和GLASS地表反照率在绝大多数时段内都能与地面观测结果保持较好的一致性,能够准确地反映地表积雪及消融过程引起地表反照率的异常变化过程,但是部分时段也与地面观测结果存在较大的差异。

3)局地积雪是导致地表反照率卫星遥感反演结果与地面观测结果产生差异的重要原因,云雪判识有助于提高GLASS地表反照率的反演精度。GLASS地表反照率反演算法在积雪条件下比MODIS地表反照率反演算法更具优势,反演结果的精度更高。

4)青藏高原地表反照率卫星遥感反演结果的精度存在不确定性:即使MODIS地表反照率和GLASS地表反照率高质量反演结果,如果未经人工分析,与地面观测结果的均方根误差高达0.0787和0.0775;经过人工分析并剔除错误的反演结果后,与地面观测结果的均方根误差可以减小至0.0280和0.0329,能够满足全球气候观测系统GCOS要求地表反照率的绝对精度不低于0.05的要求。

本文研究结果有助于进一步优化和改进MODIS地表反照率和GLASS地表反照率反演算法,提高青藏高原地表反照率卫星遥感反演结果的精度、反演质量和空间分布连续性,促进地表反照率卫星遥感反演产品在青藏高原辐射能量收支和地-气相互作用研究中的应用。

致谢

参考文献