利用Himawari-8卫星红外图像反演降雨

孙绍辉1,2 李万彪1,† 黄亦鹏1

1.北京大学物理学院大气与海洋科学系, 北京 100871; 2. 93313部队, 长春 130111; †通信作者, E-mail: lwb@pku.edu.cn

摘要 利用静止气象卫星 Himawari-8 的红外图像数据, 结合全球降雨测量计划的降水产品 2AGPROFGMI, 通过将匹配的红外亮温和亮温差数据点与近地面雨强数据点一一对应, 使二者之间建立联系。利用静止气象卫星通道多的优点, 进一步建立多个以红外亮温(BT)、亮温差(BTD)以及近地面雨强为参数的二维和三维雨强查算表。通过反演试验, 找到能够较好地识别不同等级降雨强度的查算表, 其中三维查算表(BT10.4, BTD12.4-10.4, BTD6.2-7.3)对降雨事件表现出较高的检测概率(POD=0.8817)和较低的错判率(FAR=0.4042), 可以较好地判识降雨区域。

关键词 卫星遥感; 亮温; 亮温差; 雨强; 反演

降雨的准确测量对工农业生产、水资源利用以及洪涝和干旱等自然灾害的预报有重要意义[1]。由于地面雨量计和地基测雨雷达观测网的密度不足, 很难监测和跟踪暴雨及其他灾害性天气系统。利用卫星遥感数据反演降雨, 覆盖区域广, 时空分辨率高, 可以弥补地面常规观测的不足。

卫星遥感可分为主动微波遥感、被动微波遥感和可见光/红外遥感。主动、被动微波遥感可以分别探测降雨的垂直结构和云内部的结构信息, 但目前微波遥感器搭载在极轨气象卫星上, 时间分辨率不高。静止气象卫星可见光/红外遥感可以提供高时间分辨率的降水产品, 对降雨起到很好的监测作用[23]。因此, 加强利用静止气象卫星多光谱遥感测量信息反演降雨的研究很有必要。

Ba等[4]采用温度梯度、坡度参数[5]和云顶粒子有效半径[6], 利用地球静止业务环境卫星(Geosta-tionary Operational Environmental Satellite, GOES) 0.65, 3.9, 6.7, 11 和 12 μm 5 个通道建立 GOES 多光谱反演算法(GOES multispectral rainfall algorithm, GMSRA)。Behrangi等[7]利用GOES-R[8]卫星多光谱图像, 采用PERSIANN-MSA方法估计降雨, 其中方案12对有雨和无雨的判识表现最好, 降雨检测概率(POD)为0.530, 降雨错判率(FAR)为0.355。Tao等[9]基于卫星的红外和水汽通道信息, 构建深度神经网络来估算降雨, 对研究区域有雨和无雨的判识指标POD为0.418, FAR为0.528。Zhuge等[10]运用红外亮温以及可见光反射率, 构建降雨概率判识矩阵来区分降雨云和非降雨云, 当降雨概率设定为50%时, 得到POD=0.6651, FAR=0.4685。金晓龙等[11]通过分析 TRMM (Tropical Rainfall Measuring Mission), CMORPH (Climate Prediction Center’s MORPHing technique)和GPM (Global Precipitation Measurement)卫星降水产品在山区的适用性, 发现这 3 套产品在探测弱降雨事件方面均具有较高的准确性(POD≈0.58)和较低的错判率(FAR≈0.63)。

很多现有的反演降雨研究都使用可见光、红外窗区和水汽等多个通道信息, 而Himawari-8卫星遥感资料包含红外窗区外的其他红外通道信息, 一些包含降雨信息的通道有助于改进降雨反演。此外, 采用查算表的研究方法目前比较少见。本研究利用Himawari-8多通道红外亮温数据[12]和GPM降雨资料数据[1314], 建立以匹配的红外亮温(brightness temperature, BT)、亮温差(brightness temperature difference, BTD)及近地面雨强(rain rate, RR)为参数的雨强查算表来估计降雨。

1 资料

1.1 资料介绍

本文使用的 Himawari-8 卫星数据来源于日本宇宙航空研究开发机构地球观测研究中心。数据已经过辐射定标和地理定位等预处理。本文使用的是一级网格化数据: 空间分辨率为5km, 时间分辨率为10min, 覆盖范围为60°N—60°S, 80°E—160°W。Himawari-8是新一代静止气象卫星, 搭载的辐射成像仪AHI(Advanced Himawari Imager)有16个通道, 本文使用其 10 个红外通道中中心波长为 6.2,7.3,8.6,9.6,10.4,12.4 和 13.3μm 的 7 个通道的亮温资料。全圆盘观测仅需10分钟, 具备多波段通道、时空分辨率高的优点, 有利于监测强降雨云团的发展和变化。数据取自网站http://www.eorc.jaxa.jp/ptree/ index.html。

GPM核心卫星GPMCO(GPMCoreObservato-ry)是热带降雨测量卫星TRMM的接续卫星, 覆盖范围为 65°N—65°S, 可观测整个中国地区。本文采用GPM 的二级降水产品 2AGPROFGMI, 该降雨资料由 GPROF2014 (Goddard profiling algorithm)廓线反算法得到[15]。Kummerow等[16]研究GPROF估算降雨的可靠性, 结果表明该算法对海洋和陆面的反演结果与地基雷达和地面雨量计探测结果的相关性分别达到 0.86 和 0.8。Petkovic 等[17]在利用 GPROF 研究极端降雨时, 发现该算法的反演结果与地基雷达探测结果较一致, 故采用降雨资料2AGPROFGMI的近地面雨强数据作为降雨真值比较可靠。资料取自网站https://www.gportal.jaxa.jp/ gp/top.html。

1.2 资料匹配

为了使Himawari-8资料与2AGPROFGMI资料在空间和时间上相匹配, 本文采取以下措施。

1)选取两种资料的经纬度范围为15—45°N, 90 —130°E, 满足此条件的GPMCO轨道每天有3~4条。

2)本文使用的Himawari-8资料的空间分辨率为5km, 2AGPROFGMI资料的空间分辨为8km´15km(随扫描位置变化)。在选定范围内, 以0.05°´0.05°的经纬度网格, 分别对Himawari-8和2AGPROFGMI的数据网格化, 得到选定范围内有601´801个单位网格, 每个单位网格分别对应一个亮温数据和近地面雨强数据。因此, 亮温数据点与近地面雨强数据点一一对应, 近地面雨强缺省数据为NaN。此操作使用matlab中的griddata函数。

3)Himawari-8资料和2AGPROFGMI资料中都有详细的观测时间数据, 2AGPROFGMI资料与Himawari-8资料的观测时间间隔不超过30s, 可以最大限度地避免降雨强度和结构变化导致的问题。图1(a)为研究区域2016年7月1日01:20 UTC的Himawari-8卫星红外云图(10.4μm通道亮温), 图1(b)是与卫星云图空间匹配的2AGPROFGMI资料的近地面雨强分布, 红线框内为空间匹配区域。对比图1(a)和(b)可以发现, 在空间匹配区域内, 红外云图上亮温分布与2AGPROFGMI资料的近地面雨强分布在空间上有较好的对应关系, 近地面雨强高值区对应于亮温低值区。图2是针对图1空间匹配后的两种资料在同一观测点, 二者的观测时间间隔小于30s的空间匹配图。对比图2(a)和(b), 可以更清晰地看出, 近地面雨强的高值区与红外云图亮温低值区有很好的对应关系。

width=445.8,height=174.1

(a) Himawari-8卫星红外云图(2016年7月1日01:20 UTC), 红线框内为2AGPROFGMI资料的覆盖范围; (b) 2AGPROFGMI资料的近地面雨强分布(红线框内区域)

图1 红外云图及近地面雨强分布图的空间匹配

Fig. 1 Space matching of infrared image and surface precipitation

width=467.75,height=178.55

(a) Himawari-8卫星红外云图; (b) 2AGPROFGMI资料的近地面雨强分布

图2 红外云图与近地面雨强分布图的空间和时间匹配

Fig. 2 Space and time matching of infrared image and surface precipitation

1.3 数据筛选

根据每小时降雨量, 可以划分降雨强度等级: 0.1~2.5 mm/h为小雨, 2.5~8.0 mm/h为中雨, 8.0~16.0 mm/h为大雨, ³16.0mm/h为暴雨[18]。用雨强值0.1 mm/h来区分有雨和无雨。选取2016年6—8月研究区内时间与空间相匹配, 且大雨以上(雨强³8 mm/h)样本数占总样本数3%以上的数据作为研究对象。表1为数据匹配文件, 从中选取相应的匹配数据用于研究。GPMCO卫星环绕地球一周的时间为93分钟, 用于为2AGPROFGMI数据文件命名的时间, 是绕地球一周的开始和结束时刻。

选取匹配文件的日期为 20160603, 20160621, 20160730, 20160812和20160813的数据作为训练数据集来建立反演雨强算法, 文件日期为20160701, 20160704和20160725的数据作为测试数据集来分析雨强反演精度。表2为匹配文件中不同降雨强度的样本数。

2 建立雨强等级反算法

云顶高度、云顶温度以及云的垂直结构等都与降雨有密切关系。不同波长的通道可以反映云的不同方面信息, 因此, 利用不同通道的信息, 可以较好地估计降雨。一般情况下, 降雨云的云体都比较厚, 云顶比较高, 云顶温度比较低, 所以假设冷云顶对应大值雨强。

表1 用于Himawari-8与2AGPROFGMI匹配的文件

Table 1 Matching files of Himawari-8 and 2AGPROFGMI

日期Himawari-8数据2AGPROFGMI数据 20160603NC_H08_20160603_1940_R21_FLDK.02401_02401.ncGPMCOR_GMI_1606031841_2013_012869_L2S_GL2_04A.h5 20160621NC_H08_20160621_0350_R21_FLDK.02401_02401.ncGPMCOR_GMI_1606210317_0450_013139_L2S_GL2_04A.h5 20160701NC_H08_20160701_0120_R21_FLDK.02401_02401.ncGPMCOR_GMI_1607010053_0226_013293_L2S_GL2_04A.h5 20160704NC_H08_20160704_0020_R21_FLDK.02401_02401.ncGPMCOR_GMI_1607032351_0124_013339_L2S_GL2_04A.h5 20160725NC_H08_20160725_0430_R21_FLDK.02401_02401.ncGPMCOR_GMI_1607250329_0501_013668_L2S_GL2_04A.h5 20160730NC_H08_20160730_1640_R21_FLDK.02401_02401.ncGPMCOR_GMI_1607301610_1743_013754_L2S_GL2_04A.h5 20160812NC_H08_20160812_2300_R21_FLDK.02401_02401.ncGPMCOR_GMI_1608122201_2334_013960_L2S_GL2_04A.h5 20160813NC_H08_20160813_2210_R21_FLDK.02401_02401.ncGPMCOR_GMI_1608132110_2242_013975_L2S_GL2_04A.h5

说明: 20160603表示2016年6月3日; NC_H08_20160603_1940_R21_FLDK.02401_02401.nc表示2016年6月3日世界时19:40时刻Himawari-8卫星数据文件; GPMCOR_GMI_1606031841_2013_012869_L2S_GL2_04A.h5表示GPMCO卫星自下传数据以来第12869轨数据文件, 文件名中1841_2013表示开始时刻为18:41, 结束时刻为20:13, 为跨赤道时刻。

表2 匹配文件不同等级降雨强度样本数

Table 2 Sample numbers of different grades of RR in matching files

日期样本数总样本数 无雨有雨小雨中雨大雨以上 2016060311747625740921136102918004 20160621111672446140756047913613 201607016129318523494443929314 201607045682370224396386259384 201607251793462386939525 20160730690510573034552997962 2016081210147780946951508160617956 2016081365973779261268148610376

2.1 生成反演雨强查算表的变量

本文选取研究区域(15—45°N, 90—130°E) 2016年夏季6—8月的数据, 首先找寻与降雨信息有关的通道亮温或亮温差, 并以此作为建立雨强查算表的变量。

2.1.1 单通道亮温与雨强的关系

假设冷云顶对应大值雨强, 将10.4μm通道亮温作为一个变量。由于窗区10.4μm波段相对透明, 云顶上方的水汽对窗区辐射的影响可以忽略不计, 因此可将窗区10.4μm通道的亮温视为云顶温度。利用匹配的亮温数据点与近地面雨强数据点的一一对应关系, 尝试建立亮温与雨强的关系。

图3为由上述匹配数据得到的AHI 10.4 μm通道亮温与对应雨强的散点分布图。可以看到, 尽管从趋势上看, 亮温越低, 雨强越大, 但无法获得二者之间清晰的函数关系。因此, 继续对匹配数据进行统计处理, 得到雨强值间隔1mm/h条件下, 对应AHI 10.4 μm通道亮温的统计箱形图(图4)。从图4可以得到不同雨强间隔的亮温值范围, 也得到较清晰的亮温高值对应雨强低值、亮温低值对应雨强高值的结果。

受从图4得到的结果启发, 利用匹配数据, 获得雨强与1mm/h雨强间隔内平均亮温的对应关系(图5)。图5中, 拟合得到的雨强(RR)与AHI10.4μm通道亮温(BT10.4)的关系曲线方程为

RR = 6.428´108 exp(−0.0845 BT10.4), (1)

由此建立根据单通道亮温计算雨强的函数关系。根据式(1), 可以通过 AHI10.4μm 通道亮温直接计算雨强。

width=221.15,height=158.75

图3 AHI 10.4 μm通道亮温与对应雨强的散点分布

Fig. 3 Scatter plot of AHI 10.4 μm channel brightness temperature and corresponding RR

2.1.2 卷云的影响

虽然卷云的云顶温度低, 但卷云是非降雨云, 因此估算降雨时要去除卷云的影响。Kurino[19]用地球静止气象卫星 GMS-5(Geostationary Meteorologi-cal Satellite)11μm 与 12μm 通道的亮温差(BTD11-12, 表示 BT11 − BT12, 下同)来去除非降雨卷云; Ba等[4]认为当 GOES 卫星 10.7 μm 与 12 μm 通道的亮温差大于 1K 时, 可以判定为卷云。因此, 为排除云顶温度很低的非降雨云——卷云, 本文选择 Himawari卫星 12.4μm 与 10.4μm 通道的亮温差作为查算表的一个变量。

2.1.3 与雨强关系敏感的变量的选择

Kurino[19]利用11μm与6.7μm通道的亮温差来判识能产生强降雨的深对流云; Ba等[4]利用6.9μm与10.7μm通道的亮温差来判识上冲云顶; Mecika-lski等[20]采用Meteosat-9卫星(Meteorological Sate-llite)不同红外通道的亮温差作为云厚度的指标。将这些亮温差作为查算表的变量, 本文测试了大部分红外通道亮温差组合, 结果表明: 这些变量中有对雨强敏感的, 如图6(a)所示; 有对雨强不敏感的, 如图6(b)所示。图6中, 对亮温差做了平均处理, 将在1mm/h雨强间隔内所有点的亮温差求平均, 得到该雨强间隔对应的亮温差。图6(a)和(b)中, 左下角的数据点说明亮温差对极小雨强值不敏感, 这也是影响红外遥感反演降雨精度的因素之一。

最后, 将亮温和亮温差作为变量, 去除与雨强关系不敏感的变量, 与雨强关系相似的变量只保留一个, 利用以下 13 个变量建立查算表: BT10.4, BTD12.4-10.4, BTD6.2-7.3, BTD6.2-10.4, BTD6.2-13.3, BTD7.3-9.6, BTD7.3-10.4, BTD8.6-7.3, BTD 8.6-9.6, BTD8.6-13.3, BTD9.6-10.4, BTD9.6-13.3 和BTD13.3-10.4。

2.2 建立反演雨强查算表

width=442.2,height=215.4

箱型框内横线表示对应的亮温中值, 上下两端表示中值附近一半样本的分布范围, 圆点对应亮温平均值, 加号表示异常数据

图4 雨强对应的AHI 10.4μm通道亮温统计箱形图

Fig. 4 Box plot of brightness temperature statistics of 10.4 μm channel corresponding to RR

width=214.05,height=109.65

圆点为对应于1mm/h雨强间隔内的平均亮温, 实线为圆点对应数据拟合得到的雨强与AHI 10.4 μm通道亮温(BT10.4)的关系曲线

图5 AHI 10.4 μm通道亮温(BT10.4)与雨强的关系

Fig. 5 Relationship between BT10.4 and RR

将匹配的2AGPROFGMI降雨资料的近地面雨强数据作为雨强真值, 将匹配的亮温和亮温差数据点与近地面雨强数据点一一对应, 利用训练集数据,基于BT10.4、其余12个亮温差变量中的一个以及对应的近地面雨强数据, 建立12个二维查算表, 包括以下3个步骤。

1)在训练集中读取数据: 匹配的BT10.4, 其余12个变量中的一个亮温差, 匹配的雨强真值。

2)基于BT10.4以及其余12个变量中的一个亮温差, 建立高分辨率(2K×0.2K)的二维查算表数据网格。

3)在高分辨率的二维查算表网格中, 利用mat-lab中的griddata函数和griddatan函数, 采用邻近点线性插值的方法, 利用匹配的雨强数据赋值, 生成二维查算表。

类似地, 基于BT10.4和BTD12.4-10.4、其余11个亮温差变量中的一个以及对应的近地面雨强数据, 建立11个三维查算表, 包括以下3个步骤。

1)在训练集中读取数据: 匹配的BT10.4和BTD- 12.4-10.4, 其余11个变量中的一个亮温差, 匹配的雨强真值。

2)基于BT10.4和BTD12.4-10.4以及其余11个变量中的一个亮温差, 建立高分辨率(1K×0.1K×0.1 K)的三维查算表数据网格。

3)在高分辨率的三维查算表网格中, 利用mat-lab中的griddata函数和griddatan函数, 采用邻近点线性插值的方法, 利用匹配的雨强数据赋值, 生成三维查算表。

图7显示雨强随亮温和亮温差的变化趋势。其中, 图7(a)和(b)分别是具有代表性的二维散点分布图和二维查算表格点数据分布图, 即通过训练集匹配的雨强随对应的BT10.4和BTD12.4-10.4变化的二维散点图以及生成的二维查算表格点数据分布图; 图7(c)是雨强随对应的BT10.4, BTD12.4-10.4和BTD6.2-10.4变化的三维散点分布图。由于三维查算表图像显示的信息可读性不强, 因此本文没有展示。从散点分布图(图7(a)和(c))可以发现, BTD- 12.4-10.4和BTD6.2-10.4的值越大, BT10.4越小, 对应的雨强越大; 从查算表(图7(b))也可以看到这样的趋势。

建立查算表后, 先通过静止卫星数据, 得到研究区域的BT10.4以及其他12个亮温差变量, 再选取对应的二维或三维查算表, 最后通过查算表找到对应的雨强。

3 反演试验个例

width=409.05,height=161.25

(a) BTD6.2-10.4与雨强的关系; (b) BTD7.3-13.3与雨强的关系

图6 亮温差与雨强的关系

Fig. 6 Relationship between BTD and RR

width=453.6,height=476.25

(a)雨强随BT10.4和BTD12.4-10.4的变化; (b)BT10.4和BTD12.4-10.4二维雨强查算表的格点数据分布; (c)雨强随BT10.4 (横向亮温)、BTD12.4-10.4 (横向亮温差)和BTD6.2-10.4 (纵向亮温差)的变化

图7 雨强与亮温、亮温差的二维和三维关系

Fig. 7 Two-dimensional and three-dimensional relations of RR with brightness temperature and BTD

利用测试数据集来分析雨强反演精度。为了评估反算法的效果, 引入 3 个参数作为判识指标: 检测概率(POD)、错判率(FAR)以及 Heidke 技术得分(Heidke skill score, HSS)[21], 计算公式如下:

width=58.5,height=26.85,

width=58.5,height=28.5,

width=163.35,height=30.1

计算列联表(表3)中相应指数, 可得出各指标参数值。在判断有无降雨时, 表3中A代表反演值和真值都有雨的样本数; B代表反演值无雨而真值有雨的样本数; C代表反演值有雨而真值无雨的样本数; D代表反演值和真值都无雨的样本数; 同理, 在判识不同等级降雨时也用此方法。POD和FAR值在0~1之间变化, HSS值在−1~1之间变化, 理想的情况是POD=1, FAR=0, HSS=1。

表3 雨强真值与反演值列联表

Table 3 Contingency table between observed and estimated RR

雨强反演值雨强真值 1 (是)0 (否) 1 (是)AC 0 (否)BD

图 8(a)~(d) 均采用测试集 2016 年 7 月 1 日样本数据。图 8(a)是 2016 年 7 月 1 日雨强真值分布, 即2AGPROFGMI的近地面雨强, 与图 2(a)红外云图对应; 图 8(b)~(d)分别是利用 BT10.4、二维(BT10.4, BTD12.4-10.4)查算表以及三维(BT10.4, BTD12.4-10.4, BTD6.2-10.4)查算表反演得到的雨强分布。其中, BT10.4 反演是利用式(1)获得的。从图 8 可以看到, 反演的雨强分布范围与真值雨强分布有较好的对应关系, 这种对应关系以三维查算表反演分布为最佳。三维查算表反演的雨强分布优于单变量和二维查算表反演的雨强分布。单变量反演给出较大的雨强分布范围以及范围较大的雨强大值区, 因此, 可以排除利用BT10.4进行雨强反演的方案。

现在对比和讨论二维查算表和三维查算表的反演结果。图 9 是利用测试集 2016 年 7 月 1 日样本数据, 采用两个变量(BT10.4, BTD12.4-10.4)和 3 个变量(BT10.4, BTD12.4-10.4, BTD6.2-10.4)对应查算表反演得到的雨强与雨强真值对应的散点分布图, 真值与反演值的相关系数(R)分别为 0.4750 和 0.4675, 均方根误差(RMS)分别为 3.4014 和 3.3652mm/h。可以看到, 变量个数越多的算法, 得到的结果均方根误差越小。

width=439.3,height=382.65

(a)雨强真值分布(2AGPROFGMI的近地面雨强); (b)BT10.4反演的雨强分布; (c)二维(BT10.4, BTD12.4-10.4)查算表反演的雨强分布; (d)三维(BT10.4, BTD12.4-10.4, BTD6.2-10.4)查算表反演的雨强分布

图8 不同组合的变量反演的雨强分布与雨强真值分布的比较

Fig. 8 Comparisons of actual distribution of RR with the distributions retrieved by variables of different combinations

width=459.2,height=164.4

(a)两个变量(BT10.4, BTD12.4-10.4); (b) 3个变量(BT10.4, BTD12.4-10.4, BTD6.2-10.4)

图9 不同反算法得到的雨强反演值与雨强真值对应的散点分布(2016年7月1日)

Fig. 9 Scatter plots of observed and estimated RR base on different retrieval algorithms on July 1, 2016

表 4 给出利用测试集 2016 年 7 月 1 日样本数据,采用二维查算表和三维查算表反算法, 对不同等级降雨的技术评分。可以看到, 变量个数越多的算法, 技术评分表现越好。

下面只讨论三维查算表算法的反演结果。图10(a)和(b)分别是依据测试集 2016 年 7 月 4 日的样本数据以及测试集全部样本数据, 均采用BT10.4, BTD12.4-10.4, BTD6.2-10.4对应的三维查算表, 反演得到的雨强值与真值对应的散点分布图, 真值与反演值相关系数(R)分别为0.5258和0.5018, 均方根误差(RMS)分别为3.8508和3.5808mm/h。图10显示, 在利用部分和全部测试集样本数据这两种情况下, 三维查算表算法都得到比较合理的反演精度。

现在进一步讨论三维查算表算法的技术评分。表5是分别依据测试集2016年7月4日样本数据和测试集全部样本数据, 采用三维(BT10.4, BTD12.4-10.4, BTD6.2-10.4)查算表, 对不同等级降雨反演结果的技术评分统计结果。可以看到, 该三维查算表对有雨事件有较高的监测概率(最高约为87%)和较低的错判率(最低约为19%), 能较好地判识降雨区域, 同时也能在一定程度上对降雨等级进行判识。

4 结果对比和误差讨论

虽然2.2节建立12个二维查算表和11个三维查算表, 但是测试结果显示, 变量个数越多(与降雨有关的通道信息增多), 算法的反演效果越好。正是由于通道信息少, 导致二维查算表对大雨以上降雨等级反演的技术评分普遍表现不好。总体来说, 三维查算表算法表现较好。因此, 依据测试集全部样本数据, 只对这些三维查算表算法进行统计讨论。表6列出对不同等级降雨反演技术评分表现好的两个三维查算表(BT10.4, BTD12.4-10.4, BTD9.6-10.4和 BT10.4, BTD12.4-10.4, BTD6.2-7.3)反演结果的一些统计参数。

表4 2016年7月1日不同反算法技术评分比较

Table 4 Comparison of different retrieval algorithms with skill scores on July 1, 2016

雨强/(mm·h−1)二维查算表(BT10.4, BTD12.4-10.4)三维查算表(BT10.4, BTD12.4-10.4, BTD6.2-10.4) PODFARHSSPODFARHSS <0.1 (无雨)0.62750.16770.35450.62050.15620.3560 >0.1 (有雨)0.76460.47560.35450.77900.48390.3560 0.1~2.5 (小雨)0.50600.63480.17790.53710.64390.1794 0.1~8.0 (小到中雨)0.68430.55400.26680.69960.56120.2703 >2.5 (中雨以上)0.56820.65820.35290.57010.62310.3876

width=453.6,height=170

(a)依据测试集2016年7月4日样本数据; (b)依据测试集全部样本数据

图10 三维(BT10.4, BTD12.4-10.4, BTD6.2-10.4)查算表反演雨强值与雨强真值对应的散点分布

Fig. 10 Scatter plots of observed and estimated RR base on three-dimensional (BT10.4, BTD12.4-10.4, BTD6.2-10.4) lookup table

表5 三维(BT10.4, BTD12.4-10.4, BTD6.2-10.4)查算表技术评分

Table 5 Skill scores of three-dimensional (BT10.4, BTD12.4-10.4, BTD6.2-10.4) lookup table

雨强/(mm·h−1)2016年7月4日样本数据测试集全部样本数据 PODFARHSSPODFARHSS <0.1 (无雨)0.48720.24350.37000.57840.18860.3882 >0.1 (有雨)0.87110.32620.37000.82620.39760.3882 0.1~2.5 (小雨)0.61850.53020.20880.58160.58310.2050 0.1~8 (小到中雨)0.76070.44820.23500.73370.49720.2770 >2.5 (中雨以上)0.69570.44180.51960.62960.52210.4621

表6 技术评分表现好的三维查算表

Table 6 Three-dimensional lookup table with good skill scores

查算表雨强/(mm·h−1)PODFARHSSRRMS/(mm·h−1) BT10.4, BTD12.4-10.4, BTD9.6-10.4<0.1 (无雨)0.69390.16060.45150.58633.4523 BT10.4, BTD12.4-10.4, BTD6.2-7.3>0.1 (有雨)0.88170.40420.41490.48833.6780 BT10.4, BTD12.4-10.4, BTD6.2-7.30.1~2.5 (小雨)0.62290.58090.22970.48833.6780 BT10.4, BTD12.4-10.4, BTD6.2-7.30.1~8.0 (小到中雨)0.78890.49960.30380.48833.6780 BT10.4, BTD12.4-10.4, BTD9.6-10.4>2.5 (中雨以上)0.72810.53820.49360.58633.4523 BT10.4, BTD12.4-10.4, BTD9.6-10.4>8.0 (大雨以上)0.62680.54700.49340.58633.4523

从表6可以发现, 包含BT10.4, BTD12.4-10.4, BTD6.2-7.3的三维查算表能较好地判识有雨(POD =0.8817, FAR=0.4042)、小雨(POD=0.6229, FAR= 0.5809)和小雨到中雨(POD=0.7889, FAR=0.4996)的降雨, 包含BT10.4, BTD12.4-10.4, BTD9.6-10.4的三维查算表能较好地判识无雨(POD=0.6939, FAR= 0.1606)、中雨以上(POD=0.7281, FAR=0.5382)和大雨以上(POD=0.6268, FAR=0.5470)的降雨。

Krawczyk等[22]发现, BTD6.2-7.3值约为零时, 云发展旺盛, 云顶可到达对流层顶。Elmer等[23]发现, BTD6.2-7.3和BTD9.6-10.4分别指示垂直水汽含量和云的厚度。表6中评分表现好的三维查算表的变量就包含BTD6.2-7.3和BTD9.6-10.4, 说明这两个亮温差包含较丰富的降雨信息。包含BTD9.6-10.4的查算表对大雨以上的检测概率明显提高, 说明BTD9.6-10.4非常有助于检测大降雨事件。因此, 本文根据统计结果得到表现好的三维查算表, 选择的变量很好地证实了文献[22‒23]的结果, 也说明了变量选取的合理性。

利用红外图像数据建立雨强查算表, 反演降雨的误差主要来源于以下几方面。

1)红外辐射只能看到云表面的信息, 不能穿透云、进而探测云的内部信息, 是间接估计降雨。

2)降雨与云的亮度、面积和种类有关, 它们之间的关系比较复杂。假设冷云顶对应比较大的降雨会对估计降雨带来一定的误差, 实际上有的云顶温度很低, 但对应的云降雨量很小, 甚至没有降雨。

3)将高分辨率的Himawari-8资料与低分辨率的GMI资料匹配, 将可能压低估计降雨的高值, 抬高估计降雨的低值。

4)数据的时间匹配存在误差, 红外图像和近地面雨强资料分别由两个不同的卫星平台得到, 对不同观测平台资料进行时间匹配是困难的。

5)降雨过程往往是一个天气系统造成的, 降雨区域与云的位置往往不完全一致。采用数据点一一对应的方式建立红外图像数据与近地面雨强的关系, 会造成一定的误差。

5 结论和展望

本文基于2016年6—8月Himawari-8静止卫星多通道红外图像数据和GPM降雨资料2AGPROFGMI的近地面雨强数据, 寻找红外亮温、亮温差与雨强之间的关系, 建立查算表来反演雨强。通过对算法的反演试验, 证明三维雨强查算表可以较好地判识有雨与无雨, 并且能够对不同等级的降雨有较好的识别。其中, 三维查算表(BT10.4, BTD12.4-10.4, BTD6.2-7.3)对降雨事件表现出较高的检测概率(POD=0.8817)和较低的错判率(FAR=0.4042), 能较好地判识降雨区域。基于该算法, 利用静止卫星时间分辨率高的优点, 结合地面常规观测资料, 可以连续地监测和跟踪比较大的降雨, 对防灾减灾起到很好的作用。

本文采用查算表直接估算雨强, 倘若先进行降雨检测(先判断有无雨), 再对有雨的事件进行反演, 反演精度可能会在一定程度上有所提高。云的特征与降雨密切相关, 分析降雨云的演变和纹理特征等, 对反演降雨会有很大的帮助。本文主要利用Himawari-8通道多的优点, 尝试用多种亮温差寻找新的降雨信息, 有关其他可能影响降雨的因素(如降雨云的形态、云顶纹理特征以及云发展和演变特征等), 将在后续的工作中探讨。

本文只采用一种降雨资料作为雨强真值, 如果将其他资料(如地面天气雷达测雨资料等)也作为真值, 可以进行更多的反演结果检验, 提高反演精度的可信度。另外, 我国新一代静止气象卫星FY-4A已经发射运行, 在其搭载的多通道扫描成像辐射计(Advanced Geosynchronous Radiation Imager, AGRI)上, 有类似Himawari-8静止卫星的红外通道, 本文根据Himawari-8通道亮温建立的降雨反算法, 今后可以应用到FY-4A卫星降雨反演研究中。

参考文献

[1]Behrangi A, Andreadis K, Fisher J B, et al.Satellite-based precipitation estimation and its application for streamflow prediction over mountainous western U.S. basins. Journal of Applied Meteorology & Climato-logy,2014, 53(12): 2823‒2842

[2]Prigent C. Precipitation retrieval from space: an overview. Comptes rendus Geoscience, 2010, 342(4): 380‒389

[3]刘元波, 傅巧妮, 宋平, 等. 卫星遥感反演降水研究综述. 地球科学进展, 2011, 26(11): 1162‒1172

[4]Ba M B, Gruber A. GOES multispectral rainfall al-gorithm (GMSRA). Journal of Applied Meteorology, 2001, 40(8): 1500‒1514

[5]Adler R F, Negri A J. A satellite infrared technique to estimate tropical convective and stratiform rainfall. Journal of Applied Meteorology, 1988, 27(1): 30‒51

[6]Rosenfeld D, Gutman G. Retrieving microphysical properties near the tops of potential rain clouds by multispectral analysis of AVHRR data. Atmospheric Research, 1994, 34: 259‒283

[7]Behrangi A, Hsu K L, Lmam B, et al. PERSIANN-MSA: a precipitation estimation method from satellite- based multispectral analysis. Journal of Hydrometeo-rology, 2009, 10(6): 1414‒1429

[8]Schmit T J, Gunshor M M, Menzel W P, et al. Introducing the next-generation advanced baseline imager on GOES-R. Bulletin of the American Meteo-rological Society, 2005, 86(8): 1079‒1096

[9]Tao Y, Hsu K, Gao X, et al. A two-stage deep neural network framework for precipitation estimation from bispectral satellite information.Journal of Hydrome-teorology, 2018, 19(2): 393‒408

[10]Zhuge X, Yu F, Zhang C. Rainfall retrieval and nowcasting based on multispectral satellite images. Part I: retrieval study on daytime 10-minute rain rate. Journal of Hydrometeorology, 2011, 12(6): 1255‒ 1270

[11]金晓龙, 邵华, 张弛, 等. GPM卫星降水数据在天山山区的适用性分析. 自然资源学报, 2016, 31(12): 2074‒2085

[12]Bessho K, Date K, Hayashi M, et al. An introduction to Himawari-8/9 — Japan’s new-generation geosta-tionary meteorological satellites. Journal of the Me-teorological Society of Japan, 2016, 94(2): 151‒183

[13]Hou A Y, Kakar R K, Neeck S, et al. The global precipitation measurement mission. Bulletin of the American Meteorological Society, 2014, 95(5): 701‒ 722

[14]唐国强, 万玮, 曾子悦, 等. 全球降水测量(GPM)计划及其最新进展综述. 遥感技术与应用, 2015, 30 (4): 607‒615

[15]Kummerow C D, Randel D L, Kulie M, et al. The evolution of the Goddard profiling algorithm to a fully parametric scheme. Journal of Atmospheric & Oceanic Technology, 2015, 32(12): 113‒130

[16]Kummerow C D, Hong Y, Olson W S, et al.The evolution of the Goddard profiling algorithm (GPROF) for rainfall estimation from passive micro-wave sensors. Journal of Applied Meteorology, 2001, 40(11): 1801‒1820

[17]Petkovic V, Kummerow C D. Performance of the GPM passive microwave retrieval in the Balkan flood event of 2014. Journal of Hydrometeorology, 2015, 16(6): 2501‒2518

[18]林晔, 王庆安, 顾松山, 等. 大气探测学教程. 北京: 气象出版社, 1993: 67

[19]Kurino T. A satellite infrared technique for estimating “deep/shallow” precipitation. Advances in Space Re-search, 1997, 19(3): 511‒514

[20]Mecikalski J R, Mackenzie W M J, Koenig M, et al. Cloud-top properties of growing cumulus prior to convective initiation as measured by meteosat second generation. Part I: infrared fields. Journal of Applied Meteorology & Climatology, 2010, 49(12): 521‒534

[21]Conner M D, Petty G W. Validation and intercom-parison of SSM/I rain-rate retrieval methods over the continental United States. Journal of Applied Meteo-rology, 1998, 37(7): 679‒700

[22]Krawczyk K, Jasinski J. Multispectral satellite data application to hazardous convection monitoring // In-ternational Conference on Environmental Engineering (ICEE) Selected papers. Vilnius: Vilnius Gediminas Technical University Press Technika, 2014: 1‒9

[23]Elmer N J, Berndt E, Jedlovec G J. Limb correction of MODIS and VIIRS infrared channels for the improved interpretation of RGB composites. Journal of Atmos-pheric & Oceanic Technology, 2016, 33(5): 1073‒ 1087

Retrieval of Precipitation by Using Himawari-8 Infrared Images

SUN Shaohui1,2, LI Wanbiao1,†, HUANG Yipeng1

1. Department of Atmospheric and Oceanic Sciences, School of Physics, Peking University, Beijing 100871; 2. 93313 Troop, Changchun 130111; † Corresponding author, E-mail: lwb@pku.edu.cn

Abstract Using the matched data between infrared images of the geostationary meteorological satellite Himawari-8 and the product 2AGPROFGMI of Global Precipitation Measurement (GPM), the matched infrared brightness temperature (BT), brightness temperature difference (BTD) and surface precipitation are connected by pixel to pixel. Further more, because of the advantage of the more channels of geostationary meteorological satellite, two-dimensional and three-dimensional lookup tables of rain rate (RR) are established with the matched infrared brightness temperature, brightness temperature difference and surface precipitation. The lookup tables which can identify the rain rate of different grades are found by the retrieval tests. The three-dimensional lookup tables (BT10.4, BTD12.4-10.4, BTD6.2-7.3) show higher probability of detection (POD=0.8817) and lower false alarm rate (FAR=0.4042) when detecting the rain, so it can identify the areas of rainfall.

Key words satellite remote sensing; brightness temperature; brightness temperature difference; rain rate; retrieval

doi: 10.13209/j.0479-8023.2018.101

收稿日期: 2018-03-16;

修回日期: 2018-05-18;

网络出版日期: 2018-12-27

国家科技支撑计划(2013CB430104)资助