气象学报  2013, Vol. 71 Issue (4): 769-782   PDF    
http://dx.doi.org/10.11676/qxxb2013.056
中国气象学会主办。
0

文章信息

赵亮, 朱玉祥, 杨弘, 胡娅敏, 夏明, 郭建平. 2013.
ZHAO Lung, ZHU Yuxiang, YANG Hong, HU Yamin, XIA Ming, GUO Jianping. 2013.
一种基于卫星遥感与地面测站数据融合技术的雪深动态反演方法
A dynamic approach to retrieving snow depth based on the technology of integrating satellite remote sensing and in situ data
气象学报, 71(4): 769-782
Acta Meteorologica Sinica, 71(4): 769-782.
http://dx.doi.org/10.11676/qxxb2013.056

文章历史

收稿日期:2012-10-24
改回日期:,2013-04-02
一种基于卫星遥感与地面测站数据融合技术的雪深动态反演方法
赵亮1, 朱玉祥2 , 杨弘1, 胡娅敏3, 夏明1, 郭建平4    
1. 61741部队, 北京, 100094;
2. 中国气象局气象干部培训学院, 北京, 100081;
3. 广东省气候中心, 广州, 510080;
4. 中国气象科学研究院大气成分研究所, 北京, 100081
摘要:提出了一种新的基于被动微波遥感和地面测站数据融合技术的雪深动态反演方法.这种新方法不再依赖单一的地面测站数据或卫星遥感数据,而是利用它们联合建立雪盖可信度指数,共同确定雪盖分布;然后在此基础上采用时空距离权重法设定反演系数动态参数化方案,反演雪深.这种雪深反演方法具有以下特点:针对不同时空条件下反演系数的动态差异问题,提出利用实时测站观测雪深,灵活调整雪深反演系数的解决方案,使反演系数具备随时空动态调整的能力,这是与静态反演方法最大的区别;充分利用了被动微波遥感数据时空连续性好的优势,能够在测站稀少的西部高山地区反演出空间分辨率相对较高的雪深数据,这是地面观测无法做到的.初步检验结果显示,该方法较明显地提高了中国西部高原地区和东部雪盖南缘区的反演精度,并克服了原有融合方法在中国西部雪盖面积偏小的问题,有效避免了静态反演方法在高山地区严重高估而平原地区低估雪深的问题,实现了被动微波遥感和地面观测数据的有效融合,扩大了雪深监测的有效范围
关键词积雪深度     雪盖     反演     被动微波遥感     青藏高原    
A dynamic approach to retrieving snow depth based on the technology of integrating satellite remote sensing and in situ data
ZHAO Lung1, ZHU Yuxiang2 , YANG Hong1, HU Yamin3, XIA Ming1, GUO Jianping4    
1. 61741 Troops of PLA, Beijing 100094, China;
2. CMA Training Center, Beijing 100081, China;
3. Guangdong Climate Center, Guungzhou 510080, China;
4. Institute of Atmospher Composition, Chinese Academy of Meteorological Sciences, Beijing 100081, China
Abstract:A new dynamic approach to retrieving snow depth is developed based on integration of passive microwave remote sensing and in situ data.First,the snow-cover confidence index is established by use of both the passive microwave remote sensing and in situ data to identify the snow cover; second,a new dynamic parameterized scheme (the distance weighted method) is developed based on the index.The characteristics of the snow-depth retrieval approach arc as follows:on the one hand, for the difference issue of retrieval coefficients in different spatial-temporal circumstances,a solution is proposed that retrieval coefficients arc able to be adjusted according to real-time observed snow depth,which is the biggest difference from static retricval approaches;on the other hand,the advantage of spatial-temporal continuity of the passive microwave remote sensing da-to has been fully exploited,being able to retrieve the snow depth with relative high resolution and precision in the west China where few stations arc located.The results show that the approach implements the efficient integration of passive microwave remote sensing and observed data, exerts the advantages of different source data, improves obviously the retrieval precision in the western part of China and the south marginal regions of snow cover in eastern China, and solves the question in the old integrating approach that the area of snow cover was always relatively smaller than normal in the west China, amplifying the detectable coverage area of snow depth.In contrast to the static retrieval approach,the dynamic retrieval approach avoids efficiently the question that snow depth was overestimated in mountain regions and underestimated in plain regions.
Key words: Snow depth     Snow cover     Retrieving     Passive microwavc remote sensing     the Tibetan Plateau    
1 引 言

雪可以通过改变地表反照率而影响地面对太阳辐射的吸收,可以因低的热导性而在地表和大气间形成绝热层,可以在融雪时产生吸热效应,改变大气增温率,这些效应综合起来会影响天气气候系统,最终会影响降水、气温等气候要素(罗勇,1995李培基,1996Cohen et al,1999陈乾金等,2000刘玉洁等,2003Yasunari,2007张佳华等,2008)。因此,人类估算雪盖和雪深的能力直接影响人类监测气候变化和检测气候模式的能力(Kelly et al,2003)。同时,对积雪的准确监测,可以对气候的监测和预测有重要帮助,例如,Cohen等(2003)研究显示,在雪盖期,用雪盖指数(基于夏秋雪盖)预测北美东部陆面温度比北极涛动指数更好。而且,积雪也可以成为水源,准确地估算雪水容量有助于提高淡水供给的管理能力。

然而,积雪是最动态和活跃的水文学变量之一,北半球的积雪面积可以从4.7×107 km2(1和2月)改变为4×106 km2(8月)(Kelly,2009)。传统的台站观测显然不能满足要求,准确监测积雪需要先进的现代技术。遥感方法监测积雪已有40多年历史(Frei et al,1999),其中,利用可见光和红外传感器遥感雪盖面积的技术正走向成熟,而微波遥感雪容量(雪深或雪水当量)处于发展阶段(Kongoli et al,2006Foster et al,2011)。被动微波遥感在特定频率上有能力穿过云层,几乎可以全天候工作,而且,可以提供雪水当量、雪深等信息。所以,目前,在反演雪深方面,被动微波遥感正发挥着重要作用,并已经取得一些重要进展,例如,应用多通道微波扫描辐射计(SMMR)、专用传感器微波成像仪(SSM/I)、先进的微波探测器(AMSU)和先进的微波扫描辐射计(AMSR-E)等探测数据反演雪深(Chang et al,1987Chen et al,2001Kelly et al,2003车涛等,2004陈爱军等,2005延昊等,2008)。

但是,被动微波遥感反演雪深还有很多重要的问题尚待解决。一方面,在大多数(尤其平原)地区,与可见光和红外雪盖图对比,被动微波反演经常趋向于低估雪盖面积(Chang et al,1987),雪水当量的误差也往往较大(Pulliainen et al,2001)。这可能与被动微波不能探测浅雪和湿雪有关。另一方面,在干冷的高原地区,却又往往明显高估雪深和雪面积,青藏高原和亚洲中部地区是主要的高估区。Foster等(2011)认为,这是因为尽管这里冬季气温很低,但由于大气太干燥而经常没有雪,而被动微波反演算法误认为这种冷的表面是雪盖。青藏高原地表状况复杂,冻土、寒漠等都会产生类似积雪层的散射特征(车涛等,2004李晓静等,2007李新等,2007)。因此,对于较大区域(如整个中国地区)的雪深反演,如果采用固定参数的静态反演算法,会造成较大误差。

鉴于静态被动微波反演雪深在局地和短时间尺度上存在的问题,近些年有学者尝试用动态参数化方案反演雪深(Kelly et al,2003Foster et al,2005),考虑了雪粒径的动态变化和森林覆盖度,在大多数地区和大部分时段效果较好,但在高山和高原地区效果仍然欠佳。目前解决这一问题的方法之一就是将微波遥感与其他遥感(如可见光遥感)数据结合联合反演积雪,例如,Armstrong(2001)将基于不同微波传感器(如SSM/I、AMSR-E或AMSU等)的微波雪水当量产品与基于其他辐射器(如中分辨率成像光谱仪(MODIS)等)的雪盖面积产品进行融合;美国国家海洋大气局(NOAA)国家环境卫星数据信息系统(NESDIS)使用交互式多传感器雪冰制图系统(IMS),在业务上提供北半球雪冰图作为环境预测模式的输入量(Helfrich et al,2007);杨秀春等(2009)将光学MODIS数据与被动微波AMSR-E数据协同用于监测中国草原积雪。特别值得关注的是Foster等(2005)提出了一种被动微波反演雪水当量的新算法,2007年美国空军和NASA采用该方法联合开发了融合可见光、被动微波和散射计的雪产品(Foster et al,2007),Hall等(2007)将这一产品称为AFWA-NASA(简称ANSA)融合雪盖产品,并做了初步评估。Akyurek等(2010)Hall等(2012)进一步评估了这一产品在东土耳其山区和五大湖区下游的反演准确度。

现在,可供利用的卫星遥感资料越来越多,将多种来源的数据进行有效融合,将是未来积雪反演的一个重要方向。基于天基的卫星遥感和基于地基的测站观测各有优势,结合起来使用能较好地反映积雪信息(张佳华等,2008),因此,发展空间遥感和地面观测于一体的数据融合技术,有助于提高雪深反演的准确性。就中国这样一个地形复杂的大陆尺度情况而言,虽然利用被动微波反演雪深有诸多优点,但是,在高山地区(如青藏高原地区),多种单一的被动微波反演产品,效果都不够理想。测站观测雪深虽然较准确,但由于地形复杂,测站稀少,致使其代表性差。所以,在反演雪深时,在测站稀少的地区应当以遥感数据为主,但也应适当借鉴测站数据以调整反演参数。另外,在测站相对较密的地区(中东部地区),被动微波遥感反演的积雪信息在雪盖南侧边缘区(大陆雪线附近)往往也不准确(赵亮等,2010Foster et al,2011),但这里的观测数据代表性强,所以,应充分重视和利用地面观测数据。因此,根据测站分布的密度,来调整测站数据和遥感数据的利用程度,即调整其在雪深反演中所占的权重,发展空间遥感和地面观测于一体的数据融合技术,将有可能提高雪深反演精度。目前,国际上雪产品的多源数据融合大多基于卫星遥感数据,而未将地面观测站数据加以合理融合,其不可回避的问题是在有云时段的高山地区,反演结果较差,例如,ANSA算法即存在这样的问题,而如果结合地面测站数据,则必将有利于雪深反演精度的提高。赵亮等(2010)曾对遥感-测站结合的雪深动态反演方法进行了初步探讨,但其方法在确定雪盖时没有利用遥感数据,所得雪盖不够准确,测站稀少的地区常无法获得反演系数。本文针对这些问题,试图在反演算法阶段将台站观测数据与遥感数据进行融合,使雪盖的确定更为科学准确,并在此基础上提出新的基于距离权重法的反演系数动态参数化方案以改进反演算法,提高反演准确度。2 数据和资料

亮温数据来自搭载于美国国防气象卫星计划(DMSP)卫星上的专用传感器微波成像仪(SSM/I)(Armstrong et al,2003),轨道高度为830 km,由美国国家雪冰数据中心(NSIDC)获得,时间跨度为1987年7月—2009年11月。该仪器接收到19、37和85 GHz的双极化和22 GHz 的垂直极化亮温数据。数据格式和投影方式是等面积可调整地球网格(EASE-Grids),空间分辨率为25 km,此外,85 GHz还有12.5 km的分辨率。本文进行反演前利用美国国家雪冰数据中心提供的坐标转换工具将EASE-Grid投影亮温数据转换成等经纬度投影亮温数据(0.225°分辨率),便于插值运算。观测雪深数据来自中国气象局整编的中国753个气象台站积雪日值数据集,时间跨度与专用传感器微波成像仪(SSM/I)亮温数据一致。3 存在的问题和解决方案

不同地区不同季节的积雪,其物理性质往往不同,即这种差异是时空动态变化的。所以,用固定的反演系数是不合理的。对于某一时次某一局地反演系数的确定,有重要参考意义的是该时次该格点的遥感信息和附近测站的观测信息,在技术实现过程中应使其充分发挥作用。赵亮等(2010)利用不同来源的数据,提出动态变化的反演系数代替固定反演系数的方案,很大程度上克服了积雪物理特征的区域性和季节性差异带来的反演误差,初步实现了两种不同来源资料的优势互补。其相对静态方法有诸多优势,但存在的问题是:确定雪盖时微波遥感信息没有被充分利用,造成信息浪费,也导致反演雪盖存在较大误差,即认为插值后的正的反演系数区域就是雪盖,造成了西部测站稀少地区雪盖面积偏小;另外,由于剔除了反演系数为负值的格点(常出现在东部湿度较大的长江流域和沿海地区,观测有雪但19和37 GHz水平极化亮温差小于0的区域),造成该格点无法反演,使该地区雪盖面积偏小。我们认为,前者可以通过微波遥感反演雪盖(利用积雪分类树方法得到)的方法克服(以后还可结合可见光遥感雪盖得到更准确的雪盖分布),后者则可通过保留观测有雪但亮温差小于0的格点(即允许反演系数小于0)来克服。

在解决以上两个问题的基础上,提出改进的动态反演雪深方法,其基本思想是:获取“遥感”雪盖和“观测”雪盖,在此基础上获得“融合”雪盖,建立时空动态的雪盖可信度指数,并针对不同指数值,提出“时空距离权重法”计算动态反演系数,从而得到动态反演雪深。这里,需注意两种特殊情况,在测站稀疏的地区或夏季,常出现观测无雪,但遥感反演有雪的情况,此时应该以遥感结果为主要参考依据,即认为有雪,但其反演系数确定可参考该地区一定期限内有雪时的平均反演系数。另一种情况是,当观测有雪时,不管遥感反演是否有雪,这样的区域都应予以保留,即认为该区域(格点)有雪,确定其反演系数时,可依据该格点与最近有雪站点的距离作为权重,结合该时次用插值方法得到的反演系数及其有雪时的均值来共同决定反演系数。

详细技术路线如下(图 1):

图 1 遥感-测站相结合的动态雪深反演方法技术路线流程
(设定测站数据影响半径r0=403.8 km)
Fig. 1 Technical flowchart for the temporal-spatial dynamic approach to retrieving snow depth based on integration of remote sensing and observed data
(the influencing radius r0=403.8 km)

第1步,数据预处理

(1)获取初步的观测格点雪深:利用气象学上应用较普遍的Cressman插值方法(Cressman,1959)将测站雪深插到0.225°网格点上,并去除小于0的虚假格点雪深,得到初步的观测格点雪深D0(x,y,t)。冯锦明等(2004)曾比较了多种空间插值方法,发现Cressman插值结果与原始台站观测最为接近,但其在某些区域会出现空值或负值(冯锦明等,2004),对于雪深而言,负值显然是不合理的,所以,这里需将负值剔除。

(2)获取初步的格点反演系数及其年均值:首先,需利用NSIDC提供的坐标转换工具将EASE-Grid投影的SSM/I亮温数据(25 km分辨率)转换成等经纬度投影亮温数据(0.225°分辨率);其次,获取中国范围内各通道逐日(1987年7月—2009年11月)亮温数据,缺测的点用之前时次升轨和降轨数据补齐(赵亮等,2010),并计算19和37 GHz水平极化亮温差ΔTb(x,y,t);最后,利用式(1)

获取初步的格点反演系数A0(x,y,t),并进一步得到其年内有雪时的均值A0m(x,y)。这里,需要注意的是,保留ΔTb(x,y,t)<0即A0(x,y,t)<0的格点,因为实际上亮温差小于0并不意味着一定无雪(尤其在东部地区),初步的格点反演系数暂时先保留这样的格点,后面再使用一定的方法从中筛选出有雪的点,有利于提高东部地区的反演精度。

(3)获取逐日“遥感”雪盖:利用积雪分类树方法(Grody et al,1996),获取未考虑测站实测数据的逐日纯“遥感”雪盖。

(4)获取逐日“观测”雪盖:根据某日距离某格点最近的站观测雪深是否大于0,来判断该日该格点观测是否有雪,这样,获得未考虑遥感数据的逐日纯“观测”雪盖。这样得到的“观测”雪盖在中东部地区比用空间插值方法得到的雪盖更为合理,而后者会出现许多虚假的有雪区域。而在西部得到的“观测”雪盖往往偏小,但可以根据“遥感”雪盖纠正。

(5)计算每日各格点离最近有雪站的距离t(x,y,t)。

第2步,计算雪盖可信度指数并获得“融合”雪盖

根据遥感雪盖和观测雪盖的不同组合关系(存在4种组合关系:遥感和观测都有雪,遥感有雪观测无雪,遥感无雪观测有雪以及遥感和观测都无雪),计算出雪盖可信度指数I(x,y,t),以上4种组合关系分别对应I(x,y,t)等于3、2、1和0。将雪盖可信度指数I≥1的区域(即“遥感”雪盖和“观测”雪盖),共同组成动态反演方法所需的“融合”雪盖。其中,I(x,y,t)=3有雪和I(x,y,t)=0无雪的区域可信度较高,在西部台站稀疏地区I(x,y,t)=2的区域比I(x,y,t)=1的区域可信度高,而在东部台站密集地区I(x,y,t)=1的区域比I(x,y,t)=2的区域可信度高。

第3步,针对不同的雪盖可信度指数,提出“时空距离权重法”获取动态反演所需的权重系数R(x,y,t)

初步的格点反演系数是时空动态变化的,这是其优势,但它是由测站数据空间插值得到的,而在西部地区测站稀少,当日有雪的站就更稀,致使西部遥感有雪的地区常得不到反演系数或得到的反演系数可信度低(因其距离有雪站的空间距离太远),此时可参考无时间变化的经验值(某时段的平均值)来解决这一问题,即可适当使用初步格点反演系数A0(x,y,t)年内有雪时的均值A0m(x,y)作为一定参考,因为A0m(x,y)有值的区域要比A0(x,y,t)多而且稳定。某格点距离有雪站的空间距离越远,越应更多地参考该点的“经验值”A0m(x,y),而不是“观测值”A0(x,y,t)。最终使用时空动态变化的A0(x,y,t)及无时间变化的年均值A0m(x,y)共同决定反演系数。权重系数R和1-R实际上代表在多大程度上使用这两者的结果。

A. 观测有雪(即I(x,y,t)=3或I(x,y,t)=1)时,首先,确定测站数据影响半径r0赵亮等(2010)给出了中国境内各格点到相邻最近观测站点的空间距离分布,最大值为403.8 km,为了使测站数据充分发挥作用,设定测站数据影响半径r0=403.8 km;然后,计算权重系数

其中,r(x,y,t)为之前计算的该日各格点离最近有雪站的距离。由于这里测站数据影响半径r0取为格点离站点的最大可能距离,所以,这里一定满足0≤R(x,y,t)≤1,此时其反演系数将由“观测值”A0(x,y,t)(使用了该日测站雪深插值结果)和“经验值”A0m(x,y)(初步反演系数年均值)共同决定。

B. 遥感有雪,观测无雪(I(x,y,t)=2)且年平均格点反演系数A0m(x,y)>0时,令权重系数R(x,y,t)=1,即其反演系数将由“经验值”A0m(x,y)代替。

第4步,计算动态反演系数

A. 当R(x,y,t)有值时,利用权重法计算动态反演系数

式中,存在两种极限情况:当R(x,y,t)=0时即格点与有雪站点位置重合时,完全由该日该格点初步的反演系数A0(x,y,t)作为动态反演系数A(x,y,t);而当R(x,y,t)=1时,即格点位于测站数据影响半径边缘时,完全由该位置的初步反演系数年均值A0m(x,y)作为动态反演系数A(x,y,t)。而在大多数情况下,最终的动态反演系数由初步的反演系数A0(x,y,t)和其年均值A0m(x,y)及该格点离最近有雪站的距离t(x,y,t)(决定权重系数R(x,y,t)大小)来共同决定。

B. 对于遥感有雪,观测无雪(I(x,y,t)=2),而年平均格点反演系数A0m(x,y)≤0,则动态反演系数由最近有雪格点(x′,y′)的年均值代替,即A(x,y,t)=A0m(x′,y′)。

C. 确定为无雪(雪盖可信度I(x,y,t)=0)的格点,令A(x,y,t)=0。

第5步,计算雪深

4 个例试验

以2000年1月21日为例。结合积雪分类树方法反演的当日“遥感”雪盖(图 2a)和利用最近站观测是否有雪得到的“观测”雪盖(图 2b),可明显看出,在西部高原地区,前者面积明显大于后者,这里应以遥感结果作为主要参考依据,因为这里测站稀疏,测站数据空间代表性差。但是,在长江中游及其以南地区,遥感雪盖面积又小于观测结果。由于这里测站较密,观测结果有较高的可信度,所以,两者发生矛盾时应主要以观测结果为准。可见,确定雪盖时,不能凭单一的方法确定,而应结合观测和遥感结果,针对不同地区特点,相互协同,共同决定。图 2c给出根据遥感雪盖和观测雪盖的不同组合关系获得的雪盖可信度指数I。在西部高原地区,存在大面积I=2(遥感有雪,观测无雪)的区域,鉴于西部地区两种数据的分布特点,这样的区域应予以保留,更接近实际情况;在中东部地区,存在一些I=1(观测有雪,遥感无雪)的区域,鉴于这些地区测站分布较密,观测数据可能更接近实际,这样的区域也应保留;遥感雪盖和观测雪盖都显示有(无)雪时,称之为确定有(无)雪,对应I=3(I=0)。所以,在雪盖可信度指数I≥1的区域,共同组成动态反演方法而获得“融合”雪盖。

同时,经过技术路线的第1步(数据预处理),得到当日初步格点反演系数A0(x,y,t)(0.225°分辨率,保留了ΔTb(x,y,t)<0即A0(x,y,t)<0的格点(图 2d)及其年内有雪时的均值A0m(x,y)(图 2e)。对比测站观测数据有(无)雪的站点分布(图 2b),可以看出在图 2d中,A0(x,y,t)<0的区域,实际上有相当一部分是应有雪的。所以,保留A0(x,y,t)<0的区域是合理的,当然,还应结合观测雪盖和遥感雪盖结果,将观测和遥感都无雪的区域去除。这样,中东部的反演系数精度会有所提高。但是,由于初步的反演系数A0(x,y,t)(图 2d)只用到了观测雪深,没有利用遥感雪盖(尤其没有考虑遥感有雪,观测无雪,即I=2的情况),致使遥感有雪但有雪测站较稀的地区或时次常常得不到反演系数(图 2d的西部地区),造成西部或夏季可反演的面积明显偏小,而且,即使有反演系数的地区,也可能因为离有雪站点较远,而使插值得到的初步反演系数缺乏代表性。

图 2 2000年1月21日中国地区多种要素分布(a.积雪分类树方法反演的“遥感”雪盖,b.根据最近站是否有雪得到的“观测”雪盖和有(无)雪站点分布(有雪和无雪站分别用o和+表示),c.雪盖可信度指数I(黑色实线表示“遥感”雪盖边界),d.初步的反演系数,e.2000年初步的反演系数年均值,f.动态反演权重系数,g.动态反演系数,h.动态反演雪深,i.观测雪深(Cressman插值得到),j.Chang92静态反演雪深) Fig. 2 Element distribution in China on Jan 21,2000:(a)“remote sensing” snow cover by the decision tree method;(b)“observed” snow cover and the distribution of the stations with snow(circle “o”) and without snow(plus sign “+”);(c)the snow-cover confidence index I(black solid lines denote the margin of “remote sensing” snow cover);(d)preliminary retrieval coefficients;(e)the annual mean values of preliminary retrieval coefficients in 2000;(f)the weighting coefficients of the dynamic retrieval method;(g)dynamic retrieval coefficients;(h)dynamic retrieval snow depth;(i)observed grid snow depth by the Cressman interpolation; and (j)static retrieval snow depth by the Chang92 approach

为了解决这一问题,提出用“时空距离权重法”改进反演系数,即使用时空动态变化的“观测值”A0(x,y,t)及无时间变化“经验值”A0m(x,y)共同决定反演系数。按照技术路线第2—4步,计算不同的雪盖可信度指数,将格点与最近有雪站点的距离作为计算权重的参数,计算其与测站数据的影响半径之比R(x,y,t),将这一比值与1的差值作为该日A0(x,y,t)的权重系数,这一比值作为A0m(x,y)的权重系数,两项之和最终确定为该日该格点的反演系数A(x,y,t)。经计算,从当日的权重系数R(x,y,t)(图 2f)可以看出,在中东部地区,权重系数较小,表示改进的反演系数A(x,y,t)主要由该日初步的反演系数A0(x,y,t)决定,而在西部高原地区,权重系数较大,表示改进的反演系数A(x,y,t)较大程度上由初步反演系数年均值A0m(x,y)来决定。

最终获得的该日动态反演系数A(x,y,t)(图 2g),对比之前初步的反演系数A0(x,y,t)(图 2d)可以看出,在西部高原地区,动态反演系数的面积有明显增大,台站稀疏的地区也能获得反演系数,这得益于充分利用了遥感数据,使雪盖可信度指数等于2即遥感有雪,观测无雪的区域得以保留,较大地增大了西部地区雪盖面积;在中东部地区,之前初步的反演系数由于插值方法的原因造成一些区域出现虚假的反演系数,而改进的动态反演系数A(x,y,t)有效地剔除了这些地区,这得益于确定观测雪盖时,剔除了最近站无雪的格点。仔细观察动态反演系数的分布,发现其分布特征比较复杂,空间差异很大,有些地区可以大于3,而有些地区却非常小,甚至出现负值。这暗示,若使用固定的反演系数反演雪深,必然会出现较大误差。

最后,根据式(4)反演当日雪深D(x,y,t),得到空间分辨率为0.225°的雪深分布(图 2h)。对比当日插值到格点的观测雪深(图 2i),在西部地区,动态反演雪深分布的空间连续性更好,无论是从雪盖还是雪深分布特征来看,动态反演雪深在西部地区的反演结果可能更接近实际;在中东部地区,动态反演的雪盖分布剔除了观测雪深中因插值方法造成的虚假雪盖,显然更为准确。图 2j给出利用Chang等(1992)的固定反演系数方法(简称Chang92方法)反演的当日雪深即静态反演雪深。可以看出,动态反演的雪盖面积明显大于Chang92的静态反演结果,在中东部地区更为明显,前者与实际观测更为吻合;而就雪深分布而言,在青藏高原地区,动态反演的雪深明显小于静态反演,又大于测站观测雪深。

动态反演方法的关键是反演系数的确定,中国地形复杂,西部高原地区与中东部地区地形和气候环境差异很大,不同地区的反演系数应该存在较大空间差异,动态反演方法计算的反演系数(图 2g)也反映出这一点。另外,不同的降雪过程和积雪特征的季节差异,也会导致反演系数随时间的变化,图 3给出2000年全中国、中国西部和中东部区域平均的动态反演系数逐日演变,可以看到,除了区域性差异(西部地区反演系数较低,中东部较高)外,反演系数随时间和季节的变化也十分明显,冬季偏高,夏季偏低,冬、夏季可以相差5倍,即使同一季节内,不同的降雪过程前后,反演系数也有较大变化。可以看出,用动态反演方法得到的反演系数可以较好地描述这种空间差异和时间演变。

图 3 2000年全中国、中国西部和中东部平均的逐日动态反演系数(西部和中东部以105°E为界) Fig. 3 Daily dynamic retrieval coefficients averaged over China,the western part of China and the central-eastern part of China in 2000(the boundary between the western part of and the central-eastern part of China is 105°E)
5 精度评价

为了系统地检验提出的新动态反演方法的精度,需要给出与真实场的误差,即进行精度评价。测站的观测雪深虽不能严格地代表真实值,但仍然是检验反演结果的重要依据,尤其对于中东部地区,将以观测雪深作为参考,将反演雪深与观测雪深的差作为反演误差进行精度评价。5.1 误差分析

仍以2000年1月21日为例。为了将反演结果与观测雪深比较,一种方法是将站点观测雪深插值到与反演雪深相对应的格点上,绘制反演雪深与观测雪深的散点分布,图 4给出了这两种方法计算的散点分布。比较后发现,Chang92的静态反演方法容易高估雪深,且观测有雪,反演无雪或观测无雪,反演有雪的情况较多,而动态反演方法总体上并不明显高(低)估雪深,观测有雪,反演无雪的情况也较少,并且,与观测雪深的相关也很强。

图 4 2000年1月21日反演雪深与观测雪深散点图(将站点观测雪深插值到格点)(a.动态反演,b.Chang92静态反演) Fig. 4 Scatter plots of observed vs. retrieval snow depth for each of the stations on Jan 21,2000(the station snow depth is interpolated into the grids):(a)the dynamic retrieval approach in this paper,and (b)the Chang92 static retrieval approach

下面,进一步检测误差的空间分布。将动态和静态反演雪深与插值到格点的观测雪深相减,得到误差的空间分布(图 5)。可以看出,动态反演雪深的误差和空间范围都远小于静态反演雪深。动态反演的误差大值区主要分布在青藏高原部分地区等积雪较深的区域(图 5a)。导致误差的原因,一方面可能是由于积雪越深,误差能达到的范围也就越大,另一方面是由于微波遥感反演雪深时用到了亮温差随雪深变化的线性关系,但是,在积雪很深时,它们之间并不满足简单的线性关系。Chang92方法在青藏高原地区存在较大误差(图 5b),一方面由于微波遥感在这里易高估雪深,另一方面这里的观测数据代表性差,易低估雪深,因其往往处于地势相对周围平坦和海拔较低的地区。所以,在青藏高原地区出现比观测轻微高估的结果,更为符合实际,而动态反演方法得到的雪深正是如此。此外,当天的格点场总体误差情况如图 6,动态反演的误差分布在0附近的百分率远高于Chang92方法反演结果。

图 5 2000年1月21日反演雪深误差的空间分布(a.动态反演,b.Chang92静态反演) Fig. 5 Spatial distributions of error(a)the dynamic retrieval snow depth; and (b)the Chang92 static retrieval snow depth
图 6 2000年1月21日反演雪深格点场误差百分率分布(a.动态反演,b.Chang92静态反演) Fig. 6 Histograms of total error percentage(a)the dynamic retrieval snow depth; and (b)the Chang92 static retrieval snow depth
5.2 准确率分析

为了检验新的动态反演方法在四季更替中的精度变化,采用赵亮等(2010)反演准确率的定义方法,将反演雪深插值到站点,对于有雪的站,若反演与实测结果误差在误差阈值范围内,则认为该站为有雪合理站,否则,视为有雪不合理站,据此统计出有雪合理站数占总有雪站数的百分率,称之为反演准确率Ra;另外,若再将实测和反演结果都无雪的站视为无雪合理站,则可以统计出总合理站数(有雪合理站数与无雪合理站数之和)占总站数的百分率,称之为反演准确率Rb。据此,计算得到2000年1月21日的反演准确率(表 1)。不同的误差阈值下,动态反演方法均比Chang92静态方法准确率高。为了进一步验证在其他季节动态反演方法的精度,图 7给出2000年全年逐日雪深的反演准确率RaRb。对于动态反演方法,其准确率基本处于80%上方,而Chang92方法的反演准确率在晚冬和初春较低。

表 1 动态与Chang92静态反演方法反演的2000年1月21日准确率比较(单位:%) Table 1 he comparison of accuracy rates between the dynamic retrieval approach and the Chang92 static retrieval approach(unit: %)
有(无)
雪站数
反演
准确率
-3 cm<误差<3 cm -5 cm<误差<5 cm -10 cm<误差<10 cm
动态反演Chang92动态反演Chang92动态反演Chang92
268(411)Ra96.641.198.557.898.884.0
Rb85.660.586.467.086.577.2
图 7 2000年动态和Chang92反演方法的雪深反演准确率(曲线,误差阈值:5 cm)与有雪站数(直方图)逐日演变(a.反演准确率Ra,b.反演准确率Rb) Fig. 7 Daily series of snow depth retrieval accuracy rates(error threshold: 5 cm)of the dynamic and Chang92 retrieval approaches and the number(histogram)of stations with snow in 2000:(a)retrieval accuracy rate Ra; and (b)retrieval accuracy rate Rb(see the text for further details)
5.3 交叉验证

交叉验证技术是一种客观的检验方法,将测站样本分为反演和校验两部分,用反演样本计算出校验样本的推测值,将其与实测的校验样本值比较(吴洪宝等,2010)。采用均方根反演误差可以定量衡量反演值与实测值的接近程度

其中,ZeiZoi分别为第i时次反演值和观测值,n为总时次数。

分别选取东北(50136,黑龙江漠河)、西北(51379,新疆奇台)、西部高原(55294,西藏安多)和东部雪线附近(58215,安徽寿县)4个站作为交叉验证的校验参考站。一并剔除这4个站后,再采用本文的动态反演方法反演出这4个站的雪深。将结果与实测雪深和静态反演方法进行比较(图 8)。交叉验证的结果显示,相比Chang92方法的静态反演而言,这4个站的动态反演结果都更接近观测。除58215站外,Chang92方法的结果在冬春普遍有较严重的偏高,而动态反演结果与观测偏离相对小些。其中,50136站(可以代表东北地区),冬季偏离不大,春、秋季偏离相对大些,Chang92和动态反演的均方根反演误差分别为12.7和7.6 cm;51379站(代表西北地区),1月偏小,其余季节偏离不大,Chang92和动态反演的均方根反演误差分别为12.5和5.1 cm;55294站(代表青藏高原地区),冬春反演结果比观测大一些,但没有Chang92偏离严重,Chang92和动态反演的均方根反演误差分别为10.8和3.4 cm,而且,动态反演的春季融雪速度不是很快,这与张佳华等(2008)的分析结果一致,可能更接近真实情况;58215站(代表东部雪线附近),动态反演结果很接近观测,比Chang92的静态反演效果好得多,Chang92和动态反演的均方根反演误差分别为0.58和0.19 cm,这得益于其融合了该站附近的观测结果且这里的测站分布较密。因此,本研究新的动态反演方法的效果总体来说还是较好的,但在东北地区春、秋季误差较大,青藏高原地区反演结果比观测高一些。由于篇幅所限,这里所进行的交叉验证还是比较初步的,未来还可以采用随机或逐个抽取样本的方式进行交叉验证。

图 8 4站逐日的动态反演(剔除4站后重新反演)、Chang92静态反演和观测雪深(a.黑龙江漠河(50136),b.新疆奇台(51379),c.西藏安多(55294),d.安徽寿县(58215)) Fig. 8 Snow depth daily series of the dynamical retrieval approach by use of data with the 4 stations of 50136,51379,55294 and 58215 excluded,the Chang92 static retrieval approach and the observation for the 4 stations((a),(b),(c),and (d),respectively)in 2000
6 总结和问题6.1 总 结

本研究提出了新的动态反演参数化方案(时空距离权重法),开发了利用被动微波遥感与地面观测数据联合反演雪深的新方法。该方法一方面因利用了微波遥感数据,因而能够在测站稀少的西部高山地区反演出精度相对较高的雪深数据;另一方面,因结合了实时实地的测站雪深资料,因而反演系数能够随时空灵活自动地调整,避免了固定反演系数的静态反演方法在高山地区严重高估而在平原地区低估雪深的情况。个例检验表明,动态反演系数随时空变化,使不同区域不同时间的反演方程具有不同的形式,这在一定程度上克服了积雪物理性质的时空差异带来的反演误差。精度评估表明,该动态反演方法反演的积雪空间分布连续性好,在站点稀疏区域也能得到较合理的雪深数据。与原有方法(赵亮等,2010)相比,新方法在雪盖确定方面有重要改进,不再依赖单一的地面测站数据及其空间插值结果,而是利用这两种不同来源的数据共同确定雪盖,所得雪盖分布更加合理(尤其在青藏高原地区),雪深分布准确性也有明显提高;此外,考虑了在不同时空条件下遥感和观测数据各自的可信度,设计了新的动态参数化方案,这使得反演系数的确定更加合理,为提高反演精度提供了保证。与静态反演方法相比,在华北和华中低估雪盖面积而西部高估雪深的缺点基本被克服。对比国际上较新的ANSA算法(Foster et al,2011),在ANSA算法表现较弱的有云的高山地区,新提出的被动微波结合测站数据的动态反演算法所得雪深表现较好,不存在明显高估,只是雪盖面积可能有所高估,这主要是因为目前这种算法还没有融入可见光数据的缘故。6.2 存在问题和发展方向

目前,国际上融合多种遥感数据反演雪深的研究已经取得重要进展并形成产品(Brubaker et al,2005Foster et al,20072011Hall et al,20072012),可以预见,这种多源数据融合技术将是未来遥感反演的一个重要方向。本研究发展了观测雪深与被动微波反演雪深的融合技术,但是,目前存在的不足是:还没有融入可见光、散射计等其他遥感数据,所以,在无云时雪盖的确定和融雪探测等方面还有很大提升空间,将来需尽快融入更多来源的遥感数据以提高反演准确性。另外,目前该方法没有考虑反演系数随海拔高度的变化,其插值方法只涉及二维,未来会结合地理信息,考虑在第三维即垂直方向上也能进行合理的插值化处理。不过,这种改进的遥感-测站相结合的雪深动态反演方法,在初步应用中已经表现出可喜的优越性,相信对这种方法的深入研究和完善将有利于弥补在积雪资料研究中的一些不足。

致谢:感谢国家卫星气象中心许健民院士和刘玉洁研究员的宝贵意见。

参考文献
车涛,李新,高峰.2004.青藏高原积雪深度和雪水当量的被动微波遥感反演.冰川冻土,26(3):363-368
陈爱军,刘玉洁,杜秉玉.2005.应用AMSU资料监测中国地区雪盖的初步研究.应用气象学报,16(1):35-44
陈乾金,高波,李维京等.2000.青藏高原冬季积雪异常和长江中下游主汛期旱涝及其与环流关系的研究.气象学报,58(5):582-595
冯锦明,赵天保,张英娟.2004.基于台站降水资料对不同空间内插方法的比较.气候与环境研究,6(2):261-277
李培基.1996.喜马拉雅山积雪与印度季风降水呈明显反相关关系的商榷.气象学报,54(3):379-384
李晓静,刘玉洁,朱小祥等.2007.利用SSM/I数据判识我国及周边地区雪盖.应用气象学报,18(1):12-20
李新,车涛.2007.积雪被动微波遥感研究进展.冰川冻土,29(3):487-496
刘玉洁,郑照军,王丽波.2003.我国西部地区冬季雪盖遥感和变化分析.气候与环境研究,8(1):114-123
罗勇.1995.青藏高原冬春季雪盖对东亚夏季大气环流影响的研究.高原气象,14(4):505-512
吴洪宝,吴蕾.2010.气候变率诊断和预测方法.北京:气象出版社,400pp
延昊,张佳华.2008.基于SSM/I被动微波数据的中国积雪深度遥感研究.山地学报,26(1):59-64
杨秀春,曹云刚,徐斌等.2009.多源遥感数据协同的我国草原积雪范围全天候实时监测.地理研究,28(6):1704-1712
张佳华,吴杨,姚凤梅.2008.卫星遥感藏北积雪分布及影响因子分析.地球物理学报,51(4):1013-1021
赵亮,朱玉祥,程亮等.2010.遥感-测站相结合的动态雪深反演方法初探.应用气象学报,21(6):685-697
Akyurek Z, Hall D K, Riggs U A, et al. 2010.Evaluating the utility of the ANSA blended snow cover product in the mountains of eastern Turkey. Int J Remote Sensing, 31(14):3727-3744
Armstrong R L, Brodzik M J. 2001. Recent Northern Hemisphere snow extent; A comparison of data derived Irom visible and microwave satellite sensors. Ueophys Res Lett, 28(19):3673-3676
Armstrong R L, Knowles K W, Brodzik M J,et al. 2003. DMSP SSM/I Pathfinder daily EASE-Urid Brightness temperatures. Boulder, Colorado USA; National Snow and Ice Data Center, Digital media
Brubaker K L, Pinker R, Deviatova E. 2005. Evaluation and comparison of MODIS and IMS snow cover estimates for the continental United States using station data. J Hydrometeor, 6(6) 1001-1017
Chang A T C,Foster J L, Hall D K, et al. 1992. The use of micro-wave radiometer data for characterizing snow storyge in western China. Ann Ulaciol, 16: 215-219
Chang A T G,Foster J L, Hall D K. 1987.Nimbus-7 SMMR derived global snow cover parameters. Ann Glaciol, 9: 39-44
Chen C,Nijssen B, Guo J,et al. 2001. Passive microwave remote sensing of snow constrained by hydrological simulations. IEEE Trans Geosci Remote Sensing, 39(8):1744-1756
Cohen J L, Entekhabi D. 1999. Eurasian snow cover variability and Northern Hemisphere climate predictability. Ueophys Res Lett, 26(3):345-348
Cohen J L, Saito K. 2003. Eurasian snow cover, more skillful in predicting U. S, winter climate than the NAO/AO?. Ueophys Res Lett, 30(23):doi:10. 1029/2003UL018053
Cressman U P. 1959. An operational objective analysis system. Mon Wea Rev, 87(10):367-374
Foster J L, Sun C J,Walker J P, et al. 2005. Quantifying the uncertainty in passive microwave snow water equivalent observations. Remote Sensing Environ, 94(2):187-203
Foster J L, Hall D K, Eylander J B, et al. 2007. Blended visible(MODIS),passive microwave (AMSR-E)and scatterometer(QuikSCAT) global snow products//Proceedings of the 64th Annual Eastern Snow Conlerence. St Johris Newfoundland, 27-36
Foster J L, Hall D K, Eylander J B, et al. 2011. A blended global snow product using visible, passive microwave and scatterometer satellite data. Int J Remote Sensing, 32(5):1371-1395 Frei A, Robinson D A. 1999. Northern Hemisphere snow extent Regional variability 1972-1994 .
IntJ Climatol, 19(14):1535-1560 Grody N C,Basist A N. 1996. Global identification of snowcover using SSM/I Measurements. IEEE Trans Geosci Remote Sensing, 34(1):237-249
Hall D K, Riggs U A. 2007. Accuracy assessment of the MODIS snow products. Hydrol Process, 21(12):1534-1547
Hall D K, Foster J L, Kumar S, et al. 2012. Improving the accuracy of the AFWA-NASA(ANSA) blended snow-cover product over the Lower Ureat Lakes region. Hydrol Earth Syst Sci Discuss, 9(1):1141-1161
Hellrich S R, McNamara D, Ramsay B H, et al. 2007. Enhancemenu to, and forthcoming developments in the Interactive Multisensor Snow and Ice Mapping System(IMS).Hydrol Process, 21(12):1576-1586
Kelly R E, Chang A T, Tsang L, et al. 2003. A prototype AMSRE global snow area and snow depth algorithm. IEEE Trans Geosci Remote Sensing, 41(2):230-242
Kelly R E J. 2009. The AMSR-E snow depth algorithm: Description and initial results. J Remote Sens Soc Japan (ULI/AMSR Special Issue),29(1):307-317
Kongoli C,Dean C A, Hellrich S R, et al. 2006. The retrievals of snow cover extent and snow water equivalent Irom a blended passive microwave-interactive multisensor snow product//Proceedings of the 63rd Eastern Snow Conference. Newark, 89-103
Pulliainen J,Hallikainen M. 2001. Retrieval of regional snow water equivalent Irom space-borne passive microwave observations. Remote Sensing Environ, 75(1):76-85
Yasunari T. 2007. Role of land-atmosphere interaction on Asian monsoon climate. J Meteor Soc Japan, 858; 55-75