3.1 方法原理
本次研究的能量与水平衡系统(Energy and Water Balance Monitoring System,EWBMS),用于从接收的卫星云图中生成降雨量图像,实际蒸散、降雨量和温度,最后生成的数据作为干旱监测和大尺度水文模型(Large Scale Hydro Model,LSHM)的主要输入数据,其流程如图3-1所示。静止卫星FY2C或GMS的图像每小时接收一次,可以确定云层出现频率或“云生存期”。从所有接收的每小时图像数据中,提取代表正午和子夜的VIS和TIR的合成图,然后将这些提取出来的数据处理成定量的空间连续的降雨量、辐射、显热通量、温度和蒸散的数字化图像。除了这些卫星图像外,几乎不需要外界输入数据。
图3-1 能量与水平衡系统(EWBMS)流程图
3.1.1 预处理
预处理根据卫星接收系统得到的每小时的VIS和TIR图像中,计算云的生存期、当地正午和当地子夜的合成图,云的生存期图像只用TIR图像。所观察物体在红外光谱测得的辐射率与该物体的温度直接相关。卫星观察到的物体则是地球表面或最高层的云顶。云的温度与它高于地面的高度成比例:典型的温度垂直梯度是-6.5℃/1000m。通过对云图的直方图分析,可以区分4种云层。TIR灰度值的阈值可以转换成行星温度,其相应的温度和高度见表3-1。
表3-1 云层的定义和相应温度与高度
每小时都有新的TIR图像被接收到,对于每一个像元需要确定是否有云,若有云,需要确定是哪类云。这些结果将被储存在4个文档中(每种类型的云有一个文档)。这些文档每天更新,这样可以生成每日多层云类频率文档。
能量平衡的主要输入数据是当地正午和子夜的合成图像,它们代表整个图像在正午和子夜的情形。在预处理时,这些图像是通过每小时的VIS和TIR图像进行合成的。对于图像中的每个像元,当地的正午或子夜时间是根据它的经度位置进行计算的,并用GMT制式(格林尼治标准时间)表示,然后,选择两幅由GMT时间表达的最接近当地正午或子夜的图像,由选取的两幅图像内插得到合成图像中的每一像元的值。这样,可以生成当地正午的VIS和TIR合成图像,但是子夜只有TIR的合成图像。对于某一像元,当与其最接近的当地正午或子夜两幅图像中的一幅没有时,则不需要插值。所以,当有图像遗漏时,合成图像中会有明显的线。当遗漏的用于生成合成图像的小时云图达到或多于3幅时,这一天就不能生成正午或子夜的合成图像了。
3.1.2 降雨图绘制
空间降水估计的方法基于两种信息:①从气象站得到的降水量点数据;②从风云二C静止气象卫星得到的云出现频率数据。点数据可以从世界气象组织(WMO)的全球远程通信系统(GTS)和中国气象局(CMA)的气象站获得。GTS收集了全球约1.1万个气象站中的数据,其中95%的数据可以通过GTS在6h内获得。在中国,有400个气象站通过WMO-GTS网络报告降水数据,中国气象局提供另外800个国内气象站的数据。
至今,人们已经开发出一些用气象卫星数据生成降雨量场的方法。普遍采用的有冷云生存期(CCD)技术,它将出现的非常高和“冷的”云与测量的雨量联系起来,历史的数据用于标定,CCD技术很适合于对流降雨的估计。EARS方法与此则有两点不同:首先,它采用了4种云高度和“温度阈值过余值”,而常规的只有一种云高度,所以它同时也考虑了低云和锋面降雨;其次,方法中没有根据历史数据进行标定,而是将雨量站的数据与云生存期进行近实时的结合。
EWBMS降雨处理首先是对卫星得到的冷云数据和有雨量站的像元的降水数据进行多次“局部”回归。这个“局部”回归是基于这个站与周边相邻的12个站。对站点j所得的回归方程为
式中 CDn——n层云的生存期(频率)。
然而,这些回归方程并不能很好地估计降水P。所以,对于每个站都需要确定估计值与观察值的残差,即
式(3-1)中的回归系数aj,n、bj和式(3-2)的残差Dj,用距离倒数加权法在6个雨量站间进行内插,以得到像元i的对应值。最后,空间分布降水将逐个像元进行计算,即
注意:在站点位置的降水估计值总是等于站点的报告降水值。在降雨量的点数据中,这种影响并不明显,因为测量站通常是处于相对比较低的地方。众所周知,降水量取决于地表和约11km高的对流层顶间的可降的水量。所以,可以预期,降水量跟地表与对流层顶间的大气柱的高度或质量成正比(图3-2)。
图3-2 地表的能量平衡
3.1.3 能量平衡方程
作为EWBMS的一部分,能量平衡监测系统的目的是求解地表的能量平衡方程,它可以写成
式中 LE——潜热通量,W/m2;
In——净辐射,W/m2;
H——显热通量,W/m2;
E——光合电子输运量,W/m2;
G——土壤热流通量,W/m2。
地表反照率和地表温度是主要的输入数据。用于地表水蒸发的能量(LE)等于提供给地表的净辐射能(In)减去加热空气的显热通量(H)、植被的光合电子输送量(E)和用于加热土壤的能量(G)。对于一天而言,土壤的热流通量可以认为是零(G≈0)。所以,地表能量平衡方程可以重新写成
只有正午的卫星图像才能用于处理能量的通量。通过对太阳日周期的傅里叶分析,可以从正午的辐射和显热通量得到相应的日平均值。另外,还需要相应的地理位置和天数信息。该方法假定这天大气传输的条件不变。修正的日辐射考虑了正午前后云覆盖的影响。
1.大气校正
通过VIS和TIR数据的标定,可以得到行星反照率和温度,即观察到的穿过大气的值。然而,为了计算地表能量平衡的不同的能量组成,需要地表反照率和温度,为了完成这个转换,必须执行大气校正程序。大气的吸收和散射使得从卫星可见光波段观察到(行星反照率)的反射率与地面实际反射率间的差异。大气中的吸收主要是水汽造成的,散射是由于空气分子(如N2、O2)和空气中的悬浮物质造成的。由于受水汽对红外辐射的吸收和反射的影响,行星温度低于实际的地表温度,而散射的影响较小。大气校正程序不需要大气成分和分层的信息,而使用图像相对信息,如对比度。图像对比度随着大气扰动的增强而下降。
对于可见光波段,采用Kondratyev的总辐射传输模型(1969)。这个模型用于太阳辐射的直射和漫射。对于向下和向上的总辐射通量分别建立联立微分方程。向下的是直射,而向上的是漫射(图3-3)。Kondratyev的辐射传输模型只解释了后向散射并忽略了大气对辐射的吸收。然而,大气吸收的量级为10%。在这个模型中引入了一个吸收因子(k)。在图3-3所示的微元δt中,对总辐射后向散射和吸收受到的影响作了修正。近似的微分方程为
图3-3 向下和向上辐射通量的图示
式中 α——光的后向散射系数(≈0.1);
k——光的吸收系数(≈0.03);
is——太阳的天顶角。
这个微分方程可以得到分析解,这样就得到了两个函数。一个是地表反射率与行星反射率和光学厚度的关系;另一个是通过无云大气的太阳辐射与地表反射率和光学厚度的关系。光学厚度是对大气中光学上有影响的物质数量的度量。
式中 A'——观察的行星反射率;
A——表面反照率;
τ——大气的光学厚度。
图3-4揭示的是地表的反射率和所吸收的太阳辐射,是行星反射率和光学厚度的函数。地表所吸收的太阳辐射定义为1减去地表反射率,即1-A,乘上通过大气的透射率(t)。一旦光学厚度确定,式(3-9)可将每个像元的行星反射率转换成地表反射率。当地表的反射率最小时,光学厚度的影响最大,可以从茂密的植被区中发现这种情形。为了确定日平均光学厚度,第一步是对图像中的每一像元求出旬(一旬为10d)的最小行星反射率;随后,可以得到最小行星反射率的“最黑”像元,它们是跟地表最小反射率相关的。对于图像中足够茂密的森林,典型的地表最小反射率值为7%。从最小行星反射率和最小地表反射率可以计算出光学厚度。定出光学厚度后,将它应用于整个图像,可将每一像元的行星反射率转换成地表的反射率。
图3-4 地表的反照率和所吸收太阳辐射是行星反照率和光学厚度(τ)的函数
比较传输模型和总辐射,需要对EWBMS的总辐射的输出值做一些调整,以达到二者匹配。所以,在式(3-9)中引入两个标定系数,即
这样做有两个理由:一是通过大气的透射率是在正午确定的,但从全天来讲,有效的输送要比该值小,因为大多数时间太阳是倾斜的;二是太阳辐射通过大气的透射率,平均来讲,小于通过气象卫星可见光的窗口(0.55μm~0.9μm)的透射率。采用C1=0.77和C2=0.65的标定系数使结果能很好地与观察的总辐射的值匹配。
对于热红外波段的大气校正,将采用其他方法。行星温度(T'0)与地表温度(T0)之间的关系可用式(3-11)表述,即
式中 k——大气校正系数;
im——卫星天顶角;
Ta——大气边界层顶端的空气温度,K。
在边界层顶部的空气温度(Ta)可以根据正午和午夜的像元温度的线性回归得到,如图3-5所示。对于理想传热的情形,此时的空气温度的估计值为T0,noon=T0,midnight=Ta。大气边界层厚度白天通常在1~2km范围内波动。覆盖整个区域的边界层顶空气温度图,可以将该方法应用于一个移动的200×200km的窗口而求得。为了计算校正系数,选用图像中最干的像元,并假定其相应的蒸散为0。对于每个像元,计算由式(3-12)定义的干燥指数(DI),即
图3-5 参考温度与正午和午夜行星温度散点图的偏差
式中 In——净辐射。
对于图像中最干的像元,假定其潜热通量(LE)为0,因此,净辐射等于显热通量(H)。一旦最干像元的显热通量确定,相应的地表温度可由净辐射和空气温度T0=Ta+In/α计算获得。由于方程式(3-8)的校正系数(k)是唯一的未知参数,它的值可以由计算获得。校正系数可以用于整个图像。已知校正系数和空气温度,就可以逐一像元计算地表温度。这些地表温度可以用于计算显热通量。
2.云检测
EWBMS系统能在有云或无云情况下绘制蒸散图像地图。所以,采用了云检测算法,将有云像元从无云像元中分离出来。当像元中有云或部分有云时,若错将其划分成无云的类别,蒸散会被过高估计,所以要尽量避免这种情况的发生。
在云检测中,对可见光和红外图像采用了阈测试,并定义了4种云检测准则。若符合4种准则中的一种,像元将被标为有云。这4种测试如下:
(1)A'≥A'min+阈值1。
(2)T'0,noon≤T'0,noon,max-阈值2。
(3)T'0,noon≤T'0,midnight-阈值3。
(4)T'0,noon≤Ta。
式中 A'min——旬中观察的最小行星反照率;
T'0,noon,max——正午无云行星温度,K。
第(1)种和第(2)种的云检测是根据观察的像元值和无云时的像元值二者的比较。为了得到在测试(1)中无云情况下的最小行星反射率的值,由连续10个正午可见光合成图组成最小反射率图。假定在这10d中(旬)每个像元至少有一次是无云的。测试(1)只能应用于白天时段。
测试(2)将像元行星温度与无云行星温度进行比较。这里采用了一种特殊的方法估计无云的行星温度。这种方法使用了1000km×1000km见方视窗的热红外通道数据。在这个视窗中确定那些具有最高行星温度的点,并假定哪些像元代表无云的点。然后,视窗中心像元的温度可由距离加权的平均方法计算。这种方法在欧洲由De Valk等人(1998)验证。测试(2)被用来检测正午和午夜热红外图像中的有云像元。对于第(1)种和第(2)种的测试,它们各自能测出约70%的有云像元,若两者结合,检测率增加到80%。测试(3)、(4)只是在测试(1)、(2)的基础上额外增加一些有云的像元。测试(1)、(2)中所用的阈值需由经验确定。在Metclock项目的框架下,De Valk等(1998)确定了测试(1)、(2)中的阈值。他们发现,同样的阈值可用于冬季和夏季。
采用这些测试后,发现在云的边缘的相对蒸散远大于略离云边缘处的值。高的相对蒸散意味着行星温度较低,而行星温度较低表明像元有云或部分有云。当完成能量平衡监测系统产品的处理后,一种特别的程序被用来检测那些邻接有云像元和那些其相对蒸散远高出周边无云像元的值的像元。这些像元将被标记为有云,并返回到能量平衡程序中重新计算。
3.总辐射
通过大气的总辐射透射率(t)可以由Kondratyev模型计算。正午地表的总辐射()可由式(3-13)计算,即
式中 i——正午太阳天顶角;
S——大气边缘太阳辐射的强度,1355W/m2。
下一步是将正午的总辐射()转换成总辐射的日平均值(Ig)。转换因子的值v=Ig/将通过对太阳日周期的积分来确定,它是纬度和天数的函数。当像元标志为有云时,可以用Kubelka-Munk理论推算通过云的太阳辐射(tc):
式中 Acb——积雨云的反照率(=0.92)。
为了确定每一个像元的地表的反照率(A),将合成每日的可见光图像,构作旬最小地面反射率图,并假定在这个旬中至少有一次这个像元是无云的。图3-6是在3种不同地表反射率情况下的行星反射率(有云)与tc的关系。有云时正午地表的总辐射()可以按式(3-15)估计,即
将有云时正午的总辐射转换成日平均值的方法与无云时的情形一致。
4.云的生存期
图3-6 通过云层的太阳辐射
若在能量平衡监测系统中只使用正午的图像,当正午为无云但早上和下午有云时,辐射(和其他的能量平衡通量)会被高估。同样,当正午为有云但早上和下午无云时,辐射被低估。为了得到更准确的日辐射值,在系统中将采用描述当地上午9时至下午3时的云覆盖量的小时热红外图像。生成的日云覆盖信息文件给出了上午9时至下午3时的云量信息。这些图将作为能量平衡监测系统的额外的输入数据,它能显著地改善日辐射的测量值与计算值间的关系。
红外通道的云检测采用阈值法,它是当地时间一年中的天数和纬度的函数。这是因为早上的温度较低,若采用中午的值而不针对早上的低温做调整,像元会被归类于有云。对于一年中的天数和纬度,它的情况也一样。同样也应考虑早上/下午与正午之间的太阳入射辐射的差异。用于正午合成图中的辐射数据主要取决于正午的情况,如当地11时和13时的云覆盖情况比当地时间早上9时和15时的云覆盖情况更重要。
5.净辐射
地球表面的净辐射(In)由短波(太阳)和长波(陆地、热、红外)辐射通量的差值计算而得,并可表达为日平均值。
式中 Ig——地表的日平均总太阳辐射,W/m2;
ε0La——向下所吸收的长波(热)辐射,W/m2;
L0——地表发射的长波辐射,W/m2;
ε0——地表发射率;
σ——斯蒂芬—玻尔兹曼常数,W/(m2·K4);
εa——大气发射率。
净长波辐射(Ln)可以从地表和大气温度和发射率进行计算。地表发射率(ε0)通常在0.85(沙漠)~0.95(植被)间变化。假定平均值为0.9。大气的发射率是(εa)从Brunt经验公式中导出的,采用空气比湿(Sa)作为输入量,即
式(3-17)可以转化成
这里=(T0+Ta)/2为平均温度。通常,气候净长波辐射(Lnc)在80W/m2这个量级。方程式(3-19)的第一项可以被称为辐射热流通量(Hr),它取决于地表与边界层顶端的温差。在能量平衡方程中,Hr实际上与显热通量(H)可以结合在一起。
在计算净辐射时,当像元被云覆盖时,可以假定云下面的长波的辐射可以略去(Ln=0)。然后,净辐射可以按式(3-20)计算,即6.显热通量
在无云的正午,与大气的热交换可用式(3-21)计算,即
式中 C——阻力系数,W/(m3·s·K);
Va——平均风速,m/s;
ε0——地球表面发射率;
σ——斯蒂芬—玻尔兹曼常数,W/(m2·K4);
——平均温度,K,=(T0+Ta)/2。
大气的显热换热系数(α)是对流显热换热系数(αc)和辐射显热换热系数(αr)之和。采用平均风速和地表发射率,显热通量(H)就可以从地表和大气的温差确定。
重写式(3-21),得
为了将正午的显热通量与日平均显热通量联系起来,假定Bowen比,即潜热与显热的能量比,在日周期中为一常数。随后,从式(3-25)中可以得到日平均显热通量(H),即
7.光合能耗
当有植被时,一部分的太阳能用于光合电力输运(E)固化二氧化碳。这部分能量估计为Cv——植被覆盖率。
式中 ε——日平均光合利用系数;
光合利用系数可由光合非活化模型(Rosema等,1998)估计。植被覆盖率Cv并没有预先确定。但有植被处是以高蒸散值为特征的。所以,用相对蒸散(LE/LEp)作为作物覆盖的一种表征。
8.实际蒸散
净辐射(In)、显热通量(H)和光合电子输运(E)确定后,潜热通量(即用能量单位表示的实际蒸散)可以从式(3-28)得到,即
土壤的热流通量,从日平均来讲,可以略去不计。若像元标志为有云,地表的温度未知。这样,显热通量(H)和潜热通量(LE)无法计算。在这种情况下,假定其Bowen比(β=H/LE)与前一个无云天相当。Bowen比是地表湿润情况的一个指标,假定其在被云覆盖时不变。这样,可以估计实际蒸散,即
9.1.5m高温度
高于地表1.5m处的空气日平均温度(T1.5m)可以由大气边界层顶部的空气温度(Ta)和地表温度确定(T0)。地表的日平均温度的定义是正午地表温度(T0,n)和子夜地表温度(T0,m)的平均值。通过黄河流域25个WMO-GTS站的1.5m高温度和卫星数据,加权函数加权后可得到
T1.5m温度跟10日的气温差会大于10℃,将用10d的值进行替代,以修正卫星温度测量的误差。