CN103679734A - 基于svm和pde的有眼台风二维表面风场反演方法 - Google Patents

Abstract

本发明主要针对常用算法求解红外云图有眼台风二维表面风场误差较大的特点，提出一种结合SVM和PDE的新算法对有眼台风二维表面风场进行反演方法。首先对气象局提供的红外云图结合PDE算法进行台风眼壁提取，得到眼壁各点的坐标值和灰度之后，然后利用SVM进行平均灰度和台风近中心最大风速的建模，最后对于任意给定的红外云图经过分割眼壁，提取参考点灰度作为模型输入得出参考点风速，利用线性插值反演整个二维表面风场。经过实验证明，由本发明求得的平均误差为2.9m/s，性能优于线性回归法和径向基神经网络法，为台风强度预报提供了一个有力的依据。

Description

基于SVM和PDE的有眼台风二维表面风场反演方法

技术领域

本发明属于台风强度预报应用领域。尤其涉及一种基于SVM和PDE的有眼台风二维表面风场反演方法。

背景技术

台风风场的信息相当重要，首先它是台风强度的体现，如近中心最大风速一般可以表示为台风的强度；其次它还能反映台风的结构信息，结构也是决定台风强度的因素之一；因此，提高强度预报水平的关键之处就在于怎样选取适合的技术从已知的信息中提取风场信息，并结合适当的数值预报进行强度预报。台风强度一般来说主要由近地面眼壁附近最大风速所决定，因此，低层风场的反演在整个台风强度预报过程中显的尤为重要。美国联合台风警报中心（JTWC）从1987年因经费缺乏而停止使用侦查飞机辅助后，开始使用卫星（如静止卫星、SSM/I、TRMM、QuickSCAT等）图像和遥感资料结合Dvorak客观技术(ODT)对台风低层风场进行反演。

低层风场可以分为台风内核风场和台风外围风场。目前，基于卫星微波散射计的反演法一般在低风速和低降水的环境中可以得到最好的效果，因而这类方法适合于估计远离台风眼墙（强风和强降水区）的台风外围风场。基于卫星被动微波设备，像特殊传感器微波成像仪（SSM/I）常用于估计开阔水域的表面风场，但也只能限用于台风外围风场反演。类似的，基于静止卫星云迹风也只可以反演台风外围风场，因为内核的卷云具有模糊效应。台风外围风场更多的信息是从经常飞行在西大西洋盆地上的湾流IV喷气飞机的下投式探空仪的全球定位系统（GPS）获得。但是基于安全考虑，小型高海拔喷气飞机避免进入台风内核区域。

虽然台风外围风场信息来源很多，但是台风内核风场信息目前只能通过低空飞机侦测获得，例如国家海洋和大气管理局（NOAA）的WP-3D,美国空军的WC-130等。目前对于台风内核风场的反演主要采用基于红外卫星数据和借助统计学中的多元线性回归法对飞行层的风速进行分析。虽然线性回归法算法在很多情况下简单有效，但是其拟合误差通常较大，因此本发明利用径向基神经网络算法代替线性回归法的方法，建立有眼台风内核风速和云图灰度的模型，以提高其建模精度，更加准确地刻画台风内核风场，为后续台风内核风场反演奠定基础。

发明内容

本发明的目的主要是针对常用算法求解红外云图有眼台风二维表面风场误差较大的特点，提出一种结合SVM和PDE的新算法对有眼台风二维表面风场进行反演。首先对气象局提供的红外云图结合PDE算法进行台风眼壁提取，得到眼壁各点的坐标值和灰度之后，利用SVM进行平均灰度和台风近中心最大风速的建模，最后对于任意给定的红外云图经过分割眼壁，提取参考点灰度，结合模型得到参考点风速，利用线性插值反演整个二维表面风场。

本发明的技术方案所述程序依次执行以下步骤如下：

1）在气象局提供的台风数据库中选取所有可用红外云图及其近中心最大风速的数据。

2）对步骤1中选中的图像利用PDE进行台风眼壁的分割。

台风中心位置的确定对于天气分析和台风预测具有重要的意义，但在实际工作中，台风中心位置作为台风重要特征主要依靠人工方式进行确定，因此，如何利用卫星云图，自动、准确地分割眼壁显得至关重要。基于静止卫星云图的台风眼区提取从本质上来讲属于图像分割的范畴。图像分割技术是图像处理、分析与理解、图像识别和计算机视觉领域的一项基本而又关键的技术。常见的图像分割有：

（1）基于边缘检测的图像分割方法，图像分割方法中最常用的是基于边缘检测的图像分割法。

（2）基于灰度特征的阈值分割方法。

（3）基于区域的分割方法。

（4）基于分水岭的图像分割方法。

根据台风眼的轮廓特点，本发明采用测地活动轮廓模型（GAC模型）的偏微分方程（PDE）台风云图分割法。偏微分方程图像分割技术是二十世纪八十年代产生并逐渐发展起来的一种非线性图像分割方法。其基本思想是利用动力学模型的思想，通过定义一条初始化曲线，在曲线本身的内力和图像数据所构造的外力作用下，驱使曲线向目标轮廓移动，从而获得目标物体之特征。偏微分方程图像分割的最大优点是曲线在演化的过程中始终保持连续性和光滑性，可以实现连续、封闭的目标轮廓提取，这是传统的分割方法所不能直接实现的问题。偏微分方程的图像分割技术作为一种比较新颖且有效的图像分割方法，为解决传统图像分割技术诸多困惑提供了一个很好的途径。本发明中欲提取的台风眼区轮廓一般来说是不规则的圆形，鉴于后续基于所提取的台风眼壁处的风速拟合要求所提取的眼区轮廓必须连续而且封闭，而GAC模型恰好符合本发明的要求。

1987年，Kass,Witkin和Terzopoulos从全新角度，直接用连续曲线模型来定位图像边缘，提出了Snakes活动轮廓模型（参数模型），该模型在图像某个区域的附近定义一条带能量的初始化演化曲线（显示表达），在自身内力和图像信息产生的外力的作用下不断地运动，最后收敛于目标物体的边缘。该模型首次引进了变分法，其局部极值组成可供高层视觉处理的图像分割结果，这样在寻找显著图像特征时，即可通过将图像特征推向一个适当的局部极值点来与模型进行交互。1988年，R.Osher和J.A.Sethian提出水平集方法（Level Set Method），使得活动轮廓模型进一步完善。与参数模型不同，水平集方法并不跟踪不同时刻曲线运动情况，而是在固定坐标系中更新不同时刻下的水平集函数来模拟曲线的演化。该方法大大拓宽了Snakes模型的应用范围，活动轮廓模型的理论得到了很大的发展。1997年，V.Caselles,R.Kimmel和G.Sapiro三人共同研究提出了测地线活动轮廓模型（Geodesic Active Contour,GAC）^[7],GAC模型的提出所产生的的影响是巨大的，引发了图像分割技术的新革命，对于图像处理技术的进步具有重大推动作用。

测地活动轮廓模型（GAC）简介:

活动轮廓模型是测地活动轮廓模型的基础。由Kass提出的活动轮廓模型，其主要思想是将图像分割问题转化为最小化一个封闭曲线C(p)的“能量”泛函：

$E\left[C\left(p\right)\right]=\mathrm{\α}\underset{0}{\overset{1}{\∫}}\left|C_p\left(p\right)\right|\mathrm{dp}+\mathrm{\β}\underset{0}{\overset{1}{\∫}}{\left|C_p\left(p\right)\right|}^2\mathrm{dp}-\mathrm{\λ}\underset{0}{\overset{1}{\∫}}|\▿I[C\left(p\right)\left]\right|\mathrm{dp}---\left(1\right)$

其中第一项的积分是曲线的欧几里得弧长，第二项

$\underset{0}{\overset{1}{\∫}}{\left|C_p\left(p\right)\right|}^2\mathrm{dp}=\underset{0}{\overset{1}{\∫}}{\left|\mathrm{\κ}\right|}^2{\left(\frac{\mathrm{dp}}{\mathrm{ds}}\right)}^2\mathrm{dp}---\left(2\right)$

表示曲线“振荡”的能量。那么，最小化式(1)就需要要最小化前两项和最大化第三项。最小化前两项就得使得闭合曲尽量短、尽量光滑，因为在曲线逐渐变短的过程中也自然地趋向于光滑，所以式(1)第二项忽略不计。最大化第三项即要求曲线C尽量地与图像梯度

取到极大值的位置一致。从以上原理可以发现，活动轮廓模型的原理与基于边缘的图像分割算法很接近。不同于基于边缘的分割算法的是，基于活动轮廓模型的图像分割，不会出现边缘断裂或者过分割现象。

鉴于式(1)第三项的负号给计算带来的麻烦，现引进辅助函数g(r),r∈R⁺,则式(1)变为:

$E\left[C\left(p\right)\right]={\mathrm{\α}}_1\underset{0}{\overset{1}{\∫}}\left|C_p\left(p\right)\right|\mathrm{dp}+a_2\underset{2}{\overset{1}{\∫}}g(\▿I[C\left(p\right)\left]\right)\mathrm{dp}---\left(3\right)$

观察式(3)，发现E[C(p)]主要与曲线的形状、位置和参数p的变化有关。参数p的随意性大大影响了E的结果，为了避免这种情况，将活动轮廓模型改为使用无自由参数p的测地活动轮廓模型。为使E[C(p)]最小，需要曲线通过

的局部最小值，这个思想类似于光学中费马定理，于是能量泛函式为：

$L_R\left(C\right)=\underset{0}{\overset{L\left(C\right)}{\∫}}g\left(\right|\▿I\left[C\left(s\right)\right]\left|\right)\mathrm{ds}---\left(4\right)$

其中L(C)是闭合曲线C的欧几里得弧长，L_R(C)是经过加权后的弧长。要最小化式(4)，相应的梯度下降流为：

$\frac{\∂C}{\∂t}=g\left(C\right)\mathrm{\κ}\stackrel{\→}{N}-(\▿g\·\stackrel{\→}{N})\stackrel{\→}{N}---\left(5\right)$

其中前一项中

表示的是平均曲率运动，为曲面的单位法向量。但是式(5)存在一个明显的缺陷，即它只适合分割凸的对象，而不适合用来分割凹形对象。于是，将式(5)改进如下：

$\frac{\∂C}{\∂t}=g(c+\mathrm{\κ})\stackrel{\→}{N}-(\▿g\·\stackrel{\→}{N})\stackrel{\→}{N}---\left(6\right)$

其中参数c为一可调参数，通常与演化的速度有关。

在使用数值方案之前，首先要计算出图像的灰度梯度，才能得到边缘函数

梯度算子的作用类似于微分，因此会使图像尖锐化，计算梯度的过程对噪声非常敏感。为防止噪声的干扰，可先对图像进行预处理。高斯平滑是最常用的预处理方法，公式如下：

I_σ(x,y)＝I(x,y)*g_σ(x,y)  (7)其中σ表示高斯方差。

接下来便要选取一种合适的边缘函数g了，一般可令边缘函数为

$g\left(r\right)=\frac{1}{1+{(r/k)}^p},p=\mathrm{1,2}---\left(8\right)$

式中k是待选的常数，边缘函数g的下降速度由它决定。有了边缘函数的定义后，将g(r)中的r用图像中的各个像素的梯度大小带入后，即可得到：

$g(x,y)=g\left(\right|\▿I\left|\right)----\left(9\right)$ 这样就可以开始使用数值方案来实现图像分割操作了，我们采用水平集的数值方法。

曲线演化的一般方程式为：

$\frac{\∂C}{\∂t}=\mathrm{\β}\stackrel{\→}{N}---\left(10\right)$

曲线演化水平集方法的基本方程式为：

$\frac{\∂u}{\∂t}=\mathrm{\β}|\▿u|---\left(11\right)$

对比式（6）和式（10）可知，

代入式（11）

$\frac{\∂u}{\∂t}=\mathrm{\β}|\▿u|=[g(c+\mathrm{\κ})+\▿g\·\frac{\▿u}{|\▿u|})\left]\right|\▿u|=\mathrm{cg}|\▿u|+\mathrm{div}\left(g\frac{\▿u}{|\▿u|}\right)|\▿u|---\left(12\right)$

由式(5)可知：

$\frac{\∂u}{\∂t}=\mathrm{div}\left(g\frac{\▿u}{|\▿u|}\right)|\▿u|---\left(13\right)$

由于图像以离散的形式存储，而偏微分方程描述的是连续的形式，因此在对式(12)和(13)作数值计算前，采用“半点离散化”方案离散化这两式中的散度算子。得到：

$s_{\mathrm{ij}}=\mathrm{div}{\left(g\frac{\▿u}{|\▿u|}\right)}_{\mathrm{ij}}\≈C_{1,\mathrm{ij}}{g}_{1,\mathrm{ij}}(u_{i,j+1}-u_{\mathrm{ij}})-C_{2,\mathrm{ij}}{g}_{2,\mathrm{ij}}(u_{i,j}-u_{i,j-1})+C_{3,\mathrm{ij}}{g}_{3,\mathrm{ij}}(u_{i+1}-u_{\mathrm{ij}})-C_{4,\mathrm{ij}}{g}_{4,\mathrm{ij}}(u_{i,j}-u_{i-1,j})---\left(14\right)$

最后采用一种双曲型方程的迎风差分方法^[8]，对离散化的式(12)的偏微分方程进行数值计算。

得到显示方案：

$u_{\mathrm{ij}}^{n+1}=u_{\mathrm{ij}}^n+\mathrm{\Δt}(\mathrm{cg}{\▿}_{\mathrm{ij}}^{(-)}+\max (s_{\mathrm{ij}},0){\▿}_{\mathrm{ij}}^{(-)}+\min (s_{\mathrm{ij}},0){\▿}_{\mathrm{ij}}^{+})---\left(15\right).$

3）利用SVM建立灰度风速模型

支持向量机的原理：

机器学习的任务是通过分析已知的有限样本来估计系统输入变量和输出变量之间的关系，并对其他给出的输入变量，做出输出变量的预测。现假设输入变量是X，输出变量是Y，并已知X与Y之间有某种未知的关系，可认为它们遵循某个联合概率F(X,Y)。若已给定了l个独立的样本并且集中于函数组{f(X,W)}中，求一个能使其期望风险达到最小的最优函数解{f(X,W₀)}，这样即可估计出X和Y之间具有的关系，即：

minR(W)＝L∫(Y,f(X,W))dF(X,Y)   (16)式中，W是函数广义参数,而L(Y,f(X,W))是用f(X,W)对Y进行预测所形成的误差，Y为输出变量（假设函数为单值连续函数），误差估计函数可定义为：

L(Y,f(X,W))＝(Y-f(X,W))²   (17)即采用平方和最小准则，当然也可选用别的误差估计函数，如高斯误差函数等。

在一般的应用中，样本数目不可能是无限的，于是人们提出了经验风险最小化准则，也就是由样本来定义经验风险，表达式如下：

$R_{\mathrm{emp}}\left(W\right)=\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}L(Y,f(X_i,W))---\left(18\right)$

伴随着统计学习理论的发展，出现了实际风险与经验风险不一致的情况，因此有必要定义一些函数集学习性能的指标来描述实际与经验风险之间的关系。其中最有名的是用Vapnik和Chervonenkis两人名字首字母命名的VC维概念，VC理论证明了实际风险R(W)和经验风险R_emp(W)以概率1-η满足下式：

$R\left(W\right)\≤R_{\mathrm{emp}}\left(W\right)+\sqrt{\frac{h(\ln (2n/h)+1)-\ln (\mathrm{\η}/4)}{n}}---\left(19\right)$

上式又可简化为：

(19)式中，h是函数集的VC维数；l是已知的有限样本数；η是满足0＜η＜1的参数；

称为置信区间，于是式(20)称为结构风险式。由式(20)可知，为了使结构风险获得最小值，则要求式(20)中的右边两项之和取得最小值。

实际问题中人们遇到更多的是非线性或者多维的情况，为了能解决这类负载的问题，需要引入函数映射的思想，也即在支持向量机的使用过程中加入核函数。核函数的作用是将低维非线性空间中的复杂运算映射到某个高维线性空间，在高维空间中处理可相对简单，计算完毕再返回到低维空间，在高维空间中一般是计算内积和，从而可避免维数灾难对计算过程的影响，高维线性空间中构造支持向量机的思想与一般的线性回归类似。常见的核函数有以下几种：

(1)线性核，即两个向量的内积：

K(X,X_i)＝＜X,X_i＞

(2)多项式核：K(X,X_i)＝(＜X·X_i＞+1)^d，其中d为多项式核的阶数。

(3)径向基核（RBF），传统的径向基核使用以下判定准则，即

$f\left(X\right)=\mathrm{sgn}(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i{K}_{\mathrm{\γ}}\left(\right|X-X_i\left|\right)-b))$

其中，K_γ(|X-X_i|)是一个非负的单调函数，其大小取决于两向量的距离|X-X_i|，若要使它趋于0，

则训练的样本总数要趋于无穷大。这个判定准则需要顾及众多参数，包括：γ的值、X_i的数目n、以及α_i的值。最常见判定准则是高斯函数，即：

其中σ是核函数的宽度。高斯函数与经典径向基函数的主要区别在于中心点X_i和输出权值均由支持向量机算法本身决定，而每个中心点对应着下一个支持向量。

(4)多层感知机核，满足Mercer条件的Sigmoid函数为K(X,X_i)＝tanh(γ_?＜X,X_i＞+c)，其中γ是梯度参量；c为偏离参量。采用Sigmoid函数作为内积后，则使得此时的SVM包含了一个单多层感知机，其中的隐含层节点数可以由算法自动确定。

支持向量机算法的实现：

根据已知样本空间，先选择合适的核函数

将低维非线性空间中的Xi映射到高维线性空间的

接着在高维空间的内积计算过程中，使用核函数K(Xi,Xj)来替换内积运算

从而实现非线性回归。于是问题求解的形式可转化为一个最优二次规划问题，即目标函数为：

$\min \frac{1}{2}{\left|\right|w\left|\right|}^2---\left(21\right)$

约束条件为：

引入Lagrange算子后，即可得到该问题的对偶形式：

$\min w\left(a\right)=\frac{1}{2}\underset{i,j=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)(a_j^{*}-a_j)K(X_i,X_j)-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{Y}_i(a_i^{*}-a_i)+\mathrm{\ϵ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}+a)---\left(23\right)$

约束条件为：

$\begin{array}{c}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)=0\\ a_i^{*},a_i\≥0,i=\mathrm{1,2},...,l\end{array}---\left(24\right)$

a_i≠0对应的向量即为支持向量，要得到这两项，只需求解由式(23)、(24)构成的Q-P问题。

最终我们得到了回归机的解答式：

$Y=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+b$

基于支持向量机的回归算法的一般步骤如下：

Step1：给出已知的训练样本集

其中

Step2：根据问题的需求选择合适的核函数K(X,X_i)，并设定精度范围ε；

Step3：求解受式（24）约束的式（23）的最优二次规化问题，求得

Step4：在

中选择任意一个不为0的分量，即支持向量，并计算

$b=Y_j-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+\mathrm{\ϵ},$ 当 $a_i\≠0,b=Y_j-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)-\mathrm{\ϵ},$ 当

Step5：构造非线性超平面 $Y=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+b$

拟合方法的选择：

数据获得后便可进行数据拟合，使用循环的模式，总共有多少个样本就循环多少次，每次循环只取一个样本用来作预测样本，其他的样本用作训练样本，并计算相应的误差，一次循环计算一次误差大小，每次循环所用的预测样本都要求不同，于是所有循环结束后，每个样本的预测误差都可以描绘出。

4）给定一图像及其中心经纬度和估算点经纬度。

5）对于给定图像进行分割眼壁得到参考点灰度。

6）将参考点灰度作为模型的输入值，求出参考点的风速。

7）利用线性插值求出图像估算点风速。

云图上任意两点的经纬度计算实际距离的基本原理如下：

设A、B两点的纬度分别为A和B，经度分别为la和lb，要求的是A、B两点间的大圆弧长。令A、B的中点的纬度为C，经度为lc，A、B两点间的地心夹角为θ，则根据球面三角关系，可列出如下方程：

cos(θ)＝sin(A)sin(B)+cos(A)cos(B)cos(la-lb)  (25)

cos(θ/2)＝sin(A)sin(C)+cos(A)cos(C)cos(la-lc)  (26)

cos(θ/2)＝sin(B)sin(C)+cos(B)cos(C)cos(lb-lc)  (27)由于C、lc、θ为未知量，所以式(26)，(27)为非线性方程，于是求解中点的问题转化为如何求解这个非线性方程组的问题。由于此方程组无法得到符号解，只能应用迭代法解出近似的数值解，使用计算机求解可以得到很高的精度。

设f(x)＝0为非线性方程，f为[a,b]上二阶可导函数，且f'(x)、f″(x)均不为零，不妨设f'(x)＞0，f″(x)＞0即f为[a,b]上递增凸函数，又设f(a)＜0，f(b)＞0，那么由连续函数的介值定理知，方程f(x)＝0在[a,b]内存在唯一解ζ。牛顿迭代法的基本思想是构造一点列{X_n}，使

且当n充分大时，可以将X_n作为ζ的近似值，即方程f(x)＝0的近似解。其中X₀＝b，X_n＝X_n-1/f'(X_n-1)，(n＝1,2,3…)。

Matlab为解决非线性问题提供了良好的开发环境。在Matlab中，迭代法求解非线性方程通常就是求解多元函数零点的问题。最好要先知道零点的大致位置，因为一旦零点大致位置知道，就可用数值方法搜索到精确零点。有了初始零点后，求零点的精确值，可以借助命令fsolve进行。

于是，大圆一度的弧长乘以地心夹角θ即可得到大圆距离D：

D＝(2πR_z/360o)·θ＝111·θ  (28)式中R_z为地球半径。在方程(25)、(26)、(27)以及式(28)的计算中，北纬和东经取正值，南纬和西经取负值。

线性插值的常规算法：

图1所示为一平面上任意两点(x₁,y₁)、(x₂,y₂)，可代表二维图像上的任意两个像素，并用直线段将这两点连起来，(x,y)是直线段上的任意一点，字母S表示该点具有的数值（可以是灰度值或其他数据），本发明使用的是风速，于是S(x₁,y₁)表示点(x₁,y₁)处的风速大小，S(x₂,y₂)表示点(x₂,y₂)处的风速大小，S(x,y)表示点(x,y)处的风速大小。线性插值的观点是把两点之间线段上所有点所的风速值看作是线性均匀变化的。因此，当已知S(x₁,y₁)和S(x₂,y₂)时，则S(x,y)的大小为：

$S(x,y)=\frac{\sqrt{{(x-x_1)}^2+{(y-y_1)}^2}}{\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2}}(S(x_2,y_2)-S(x_1,y_1))+S(x_1,y_1)---\left(29\right)$

估算点风速求解：

图2所示为一台风眼壁的简化示意图，涂红色的部分是已经被分割出来的眼壁，涂绿色的部分是分割后的眼壁的中心位置，由于实际台风的眼壁并不呈正圆形，于是简化示意图也不画成圆形。台风年鉴上不仅提供了最大风速，还提供了中心的经纬度，中心经纬度可通过图像的坐标与经纬度的对应关系转换为坐标。现在以台风中心处的像素点为原点O'，水平向右为X'轴，竖直向下为Y'轴建立一个新的坐标系。当然首先要记录新坐标系与原坐标系之间的转换关系：

$\begin{array}{c}{x}^{\text{'}}=f_1\left(x\right)\\ y^{\text{'}}=f_2\left(y\right)\end{array}---\left(30\right)$ 式中，x,y分别为原坐标系中像素点的坐标，x',y'分别为新坐标系中的坐标，f₁,f₂为对应关系。将已经被分割好的图像的所有像素从原图像坐标系映射到新坐标系，这样，台风眼壁就被新的坐标系划分为四个区，对于图1所示简化示意图，规定右下角的区域为第一象限，右上角的区域为第二象限，左上角的区域为第三象限，左下角的区域为第四象限。眼壁处各点的风速可通过步骤6）中叙述的灰度风速模型进行预测得到。对于图像中任意给定一点的坐标，要求出对应的风速，则要先找出所需要的参考点，也即公式(29)中的两点(x₁,y₁)和(x₂,y₂)。其中一点即来自眼壁，我们规定，以图2建立的新坐标系为准，从Y'轴的正方向开始，沿着逆时针方向就可求出眼壁处各个点相对于Y'轴的夹角，对带估计风速的像素点也采用同样的规定得到与Y'轴的夹角，然后与眼壁各像素点的夹角进行比较，最接近的点便是参考点。

另外一个参考点可根据情况来确定选取台风中心还是台风外的一点。根据常识我们知道，在台风来临时虽然伴随着大风大雨，但是在台风中心处却是无风或微风，于是我们可认为中心处的风速为零，另外，一般的台风大小在半径1000公里，以外部分已不受台风影响，于是，台风外的一点也可作为风速为零的参考点，本发明的实验中取半径1000km处的风速为零。

计算前先要判断待估算的点是在眼壁内还是在眼壁外，只要计算待估算点到中心的距离是大于还是小于参考点到中心的距离即可，由于计算的结果跟实际距离有关，因此这里的坐标距离都应转化为实际距离，可通过坐标与经纬度的对应关系，计算实际距离。若待估算点在眼壁内，则在新坐标系中风速计算公式(29)改造为：

式中，(x₂,y₂)为参考点坐标，(x,y)为待估算点坐标，d₁为(x,y)到中心的实际距离，d为(x₂,y₂)到中心的实际距离，单位为km，S(x,y)为待估算点的风速，S(x₂,y₂)为参考点的风速。若待估算点在眼壁外，则取1000km处的风速为零，(29)式改造为：

$S(x,y)=\frac{1000-d_1}{1000-d}S(x_2,y_2)---\left(32\right)$

式中各个量的含义同式(29)。

由于采用了本发明所述的技术方案，得到的台风二维风场信息，平均误差和运行时间均比线性回归法和径向基神经网络小，对于台风强度的预测具有重要的参考价值。

附图说明

图1为现有技术中线性插值常规算法；

图2为眼壁的简化示意图。

图3为本发明基于SVM和PDE的有眼台风二维表面风场反演方法的算法流程图。

图4线性回归法拟合的误差散点图。

图5径向基神经网络法拟合的误差散点图。

图6SVM拟合的误差散点图。

图7线性回归法拟合的误差柱状图。

图8径向基神经网络法拟合的误差柱状图。

图9SVM拟合的误差柱状图。

图102007年9月17日20时台风云图。

图11初始化2007年9月17日20时的台风云图。

图122007年9月17日20时台风云图分割眼壁后的效果图。

图13线性回归法求得的风速的误差散点图。

图14径向基神经网络法求得的风速的误差散点图。

图15SVM求得的风速的误差散点图。

图16线性回归法求得的风速的误差柱状图。

图17径向基神经网络法求得的风速的误差柱状图。

图18SVM求得的风速的误差柱状图。

具体实施

如图1-图18为本发明一种基于SVM和PDE的有眼台风二维表面风场反演方法。下面结合说明书附图来进一步说明本发明，图3所示为本发明的流程图。

具体步骤如下：

1）在气象局提供的台风数据库中选取2005年到2008年的云图及其近中心最大风速的数据。

2）对步骤1中选中的图像利用PDE进行台风眼壁的分割，得到眼壁灰度值，一般来说，由于分割算法的误差和灰度图像本身拍摄的误差，使得眼壁处各像素点的灰度会有波动，并且一般中国气象局上海台风研究所台风年鉴中每一时刻只提供一个最大风速值，于是了为了减少随机误差，本发明将不使用单个像素点作为变量，而是对分割出来的整个眼壁的所有像素灰度取平均值，作为其中一个变量即：

$\stackrel{\‾}{Y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{Y}_i---\left(33\right)$

是眼壁处像素的平均灰度，Y_i是眼壁处第i个像素的灰度值，n是眼壁处的像素个数。得到平均灰度值后，将其与对应的近中心最大风速列于表1所示，

其中c为式（6）的可调参数：

表1云图分割后的台风眼壁灰度与最大风速数据

3）利用SVM建立灰度风速模型，并对建模效果进行拟合，同时将本发明的方法和线性回归法、径向基神经网络法进行对比。将误差的散点图和柱状图列于图4-图9，将误差值和运行时间列于表2：

表2各拟合算法的平均误差和运行时间

4）选定2007年9月17日20时（世界时）的韦帕云图作为本发明的测试图像，计算该时刻台风各处的风速值，此时的台风中心为东经125.7°，北纬23.4°。

5）对于给定图像进行分割眼壁得到参考点灰度，如表3所示，眼壁分割效果如图10-图12所示：

表3眼壁各点坐标值及其灰度值

6）将参考点灰度作为模型的输入值，求出参考点的风速如表4所示：

表4眼壁各点灰度值、风速

7）利用线性插值求出给定图像的30个估算点风速如表5所示，同时将利用线性回归法和径向基神经网络法作为对比，求的其各自的估算点风速分别如表6和表7所示。由本发明求的风速平均误差的散点图和柱状图，线性回归法和径向基神经网络法求的风速平均误差的散点图和柱状图如图13-图18所示。

表5预测风速与实际风速比较（SVM）

表6预测风速与实际风速比较（径向基神经网络）

表7预测风速与实际风速比较（回归法）

从表5-表7和图13-图18中可以看出，本发明的计算结果与实际测量结果还是比较接近的，除了有个别点的误差较大，如最大的误差达到了6.73m/s，大部分点的误差都在4m/s以内，平均误差为2.9m/s，而径向基神经网络比SVM平均误差稍微大一点，达到3.08，但是回归法平均误差为4.61，大大高于SVM和径向基神经网络，这证明了本发明提出的算法在台风风速反演上的有效性。

Claims (5)

1.基于SVM和PDE的有眼台风二维表面风场反演方法，其特征在于：首先对气象局提供的红外云图结合PDE算法进行图像分割，提取台风眼壁，得到眼壁各点的坐标值和灰度之后，再利用SVM进行平均灰度和台风近中心最大风速的建模，最后对于任意给定的红外云图进行分割眼壁，提取参考点灰度，结合灰度风速模型，求出参考点风速，利用线性插值反演整个二维表面风场，

具体步骤如下：

步骤1：在气象局提供的台风数据库中选取所有可用的红外云图及其近中心最大风速数据；

步骤2：对步骤2中选中的图像分割其台风眼壁；

步骤3：建立灰度风速模型；

步骤4：给定一图像及其中心经纬度和估算点经纬度；

步骤5：对于给定图像进行分割眼壁得到参考点灰度；

步骤6：将参考点灰度作为步骤3建立的模型的输入值，求出参考点的风速；

步骤7：利用线性插值法求出图像估算点风速。

2.如权利要求1所述的基于SVM和PDE的有眼台风二维表面风场反演方法，其特征在于，所述步骤2中的分割方法采用PDE方法为：根据台风眼的轮廓特点，采用活动轮廓模型，即将图像分割问题转化为最小化一个封闭曲线C(p)的“能量”泛函：

$E\left[C\left(p\right)\right]=\mathrm{\α}\underset{0}{\overset{1}{\∫}}\left|C_p\left(p\right)\right|\mathrm{dp}+\mathrm{\β}\underset{0}{\overset{1}{\∫}}{\left|C_p\left(p\right)\right|}^2\mathrm{dp}-\mathrm{\λ}\underset{0}{\overset{1}{\∫}}|\▿I[C\left(p\right)\left]\right|\mathrm{dp}---\left(1\right)$

其中第一项的积分是曲线的欧几里得弧长，第二项

$\underset{0}{\overset{1}{\∫}}{\left|C_p\left(p\right)\right|}^2\mathrm{dp}=\underset{0}{\overset{1}{\∫}}{\left|\mathrm{\κ}\right|}^2{\left(\frac{\mathrm{dp}}{\mathrm{ds}}\right)}^2\mathrm{dp}---\left(2\right)$

表示曲线“振荡”的能量，那么，最小化式(1)就需要要最小化前两项和最大化第三项，最小化前两项就得使得闭合曲尽量短、尽量光滑，因为在曲线逐渐变短的过程中也自然地趋向于光滑，所以式(1)第二项忽略不计，最大化第三项即要求曲线C尽量地与图像梯度

取到极大值的位置一致，

鉴于式(1)第三项的负号给计算带来的麻烦，现引进辅助函数g(r),r∈R⁺,则式(1)变为:

$E\left[C\left(p\right)\right]={\mathrm{\α}}_1\underset{0}{\overset{1}{\∫}}\left|C_p\left(p\right)\right|\mathrm{dp}+a_2\underset{2}{\overset{1}{\∫}}g(\▿I[C\left(p\right)\left]\right)\mathrm{dp}---\left(3\right)$

观察式(3)，发现E[C(p)]主要与曲线的形状、位置和参数p的变化有关，参数p的随意性大大影响了E的结果，为了避免这种情况，将活动轮廓模型改为使用无自由参数p的测地活动轮廓模型，为使E[C(p)]最小，需要曲线通过

的局部最小值，这个思想类似于光学中费马定理，于是能量泛函式为：

$L_R\left(C\right)=\underset{0}{\overset{L\left(C\right)}{\∫}}g\left(\right|\▿I\left[C\left(s\right)\right]\left|\right)\mathrm{ds}---\left(4\right)$

其中L(C)是闭合曲线C的欧几里得弧长，L_R(C)是经过加权后的弧长，要最小化式(4)，相应的梯度下降流为：

$\frac{\∂C}{\∂t}=g\left(C\right)\mathrm{\κ}\stackrel{\→}{N}-(\▿g\·\stackrel{\→}{N})\stackrel{\→}{N}---\left(5\right)$

其中前一项中

表示的是平均曲率运动，为曲面的单位法向量，但是式(5)存在一个明显的缺陷，即它只适合分割凸的对象，而不适合用来分割凹形对象，于是，将式(5)改进如下：

$\frac{\∂C}{\∂t}=g(c+\mathrm{\κ})\stackrel{\→}{N}-(\▿g\·\stackrel{\→}{N})\stackrel{\→}{N}---\left(6\right)$

其中参数c为一可调参数，通常与演化的速度有关，

在使用数值方案之前，首先要计算出图像的灰度梯度，才能得到边缘函数

梯度算子的作用类似于微分，因此会使图像尖锐化，计算梯度的过程对噪声非常敏感，为防止噪声的干扰，可先对图像进行预处理，高斯平滑是最常用的预处理方法，公式如下：

I_σ(x,y)＝I(x,y)*g_σ(x,y)  (7)

其中σ表示高斯方差，

接下来便要选取一种合适的边缘函数g了，一般可令边缘函数为

$g\left(r\right)=\frac{1}{1+{(r/k)}^p},p=\mathrm{1,2}---\left(8\right)$

式中k是待选的常数，边缘函数g的下降速度由它决定，有了边缘函数的定义后，将g(r)中的r用图像中的各个像素的梯度大小带入后，即可得到：

$g(x,y)=g\left(\right|\▿I\left|\right)---\left(9\right)$

这样就可以开始使用数值方案来实现图像分割操作了，我们采用水平集的数值方法，

曲线演化的一般方程式为：

$\frac{\∂C}{\∂t}=\mathrm{\β}\stackrel{\→}{N}---\left(10\right)$

曲线演化水平集方法的基本方程式为：

$\frac{\∂u}{\∂t}=\mathrm{\β}|\▿u|---\left(11\right)$

对比式（6）和式（10）可知，代入式（11）

$\frac{\∂u}{\∂t}=\mathrm{\β}|\▿u|=[g(c+\mathrm{\κ})+\▿g\·\frac{\▿u}{|\▿u|})\left]\right|\▿u|=\mathrm{cg}|\▿u|+\mathrm{div}\left(g\frac{\▿u}{|\▿u|}\right)|\▿u|---\left(12\right)$

由式(5)可知：

$\frac{\∂u}{\∂t}=\mathrm{div}\left(g\frac{\▿u}{|\▿u|}\right)|\▿u|---\left(13\right)$

由于图像以离散的形式存储，而偏微分方程描述的是连续的形式，因此在对式(12)和(13)作数值计算前，采用“半点离散化”方案离散化这两式中的散度算子，得到：

$s_{\mathrm{ij}}=\mathrm{div}{\left(g\frac{\▿u}{|\▿u|}\right)}_{\mathrm{ij}}\≈C_{1,\mathrm{ij}}{g}_{1,\mathrm{ij}}(u_{i,j+1}-u_{\mathrm{ij}})-C_{2,\mathrm{ij}}{g}_{2,\mathrm{ij}}(u_{i,j}-u_{i,j-1})+C_{3,\mathrm{ij}}{g}_{3,\mathrm{ij}}(u_{i+1}-u_{\mathrm{ij}})-C_{4,\mathrm{ij}}{g}_{4,\mathrm{ij}}(u_{i,j}-u_{i-1,j})---\left(14\right)$ 最后采用一种双曲型方程的迎风差分方法^[8]，对离散化的式(12)的偏微分方程进行数值计算；

得到显示方案：

$u_{\mathrm{ij}}^{n+1}=u_{\mathrm{ij}}^n+\mathrm{\Δt}(\mathrm{cg}{\▿}_{\mathrm{ij}}^{(-)}+\max (s_{\mathrm{ij}},0){\▿}_{\mathrm{ij}}^{(-)}+\min (s_{\mathrm{ij}},0){\▿}_{\mathrm{ij}}^{+})---\left(15\right).$

3.如权利要求1所述的基于SVM和PDE的有眼台风二维表面风场反演方法，其特征在于，所述步骤3中的建模方法采用SVM方法，假设输入变量是X，输出变量是Y，并已知X与Y之间有某种未知的关系，可认为它们遵循某个联合概率F(X,Y)，若已给定了l个独立的样本

并且集中于函数组{f(X,W)}中，求一个能使其期望风险达到最小的最优函数解{f(X,W₀)}，这样即可估计出X和Y之间具有的关系，即：

minR(W)＝L∫(Y，f(X，W))dF(X，Y)               (16)

式中，W是函数广义参数,而L(Y,f(X,W))是用f(X,W)对Y进行预测所形成的误差，Y为输出变量（假设函数为单值连续函数），误差估计函数可定义为：

L(Y,f(X,W))＝(Y-f(X,W))²  (17)

即采用平方和最小准则，

在一般的应用中，样本数目不可能是无限的，于是人们提出了经验风险最小化准则，也就是由样本来定义经验风险，表达式如下：

$R_{\mathrm{emp}}\left(W\right)=\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}L(Y,f(X_i,W))---\left(18\right)$

通过VC理论证明了实际风险R(W)和经验风险R_emp(W)以概率1-η满足下式：

$R\left(W\right)\≤R_{\mathrm{emp}}\left(W\right)+\sqrt{\frac{h(\ln (2n/h)+1)-\ln (\mathrm{\η}/4)}{n}}---\left(19\right)$

上式又可简化为：

(19)式中，h是函数集的VC维数；l是已知的有限样本数；η是满足0＜η＜1的参数；

称为置信区间，于是式(20)称为结构风险式，由式(20)可知，为了使结构风险获得最小值，则要求式(20)中的右边两项之和取得最小值，

此时在支持向量机的使用过程中加入核函数，核函数的作用是将低维非线性空间中的复杂运算映射到某个高维线性空间，在高维空间中处理可相对简单，计算完毕再返回到低维空间，在高维空间中一般是计算内积和，从而可避免维数灾难对计算过程的影响，高维线性空间中构造支持向量机的思想与一般的线性回归类似，常见的核函数有以下几种：

(1)线性核，即两个向量的内积：

K(X,X_i)＝＜X,X_i＞

(2)多项式核：K(X,X_i)＝(＜X?X_i＞+1)^d，其中d为多项式核的阶数，

(3)径向基核（RBF），传统的径向基核使用以下判定准则，即

$f\left(X\right)=\mathrm{sgn}(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\α}}_i{K}_{\mathrm{\γ}}\left(\right|X-X_i\left|\right)-b))$

其中，K_γ(|X-X_i|)是一个非负的单调函数，其大小取决于两向量的距离|X-X_i|，若要使它趋于0，

则训练的样本总数要趋于无穷大，这个判定准则需要顾及众多参数，包括：γ的值、X_i的数目n、以及α_i的值，最常见判定准则是高斯函数，即：

其中σ是核函数的宽度，高斯函数与经典径向基函数的主要区别在于中心点X_i和输出权值均由支持向量机算法本身决定，而每个中心点对应着下一个支持向量，

(4)多层感知机核，满足Mercer条件的Sigmoid函数为K(X,X_i)＝tanh(γ·＜X,X_i＞+c)，其中γ是梯度参量；c为偏离参量，采用Sigmoid函数作为内积后，则使得此时的SVM包含了一个单多层感知机，其中的隐含层节点数可以由算法自动确定，

支持向量机算法的实现：

根据已知样本空间，先选择合适的核函数

将低维非线性空间中的X_i映射到高维线性空间的

接着在高维空间的内积计算过程中，使用核函数K(X_i,X_j)来替换内积运算

从而实现非线性回归，于是问题求解的形式可转化为一个最优二次规划问题，即目标函数为：

$\min \frac{1}{2}{\left|\right|w\left|\right|}^2---\left(21\right)$

约束条件为：

引入Lagrange算子后，即可得到该问题的对偶形式：

$\min w\left(a\right)=\frac{1}{2}\underset{i,j=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)(a_j^{*}-a_j)K(X_i,X_j)-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{Y}_i(a_i^{*}-a_i)+\mathrm{\ϵ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}+a)---\left(23\right)$

约束条件为：

$\begin{array}{c}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)=0\\ a_i^{*},a_i\≥0,i=\mathrm{1,2},...,l\end{array}---\left(24\right)$

a_i≠0对应的向量即为支持向量，要得到这两项，只需求解由式(23)、(24)构成的Q-P问题，

最终我们得到了回归机的解答式：

$Y=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+b$

基于支持向量机的回归算法的一般步骤如下：

步骤1：给出已知的训练样本集

其中

步骤2：根据问题的需求选择合适的核函数K(X,X_i)，并设定精度范围ε；

步骤3：求解受式（24）约束的式（23）的最优二次规化问题，求得

步骤4：在

中选择任意一个不为0的分量，即支持向量，并计算

$b=Y_j-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+\mathrm{\ϵ},$ 当 $a_i\≠0,b=Y_j-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)-\mathrm{\ϵ},$ 当

步骤5：构造非线性超平面 $Y=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}(a_i^{*}-a_i)K(X_i,X_j)+b$

拟合方法的选择：

数据获得后便可进行数据拟合，使用循环的模式，总共有多少个样本就循环多少次，每次循环只取一个样本用来作预测样本，其他的样本用作训练样本，并计算相应的误差，一次循环计算一次误差大小，每次循环所用的预测样本都要求不同，于是所有循环结束后，每个样本的预测误差都可以描绘出。

4.如权利要求1所述的基于SVM和PDE的有眼台风二维表面风场反演方法，其特征在于，所述步骤7中的云图上任意两点的经纬度计算实际距离方法为：设A、B两点的纬度分别为A和B，经度分别为la和lb，要求的是A、B两点间的大圆弧长，令A、B的中点的纬度为C，经度为lc，A、B两点间的地心夹角为θ，则根据球面三角关系，可列出如下方程：

cos(θ)＝sin(A)sin(B)+cos(A)cos(B)cos(la-lb)  (25)

cos(θ/2)＝sin(A)sin(C)+cos(A)cos(C)cos(la-lc)  (26)

cos(θ/2)＝sin(B)sin(C)+cos(B)cos(C)cos(lb-lc)  (27)由于C、lc、θ为未知量，所以式(26)，(27)为非线性方程，应用迭代法解出近似的数值解，

设f(x)＝0为非线性方程，f为[a,b]上二阶可导函数，且f'(x)、f″(x)均不为零，不妨设f'(x)＞0，f″(x)＞0即f为[a,b]上递增凸函数，又设f(a)＜0，f(b)＞0，那么由连续函数的介值定理知，方程f(x)＝0在[a,b]内存在唯一解ζ，牛顿迭代法的基本思想是构造一点列{X_n}，使

且当n充分大时，可以将X_n作为ζ的近似值，即方程f(x)＝0的近似解，其中X₀＝b，X_n＝X_n-1/f'(X_n-1)，(n＝1,2,3…)，

于是，大圆一度的弧长乘以地心夹角θ即可得到大圆距离D：

D＝(2πR_z/360^o)·θ＝111·θ  (28)式中R_z为地球半径。

5.如权利要求4所述的基于SVM和PDE的有眼台风二维表面风场反演方法，其特征在于，在方程(25)、(26)、(27)以及式(28)的计算中，北纬和东经取正值，南纬和西经取负值。

Images