CN111723464A - 一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 - Google Patents
一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 Download PDFInfo
- Publication number
- CN111723464A CN111723464A CN202010411137.1A CN202010411137A CN111723464A CN 111723464 A CN111723464 A CN 111723464A CN 202010411137 A CN202010411137 A CN 202010411137A CN 111723464 A CN111723464 A CN 111723464A
- Authority
- CN
- China
- Prior art keywords
- wind
- typhoon
- image
- cloud
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 58
- 238000004088 simulation Methods 0.000 title claims abstract description 31
- 230000008569 process Effects 0.000 claims abstract description 4
- 239000013598 vector Substances 0.000 claims description 154
- 238000004364 calculation method Methods 0.000 claims description 38
- 238000012545 processing Methods 0.000 claims description 24
- 239000006185 dispersion Substances 0.000 claims description 23
- 238000003908 quality control method Methods 0.000 claims description 21
- 238000001514 detection method Methods 0.000 claims description 17
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 16
- 238000004422 calculation algorithm Methods 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 12
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000010586 diagram Methods 0.000 claims description 8
- 238000005516 engineering process Methods 0.000 claims description 7
- 238000003384 imaging method Methods 0.000 claims description 7
- 238000012935 Averaging Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 6
- 239000002131 composite material Substances 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 230000002093 peripheral effect Effects 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 230000003746 surface roughness Effects 0.000 claims description 3
- 238000005457 optimization Methods 0.000 description 9
- 230000011218 segmentation Effects 0.000 description 6
- 238000000605 extraction Methods 0.000 description 3
- 240000006053 Garcinia mangostana Species 0.000 description 2
- 235000017048 Garcinia mangostana Nutrition 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241000345998 Calamus manan Species 0.000 description 1
- 240000000161 Lagerstroemia indica Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 239000009438 liyan Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 235000012950 rattan cane Nutrition 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 238000012502 risk assessment Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，涉及卫星气象遥感信息技术应用领域，通过遥感图像分析手段，结合最小二乘法椭圆拟合，模拟得到海面椭圆形台风风场的椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与x轴正方向的夹角参数，给出海面风场在不同方位的最大风速空间分布范围。该方法具有操作简单，实时性好的特点，能够模拟海面风场最大风速半径的时空异质性分布特征，得到的最大风速点拟合椭圆可以有效描述海面风场能量的时空分布范围与方位，能为风暴潮、海浪等海洋动力过程模拟提供更精确的外部驱动。
Description
技术领域
本发明涉及卫星气象遥感信息技术应用领域，特别是一种基于遥感影像特征的台风椭圆型风场参数化模拟方法。
背景技术
台风是一种发生频率高、影响范围广、突发性强、成灾危害大的自然灾害。台风灾害的致灾因子有大风、暴雨、风暴潮、海浪等，其中大风是风灾的重要组成部分，也是风暴潮、海浪等致灾因子的重要驱动因素(见：“方伟华，林伟.面向灾害风险评估的台风风场模型研究综述.地理科学进展,2013,32(06):852-867.”)。
无论是计算台风大风还是构造人造台风，最大风速半径Rmax都是一个重要而复杂的因子，对预报台风引起的大风、海浪和风暴恶劣天气非常关键。目前，台风风场常用的计算方法有藤田(1952)提出的藤田公式(见：“Fujita,T.,1952.Pressure DistributionWithin Typhoon.Geophys.Mag.,23:437-451.”)和美国气象局水文气象部(1954)提出的Myers公式等(见：“陈孔沫.台风气压场和风场模式.海洋学报,1981,3(1):44-55.”)。这些圆对称气压模式的应用效果表明它们能基本反应热带气旋域内实际气压的分布。但是，许多观测事实和研究成果都表明台风的结构具有明显的非对称特征，而包含了轴对称台风的数值模式常常带有某种系统性偏差。胡邦辉等在“热带气旋海面最大风速半径的计算”一文中提出了一种适合海面非对称热带气旋最大风速半径的计算公式，通过给出台风中心气压p0、外围环境气压p∞、台风最大风速Vmax、摩擦系数k、环境温度T、摩擦阻力偏离实际风矢量反方向的夹角φ、纬度等众多参数，利用理论模型方式推导出Rmax(见：“胡邦辉,谭言科，王举.热带气旋海面最大风速半径的计算.应用气象学报,2004(04):427-435.”)。其它的台风最大风速半径Rmax计算方法主要包括：(1)基于历史气象数据对Rmax与台风中心气压、最大风速、纬度等参数进行直接拟合来推算Rmax；(2)基于经验风场模型，建立其他等级风圈半径与Rmax的转换关系(见：“雷小途,陈联寿.热带气旋风场模型构造及特征参数估算.地球物理学报,2005,48(1):25-31.”)；(3)基于台风地面实测数据及经验模型，采用误差分析法确定最优Rmax。
近年来，卫星资料被广泛应用于台风风场反演。根据采用的卫星资料的不同，主要方法有两类：(1)基于可见光和红外遥感数据的云迹风与水汽导风反演；(2)基于微波遥感数据的海面风场反演。云迹风与水汽导风方法利用连续观测的卫星云图中云和水汽特征来反演风场，得到风矢量；微波海面风场反演是利用卫星散射计、高度计、微波辐射计等进行测风，存在的问题是获得的海面风场的时间和空间分辨率较低(见：“李艳兵，黄思训，翟景秋.卫星反演风场进展概述.气象科学,2009,29(2):277-284.”)。近年来，合成孔径雷达(SAR)作为微波雷达也能获得高空间分辨率的海面风场信息，但影像覆盖范围小，资料成本高，其时空覆盖率远远不能满足中大尺度海面风场动态监测的需要(见：“周旋,杨晓峰,李紫薇,于暘,毕海波,马胜,李晓峰.基于星载SAR数据的台风参数估计及风场构建.中国科学:地球科学,2014,44(02):355-366.”)。利用QuikSCAT卫星遥感风场数据，陈德文等结合JTWC台风参数资料，将遥感平均风剖面与Holland台风模型进行最小二乘拟合，反演了台风最大风速半径，证明了卫星遥感风场可用于反演台风最大风速半径(见：“陈德文，董剑，袁方超.基于QuikSCAT卫星遥感风场的台风最大风速半径反演及个例分析.海洋通报,2012.31(04):376-383.”)。
综上所述，由于实测气象数据的不足，Rmax常通过其他台风参数结合相关统计或物理模型推算获得。无论是基于经验关系还是理论模型，台风最大风速半径计算相对复杂，计算结果的不确定性大，且计算量较大。而基于微波遥感数据获得的海面风场的时间和空间分辨率较低，不适合用于高精度、实时台风风场反演和台风关键参数提取。
发明内容
本发明所要解决的技术问题是克服现有技术的不足而提供一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，本发明方法通过遥感图像处理技术实现，不仅实时性好，时空分辨率高，而且得到的是台风引起的海面风场在不同方向的最大风速半径结果，体现了海面风场最大风速分布的各向异性特征，更符合实际的风场能量空间分布规律，能为风暴潮、海浪等海洋动力过程模拟提供更精确的外部驱动。
本发明为解决上述技术问题采用以下技术方案：
根据本发明提出的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，包括以下步骤：
步骤S1、对包含台风结构的气象卫星遥感影像的红外波段进行分段线性拉伸增强、图像二值化、台风区域裁剪预处理；
步骤S2、基于红外波段的像元灰度特征，自动设定可变尺寸搜索窗口，结合移动搜索和阈值控制技术，确定台风中心位置；
步骤S3、由三个相邻时相的红外波段数据，基于模板相关匹配技术两两匹配，获得两个时次的原始云导风矢量，对其进行时间连续性质量控制，得到中间时相的台风实时风场；
步骤S4、建立气象卫星遥感影像水汽波段与亮温之间的拟合关系，结合大气廓线进行查找或推算，指定风矢量的气压高度，进行云压检测和空间连续性质量控制，得到中间时相的台风立体风场；
步骤S5、利用三维风场的风速插值计算，由台风立体风场推算海面10m高度风场；利用海面粗糙度和海面平均波高进行风速修正，得到海面风场；
步骤S6、利用D∞算法，推算海面风场以台风中心为中心、在不同方向上的最大风速点，过滤极值风速点后，基于最小二乘法进行椭圆拟合，模拟最大风速点的椭圆分布，输出椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与平面直角坐标系x轴正方向的夹角，得到台风引起的海面风场在不同方向的最大风速分布范围。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，所述步骤S1的具体操作如下：
步骤S11、对气象卫星遥感影像的红外波段进行像元值直方图统计，选择需要拉伸的下限像元值Flow和上限像元值Fhig，确定拉伸区间[Flow,Fhigh]；将位于拉伸区间[Flow,Fhigh]内的像元值映射至[0,255]，得到台风特征增强处理后的灰度影像；
步骤S12、对步骤S11所得到的台风特征增强处理后的灰度影像进行统计，确定二值化处理的下限像元值Blow和上限像元值Bhigh，进行灰度影像的二值化处理：当像元值在[Blow,Bhigh]之间时，将该像元赋值1，对应增强后的台风风区像元；否则赋值0，对应增强后的台风眼或者台风风区外围区域像元，生成台风二值影像；
步骤S13、给定左上角和右下角行列号信息，对步骤S12生成的台风二值影像进行裁剪，得到包含台风中心的限定范围二值影像。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，所述的步骤S2的具体操作如下：
步骤S21、假定台风眼区域的像元均匀分布在一个正方形搜索窗口内，计算步骤S13得到的限定范围二值影像的像元值之和SP，取正方形搜索窗口边长
取判别台风区的阈值TPthres＝INT(0.9×SP)；
步骤S22、利用设定的正方形搜索窗口在步骤S13得到的限定范围二值影像中按行、列号进行移动搜索，统计位于搜索窗口范围内的二值影像像元值之和BSP；如果BSP>TPthres，记下搜索窗口中心在二值影像中的行列号，换算成对应的地理坐标x,y；当整个影像搜索完成后，对所有记录的地理坐标x,y求平均，获得台风中心位置；
步骤S23、如果没有搜索到台风中心，取步长Lstep＝TPthre×sr％，sr为台风区阈值缩减因子，令TPthresh＝TPthres-Lstep，重新按步骤S22进行搜索，并输出台风中心位置；如果TPthres一直减小到值为负还没有搜索到台风中心，则认为当前的限定范围二值影像中不存在台风中心，搜索过程结束。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，步骤S23中，sr的取值范围在5～30之间。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，所述的步骤S3的具体操作如下：
步骤S31、选择相邻时相的两幅气象卫星遥感影像红外波段数据，记为云图a和云图b，两幅云图的成像时间间隔为△t；
步骤S32、根据台风天气系统的真实尺度和卫星影像的空间分辨率，这两者相除并取整，确定示踪云大小，示踪云大小取为M×M，M代表在云图a的数据矩阵中截取的示踪云子图像的行、列数；取N×N作为搜索区的大小，N代表在云图b的数据矩阵中截取的搜索区子图像的行、列数；
步骤S33、在云图b上截取搜索区子图像；给定搜索区的像元灰度平均值条件阈值SDThresh，如果搜索区的像元灰度平均值小于SDThres，表示当前截取的搜索区不在台风风区内，不参与台风实时风矢量计算；否则，表示当前截取的搜索区是有效搜索区；在云图a上对应于该搜索区的范围内，以M为行、列数，截取示踪云子图像，子图像的中心与搜索区的中心一致；
步骤S34、计算示踪云子图像在搜索区子图像内的相关系数矩阵；相关系数R(m,n)计算公式如下：
式中，i,j∈[1,M]，i,j分别为示踪云子图像当中的像元行、列号；m,n∈[1,N-M]，m,n分别为示踪云子图像在搜索区子图像内移动的行、列方向的步长改变量；f(i,j)为示踪云子图像的像元灰度值；g(i+m,j+n)为搜索区子图像内目标匹配子图像的像元灰度值；
为示踪云子图像的平均灰度；
为目标匹配子图像的平均灰度；如果相关系数R(m,n)大于预设的子图像匹配度相关系数阈值RThres，认为是有效匹配；
步骤S35、对于有效匹配，记录对应的目标匹配子图像的中心行、列号，换算成对应的地理坐标x,y；按步骤S34进行循环，完成示踪云子图像在搜索区子图像内的相关系数矩阵计算后，对所有有效匹配记录对应的目标匹配子图像的地理坐标x,y求平均，得到云图a中的示踪云子图像在云图b中的搜索区子图像内的匹配子图像的中心平均地理坐标(xm,ym)；
步骤S36、计算搜索区子图像中心与匹配云子图像的平均中心的距离D，除以两幅云图的成像时间间隔Δt，得到该搜索区中心的台风实时风矢量的风速大小V；风矢量的位置为搜索区子图像的中心位置，方向由搜索区子图像中心指向匹配云子图像的平均中心，用方位角θ表示，计算公式如下：
式中，(xs,ys)为搜索区子图像中心的地理坐标；(xm,ym)为匹配云子图像的平均中心的地理坐标；
步骤S37、对云图b按行、列顺序以M为步长进行循环，截取相应的搜索区子图像；按照步骤S33在云图a上截取示踪云子图像，按照步骤S34、S35进行相关系数计算和图像匹配，按照步骤S36进行计算，得到每个搜索区中心位置的台风实时风矢量；
步骤S38、按照步骤S31～S37操作，对三个相邻时相时相的红外波段数据进行两两匹配，获得连续两个时次的台风风场的原始云导风矢量；
步骤S39、以前一时次的原始云导风矢量作为后一时次原始云导风矢量的参照，计算出连续两个时次云导风的矢量差；如果矢量差大于两个时次平均标量风，则认为该风矢量不具有时间连续性，将后一时次原始云导风矢量中位于该位置的风矢量剔除，得到经过时间连续性质量控制的中间时相的台风实时风场。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，步骤S32中，N等于3倍的M。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，步骤S4的具体操作如下：
步骤S41、选择气象卫星遥感影像的水汽波段，建立水汽波段灰度值与亮温值t之间的拟合关系，得到台风期间亮温的空间分布；
步骤S42、根据S38计算得到的原始云导风矢量中的台风实时风矢量位置，提取该位置对应的亮温值t；利用亮温值t在标准大气温度垂直廓线表中查找该亮温值对应的气压高度p；如果能够直接查找，则直接完成对该台风实时风矢量的气压高度指定；
步骤S43、如果亮温值t不在标准大气温度垂直廓线表的亮温范围内，那么首先根据廓线表中的亮温最大值t1和最小值t2，分别查出其对应的气压高度p1和p2；然后对(t1，p1)和(t2，p2)进行对数线性内插，求得亮温值t对应的气压高度p；方法如下：
t1＝a+b×ln(p1)
t2＝a+b×ln(p2)
通过上述两公式计算系数a和b：
则气压高度p为：
其中，e为自然底数；
根据计算的气压高度p，实现该台风实时风矢量的高度指定；
步骤S44、重复步骤S42和S43，完成S38步骤得到的两个时次原始云导风矢量中的所有台风实时风矢量的高度指定；
步骤S45、对完成高度指定的连续两个时次云导风矢量图层，以前一时次云导风矢量作为后一时次云导风矢量的参照，以要判定的风矢量所在位置为中心，计算其A×A邻域内云压变化的标准差；设定云压变化标准差的阈值，如果计算得到的云压变化的标准差差值大于该阈值，则认定该位置的实时风矢量位于复合云层，是云压检测不合格风矢量，将后一时次云导风矢量中位于该位置的风矢量剔除，得到经过云压检测控制的中间时相的台风立体风场；其中，A为正奇数；
步骤S46、以经过云压检测控制的中间时相的台风立体风场作为参照风场，对经过时间连续性质量控制的中间时相的台风实时风场中的风矢量进行判别和高度赋值处理：如果该风矢量在参照风场中有对应风矢量，则把对应风矢量的高度赋给该风矢量；否则，则认为该风矢量没有通过云压检测，予以剔除；从而得到经过时间连续性质量控制且满足云压检测的中间时相的台风立体风场；
步骤S47、进一步对步骤S46得到的中间时相的台风立体风场的风矢量进行空间连续性质量控制，分为风向控制和风速控制；每个风矢量均对其A×A邻域内的风矢量计算风向离差均方根和风速离差均方根，计算公式如下：
式中，Ds、Dd为目标风矢量的风速离差均方根和风向离差均方根；Vc、θc为目标风矢量的风速和风向；o＝1，2，……，B；Vo、θo为第o个邻域风矢量的风速和风向；B为邻域内的风矢量数；
分别设定风向离差均方根和风速离差均方根阈值，若计算得到的离差均方根中的任意一个超过对应的阈值，则认为该风矢量为不合格风矢量，予以剔除；对所有的风矢量进行判别，最终得到满足时空质量控制要求的中间时相的台风立体风场。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，A为5。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，所述的步骤S5的具体操作如下：
步骤S51、利用三维风场的风速插值公式对时空质量控制处理后的台风立体风场进行插值计算，得到海面10m高度风场；风速V10计算公式如下：
式中，V为台风立体风场实时风矢量的风速值；c为经验常数；pb为基准气压高度值；p为台风立体风场实时风矢量的气压高度；
步骤S52、根据海面粗糙度和海面平均高度，利用海面10m高度风场的风速V10换算海面风速Vs，计算得到台风期间的海面风场；计算公式如下：
式中，Vs为海面风速；Z0为海面粗糙度；Zs为海面平均高度。
作为本发明所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法进一步优化方案，所述步骤S6的具体操作如下：
步骤S61、在计算矢量格式的海面风场的同时，同步输出GeoTiff格式的重采样风速图像，图像的空间分辨率为原始影像空间分辨率的M倍，图像的像元值为海面风速值，图像的地理坐标和海面风场的地理坐标一致；
步骤S62、利用D∞算法，推算海面风场在不同方向上的最大风速半径分布；具体实现方法是以台风中心为中心，向周围任意方向作射线，识别与该射线相交的重采样风速图像上的所有像元，记录像元值最大的像元位置，为该方向的最大风速点，其与台风中心的距离即为该方向上的海面最大风速半径；
步骤S63、以台风中心为中心，以平面直角坐标系x轴正方向为起始方向，按等角度间隔，按照步骤S62的方法，逆时针在360度范围内计算各角度方向的最大风速点，得到海面最大风速半径在各方向的长度；对海面最大风速点进行极值过滤，去掉偏离台风中心的最大风速点；
步骤S64、采用基于最小二乘的椭圆拟合算法，拟合最大风速半径相关的椭圆，输出椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与x轴正方向的夹角，给出台风引起的海面风场的最大风速分布范围。
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
(1)利用高时空分辨率的气象卫星遥感影像红外波段和水汽波段数据(例如Himawari-8/9日本气象卫星：时间分辨率10分钟，空间分辨率1km)，通过遥感图像解译和模拟技术，获取台风中心位置和台风立体风场，进而实时推算海面风场的最大风速的各向异性分布；
(2)与传统的台风风场数值模拟给出单一最大风速半径结果相比，本方法得到的海面风场最大风速分布结果不仅时空分辨率高，实时性强，而且可以得到海面风场最大风速结果在不同方向的异质性分布，以椭圆的形式体现海面风场能量的空间分布形态，是对单一台风最大风速半径估算结果的重要改进。
附图说明
图1为本发明一种基于遥感影像特征的台风椭圆型风场参数化模拟方法的步骤流程图。
图2为本发明方法一个具体实施例的气象卫星红外波段预处理结果图；其中，(a)是气象卫星红外波段原图，(b)是线性分段拉伸处理结果图；(c)是台风区二值影像图；(d)是台风中心区域裁剪结果图。
图3为本发明方法一个具体实施例得到的台风中心提取结果图。
图4为本发明方法一个具体实施例得到的台风实时风场图。
图5为本发明方法一个具体实施例得到的经质量控制后的台风立体风场图。
图6为本发明方法一个具体实施例得到的海面风场矢量图。
图7为本发明方法一个具体实施例得到的海面风场栅格图。
图8为本发明方法一个具体实施例得到的海面最大风速点拟合椭圆分布图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图及具体实施例对本发明进行详细描述。
本具体实施例选择2018年台风“山竹”三个时相的Himawari-8气象卫星影像作为遥感数据源，影像成像时刻为2018年9月15日23：50分、2018年9月16日00:00分和00:10分，采用ENVI IDL5.3编写遥感图像处理算法、matlab编写椭圆拟合算法、ArcGIS10.2作为矢量数据处理工具，采用C#开发可执行程序框架，进行算法集成，完成海面风场各向异性的最大风速范围遥感模拟。
一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，如图1所示，具体包括以下步骤：
S1、选取2018年9月16日00:00分的Himawari-8气象卫星遥感影像红外波段B13(IR1)为例，作为台风中心提取的原始数据，如图2中的(a)所示，具体操作如下：
S11、对台风影像的红外波段B13(IR1)进行直方图统计，可知该波段的像元值分布范围为384～975，影像中台风云区的灰度值对比度较低。台风中心大致位于台风眼的中心位置，为了突出台风眼区域与台风螺旋结构区域的对比度，增强台风眼的可识别特征，采用分段线性拉伸方法对红外波段影像进行增强处理，公式如下：
DNen＝0；DN≤Flow (15)
DNen＝255；DN>Fhigh (17)
式中，DN和DNen分别为红外波段原始像元值和分段线性拉伸增强处理后的像元值；Flow和Fhigh分别为需要拉伸的下限像元值和上限像元值，在本具体实施例中取500和800。通过公式(15)～(17)变换，可将红外波段影像位于[Flow,Fhigh]内的像元值映射至[0,255]，将所有小于Flow的像元值压缩为0，将所有大于Fhigh的像元值压缩为255。具体增强处理结果图如图2中的(b)所示。
S12、设定阈值区间[Blow,Bhigh]，对步骤S11所得到的台风特征增强处理后的灰度影像进行统计，确定二值化处理的下限像元值Blow和上限像元值Bhigh，进行灰度影像的二值化处理。具体处理方法是：对影像像元按行、列号进行循环，当像元值在[Blow,Bhigh]之间时，将该像元赋值1，对应增强后的台风风区像元；否则赋值0，对应增强后的台风眼或者台风风区外围区域像元，生成台风二值影像。在本具体实施例中，取Blow＝10，取Bhigh＝200，得到的台风二值影像如图2中的(c)所示。
S13、给定左上角和右下角行列号信息，对S12处理的台风二值影像进行裁剪，缩小台风中心的搜索范围，得到包含台风中心的限定范围二值影像，并给裁剪后的影像赋予正确的投影和坐标信息。结果如图2中的(d)所示。
S2、基于红外波段的像元灰度特征，自动设定可变尺寸搜索窗口，结合移动搜索和阈值控制技术，确定台风中心位置。具体操作如下：
S21、假定台风眼区域的像元均匀分布在一个正方形搜索窗口内，计算步骤S13得到的限定范围二值影像的像元值之和SP，取正方形搜索窗口边长
取判别台风区的阈值TPthresh＝INT(0.9×SP)，其含义为如果搜索窗口内的属于台风眼区域的像元值之和达到一定大小，说明台风中心位于该区域内；
S22、利用设定的正方形搜索窗口在S13得到的限定范围二值影像中按行、列号进行移动搜索，统计位于搜索窗口范围内的二值影像像元值之和BSP。如果BSP>TPthresh，记下搜索窗口中心在二值影像中的行列号，换算成对应的地理坐标x,y。当整个影像搜索完成后，对所有记录的地理坐标x,y求平均，获得台风中心位置。
S23、如果没有搜索到台风中心，取台风区阈值缩减因子sr＝10，步长Lstep＝TPthresh×10％，取新的判别台风区阈值TPthres＝TPthres-Lstep，重新按S22步骤进行搜索，并输出台风中心位置。如果TPthres一直减小到值为负还没有搜索到台风中心，则认为当前的限定范围二值影像中不存在台风中心，搜索过程结束。
图3显示了本具体实施例的台风中心提取结果图。结果显示，本实施例反演的台风中心位于中央气象台资料台风中心的北偏西55.66度，距离12.32km。
S3、分别选择2018年9月15日23:50分、2018年9月16日00:00分和00:10分的三个相邻时相的Himawari-8气象卫星遥感影像红外波段B13(IR1)，作为台风风场提取的原始数据，基于模板相关匹配技术两两匹配，获得两个时次台风风场的原始云导风矢量，估算台风实时风速。下面以2018年9月16日00:00分和00:10分两个时相的遥感影像波段数据为例，具体操作如下：
S31、将00:00分的红外波段数据记为云图a，将00:10分的红外波段数据记为云图b，两幅云图的成像时间间隔10分钟，记为△t；
S32、根据台风天气系统的真实尺度和卫星影像的空间分辨率，两者相除并取整，确定示踪云大小，取为M×M，M代表在云图a的数据矩阵中截取的示踪云子图像的行、列数；取N×N作为搜索区的大小，N代表在云图b的数据矩阵中截取的搜索区子图像的行、列数，取N等于3倍的M；
S33、在云图b上截取搜索区子图像。给定搜索区的像元灰度平均值条件阈值SDThresh，如果搜索区的像元灰度平均值小于SDThresh，表示当前截取的搜索区不在台风风区内，不参与台风实时风矢量计算；否则，表示当前截取的搜索区是有效搜索区。在本具体实施例中，取SDThresh＝600。在云图a上对应于该搜索区的范围内，以M为行、列数，截取示踪云子图像，子图像的中心与搜索区的中心一致；
S34、计算示踪云子图像在搜索区子图像内的相关系数矩阵。相关系数R(m,n)计算公式如下：
式中，i,j∈[1,M]分别为示踪云子图像当中的像元行、列号；m,n∈[1,N-M]分别为示踪云子图像在搜索区子图像内移动的行、列方向步长改变量；f(i,j)为示踪云子图像的像元灰度值；g(i+m,j+n)为搜索区子图像内目标匹配子图像的像元灰度值；
为示踪云子图像的平均灰度；
为目标匹配子图像的平均灰度。取子图像匹配度相关系数阈值RThresh＝0.9，如果相关系数R(m,n)大于RThresh，认为是有效匹配；
S35、对于有效匹配，记录对应的目标匹配子图像的中心行、列号，换算成对应的地理坐标x,y。按步骤S34进行循环，完成示踪云子图像在搜索区子图像内的相关系数矩阵计算后，对所有有效匹配记录对应的目标匹配子图像的地理坐标x,y求平均，得到云图a中的示踪云子图像在云图b的搜索区子图像内的匹配子图像的中心平均地理坐标(xm,ym)；
S36、计算搜索区子图像中心(xs,ys)与匹配子图像的平均中心的距离D，除以两幅云图的成像时间间隔Δt，得到该搜索区中心的台风实时风矢量的风速大小V。风矢量的位置为搜索区子图像的中心位置，方向由搜索区子图像中心指向匹配云子图像的平均中心，用方位角θ表示，计算公式如下：
式中，(xs,ys)为搜索区子图像中心的地理坐标；(xm,ym)为匹配云子图像的平均中心的地理坐标；
S37、对云图b按行、列顺序以M为步长进行循环，截取相应的搜索区子图像；按照步骤S33在云图a上截取示踪云子图像，按照步骤S34、S35进行相关系数计算和图像匹配，按照步骤S36计算得到每个搜索区中心位置的台风实时风矢量。图4显示了本具体实施例模拟得到的逆时针螺旋结构的台风实时风场；
S38、按照步骤S31～S37操作，对三个相邻时相的红外波段数据进行两两匹配，获得2018年9月15日23:50分和2018年9月16日00:00分两个时次的台风风场的原始云导风矢量；
S39、以前一时次的原始云导风矢量作为后一时次原始云导风矢量的参照，计算出连续两个时次云导风的矢量差。如果矢量差大于两个时次平均标量风，则认为该风矢量不具有时间连续性，将后一时次原始云导风矢量中位于该位置的风矢量剔除，得到经过时间连续性质量控制后的中间时相(2018年9月16日00:00分)的台风实时风场。
S4、假定台风云区主要为不透明高云，水汽通道卫星探测的辐射值与云下的大气背景无关。建立气象卫星遥感影像水汽波段与亮温之间的拟合关系，由亮温进行查找或推算，确定台风实时风矢量的气压高度，得到台风立体风场。具体操作如下：
S41、查询Himawari-8气象卫星遥感影像的数据校正表，选择Himawari-8气象卫星遥感影像的水汽波段B8(IR3)，建立水汽波段灰度值与亮温值t之间的拟合关系，得到台风期间亮温的空间分布。在本具体实施例中，拟合方程如下：
BTWV＝-0.083×DNWV+311.11 (23)
式中，BTWV为水汽波段的亮温值，单位为K；DNWV为水汽波段的灰度值。
S42、根据S38计算得到的原始云导风矢量中的台风实时风矢量位置，提取该位置对应的亮温值t。利用亮温值t在标准大气温度垂直廓线表中查找该亮温值对应的气压高度p。如果能够直接查找，则直接完成对该台风实时风矢量的气压高度指定；
S43、如果亮温值t不在标准大气温度垂直廓线表的亮温范围内，那么首先根据廓线表中的亮温最大值t1和最小值t2，分别查出其对应的气压高度p1和p2；然后对(t1，p1)和(t2，p2)进行对数线性内插，求得亮温值t对应的气压高度p。气压外推方法如下：
t1＝a+b×ln(p1) (24)
t2＝a+b×ln(p2) (25)
通过上述两公式计算系数a和b：
则气压高度p为：
其中，e为自然底数；
根据计算的气压高度p，实现该台风实时风矢量的高度指定。
S44、重复步骤S42和S43，完成S38步骤得到的2018年9月15日23:50分和2018年9月16日00:00分两个时次原始云导风矢量中的所有台风实时风矢量的高度指定；
S45、对完成高度指定的连续两个时次云导风矢量图层，以前一时次(2018年9月15日23:50分)云导风矢量作为后一时次(2018年9月16日00:00分)云导风矢量的参照，以要判定的风矢量所在位置为中心，计算其5×5邻域内云压变化的标准差。设定云压变化标准差的阈值(依据经验可设为100)，如果计算得到的云压变化的标准差差值大于该阈值，则认定该位置的实时风矢量位于复合云层，为云压检测不合格风矢量，将后一时次云导风矢量中位于该位置的风矢量剔除，得到经过云压检测控制的中间时相(2018年9月16日00:00分)的台风立体风场。
S46、以经过云压检测控制的中间时相的台风立体风场作为参照风场，对经过时间连续性质量控制的中间时相台风实时风场中的实时风矢量进行判别和高度赋值处理：如果该风矢量在参照风场中有对应风矢量，则把对应风矢量的高度赋给该矢量；否则，则认为该矢量没有通过云压检测，予以剔除。从而得到经过时间连续性质量控制且满足云压检测的中间时相(2018年9月16日00:00分)的台风立体风场；
S47、进一步对S46步骤得到的中间时相台风立体风场的风矢量进行空间连续性质量控制，分为风向控制和风速控制。每个风矢量均对其5×5邻域内的风矢量计算风向离差均方根和风速离差均方根，计算公式如下：
式中，Ds、Dd为目标风矢量的风速离差均方根和风向离差均方根；Vc、θc为目标风矢量的风速和风向；o＝1，2，……，B；Vo、θo为第o个邻域风矢量的风速和风向；B为邻域内的风矢量数，邻域取台风风场中以目标风矢量为中心的5×5区域，因此B＝25。
分别设定风向离差均方根和风速离差均方根阈值，本具体实施例中，风矢量方向以正东方向为0°，逆时针把平面分为360等分，风向离差均方根阈值取120°，风速离差均方根阈值取20m/s，两者均为经验值。若计算得到的离差均方根中的任意一个超过所设阈值，则认为该风矢量无效，予以剔除。对所有风矢量都经过判别后，最终得到满足时空质量控制要求的时空后中间时相(2018年9月16日00:00分)的台风立体风场。图5显示了本具体实施例模拟得到的经过时空质量控制后的台风立体风场。
S5、利用三维风场的风速插值公式对时空质量控制处理后的台风立体风场进行插值计算，得到海面10m高度风场；利用海面粗糙度和海面平均波高进行风速修正，得到海面风场。具体操作如下：
S51、利用三维风场的风速插值公式对台风立体风场进行插值计算，得到海面10m高度风场。风速V10计算公式如下：
式中，V为台风立体风场实时风矢量的风速值；c为经验常数，在本具体实施例中取7.0；p为前述计算得到的台风立体风场实时风矢量的气压高度；pb为基准气压高度值。在标准大气条件下，海面的气压约为1013hpa，在海拔3000m以内，每升高10m，大气压减小100pa，由此假定海平面10米高度的气压值为1012hpa。据此在本具体实施例中，基准气压高度值pb取1012hpa。
S52、根据海面粗糙度和海面平均高度，利用海面10m高度风场的风速V10换算海面风速Vs，计算得到台风期间的海面风场。计算公式如下：
式中，Vs为海面风速；Z0为海面粗糙度，取值范围为0.001～0.01m，在本具体实施例中取0.008；Zs为海面平均高度，可取台风期间海面的平均波高，在本具体实施例中取2m。
图6显示了由本具体实施例模拟得到的逆时针螺旋结构的台风立体风场修正后的海面风场。
S6、利用D∞算法，推算海面风场以台风中心为中心、在不同方向上的最大风速点，过滤极值风速点后，基于最小二乘法进行椭圆拟合，得到最大风速点的椭圆分布，模拟台风引起的海面风场的最大风速分布范围。具体操作如下：
S61、在计算矢量格式的海面风场的同时，同步输出GeoTiff格式的重采样风速图像，图像的空间分辨率为原始影像空间分辨率的M倍，图像的像元值为海面风速值，图像的地理坐标和海面风场的地理坐标一致。图7显示了本具体实施例得到的海面风场矢量图进行矢量～栅格转换后得到的海面风场栅格影像图；
S62、利用D∞算法，推算海面风场在不同方向上的最大风速半径分布。具体实现方法是以台风中心为中心，向周围任意方向作射线，识别与该射线相交的重采样风速图像上的所有像元，记录像元值最大(即该方向上海面风速值最大)的像元位置，为该方向的最大风速点，其与台风中心的距离即为该方向上的海面最大风速半径。具体计算时，需要分象限来处理，然后在每个象限里对射线的不同角度分别处理。以第一象限为例：
1)首先计算台风中心点至风速图像右上角的长度l和角度θ；
2)对当前射线的角度α(0<α≤90)进行判断，当射线角度小于θ时，将长度l按列分段，即用长度l除以列方向的分辨率，得到段数；当射线角度大于该角度θ时，将长度l按行分段，即用长度l除以行方向的分辨率，得到段数；
3)按段数对各分段从0开始累加，计算各分段点的坐标。然后通过分段点的坐标反算各分段点的行、列号，提取各分段点对应的风速图像像元的风速值；
4)比较所有分段点的风速值的大小，取得最大风速值所在的分段点的行、列号，推算出对应的坐标值；
5)计算该坐标值与台风中心之间的距离，得到该方向上的海面最大风速半径；
6)如果射线位于第二、三、四象限，则分别计算台风中心点至风速图像左上角、左下角和右下角的长度l和角度θ，然后按照2)～5)的步骤计算，得到该方向上的海面最大风速半径。
S63、以台风中心为中心，以平面直角坐标系x轴正方向为起始方向，按等角度间隔，按照S62的步骤，逆时针在360度范围内对4个象限内的射线分别处理，计算各角度方向的海面最大风速点，得到海面最大风速半径在各方向的长度。图8中的海面风场最大风速点位显示的是以台风中心为中心，以15度为间隔，计算得到的24个方位的最大风速点位。对海面最大风速点进行极值过滤，去掉偏离台风中心的最大风速点。在本具体实施例中，对24个方位的海面最大风速半径值从大到小排列，去掉前60％的最大风速点，结果如图8中的椭圆拟合的最大风速点位所示。
S64、参考Radim 
和Jan Flusser在“Numerically Stable Direct LeastSquares Fitting of Ellipses”一文中提出的椭圆拟合方法，对五个及五个以上的最大风速点位采用基于最小二乘的椭圆拟合，拟合出最大风速半径相关的椭圆，输出四个椭圆参数：椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与x轴正方向的夹角，给出台风引起的海面风场的最大风速分布范围，结果如图8中的海面最大风速点拟合椭圆所示。针对本具体实施例，最终得到的台风“山竹”在2019年9月16日00:00分引起的海面椭圆台风风场参数为：椭圆中心经纬度坐标(东经117.1°，北纬19.7°)，椭圆长轴长度203.21km，短轴长度163.61km，椭圆长轴与平面直角坐标系x轴正方向的夹角127.88度。由台风引起的海面风场的能量分布主要集中在该区域范围内。
上面结合附图和附表对本发明的实施方式作了详细地说明，但是本发明并不局限于上述实施方式所提到的算法和所使用的气象卫星遥感数据类型。在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化，进一步提高海面椭圆台风风场参数的精度。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围内。
Claims (10)
1.一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，包括以下步骤：
步骤S1、对包含台风结构的气象卫星遥感影像的红外波段进行分段线性拉伸增强、图像二值化、台风区域裁剪预处理；
步骤S2、基于红外波段的像元灰度特征，自动设定可变尺寸搜索窗口，结合移动搜索和阈值控制技术，确定台风中心位置；
步骤S3、由三个相邻时相的红外波段数据，基于模板相关匹配技术两两匹配，获得两个时次的原始云导风矢量，对其进行时间连续性质量控制，得到中间时相的台风实时风场；
步骤S4、建立气象卫星遥感影像水汽波段与亮温之间的拟合关系，结合大气廓线进行查找或推算，指定风矢量的气压高度，进行云压检测和空间连续性质量控制，得到中间时相的台风立体风场；
步骤S5、利用三维风场的风速插值计算，由台风立体风场推算海面10m高度风场；利用海面粗糙度和海面平均波高进行风速修正，得到海面风场；
步骤S6、利用D∞算法，推算海面风场以台风中心为中心、在不同方向上的最大风速点，过滤极值风速点后，基于最小二乘法进行椭圆拟合，模拟最大风速点的椭圆分布，输出椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与平面直角坐标系x轴正方向的夹角，得到台风引起的海面风场在不同方向的最大风速分布范围。
2.根据权利要求1所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，所述步骤S1的具体操作如下：
步骤S11、对气象卫星遥感影像的红外波段进行像元值直方图统计，选择需要拉伸的下限像元值Flow和上限像元值Fhigh，确定拉伸区间[Flow,Fhigh]；将位于拉伸区间[Flow,Fhigh]内的像元值映射至[0,255]，得到台风特征增强处理后的灰度影像；
步骤S12、对步骤S11所得到的台风特征增强处理后的灰度影像进行统计，确定二值化处理的下限像元值Blow和上限像元值Bhig，进行灰度影像的二值化处理：当像元值在[Blow,Bhigh]之间时，将该像元赋值1，对应增强后的台风风区像元；否则赋值0，对应增强后的台风眼或者台风风区外围区域像元，生成台风二值影像；
步骤S13、给定左上角和右下角行列号信息，对步骤S12生成的台风二值影像进行裁剪，得到包含台风中心的限定范围二值影像。
3.根据权利要求2所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，所述的步骤S2的具体操作如下：
步骤S21、假定台风眼区域的像元均匀分布在一个正方形搜索窗口内，计算步骤S13得到的限定范围二值影像的像元值之和SP，取正方形搜索窗口边长
取判别台风区的阈值TPthresh＝INT(0.9×SP)；
步骤S22、利用设定的正方形搜索窗口在步骤S13得到的限定范围二值影像中按行、列号进行移动搜索，统计位于搜索窗口范围内的二值影像像元值之和BSP；如果BSP>TPthresh，记下搜索窗口中心在二值影像中的行列号，换算成对应的地理坐标x,y；当整个影像搜索完成后，对所有记录的地理坐标x,y求平均，获得台风中心位置；
步骤S23、如果没有搜索到台风中心，取步长Lstep＝TPthresh×sr％，sr为台风区阈值缩减因子，令TPthres＝TPthresh-Lstep，重新按步骤S22进行搜索，并输出台风中心位置；如果TPthres一直减小到值为负还没有搜索到台风中心，则认为当前的限定范围二值影像中不存在台风中心，搜索过程结束。
4.根据权利要求3所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，步骤S23中，sr的取值范围在5～30之间。
5.根据权利要求1所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，所述的步骤S3的具体操作如下：
步骤S31、选择相邻时相的两幅气象卫星遥感影像红外波段数据，记为云图a和云图b，两幅云图的成像时间间隔为△t；
步骤S32、根据台风天气系统的真实尺度和卫星影像的空间分辨率，这两者相除并取整，确定示踪云大小，示踪云大小取为M×M，M代表在云图a的数据矩阵中截取的示踪云子图像的行、列数；取N×N作为搜索区的大小，N代表在云图b的数据矩阵中截取的搜索区子图像的行、列数；
步骤S33、在云图b上截取搜索区子图像；给定搜索区的像元灰度平均值条件阈值SDThres，如果搜索区的像元灰度平均值小于SDThresh，表示当前截取的搜索区不在台风风区内，不参与台风实时风矢量计算；否则，表示当前截取的搜索区是有效搜索区；在云图a上对应于该搜索区的范围内，以M为行、列数，截取示踪云子图像，子图像的中心与搜索区的中心一致；
步骤S34、计算示踪云子图像在搜索区子图像内的相关系数矩阵；相关系数R(m,n)计算公式如下：
式中，i,j∈[1,M]，i,j分别为示踪云子图像当中的像元行、列号；m,n∈[1,N-M]，m,n分别为示踪云子图像在搜索区子图像内移动的行、列方向的步长改变量；f(i,j)为示踪云子图像的像元灰度值；g(i+m,j+n)为搜索区子图像内目标匹配子图像的像元灰度值；
为示踪云子图像的平均灰度；
为目标匹配子图像的平均灰度；如果相关系数R(m,n)大于预设的子图像匹配度相关系数阈值RThresh，认为是有效匹配；
步骤S35、对于有效匹配，记录对应的目标匹配子图像的中心行、列号，换算成对应的地理坐标x,y；按步骤S34进行循环，完成示踪云子图像在搜索区子图像内的相关系数矩阵计算后，对所有有效匹配记录对应的目标匹配子图像的地理坐标x,y求平均，得到云图a中的示踪云子图像在云图b中的搜索区子图像内的匹配子图像的中心平均地理坐标(xm,ym)；
步骤S36、计算搜索区子图像中心与匹配云子图像的平均中心的距离D，除以两幅云图的成像时间间隔Δt，得到该搜索区中心的台风实时风矢量的风速大小V；风矢量的位置为搜索区子图像的中心位置，方向由搜索区子图像中心指向匹配云子图像的平均中心，用方位角θ表示，计算公式如下：
式中，(xs,ys)为搜索区子图像中心的地理坐标；(xm,ym)为匹配云子图像的平均中心的地理坐标；
步骤S37、对云图b按行、列顺序以M为步长进行循环，截取相应的搜索区子图像；按照步骤S33在云图a上截取示踪云子图像，按照步骤S34、S35进行相关系数计算和图像匹配，按照步骤S36进行计算，得到每个搜索区中心位置的台风实时风矢量；
步骤S38、按照步骤S31～S37操作，对三个相邻时相时相的红外波段数据进行两两匹配，获得连续两个时次的台风风场的原始云导风矢量；
步骤S39、以前一时次的原始云导风矢量作为后一时次原始云导风矢量的参照，计算出连续两个时次云导风的矢量差；如果矢量差大于两个时次平均标量风，则认为该风矢量不具有时间连续性，将后一时次原始云导风矢量中位于该位置的风矢量剔除，得到经过时间连续性质量控制的中间时相的台风实时风场。
6.根据权利要求5所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，步骤S32中，N等于3倍的M。
7.根据权利要求5所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，步骤S4的具体操作如下：
步骤S41、选择气象卫星遥感影像的水汽波段，建立水汽波段灰度值与亮温值t之间的拟合关系，得到台风期间亮温的空间分布；
步骤S42、根据S38计算得到的原始云导风矢量中的台风实时风矢量位置，提取该位置对应的亮温值t；利用亮温值t在标准大气温度垂直廓线表中查找该亮温值对应的气压高度p；如果能够直接查找，则直接完成对该台风实时风矢量的气压高度指定；
步骤S43、如果亮温值t不在标准大气温度垂直廓线表的亮温范围内，那么首先根据廓线表中的亮温最大值t1和最小值t2，分别查出其对应的气压高度p1和p2；然后对(t1，p1)和(t2，p2)进行对数线性内插，求得亮温值t对应的气压高度p；方法如下：
t1＝a+b×ln(p1)
t2＝a+b×ln(p2)
通过上述两公式计算系数a和b：
则气压高度p为：
其中，e为自然底数；
根据计算的气压高度p，实现该台风实时风矢量的高度指定；
步骤S44、重复步骤S42和S43，完成S38步骤得到的两个时次原始云导风矢量中的所有台风实时风矢量的高度指定；
步骤S45、对完成高度指定的连续两个时次云导风矢量图层，以前一时次云导风矢量作为后一时次云导风矢量的参照，以要判定的风矢量所在位置为中心，计算其A×A邻域内云压变化的标准差；设定云压变化标准差的阈值，如果计算得到的云压变化的标准差差值大于该阈值，则认定该位置的实时风矢量位于复合云层，是云压检测不合格风矢量，将后一时次云导风矢量中位于该位置的风矢量剔除，得到经过云压检测控制的中间时相的台风立体风场；其中，A为正奇数；
步骤S46、以经过云压检测控制的中间时相的台风立体风场作为参照风场，对经过时间连续性质量控制的中间时相的台风实时风场中的风矢量进行判别和高度赋值处理：如果该风矢量在参照风场中有对应风矢量，则把对应风矢量的高度赋给该风矢量；否则，则认为该风矢量没有通过云压检测，予以剔除；从而得到经过时间连续性质量控制且满足云压检测的中间时相的台风立体风场；
步骤S47、进一步对步骤S46得到的中间时相的台风立体风场的风矢量进行空间连续性质量控制，分为风向控制和风速控制；每个风矢量均对其A×A邻域内的风矢量计算风向离差均方根和风速离差均方根，计算公式如下：
式中，Ds、Dd为目标风矢量的风速离差均方根和风向离差均方根；Vc、θc为目标风矢量的风速和风向；o＝1，2，……，B；Vo、θo为第o个邻域风矢量的风速和风向；B为邻域内的风矢量数；
分别设定风向离差均方根和风速离差均方根阈值，若计算得到的离差均方根中的任意一个超过对应的阈值，则认为该风矢量为不合格风矢量，予以剔除；对所有的风矢量进行判别，最终得到满足时空质量控制要求的中间时相的台风立体风场。
8.根据权利要求7所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，A为5。
9.根据权利要求1所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，所述的步骤S5的具体操作如下：
步骤S51、利用三维风场的风速插值公式对时空质量控制处理后的台风立体风场进行插值计算，得到海面10m高度风场；风速V10计算公式如下：
式中，V为台风立体风场实时风矢量的风速值；c为经验常数；pb为基准气压高度值；p为台风立体风场实时风矢量的气压高度；
步骤S52、根据海面粗糙度和海面平均高度，利用海面10m高度风场的风速V10换算海面风速Vs，计算得到台风期间的海面风场；计算公式如下：
式中，Vs为海面风速；Z0为海面粗糙度；Zs为海面平均高度。
10.根据权利要求1所述的一种基于遥感影像特征的台风椭圆型风场参数化模拟方法，其特征在于，所述步骤S6的具体操作如下：
步骤S61、在计算矢量格式的海面风场的同时，同步输出GeoTiff格式的重采样风速图像，图像的空间分辨率为原始影像空间分辨率的M倍，图像的像元值为海面风速值，图像的地理坐标和海面风场的地理坐标一致；
步骤S62、利用D∞算法，推算海面风场在不同方向上的最大风速半径分布；具体实现方法是以台风中心为中心，向周围任意方向作射线，识别与该射线相交的重采样风速图像上的所有像元，记录像元值最大的像元位置，为该方向的最大风速点，其与台风中心的距离即为该方向上的海面最大风速半径；
步骤S63、以台风中心为中心，以平面直角坐标系x轴正方向为起始方向，按等角度间隔，按照步骤S62的方法，逆时针在360度范围内计算各角度方向的最大风速点，得到海面最大风速半径在各方向的长度；对海面最大风速点进行极值过滤，去掉偏离台风中心的最大风速点；
步骤S64、采用基于最小二乘的椭圆拟合算法，拟合最大风速半径相关的椭圆，输出椭圆中心坐标、椭圆长轴和短轴长度、椭圆长轴与x轴正方向的夹角，给出台风引起的海面风场的最大风速分布范围。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010411137.1A CN111723464B (zh) | 2020-05-15 | 2020-05-15 | 一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010411137.1A CN111723464B (zh) | 2020-05-15 | 2020-05-15 | 一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111723464A true CN111723464A (zh) | 2020-09-29 |
| CN111723464B CN111723464B (zh) | 2024-04-09 |
Family
ID=72564482
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010411137.1A Active CN111723464B (zh) | 2020-05-15 | 2020-05-15 | 一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111723464B (zh) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113191299A (zh) * | 2021-05-14 | 2021-07-30 | 中国科学院大气物理研究所 | 一种涡旋的识别方法、装置、存储介质及电子设备 |
| CN113505107A (zh) * | 2021-05-26 | 2021-10-15 | 中国再保险(集团)股份有限公司 | 一种台风文件压缩方法及压缩系统 |
| CN113935533A (zh) * | 2021-10-20 | 2022-01-14 | 山东省气象科学研究所(山东省海洋气象科学研究所、山东省气象局培训中心) | 一种黄渤海海区大风推算方法 |
| CN114324973A (zh) * | 2022-03-17 | 2022-04-12 | 南方海洋科学与工程广东省实验室(广州) | 台风风速反演方法、装置、电子设备及存储介质 |
| CN114740550A (zh) * | 2022-06-14 | 2022-07-12 | 广东海洋大学 | 一种连续风暴事件智能识别预警方法及系统 |
| CN115204712A (zh) * | 2022-07-26 | 2022-10-18 | 中国气象局上海台风研究所(上海市气象科学研究所) | 一种海上和沿海风电场选址评估方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101933448A (zh) * | 2010-07-26 | 2011-01-05 | 北京师范大学 | 一种制作热带气旋风带的方法 |
| CN102132662A (zh) * | 2011-01-10 | 2011-07-27 | 北京师范大学 | 改进的制作热带气旋风带的方法 |
| CN106443830A (zh) * | 2016-06-16 | 2017-02-22 | 杭州师范大学 | 一种基于多源卫星数据的台风监测及评价监测精度的方法 |
| CN107526852A (zh) * | 2016-06-21 | 2017-12-29 | 中国辐射防护研究院 | 一种核设施事故场外后果实时在线评价方法及系统 |
| KR20200018127A (ko) * | 2018-08-10 | 2020-02-19 | 한국전자통신연구원 | 정지 궤도 위성 및 극 궤도 위성을 이용하여 태풍 주변 3차원 공간에 대한 수평 바람장 산출 방법 |
-
2020
- 2020-05-15 CN CN202010411137.1A patent/CN111723464B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101933448A (zh) * | 2010-07-26 | 2011-01-05 | 北京师范大学 | 一种制作热带气旋风带的方法 |
| CN102132662A (zh) * | 2011-01-10 | 2011-07-27 | 北京师范大学 | 改进的制作热带气旋风带的方法 |
| CN106443830A (zh) * | 2016-06-16 | 2017-02-22 | 杭州师范大学 | 一种基于多源卫星数据的台风监测及评价监测精度的方法 |
| CN107526852A (zh) * | 2016-06-21 | 2017-12-29 | 中国辐射防护研究院 | 一种核设施事故场外后果实时在线评价方法及系统 |
| KR20200018127A (ko) * | 2018-08-10 | 2020-02-19 | 한국전자통신연구원 | 정지 궤도 위성 및 극 궤도 위성을 이용하여 태풍 주변 3차원 공간에 대한 수평 바람장 산출 방법 |
Non-Patent Citations (1)
| Title |
|---|
| WANG XIUQIN 等: "An elliptical wind field model of typhoons", JOURNAL OF OCEAN UNIVERSITY OF CHINA, 30 April 2004 (2004-04-30), pages 34 - 39 * |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113191299A (zh) * | 2021-05-14 | 2021-07-30 | 中国科学院大气物理研究所 | 一种涡旋的识别方法、装置、存储介质及电子设备 |
| CN113505107A (zh) * | 2021-05-26 | 2021-10-15 | 中国再保险(集团)股份有限公司 | 一种台风文件压缩方法及压缩系统 |
| CN113505107B (zh) * | 2021-05-26 | 2023-11-10 | 中国再保险(集团)股份有限公司 | 一种台风文件压缩方法及压缩系统 |
| CN113935533A (zh) * | 2021-10-20 | 2022-01-14 | 山东省气象科学研究所(山东省海洋气象科学研究所、山东省气象局培训中心) | 一种黄渤海海区大风推算方法 |
| CN113935533B (zh) * | 2021-10-20 | 2022-09-02 | 山东省气象科学研究所(山东省海洋气象科学研究所、山东省气象局培训中心) | 一种黄渤海海区大风推算方法 |
| CN114324973A (zh) * | 2022-03-17 | 2022-04-12 | 南方海洋科学与工程广东省实验室(广州) | 台风风速反演方法、装置、电子设备及存储介质 |
| CN114740550A (zh) * | 2022-06-14 | 2022-07-12 | 广东海洋大学 | 一种连续风暴事件智能识别预警方法及系统 |
| CN115204712A (zh) * | 2022-07-26 | 2022-10-18 | 中国气象局上海台风研究所(上海市气象科学研究所) | 一种海上和沿海风电场选址评估方法 |
| CN115204712B (zh) * | 2022-07-26 | 2023-02-03 | 中国气象局上海台风研究所(上海市气象科学研究所) | 一种海上和沿海风电场选址评估方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111723464B (zh) | 2024-04-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111723464A (zh) | 一种基于遥感影像特征的台风椭圆型风场参数化模拟方法 | |
| CN114758252B (zh) | 基于图像的分布式光伏屋顶资源分割与提取方法及系统 | |
| CN111767865A (zh) | 一种利用航拍影像和激光数据反演红树林生物量的方法 | |
| CN108428220B (zh) | 静止轨道卫星序列遥感影像海岛礁区域自动几何校正方法 | |
| CN110555841B (zh) | 基于自注意图像融合和dec的sar图像变化检测方法 | |
| CN110569797B (zh) | 地球静止轨道卫星影像山火检测方法、系统及其存储介质 | |
| CN108195736B (zh) | 一种三维激光点云提取植被冠层间隙率的方法 | |
| CN112712535A (zh) | 基于模拟困难样本的Mask-RCNN滑坡分割方法 | |
| CN113066112B (zh) | 一种基于三维模型数据的室内外融合方法及装置 | |
| CN113033315A (zh) | 一种稀土开采高分影像识别与定位方法 | |
| CN116182805A (zh) | 基于遥感图像的国土测绘方法 | |
| CN108364279A (zh) | 确定静止轨道遥感卫星指向偏差的方法 | |
| CN111563408A (zh) | 多层次感知特征渐进自学习的高分辨率影像滑坡自动检测方法 | |
| CN114998545A (zh) | 一种基于深度学习的三维建模阴影识别系统 | |
| CN115439753A (zh) | 一种基于dem的陡峭河岸识别方法及系统 | |
| CN115272876A (zh) | 一种基于深度学习的遥感图像船舶目标检测方法 | |
| CN113902792A (zh) | 基于改进RetinaNet网络的建筑物高度检测方法、系统和电子设备 | |
| CN112906689B (zh) | 一种基于缺陷检测与分割深度卷积神经网络的图像检测方法 | |
| CN112166688B (zh) | 基于小卫星的沙漠与沙漠化土地监测方法 | |
| CN116663739A (zh) | 一种复杂地形风机出力预测方法、系统、设备和储存介质 | |
| CN114863258B (zh) | 海天线场景中基于视角转换检测小目标的方法 | |
| CN112668615B (zh) | 一种基于深度跨尺度外推融合的卫星云图预测方法 | |
| CN108563674A (zh) | 基于rs和gis的海域地理要素测量方法、系统及装置 | |
| CN114972997B (zh) | 一种基于全天空图像3d云层重建的跟踪式光伏发电优化方法 | |
| CN110176003A (zh) | 基于交比图像量测技术检测建筑物表面损伤面积的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |