张路,胡潭高,张毅,张登荣,李瑶,沈黎达
( 1. 杭州师范大学 遥感与地球科学研究院,浙江 杭州 311121;2. 浙江省城市湿地与区域变化研究重点实验室,浙江 杭州311121;3. 广西壮族自治区环境保护科学研究院,广西 南宁 530022;4. 自然资源部空间海洋遥感与应用研究重点实验室,北京 100081)
台风是指发生在热带西北太平洋的具有气旋结构的大型天气系统,我国是全球遭受台风破坏影响最显著的国家之一,每年台风在我国造成的经济损失都高达上百亿元[1-2]。台风风圈是描述台风的一个重要参数,可以在一定程度上表达台风的结构特征,决定一个台风潜在的破坏力和可能影响范围[3]。在台风预报中,台风风圈大小对模型计算、强度、路径的预测结果具有重要意义[4-5],有助于沿海地区采取有针对性的防台措施[6]。另外,风圈对于各种基于热带气旋的应用也至关重要[7-8],美国国家飓风中心主要利用风圈大小来决定警报的强度,并且预估台风的破坏性。台风风圈是对台风具体影响范围与程度的直观体现,当前台风风圈的确定方法主要依赖于台风水平方向的风剖面信息。台风水平方向风剖面信息可以以风速-半径关系曲线的方式展现,即台风风剖面曲线[9]。台风风剖面曲线是一条台风结构内与台风中心不同距离的各点平均风速变化的曲线,因此通过提取风剖面信息,能够快速直观地展示不同等级风圈的半径。
风圈研究一直以来都是台风相关研究中的热点之一,受到国内外许多学者的关注。早期的风圈研究主要基于海岸气象站、船舶、浮标等的地面/海面实测数据。但是,由于测量站点记录的实测数据多为点源数据,缺乏对台风的整体观测能力,因此现今的风圈研究中地面/海面站点观测数据多作为辅助数据使用[10-11]。自20 世纪60 年代开始,随着航空技术的发展,基于航空平台的对热带气旋观测数据开始增多。航空观测的方法主要为通过飞机等航空器飞越台风上空或内部云系结构通过遥感手段对台风进行观测,航空遥感能够对台风外层与内层进行观测,获取高精度的台风参数[12-13]。但由于其危险性较高,成本昂贵[14],因此在非必要的时候也少有进行对台风的航空观测。20 世纪70 年代后,随着卫星遥感技术的发展,卫星遥感监测迅速成为台风监测研究的主要手段,它具备了对台风整体结构进行观测的能力以及时间分辨率相对较高的优点[15]。
近年来,越来越多的国内外学者开始尝试利用多源卫星遥感数据进行台风的相关研究。冯倩[16]进行了搭载多传感器卫星对海面风场的遥感研究;杨亮[17]研究了基于遥感数据的西北太平洋海面风场时空特征;窦芳丽等[18]基于FY-3 双频散射计数据反演了台风降水区风场;吴俞等[19]基于卫星资料研究了台风“海燕”对海南岛风雨分布的影响。上述研究说明了基于遥感数据研究台风信息是一种可行且主流的方式,因此相关学者进行了多种基于卫星散射计风场产品的台风风圈研究,如陈德文等[9]基于QuikSCAT 散射计风场数据对台风最大风速半径反演及个例进行了研究;Park 等[20]基于微波散射计数据,韩国气象局(KMA)与东京台风中心(TTC)的记录数据对西北太平洋热带气旋风圈半径进行了对比分析;Chavas 和Emanuel[21]利用QuikSCAT 散射计数据进行数学建模,分析了多个海域的热带气旋的结构与风圈半径;Said 和Long[22]将QuikSCAT 风力产品进行空间分辨率增强,估算了热带气旋的特征参数,并进行了精度评价与误差分析。其中,陈德文等[9]与Said 和Long[22]都利用了Holland 风场模型对风圈进行了反演评估,Holland 风场模型是一个台风气压-风速数值模型,是目前主流的台风风场模型之一。
海洋二号A 星(HY-2A)是我国第一颗海洋动力环境监测卫星,主要搭载了散射计、辐射计和高度计3 种传感器。其中,散射计能够提供探测海域的海面风速、海面风向等产品,已有相关学者利用散射计风场资料开展了台风监测研究,如刘晓燕等[23]基于HY-2A 散射计风场资料在台风“菲特”预报中的应用;赵勇等[24]基于HY-2A 散射计风场数据研究了台风“苏力”的海表面风场结构;兰友国等[25]综述了HY-2A 微波散射计在台风遥感监测中的应用;邹巨洪等[26]基于HY-2A 风矢量产品进行了热带气旋自动识别的研究;胡潭高等[27]通过将HY-2 散射计数据与FY-2 微波数据结合,提出了一种有效提高台风监测能力的方法;杨典等[28]基于HY-2A 散射计数据提出了一种台风强度的诊断方法。
本文将尝试基于HY-2A 散射计提供的海面风场资料,并结合Holland 风场模型,研究台风风圈信息的提取方法。首先,通过分析Holland 风场模型中的两种台风形状系数对风剖面信息提取精度的影响,提出一种基于微波散射计产品的较为稳定可靠的台风风剖面信息提取方法;然后,在此基础上确定不同等级的风圈范围;最后,通过多期台风应用实例,验证风圈提取精度。
2.1.1 HY-2A 散射计海面风场数据
本研究选用的散射计风场数据是HY-2A 卫星地面数据处理系统业务化处理生产的微波散射计L2B级产品,采用了最大似然估计及圆中数滤波方法,数据版本为业务处理版本,采用基于二维变分分析方法得到风场产品精度有所提高,但仍未对业务化处理版本的数据进行批量重处理,故在本研究中采用[29]。微波散射计风场产品文件中,记录了各观测点的经度、纬度、监测时间、风速、风向等海洋动力环境参数。产品记录了观测点海面风场情况,产品风速精度为2 m/s(风速≤24 m/s)与风速的10%(风速>24 m/s),最大记录风速值为38.5 m/s;记录风向值为0°~360°,精度为±20°,每一景数据的刈幅宽度约为1 350 km[30]。产品的空间分辨率为25 km,时间分辨率为12 h。
2.1.2 最佳路径数据集
美国联合台风预警中心(Joint Typhoon Warning Center,JTWC)记录了从20 世纪50 年代至今的基本所有全球洋面上热带气旋的各项主要参数,根据海域划分整合了不同海域上热带气旋的最佳路径数据集。最佳路径数据集提供了时间分辨率为6 h 的各热带气旋的气旋中心经纬度、近中心最大风速、中心最低气压、各级风圈半径等主要的热带气旋科学参数。
本研究基于HY-2A 卫星的散射计数据,首先对散射计海面风场数据进行预处理,获得散射计风场图,通过风场图提取台风主要参数;然后,利用散射计风场资料,改进基于Holland 风场模型,反演风速剖面信息;最后,以2012-2017 年间16 次台风为例,进行精度验证。总体技术路线图1 所示。
2.2.1 散射计风场资料预处理
从国家卫星海洋应用中心网站下载覆盖对应热带气旋的HY-2A 散射计L2B 级风场产品,确认对应时刻热带气旋的中心后截取距离气旋中心经纬度±5°的矩形区域内的所有记录点的风速值与风向值,并基于记录值对台风结构区及其周边区域进行插值重采样处理,生成对应区域内的风场图。基于生成的风场图,提取该景散射计L2B 级风场产品记录的热带气旋中心经纬度、最大风速、最大风速半径(Radius of Maximum Wind,RMW)3 个热带气旋结构参数[31-32]。提取方法如下:台风中心为风场图中,以低风速值区域中的风向涡旋中心为台风中心,提取台风中心经纬度;以台风中心为原点,半径100 km 内风速最大插值结果为台风最大风速;最大风速像元与台风中心的距离作为台风最大风速半径。
2.2.2 基于散射计和Holland 风场模型的风速剖面反演
Holland[33]基于Schloemer 指数气压分布模型,于1980 年提出一种描述台风内气压径向分布的模型
式中,r为与台风中心的距离;Pr为距离台风中心r处的气压;Pc为台风中心气压;Pn为台风外围气压;RMW为台风最大风速半径;B为台风形状系数。
由于散射计没有气压数据,需要利用其他相关参数,而台风的气压参数与风速参数存在联系,因此陈德文等[9]在Holland 气压模型的基础上,通过风速与气压的拟合得到与台风最大风速相关的台风风速剖面模型
式中,Vr为距离台风中心r处的风速;Vmax为台风最大风速;r为与台风中心的距离;RMW为台风最大风速半径;B为台风形状系数。
Holland 气压模型与Schloemer 气压模型最大的改进意义在于引入了HollandB参数,其与台风最大风速半径均为Holland 模型中重要参数,B的合理取值范围为1.0~2.5。Willoughby 和Rahn[34]及杨万康等[35]通过航空探测数据对台风模型中的各参数进行拟合,获得了对参数B的拟合计算公式
图 1 总体技术路线Fig. 1 Overall technical roadmap
式中,Vmax为台风最大风速;RMW为台风最大风速半径;L为台风中心纬度。
Vickery 等[35]及贾永君等[36]将航空探测数据与HRD 观测数据进行结合,对相关台风特征参数进行了统计分析,获得了台风形状系数B的计算公式
式中,RMW为台风最大风速半径;L为台风中心纬度。
利用式(3)与式(4),可以将计算出的B值分别代入式(2)进行台风剖面风速曲线拟合,两种B值的计算方法(下文称为WilloughbyB值计算法与VickeryB值计算法)对台风剖面风速曲线造成重要影响,不同的B值会使每一景散射计数据反演出2 条不一样的台风剖面风速曲线。
2.2.3 模型评价与验证
(1)数据匹配
JTWC 发布的台风最佳路径数据集中,记录点的时间坐标为台风活动期间内每日的0 时、6 时、12 时与18 时。散射计风场数据产品中,西北太平洋数据产品的时间坐标为每日的7-11 时与19-23 时。由于JTWC 记录台风的时间坐标与散射计风场数据记录台风的时间坐标不一致,不能直接使用JTWC 的记录参数作为参考真值,因此本文使用三次样条插值法对JTWC 记录值进行插值,获取与散射计风场数据记录时间相同的JTWC 最佳路径数据集台风参数。
在区间[a,b]上,以xi(i=1,2,3,···,n)为分段点的函数F(x),其中a=x1<x2<x3<···<xn=b,其三次样条函数S(xi)可表示为
①通过所有分段点xi,即满足
②存在二阶连续性
即
事实上对于每一个区间内都为三次多项式。
(2)精度评价
通过式(3)与式(4)分别计算每一景散射计海面风场数据的两个台风形状系数,然后利用Holland 风场模型提取对应的风圈半径,以JTWC 的最佳路径数据集风圈插值结果作为参考真值,通过计算Holland模型风圈半径结果与JTWC 风圈插值结果两者之间的均方根误差与平均绝对误差,进行精度评价,确定一种精度较高、适宜反演西北太平洋热带气旋风圈的台风形状系数B计算方法。均方根误差(RMSE)的计算公式如下
平均绝对误差(MAE)的计算公式如下
本文选择2012 年02 号台风“珊瑚”与2013 年23 号台风“菲特”这两期典型的台风作为应用实例(图2),以验证本文方法的精度,选取原因说明如下:(1)台风“珊瑚”的活动时间为5 月底,属于生成时间较早的热带气旋,而台风“菲特”活动时间为9 月底10 月初,属于生成时间为中后期的热带气旋,两个样本在时间范围上相对具有较好的代表性;(2)台风“珊瑚”属于远洋型热带气旋,即整个生命周期内都在远洋海域活动,而台风“菲特”属于登陆型热带气旋,生成后朝着陆地方向移动并登陆,两个样本的空间范围与活动性质不同,相对具有较好的代表性;(3)两个台风样本的最大风速值在HY-2A 散射计有效记录最大风速值38.5 m/s 以下[30,37]的时间相对较长,有效样本数量较多。另外,为了验证本文方法的适用性,还选择了2012-2017 年间,对满足散射计海面风场数据质量较好,对台风结构有完整覆盖,或覆盖了台风大部结构区域的14 个台风进行风圈提取验证。
图 2 台风“珊瑚”与“菲特”路径Fig. 2 Typhoon Sanvu and Typhoon Fitow road map
台风“珊瑚”于2012 年5 月21 日02 时生成于关岛东南方附近洋面,生成后朝西北方向移动,强度逐渐加强,移动方向逐渐朝北偏移。5 月24 日于21°N,139°E 附近加强为台风,并转为朝东北方向移动,5 月26 日减弱为热带风暴并继续朝东北方向移动,强度逐渐减弱,5 月28 日变性为温带气旋后于当日20 时于日本以东洋面上消散。
在台风“珊瑚”的完整生命周期内,我们选取6 景能够较好反映台风结构的微波散射计风场资料用于风剖面信息提取,具体为北京时间2012 年5 月23 日09 时、5 月23 日21 时、5 月24 日09 时、5 月24 日21 时,5 月25 日09 时,5 月26 日20 时。上述6 景数据均为台风“菲特”距离海岛或陆地较远的海域上,散射计数据受海岛陆地影响较小。从图3 可以看出5 月23 日09 时“珊瑚”已经具有较为明显的环流结构,形状结构也趋于对称,具备较好的拟合条件。
图 3 台风“珊瑚”风场Fig. 3 Wind field map of Typhoon Sanvu
对6 次微波散射计数据进行可视化,生成每一景散射计数据对应的风场图。基于风场图提取每一景散射计数据记录时,台风“珊瑚”的中心经纬度、最大风速、最大风速半径4 个重要台风结构参数(表1),用于Holland 模型进行风速剖面曲线拟合。
分别对6 景散射计数据进行模型运算,生成对应的风速剖面拟合曲线(图4,图5)。从表1 得知,除了5 月23 日09 时散射计提取的最大风速与JTWC 记录的最大风速有3.5 m/s 的较大误差外,6 景散射计数据提取的台风经纬度、最大风速、最大风速半径4 个参数与JTWC 记录参数的偏差较小。表2 表述了6 景散射计数据反演得出的风圈值与JTWC 记录值,在中心风力较低的情况下,即5 月23 日09 时与5 月26日20 时2 景中,基于VickeryB值计算法的风圈反演精度优于基于WilloughbyB值计算法的风圈半径反演精度,另外4 景时“珊瑚”的中心风力较大,基于WilloughbyB值计算法的风圈反演精度则优于VickeryB值计算法。
表 1 台风“珊瑚”结构参数对比Table 1 Comparison of structural parameters of Typhoon Sanvu
从图4 与图5 可以看出对于台风“珊瑚”,WilloughbyB值计算法的风速剖面拟合曲线对比VickeryB值计算法更为尖锐,即在最大风速半径外,风速随着半径增加而下降的速率更高。同一景数据中,最大风速半径外区域的同距离点,VickeryB值计算法对应的风速拟合值均比WilloughbyB值计算法更高。6 景WilloughbyB值计算法的风速拟合曲线在距中心远距离处均趋于同一值,而VickeryB值计算法的拟合曲线之间则趋于固定的差值。对两种方法反演拟合的风圈值与JTWC 最佳路径数据集记录的风圈值进行比较,台风“珊瑚”34 kt 风圈中,基于WilloughbyB值计算法的风圈半径RMSE 为23.6 km,MAE 为19.3 km;基于VickeryB值计算法风圈半径的RMSE 为36.7 km,MAE 为28.3 km。台风“珊瑚”50 kt风圈中,WilloughbyB值计算法的风圈半径RMSE 为18.0 km,MAE 为15.3 km;基于VickeryB值计算法风圈半径的RMSE 为27.4 km,MAE 为22.5 km。
图 4 台风“珊瑚”Willoughby B 值计算法风速剖面拟合曲线Fig. 4 Wind speed profile fitting curve of Typhoon Sanvu with value B calculated by Willoughby's method
图 5 台风“珊瑚”Vickery B 值计算法风速剖面拟合曲线Fig. 5 Wind speed profile fitting curve of Typhoon Sanvu with value B calculated by Vickery's method
台风“菲特”于2013 年9 月29 日20 时生成于帕劳以北洋面,后朝西北偏北方向移动,强度逐渐加强,10 月04 日于冲绳以南洋面加强为强台风,后转向西北偏西方向移动,穿过钓鱼岛后10 月07 日凌晨于福建省福鼎市登陆,登陆时最大风速为42 m/s,登陆后强度迅速减弱,同日下午于福建省南平市以西附近消散。
在台风“菲特”的完整生命周期内,我们选取6 景能够较好反映台风结构的微波散射计风场资料用于风剖面信息提取,具体为北京时间2013 年10 月01 日09 时、10 月01 日21 时、10 月02 日09 时、10 月02 日21 时,10 月03 日09 时,10 月03 日21 时。上述6 景数据均为台风“菲特”距离海岛或陆地较远的海域上,散射计数据受海岛陆地影响较小。从图6 可以看出10 月1 日时“菲特”已经具有较为明显的环流结构,形状结构也趋于对称,具备较好的拟合条件。
对6 次微波散射计数据进行可视化,生成每一景散射计数据对应的风场图。基于风场图提取每一景散射计数据记录时,台风“菲特”的中心经纬度、最大风速、最大风速半径4 个重要台风结构参数(表3),用于Holland 模型进行风速剖面曲线拟合。
分别对6 景散射计数据进行模型运算,生成对应的风速剖面拟合曲线(图7,图8)。从表3 得知,10 月01 日09 时至10 月03 日09 时的5 景散 射计数据提取的台风经纬度、最大风速、最大风速半径4 个参数与JTWC 记录参数误差较小,而从10 月03 日21 时散射计数据提取的台风参数与JTWC 记录值误差较大;表4 为6 景散射计数据反演得出的风圈值与JTWC 记录值,可以看出10 月01 日09 时至10 月03 日09 时这5 景数据中,基于WilloughbyB值计算法的风圈反演精度优于VickeryB值计算法,而10 月03 日21 时这一景数据基于VickeryB值计算法反演精度优于WilloughbyB值计算法。
图 6 台风“菲特”风场图Fig. 6 Wind field map of Typhoon Fitow
表 3 台风“菲特”结构参数对比Table 3 Comparison of structural parameters of Typhoon Fitow
从图7 与图8 可以看出同台风“珊瑚”一样,对于台风“菲特”,基于WilloughbyB值计算法的风速剖面拟合曲线对比VickeryB值计算法也更为尖锐,对于同一景数据中最大风速半径外区域的同距离点,基于VickeryB值计算法的风速拟合值均比WilloughbyB值计算法更高。6 景基于WilloughbyB值计算法的风速拟合曲线在距中心远距离处也趋于同一值,而基于VickeryB值计算法的拟合曲线之间也趋于固定的差值。对两种方法反演拟合的风圈值与JTWC 最佳路径数据集记录的风圈值进行比较,台风“菲特”34 kt风圈中,基于WilloughbyB值计算法风圈半径的RMSE 为46.7 km,MAE 为38.6 km;基于VickeryB值计算法风圈半径的RMSE 为65.6 km,MAE 为52.0 km。台风“菲特”50 kt 风圈中,基于WilloughbyB值计算法风圈半径的RMSE 为21.4 km,MAE 为17.0 km;基于Vickery B 值计算法风圈半径的RMSE 为31.8 km,MAE 为24.3 km。
为了进一步评价基于HY-2A 微波散射计风场产品在多个西太平洋台风风圈提取中的适用性与精度,本文还对2012-2017 年间共14 个台风进行了风圈提取与精度评价。14 个台风的RMSE 与MAE 分别如表5 和表6 所示。
由表5 和表6 可知,对于上述14 个台风的34 kt风圈中,基于WilloughbyB值计算法风圈半径的平均均方根误差为25.6 km,平均偏差为22.2 km;基于VickeryB值计算法风圈半径的均方根误差为48.0 km,平均偏差为42.1 km。50 kt 风圈中,基于WilloughbyB值计算法风圈半径的均方根误差为13.3 km,平均偏差为11.4 km;基于VickeryB值计算法风圈半径的均方根误差为21.6 km,平均偏差为18.9 km。
图 7 台风“菲特”Willoughby B 值计算法风速剖面拟合曲线Fig. 7 Wind speed profile fitting curve of Typhoon Fitow with value B calculated by Willoughby's method
图 8 台风“菲特”Vickery B 值计算法风速剖面拟合曲线Fig. 8 Wind speed profile fitting curve of Typhoon Fitow with value B calculated by Vickery's method
表 4 台风“菲特”风圈对比Table 4 Comparison of wind radii of Typhoon Fitow
表 5 2012-2017 年间14 个台风风圈的RMSETable 5 RMSE of 14 typhoons’ wind radii between 2012 and 2017
以上应用实例的结果表明:对比WilloughbyB值计算法和与VickeryB值计算法对应的两种风速剖面拟合曲线,可以看出WilloughbyB值计算法的拟合曲线更加贴近于台风水平风速剖面的实际情况,即台风外围区域风速趋于同一值。另外,对2012-2017 年期间16 期台风风圈半径的反演精度表明:基于WilloughbyB值计算法与HY-2A 海面风场资料提取的34 kt 与50 kt 风圈半径的平均均方根误差为37.6 km与18.3 km,平均偏差为32.4 km 与15.8 km;基于VickeryB值计算法与HY-2A 海面风场资料提取的34 kt 与50 kt 风圈半径的平均均方根误差为37.6 km与18.3 km,平均偏差为32.4 km 与15.8 km。因此,对于西北太平洋区域的热带气旋,Willoughby 等提出的台风形状系数计算法对比Vickery 等的计算法更为准确。综上所述,HY-2A 海面风场资料结合W 法可以有效提取西北太平洋地区的台风风剖面信息。
基于散射计海面风场数据提取台风风圈半径的误差主要来自于以下4 个方面:(1)台风的结构实际并非完全对称结构,本文中通过散射计海面风场数据生成的台风海面风场图中也能反映出这一台风特征;而Holland 模型在结构上属于完全对称结构模型,与台风的真实状态之间存在一定的偏差。(2)HY-2A 散射计数据受海陆影响较大,近海岛与近岸区域的风场数据精度误差增加,最终导致风圈半径提取的精度下降。(3)JTWC 最佳路径数据集记录数据时间为6 小时,每日从0 时(UTC 时)开始记录,而HY-2A 卫星获取西北太平洋区域数据的时间为每日的08 时至11时与20 时至23 时(UTC 时)之间,因此最佳路径数据集数据需要进行插值,造成误差。(4)原始Holland 风场模型属于气压模型,而本文使用的风场模型为风速模型,需要基于台风内部气压与风速之间的关系进行气压与风速的转换,转换结果就有导致误差的可能。
本研究还存在一定的局限性,有待于今后进一步开展相关研究:(1)本研究的方法比较适用于结构完整的热带气旋,对于早期的热带低压与热带风暴级别的气旋反演效果不是特别理想,如何提取尚未成熟期的台风风圈半径是今后研究的方向之一。(2)根据HY-2A 散射计海面风场资料的说明,对于强台风级别(中心最大风速40 m/s)以上的热带气旋其海面风场反演结果与实际风速之间存在较大的偏差,这将导致本文方法在该风速范围内提取得到的风圈半径存在较大偏差。
表 6 2012-2017 年间14 个台风风圈的MAETable 6 MAE of 14 typhoons’ wind radii between 2012 and 2017