基于TV_l1光流法的风云四号云导风估计方法及系统与流程

基于tv_l1光流法的风云四号云导风估计方法及系统
技术领域
1.本发明涉及卫星遥感技术领域，尤其涉及基于tv_l1光流法的风云四号云导风估计方法及系统。


背景技术：

2.云导风又叫云迹风，即通过云的运动推测出的风速。云导风算法是根据静止气象卫星拍摄的连续两幅或者三幅卫星云图，追踪图像上的同一云块(目标云或靶云)的运动，并计算这一特征云图的位移，同时计算出这一目标云所在的高度层次，从而得到这一云高层致使目标云移动的非常规大气风。由整个云图上所有的大气风组成的风矢量场，被称为“云迹风”。
3.研究表明，在天气预报、台风、暴雨等自然灾害预警等方面中，风场资料起着至关重要的作用。虽然目前获得风场的手段很多，无线电探空仪、风廓线仪、多普勒雷达等。但对于海洋、高原、沙漠等测站稀少或无测站地区，风场资料依旧匮乏，由于气象卫星覆盖面广，利用卫星观测和反演的云导风是这些地区主要的或者唯一的风场资料，它有效的弥补了全球陆地风场资料。云导风目前已经成为一种很重要的卫星产品，已有结果证明它在数值预报中的使用成果。早在1999年国外学者tomassini，kelly等指出在ecmwf模式下，云导风资料的引用提高了风场的准确度，改善了热带以及南半球的天气预报效果。庄照荣、薛纪善把云导风资料放到grapes三维变分同化分析场中不仅提高了风场质量而且对台风的路径和降水的预报也更准确。美国、日本气象卫星中心的云导风产品定期公开发行。
4.自世界第一颗静止气象卫星发射成功，就有人利用卫星云图反演云导风。最早是在上世纪70年代美国学者lzawa和fujita利用电影动画技术人为直观在可见光画面上判断云的移动，从而获得云导风场。但是这种人工手动识别算法完全依靠人的主观判断，不仅耗时费力，而且准确率很低。1971年endlich使用模式识别技术，通过提取云图中云团的特征量对云团进行匹配实现云导风的计算。同年，leese使用红外云图通过最大相关系数法计算，该方法主要是通过对连续的云图中大小相同的像素块进行匹配从而计算云导风。从此开启云导风的自动计算技术，由于相关系数法在搜索整个云图区域时，计算时间较长，leese后来在对云图预处理时引入快速傅里叶变换，提高了相关系数法的计算速度。经过数十年的发展，云导风的算法取得很大进展，相关系数法一直沿用至今而且是现在云导风算法中的主流算法。purdom j f w在研究快速云图导风时，发现运用相关系数法会产生“亚像素尺度位移问题”。针对这个问题，王振会开发的云导风客观导出系统(cwis)，创新地将傅里叶相位分析技术和相关系数分析法结合到云导风的计算中，而且提出“盲人摸路”的搜索方式。经过数年的发展，许建民，王振会，孙林等改进了这种方法提出tcfm系统提高了算法精确度，朱平将这种方法应用于半小时间隔卫星云图导风的计算中，结果证明傅里叶相位分析法和相关系数法相结合的方法很有实际应用价值，对示踪云的追踪精度有了进一步的提高。盲人摸路的搜索方式，是一种不同于其他国家的搜索技术，提升了计算机处理效率和风矢精度。龙智勇、石汉青等在传统的相关系数法基础上增加了图像灰度梯度特征的考虑。
王昌帅提出并行算法大大提高了云导风反演算法的速率。马侠霖、罗鹏针对灰度变化不明显问题，提出基于尺度不变特征变换，对图像灰度进行拉伸，再用相关系数法进行匹配，提高了匹配效果，但是同时也提高了计算复杂度。寿亦萱在研究模式识别在气象研究中的应用时，运用bp神经网络对云图的纹理进行分类，试验证明神经网络可以清楚的观察到一些利用常规资料不易发现的中尺度天气系统的结构特征，对云团的结构和属性有更好的识别，在此基础上运用相关系数求云迹风得到的风场。刘闵在研究中尺度云团的跟踪和预测中发现引入金字塔模型的光流法能更好的捕捉到卫星云图中像素的运动，比传统相关系数法具有更好的准确性，计算效率也更高。
5.目前的风云四号气象卫星的业务云导风产品主要面临三个问题：
6.(1)业务云导风产品采用5
×
5或者7
×
7的滑动窗区域作为最小搜素体，空间分辨率比较低(大约为20km
×
20km)；
7.(2)需要三幅图像计算一帧导风产品，时间分辨率与计算时效滞后；
8.(3)没有一套完备的产品质控方案。


技术实现要素：

9.本发明的目的在于克服现有技术缺陷，提出了基于tv_l1光流法的风云四号云导风估计方法及系统。
10.为了实现上述目的，本发明提出了一种基于tv_l1光流法的风云四号云导风估计方法，所述方法包括：
11.对风云四号气象卫星连续采集的两帧卫星云图进行预处理；
12.采用rof图像细节提取模型，对预处理后的每帧图像进行分解，分别得到对应的基本图像与细节图像；
13.对每帧细节图像建立金字塔结构，基于tv_l1光流法按照金字塔结构逐层迭代，计算得到两帧细节图像间的光流场；
14.由光流场结合地球投影模型计算风矢量；
15.根据标准大气温低垂直廓线表由云顶亮温值计算得到云顶高度；
16.结合云顶高度对风矢量进行质量检测，剔除不合格的风矢量，由卫星云图上所有的合格风矢量得到风矢量场。
17.作为上述方法的一种改进，所述对风云四号气象卫星连续采集的两帧卫星云图进行预处理；具体包括：
18.对每帧卫星云图进行通道分离处理，将卫星云图中不同通道的测值分离，每一个通道形成一幅圆盘图；再进行定位和定标处理得到该卫星云图预处理后的图像。
19.作为上述方法的一种改进，所述采用rof图像细节提取模型，对预处理后的每帧图像进行分解，分别得到对应的基本图像与细节图像；具体包括：
20.采用rof图像细节分割的模型，对预处理后的图像i进行分解，得到该图像的基本图像i
s
为：
[0021][0022]
其中，θ为太阳天顶角，为对基本图像i
s
进行梯度运算；
[0023]
由下式得到该图像的细节图像i
t
为：
[0024]
i
t
＝i
‑
0.95i
s
。
[0025]
作为上述方法的一种改进，所述对每帧细节图像建立金字塔结构，基于tv_l1光流法按照金字塔结构逐层迭代，计算得到两帧细节图像间的光流场；具体包括：
[0026]
根据尺度不变原理对每帧细节图像分别构造金字塔结构，包括l层，从粗尺度的第1层逐层迭代，直至到精细层第l层，迭代时只更新每层之间的运动场差异；得到每帧图像的光流场；
[0027]
基于tv_l1光流法计算前后两帧细节图像间的光流场。
[0028]
作为上述方法的一种改进，所述由光流场结合地球投影模型计算风矢量；具体包括：
[0029]
根据下式计算光流场矢量起点所在纬度为的地球半径r为：
[0030][0031]
其中，r
p
为地球的极地半径，ε为地球的扁率；
[0032]
根据下式计算光流场矢量起始位置和终点位置之间的地心角γ为：
[0033][0034]
其中，分别为光流场矢量起点和终点的纬度；δλ为光流场矢量起点和终点的经度λ0与λ1的差值；
[0035]
根据下式计算风速f为：
[0036]
f＝γ*r/δt
[0037]
其中，δt为两帧卫星云图之间的时间差；
[0038]
根据下式计算风向的大小d为：
[0039][0040]
根据λ0，λ1的大小判断风向，
[0041]
若λ1>λ0，风向dd＝d；若λ1≤λ0，风向dd＝360
°‑
d。
[0042]
作为上述方法的一种改进，所述根据标准大气温低垂直廓线表由云顶亮温值计算得到云顶高度；具体包括：
[0043]
根据云顶亮温值t，单位为k，在标准大气温低垂直廓线表中查找云顶亮温值t对应的云顶高度h，单位为hpa；
[0044]
如果查找不到，则判断云顶亮温值t的范围；
[0045]
如果t∈(216.14，288.14)，则从标准大气温低垂直廓线表中找到与t最相近的两个温度值t1，t2，并且t1>t，t2<t，由t1，t2分别得到对应的云顶高度h1和h2；
[0046]
如果t<216.64，则令t1＝216.64，t2＝228.43，从标准大气温低垂直廓线表中找到对应的云顶高度h1和h2；
[0047]
如果t>288.14，则令t1＝287.49，t2＝288.14；从标准大气温低垂直廓线表中找到对应的云顶高度h1和h2；
[0048]
根据温度值t1，t2和对应的云顶高度h1和h2，采用线性插值法计算得到云顶亮温值t对应的云顶高度h：
[0049][0050]
作为上述方法的一种改进，所述结合云顶高度对风矢量进行质量检测，剔除不合格的风矢量，由卫星云图上所有的合格风矢量得到风矢量场；具体包括：
[0051]
对于云图边界区求出的风矢量或者风矢量所在位置加上风矢量后超出云图拍摄的范围，则该风矢量不合格；
[0052]
当云顶气压高度低于100hpa或者高于1000hpa对应的风矢量不合格；
[0053]
当风矢量所在位置的云量为0，则该风矢量不合格；
[0054]
当亮温值缺失或者亮温值超出预设范围时，对应的风矢量不合格；
[0055]
以风矢量所在位置为中心，扩大到周围15*15像素块，如果15*15像素块中存在复合云层，则该风矢量不合格；
[0056]
当相邻两时刻风矢量水平分量加速度超过阈值时，该风矢量不合格；
[0057]
当相邻两时刻风矢量垂直分量加速度超过阈值时，该风矢量不合格；
[0058]
当风矢量的风速小于3m/s，该风矢量不合格；
[0059]
当相邻两时刻风矢量的风向改变超过30度，该风矢量不合格；
[0060]
当风矢量所处位置的云高不满足通道云高的范围，该风矢量不合格；
[0061]
当以风矢量所在位置为中心，扩大到周围5*5像素中的相邻两时刻平均云压分别为p1和p2，如果该区域相邻两时刻平均云压差的绝对值超过100hpa，该风矢量不合格；
[0062]
当太阳天顶角不属于预设的天顶角范围，该风矢量不合格；
[0063]
删除不合格的风矢量得到合格风矢量；由卫星云图上所有的合格风矢量得到风矢量场。
[0064]
一种基于tv_l1光流法的风云四号云导风估计系统，所述系统包括：预处理模块、图像分解模块、光流场计算模块、风矢量计算模块、云顶高度计算模块和质量检测模块；其中，
[0065]
所述预处理模块，用于对风云四号气象卫星连续采集的两帧卫星云图进行预处理；
[0066]
所述图像分解模块，用于采用rof图像细节提取模型，对预处理后的每帧图像进行分解，分别得到对应的基本图像与细节图像；
[0067]
所述光流场计算模块，用于对每帧细节图像建立金字塔结构，基于tv_l1光流法按照金字塔结构逐层迭代，计算得到两帧细节图像间的光流场；
[0068]
所述风矢量计算模块，用于由光流场结合地球投影模型计算风矢量；
[0069]
所述云顶高度计算模块，用于根据标准大气温低垂直廓线表由云顶亮温值计算得到云顶高度；
[0070]
所述质量检测模块，用于结合云顶高度对风矢量进行质量检测，剔除不合格的风矢量，由卫星云图上所有的合格风矢量得到风矢量场。
[0071]
与现有技术相比，本发明的优势在于：
[0072]
本发明采用tv_l1光流法进行风云四号云导风的估计，可以克服目前业务最大相关法云导风产品时间与空间分辨率差的缺点，提升产品时空分辨率的同时提高了云导风产品的精度。
附图说明
[0073]
图1是本发明的基于tv_l1光流的云导风反演流程图；
[0074]
图2(a)是风云四号气象卫星图像分解的图像基本部分；
[0075]
图2(b)风云四号气象卫星图像分解的细节部分；
[0076]
图3是风云四号云导风与基于tv_l1光流的云导风结果对比，其中，图3(a)是目前业务风云四号云导风产品；图3(b)是基于tv_l1光流的云导风产品。
具体实施方式
[0077]
下面结合附图和实施例对本发明的技术方案进行详细的说明。
[0078]
实施例1
[0079]
如图1所示，本发明的实施例1提供了基于tv_l1光流法的风云四号云导风估计方法。
[0080]
风场估算从概念上说是简单的，但是在实践中不容易做准确。为了把风矢量计算好，另外两个环节的工作也十分重要，这两个环节就是数据预处理和质量控制。因此卫星风的基本算法由五部分组成：数据预处理、光流计算、光流到风场的转换、示踪图像块的高度指定和质量控制。
[0081]
步骤1)数据预处理
[0082]
数据预处理进行敏感通道选取，定位、定标等文件数据提取，质量检验和容错处理等工作。敏感通道选取是从原始数据中选取对导风敏感的数据光谱通道，便于在后续计算中数据的反复调用。定位、定标工作的质量对卫星风的质量有非常大的影响。没有高质量的定位、定标预处理工作，卫星风的质量是无法保证的。
[0083]
步骤2)基于tv_l1的光流计算
[0084]
首先，光流法基于灰度不变原理，认为图像灰度不会发生变化，只会在空间位置上发生移动，如式(1)所示。
[0085]
i0(x,y,t)＝i0(x δx,y δy,t δt)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0086]
式中，i0表示图像，x，y表示图像水平或者垂直方向坐标，t表示时间，δx，δy分别为图像水平或者垂直方向坐标偏移量，δt为时间间隔。
[0087]
由于式(1)是欠定方程，需要正则化才能求解，式(2)在式(1)的基础上加入全变差域正则化项(tv_l1)。
[0088][0089]
式中，u
h
，u
v
分别为图像水平或者垂直方向坐标偏移量,为梯度算子，λ为数据惩罚项与正则化平衡因子。
[0090]
由于实际遥感图像并不能满足灰度恒定不变这样的假设，如可见光图像随着太阳高度强度有变化。针对这个问题本技术采用了rof图像细节分割的模型，将图像进行图像基本与图像细节分解。计算方法如式(3)与式(4)所示，见图2(a)、(b)。
[0091][0092]
式(3)中，i为原始图像，i
s
为图像基本部分。
[0093]
i
t
(x)＝i(x)
‑
0.95i
s
(x)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0094]
式(4)中，t为图像细节部分。
[0095]
另外，求解式(2)还有个缺陷，解会集中在图像的小的移动，对于大的移动会产生偏差。本技术根据尺度不变原理(scale invariant feature transform,sift)构造出金字塔结构，从粗尺度逐层迭代到精细层，迭代时候只需更新每层之间的运动场差异。这样解决了式(2)对于大尺度运动估计误差的问题。
[0096]
基于金字塔结构的tv_l1光流计算如下：
[0097]
算法1:基于金字塔结构的tv_l1光流数值计算
[0098]
输入:前后时刻结果rof分解的两帧细节图像i0以及i1[0099]
输出：从i0到i1的光流场u
[0100][0101][0102]
涉及到的公式与矩阵如下所示：
[0103]
水平方向微分滤波与垂直方向微分滤波分别为：
[0104][0105]
图像经过水平方向微分滤波与垂直方向微分滤波分别为(p＝0,1)：
[0106][0107]
图像的变形根据的偏微分方程：
[0108][0109]
式中，α为平衡因子，一般取值为0.5。
[0110][0111]
式中i
t
＝bi(i0,u)
‑
i1，bi为双线性插值；f
x
以及g
x
分别定义为f
x
＝|i
t
i
h
du
h
i
v
du
v
|2以及令i
h
＝diag(i
h
),i
v
＝diag(i
v
),ψ
′
＝φ’d
(f
x
),以及φ
′
＝φ’d
(f
x
)为对角矩阵,
[0112][0113]
式中，式中，以及
[0114]
步骤3)光流到风场的转换
[0115]
根据图像位移以及地球投影模型计算风矢量的公式如下：设λ0，λ1分别为起点和终点时示踪图像块的纬度和经度，δλ和δt分别为经度差和时间差，r
p
为地球的极地半径，ε为地球的扁率，为图像块所在纬度。
[0116][0117]
为图像块所在纬度地球的半径，那么
[0118][0119]
就是图像块起始位置和终点位置之间的地心角。风速f为
[0120]
f＝γ*r/δt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0121]
风向的大小d为
[0122][0123]
若λ1>λ0，风向dd＝d；若λ1≤λ0，风向dd＝360
°‑
d。
[0124]
步骤4)示踪图像块的高度指定
[0125]
原理上，云导风的计算是通过反演云顶亮温值得到的，所以这里指的云高就是云顶的高度。目前常见的方法是esa提出的利用水汽和红外双通道定高法，二氧化碳多通道法等。选择11.2um红外云图的云顶亮温值根据大气温度垂直廓线来推算云顶高度。
[0126]
设云顶亮温值t(k)，对应的云高为h(hpa)，云顶高度计算如下：
[0127]
(1)如果云顶亮温值t，可以在表中找到与之相等的温度值，那么表中与亮温值相对应的就是所求的云高h；
[0128]
(2)如果在表中找不到相等的温度值，首先确认这个亮温值是否超过了表中所给
温度值的最值，也就是说是216.64和288.14之间：
[0129]
a.如果t∈(216.14，288.14)，那么找到表中与t最相近的两个温度值，并且这两个亮温值一个要大于t，一个要小于t。假设t1>t，t2<t，相对应的云高压为h1和h2；
[0130]
b.如果t<216.64，那么令t1＝216.64，t2＝228.43；
[0131]
c.如果t>288.14那么令t1＝287.49，t2＝288.14；
[0132]
利用线性插值法：
[0133][0134]
求得t对应的云高压值。
[0135]
表1标准大气温低垂直廓线一览表
[0136][0137][0138]
步骤5)质量监控
[0139]
质量监控是云导风计算中的关键部分，质量监控是提高云导风准确率的重要因素之一。由于云区的时空变化以及在计算过程中可能存在的错误或者限制，使部分风矢不合理。利用客观规律或者主观判断来调整或者删除认为不合理的风矢，就是云导风的质量监控。
[0140]
为了输出高质量的云导风，通过质量监控去掉不合理的风矢，例如，同一经纬度的风矢在很短的时间间隔其大小和方向变化应该在一定范围之内；某一区域内风矢的变化也
应该在一定范围之内；若某一区域有复合云层，那么这一区域的云导风不具有参考价值，所以需要去掉这些风矢。具体的标识码如下表。
[0141]
表2云导风质量标识码
[0142][0143][0144]
设上一时刻的风场作为这一时刻风场质量监控的参照，我们假设上一时刻风速speed1,风向dir1,云压p1，这一时刻的风速speed2,风向dir2，云压p2,如果满足表中任意一个条件，本文就认为这个风矢不合理，删除并做相应的标记。下面是对表4.2质量标识码的详细解释：
[0145]
(1)质量标识码1：云图边界区求出的风矢或者风矢所在位置加上风矢后超出云图拍摄的范围，那么认为此风不合理，删除并令flag＝1；
[0146]
(2)质量标识码2：当云压低于100hpa或者高于1000hpa都认为是无云区，所以这个风速认为不合理(p2<100||p2>1000)，删除并令flag＝2；
[0147]
(3)质量标识码3：如果风矢所在位置的云量为0，那么flag＝3；
[0148]
(4)质量标识码4：此处亮温值缺失或者亮温值超出可行范围，如可见光的亮温范围1到200，红外范围150到340，flag＝4；
[0149]
(5)质量标识码5：以风矢所在位置为中心，扩大到周围15*15像素块，如果这以像素块中存在复合云层，那么认为风矢不合理，删除并令flag＝5；
[0150]
参照美国goes
‑
r云导风系统对复合云层检测的方法，以15*15像素块亮温值为研究单位，在这个区域选择3*3子区域，分别计算这些子区域的亮温均值和方差，采用聚类分析的方法识别某一亮温均值 /
‑
3倍的平均方差范围内包含点的个数，如果没有复合云层，那么其中会有两个最大聚类总包含所有子区域80％以上的点。如果出现第三个峰值，就含有复合云层。
[0151]
(6)质量标识码6：当前后两个时刻云导风水平分量加速度超过阈值时，认为风速不合理，在这里可见光通道阈值为5m/s，红外通道阈值10m/s，flag＝6；
[0152]
(7)质量标识码7：当前后两个时刻云导风垂直分量加速度超过阈值时，认为风速不合理，在这里可见光通道阈值为5m/s，红外通道阈值10m/s，flag＝7；
[0153]
(8)质量标识码8：即满足6又满足7；
[0154]
(9)质量标识码9：剔除风速太小的风矢，如果所求风速小于3m/s，flag＝8；
[0155]
(10)质量标识码10：当前后两个时刻的云导风风向改变超过30度时flag＝9；
[0156]
(11)质量标识码11：由于针对时间所处的状态白天还是晚上，云压高于600hpa或者低于600hpa，都设置不同的通道来计算云导风，所以不同通道计算的云高有一定范围，如果云导风所处位置的云高不满足通道云高的范围，就删除风矢，并做标记。
[0157]
(12)质量标识码12：当以风矢量所在位置为中心，扩大到周围5*5像素中的相邻两时刻平均云压分别为p1和p2，如果该区域相邻两时刻平均云压差的绝对值超过100hpa，标记flag＝11；
[0158]
(13)质量标识码13：当太阳天顶角不属于某一范围，即不属于白天也不属于晚上时，认为风矢不合理，删除并做相应的标记；
[0159]
(14)当质量标识码标记为0时，即风矢经过了所有质量检测，属于合格风矢。
[0160]
效果对比
[0161]
图3显示了业务风云四号云导风与基于tv_l1光流的云导风结果对比图。可以看出业务风云四号云导风分辨率粗糙，并且在台风眼区旋转较快的区域云导风产品误差较大。而基于tv_l1光流的云导风产品分辨率高，可以清晰找到台风眼区的旋转中心。
[0162]
实施例2
[0163]
本发明的实施例2提出了一种基于tv_l1光流法的风云四号云导风估计系统，该系统包括：预处理模块、图像分解模块、光流场计算模块、风矢量计算模块、云顶高度计算模块和质量检测模块；具体处理同实施例1，其中，
[0164]
所述预处理模块，用于对风云四号气象卫星连续采集的两帧卫星云图进行预处理；
[0165]
所述图像分解模块，用于采用rof图像细节提取模型，对预处理后的每帧图像进行
分解，分别得到对应的基本图像与细节图像；
[0166]
所述光流场计算模块，用于对每帧细节图像建立金字塔结构，基于tv_l1光流法按照金字塔结构逐层迭代，计算得到两帧细节图像间的光流场；
[0167]
所述风矢量计算模块，用于由光流场结合地球投影模型计算风矢量；
[0168]
所述云顶高度计算模块，用于根据标准大气温低垂直廓线表由云顶亮温值计算得到云顶高度；
[0169]
所述质量检测模块，用于结合云顶高度对风矢量进行质量检测，剔除不合格的风矢量，由卫星云图上所有的合格风矢量得到风矢量场。
[0170]
创新点：
[0171]
tv_l1光流法进行风云四号云导风的估计，解决目前业务最大相关法的时空分辨率粗糙的问题，本身就是创新点。
[0172]
针对实际tv_l1光流法估计的灰度不变假设实际可能不符合的问题，创新地使用图像的细节作为输入，有效地解决太阳高度角对风云四号图像的影响。另外，针对大尺度运动估计偏差大的难题，根据尺度不变原理采用金字塔结构，逐层估计，有效地解决了该问题。
[0173]
需要说明的是，本发明不仅可以应用于风云四号气象卫星云导风估计，还可以应用于其他新一代静止轨道气象卫星的云导风估计。例如日本的葵花8号，葵花8号成像仪数据空间分辨率为2公里，时间分辨率为10分钟，其他步骤同实施例1。
[0174]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。