台风预测与模拟

1、对台风的模拟与预测的研究摘要台风是我国东南沿海每年遭受的严重自然灾害之一，为了减少人民生命财产损失，准确有效的台风预报显得尤为必要。本文围绕台风相关数据的预测、台风的路径图以及预测福州10100年内可能遭遇的台风的最大风力问题进行了讨论。对台风的预测建立了神经网络预测模型，并对求解结果进行了详细的分析。针对问题一，首先把台风B中心的经纬度以半小时为单位的形式算出来，然后分别求出台风B中心和观察点的距离；然后利用拟合工具箱，分别得出距离与气压的关系式：和距离与风速的关系式：；最后算出台风A中心到福州的距离大致为760公里，把该距离代入所得关系式，即可得出福州此时的气压为（1.0016e+03）百

2、帕，风速为4.5195m/s。针对问题二，为了预测该台风此后72小时内的中心气压和最大风速，需要该台风以往的数据，因此，我们首先收集了该台风相关的数据，然后运用广义回归神经网络模型，对此后72小时内的中心气压和最大风速进行预测，得出结果，并作出该台风的路径图。并且，为了更直接地看出预测数据与真实数据的误差，我们分别作出时间与气压、时间与风速的关系图，从图中可看出，预测值与真实值出入不大，且路径走向一致。针对问题三，我们建立了广义回归神经网络模型，对福州10-100年内可能遭遇的最大风力进行预测。首先把收集到的数据作为训练样本，利用广义回归神经网络模型求解，最终预测出福州10-100年内可能遭遇

3、的最大风速为：44.6721m/s，再根据风力等级划分表确定出最大风力为：14级。关键词一、问题重述台风是热带气旋的一个类别;按世界气象组织定义:热带气旋中心持续风速达到 12 级(即每秒 32.7 米或以上)称为飓风(hurricane),飓风的名称使用在北大西洋及东太平洋;而北太平洋西部(赤道以北,国际日期线以西,东经 100 度以东)使用的是台风(typhoon)。台风是我国东南沿海每年遭受的严重自然灾害之一，台风水平结构分为台风眼区域、最大风雨区以及外围区;垂直结构为下层流入层、中层过渡层及上层流出层。通常在最大风雨层发生强风、强降雨;在下层区域吸收能量,形成低气压。台风会带来巨大降水

4、,同时也会带来巨大灾难,其形成的自然灾害种类包括风灾,潮灾以及水灾,其中以潮灾造成的损失最为巨大。据统计,风暴潮造成损失居全世界之首1。为了减少人民生命财产损失，准确有效的台风预报显得尤为必要。请收集相关数据，建立数学模型，完成下面问题。问题1 请结合附件1和2，根据气象学和空气动力学原理，建立数学模型，给出此时福州台风相关数据预测。问题2 请收集相关数据，根据气象学和空气动力学原理，建立数学模型，给出此后72小时内的该台风相关预报数据并画出路径图问题3 台风对沿海建筑的破坏尤为明显，和抗震等级设计一样，为了设计高层建筑的抗风能力，需要估算建筑物设计年限内可能遭遇的最大台风风力。请收集相关数据

5、，根据气象学和空气动力学原理，建立数学模型，给出福州10-100年内可能遭遇的最大风力。二、问题分析我们共需要解决三个问题：（1）给出此时福州台风相关数据预测；（2）给出此后72小时内的该台风相关预报数据并画出路径图（3）给出福州10-100年内可能遭遇的最大风力对于问题一，根据附件1和附件2，可知它们的联系是经纬度，即位置。但是，附件1和附件2的时间分隔并不一致。因此，首先把台风B中心的经纬度以半小时为单位的形式算出来，然后分别求出台风B中心和观察点的距离；然后利用拟合工具箱，分别得出距离与气压、距离与风速的关系式；最后算出台风A中心到福州的距离，把该距离代入所得关系式，即可得出福州此时的气

6、压和风速。对于问题二，欲预测福州11时后72小时内的大气压和风速，并且画出路径图，这需要该台风以往的数据。我们发现该台风在11时的参数与台风“龙王”11时的路径参数一致，因此，我们可以收集“龙王”的数据，以此作为基准，利用广义回归神经网络方法，对福州11时后72小时内的大气压和风速进行预测。对于问题三，主要在于预测福州10-100年内可能遭遇的最大风力，首先收集福州近年遭遇台风的具体情况，把所得数据作为训练样本，再利用广义神经网络得出预测样本，即福州10-100年内可能遭遇的最大风力。三、模型的假设（1）附件所给数据和我们查找的数据准确无误；（2）假设台风的行走路径没有障碍；（3）排除其他突发

7、性情况，例如气流漩涡；四、符号说明符号意义大气压风速台风A中心与福州的距离五、模型的建立与求解5.1问题一5.1.1数据处理因为附件1和2的时间分隔不一致，附件1是以一小时为单位，而附件2是以半小时为单位，因此，首先用分段线性插值把经纬度以半小时为单位的形式算出来。结果见附录1。然后，根据台风B中心与观察点B的经纬度，分别求出相对应的两地之间的距离。结果见附录2（包括距离、气压和风速）。5.1.2模型求解利用拟合工具箱分别得到气压与距离、风速与距离的关系式。详细结果如下：（1）气压与距离关系的结果General model Power2: f(x) = a*xb+cCoefficients (

8、with 95% confidence bounds): a = -203.5 (-543.4, 136.4) b = -0.615 (-1.226, -0.004049) c = 1005 (991.7, 1018)Goodness of fit: SSE: 99.83 R-square: 0.8642 Adjusted R-square: 0.8552 RMSE: 1.824图1 距离与气压关系拟合图（2）距离与风速关系式的结果： General model Power2: f(x) = a*xb+c Coefficients (with 95% confidence bounds): a

9、 = 186.5 (-83.6, 456.6) b = -0.5076 (-1.107, 0.09152) c = -1.913 (-24.99, 21.17) Goodness of fit: SSE: 142 R-square: 0.8672 Adjusted R-square: 0.8584 RMSE: 2.176图2 距离与风速关系拟合图5.1.2模型检验（1）检验理论线性回归预测是通过一组统计观测数据确定最优拟合线性关系，但我们需要对这种关系拟合的效果好坏进行评判，这种评判通常称为模型检验2。评判的结论将直接影响人们对线性回归模型的信任程度，从而也影响对预测结果的信任程度。评判的标准

10、主要是拟合的误差，如果拟合误差比较小，拟合效果就好，预测结果的信任程度就高。如果拟合误差较大，拟合效果就不太好，严重时还必须重新考察历史数据、选择变量，再重新拟合。为了评判误差产生程度，我们介绍两个基本定量。1R-square(确定系数)如果因变量的一组统计观测数据的平均值为，所有统计观测数据值都分布在这个均值的上下，我们可以求出其总的误差平方和。计算公式为： (1)如果通过线性回归模型拟合的值为，那么可以求出回归的误差平方和。计算公式为： (2)可以认为回归的误差平方和是的一部分，也就是说回归模型部分解释了实际观测值对均值的偏离，而剩余部分为，即 (3)因此可以解释同时也可以写成 (4)显然

11、，回归模型拟合较好，则总的误差平方和越能够用回归的误差平方和来表示，模型所描述的线性关系就越准确。所以，我们定义确定性系数为回归的误差平方和占总误差平方和的比例，即 (5)可以看出，的值在之间，如果的值接近1，说明实际数据对均值的绝大部分都可以由回归明显来解释，模型的拟合效果就越好；如果的值接近零，说明实际数据对均值的绝大部分都不能由回归明显来解释，即模型拟合得不好。2.(均方根)均方根误差亦称标准误差，其定义为,在有限测量次数中，均方根误差常用下式表示： (6)式中：为测量次数；为一组测量值与真值的偏差。可以看出越小越好。（2）检验结果1、在气压与距离的关系拟合结果中，R-square= 0

12、.8642，RMSE=1.824，所以该模型基本通过检验2、在距离与风速的关系拟合结果中，R-square= 0.8672，RMSE=2.176，所以该模型基本通过检验5.1.3模型求解根据台风中心A和福州的经纬度（119.28，26.08），得出两地之间的距离大约为760公里。根据上面拟合的结果，可得气压与距离的关系式为 (7)把距离代入得到福州此时的大气压为（1.0016e+03）百帕。同理，风速与距离的关系式为 (8)把距离代入得到福州此时的风速为4.5195m/s。根据风力等级表3，可知，此时正为3级风。5.2问题二5.2.1模型理论1、广义回归神经网络5（GRNN） 广义回归神经网络

13、是美国学者Don-ald F.Spencht在1991年提出的，它是径向基神经网络的一种。GRNN具有很强的非线性映射能力和柔性网络结构以及高度的容错性和鲁棒性，适用于解决非线性问题。GRNN在逼近能力和学习速度上较RBF网络有更强的优势，网络最后收敛于样本量积聚较多的优化回归面，并且在样本数据较少时，预测效果也较好。此外，网络还可以处理不稳定的数据。因此，GRNN在信号过程、结构分析、教育产业、能源、食品科学、控制决策系统、药物设计、金融领域、生物工程等各个领域得到了广泛的应用。2、GRNN的网络结构在结构上与网络结构相似。它是由四层构成，如图所示，分别为输入层、模式层、求和层和输出层。对应

14、网络输入,其输出为。图1 广义回归网络结构图（1）输入层输入层神经元的数目等于学习样本中输入向量的维数，各神经元是简单的分布单元，直接将输入变量传递给模式层。（2）模式层模式层神经元数目等于学习样本的数目，各神经元对应不同的样本，模式层神经元传递函数为 (9)神经元的输出为输入变量与其对应的样本之间的距离平方的指数形式。式中，为网络输入变量；为第个神经元对应的学习样本。（3）求和层求和层中使用两种类型神经元进行求和。一类的计算公式为,它对所有模式层神经元的输出进行算术求和，其模式层与各神经元的连接权值为1，传递函数为 (10)另一类计算公式为，它对所有模式层神经元进行加权求和，模式层中第个神经

15、元与求和层中第个分子求和神经元之间的连接权值为第个输出样本中的第个元素，传递函数为 (11)（4）输出层输出层中的神经元数目等于学习样本中输出向量的维数，各神经元将求和层的输出相除，神经元的输出对应估计结果的第个元素，即 (12)3. 理论基础广义回归神经网络的理论基础是非线性回归分析，非独立变量相对于独立变量的回归分析实际上是计算具有与最大概率值为，设随机变量和随机变量的联合概率密度函数为，已知的观测值为，则相当于的回归，也即条件均值为： (13)即为在输入为条件下，的预测输出。应用Parzen非参数估计，可由样本数据集，估算密度函数。式中，为随机变量和的样本观测值；为样本容量；为随机变量的

16、维数；为高斯函数的宽度系数，在此称为光滑因子。用代替代入式（第一个），并交换积分与加和的 顺序： 由于，对两个积分进行计算后可得网络的输出为： 估计值为所有样本观测值的加权平均，每个观测值的权重因子为相应的样本与之间Euclid距离平方的指数。当光滑因子非常大的时候，近似于所有样本因变量的均值。相反，当光滑因子趋向于0的时候，和训练样本非常接近，当需预测的点被包含在训练样本集中，公式求出的因变量的预测值会和样本中对应的因变量非常接近，而一旦碰到样本中未能包含进去的点，有可能预测效果会非常差，这种现象说明网络的泛化能力差。当取值适中，求预测值时，所有训练样本的因变量都被考虑了进去，与预测点距离近

17、的样本点对应的因变量被加了更大的权。5.2.2模型的求解我们发现台风C和台风“龙王”11时的各个参数近似，利用两个台风之间的相似性，于是，我们收集“龙王”的路径数据（见附录3），把11时以前的数据作为训练样本。训练结果如下表和下图所示：图3 福州11时后的72小时路径预测图根据上图，绿色线段反映“龙王”真实路径，红色线段是预测出的“龙王”11点以前的路径，两条路径出入很小，可以看做预测路径与真实路径吻合。表1 11时以后72小时内的预测气压与风速表时间纬度经度气压风速1132.0122.12933.3851.104131.7722.10934.2550.517131.9422.08934.53

18、49.5611132.0522.05934.7750.1514131.3622.07934.6251.0017131.1822.05934.9050.0720131.0722.05935.2049.1323130.7522.08934.9550.0326130.6422.08934.9450.0429130.5022.11935.1349.9232130.3222.14934.9950.0135130.0722.16935.0049.2638129.9522.16935.0049.2641129.6922.16935.0050.0044129.8222.22933.4151.5947129.4

19、322.23933.0752.8150129.0922.24932.4252.2753128.7222.28929.5654.4356128.4722.31927.5354.7459128.4122.31925.2856.8662127.9122.37923.5857.9165127.7622.40923.4657.9768127.7522.40923.2358.1471127.2422.45921.0659.4074127.3722.39920.7058.6477126.8222.51920.4559.7580126.6322.55920.6658.6883126.3822.58920.31

20、58.8586126.1922.66920.2059.8989125.7822.71920.1260.9392125.5022.75920.0859.9395125.3422.75920.1059.9698125.0022.86920.0459.98101124.8722.96920.0460.98104124.4122.97920.0259.99107124.4622.97921.2960.00110124.0923.05920.9360.00113123.7923.12923.5759.02116123.9423.14922.5860.00119123.4223.19923.8459.14

21、121123.1423.26925.8058.41124123.1723.26927.1357.46127122.6323.41927.8457.18130122.3523.50929.7856.26133122.4723.43932.6254.06136121.9623.62935.6951.66139121.4823.76935.2052.76142121.0623.72939.0349.73145121.1123.72940.3050.75148120.6823.66942.5649.55151120.7423.78943.6548.35154120.3523.77945.2447.47

22、157120.1923.89946.0947.03160120.0423.90948.7945.64163119.8723.87953.5743.38166119.6023.99954.4442.96169119.5124.03956.3742.03172119.3524.11957.9140.77175119.2424.16959.6040.43178119.2424.23962.1939.12181119.0524.27964.3538.49184118.8124.29966.2337.05187118.5324.37970.6534.34190118.5924.34971.3434.35

23、193118.2424.43975.5231.87196118.1024.46979.3529.90199117.9324.52982.1028.10202117.8624.56982.9527.42205117.6624.64986.2225.25为了更直观地表现出预测值与真实数据的误差，做出如下两个关于时间和气压、时间和风速的图形：图 4 时间与气压图其中蓝色线段代表真实值，红色线段代表预测值。从上图可以看出，预测气压值与真实气压值相差不多，并且两条曲线的走向一致，因此，该神经网络模型预测效果较好。图 5 时间与风速图其中蓝色线段代表真实值，红色线段代表预测值。从上图可以看出，预测风速值与

24、真实风速值相差不多，并且两条曲线的走向一致，因此，该神经网络模型预测效果较好。5.3问题三为了预测福州10-100年内可能遭遇的最大风力。我们通过网上查找台风年鉴4关于福州近40年内遭遇台风的具体情况，统计该年遭遇台风的最大风力以及最大风速，得到数据统计如下表：表 2 福州近40年内遭遇台风情况时间最大风力(级)最大风速（m/s）1980123519821120198340198545199092419928201998321999352001123820051129.5200611302009102620101235201211332013134020141235根据以上数据作为训练样本，使

25、用广义回归神经网络模型预测得出福州10-100年内可能遭遇的最大风力。通过MATLAB编程（具体程序见附录）实现，下图为福州10-100年内可能遭遇的最大风力：图6 福州10-100年内可能遭遇的最大风力运行结果为： mg = 44.6721综上所述：福州10-100年内可能遭遇的最大风速为：44.6721m/s，再根据风速和风力等级的关系划分表：表 3 风力等级表 风力等级表-风力等级划分风力等级陆地地面物体征象海面波浪/浪高（米）相当公里/时风速（米/秒）0静，烟直上平静/0.0小于100.21烟能表示风向微波峰无飞沫/0.1150.31.5.1313414937.041.41415016

26、641.546.11516718346.250.9根据上表得出最大风力为：14级。6、 模型的优缺点参考文献1刘清容，于建生，韩笑.风暴潮研究综述及防灾减灾对策.应用科技：2009.226-2272 .2014.08.24.3 .2014.08.25.4 .2014.08.25.5司守奎，孙玺菁.数学建模算法与应用M.国防工业出版社.2013.02.6姜启源，谢金星，叶俊.数学模型M.高等教育出版社.2011.01.附录附录1 用分段线性插值所求以半小时为单位的经纬度纬度经度21.7117.421.7117.321.7117.221.75117.1521.8117.121.8511721.91

27、16.921.95116.822116.722.05116.622.1116.522.15116.3522.2116.222.2116.122.211622.2115.8522.2115.722.25115.622.3115.522.35115.322.4115.122.45114.9522.5114.822.5114.6522.5114.522.55114.4522.6114.422.65114.322.7114.222.7114.122.711422.7113.8522.7113.7附录2 台风中心B各个位置与观测点B之间的距离距离 气压 风速214.534 997.9 6.0205.47

28、5 996.8 9.0196.537 996.4 10.7189.252 995.8 10.5181.995 995.2 11.5170.261 994.2 12.4158.530 993.3 12.5146.802 993.0 13.9135.077 993.0 16.3123.358 992.7 15.5111.642 992.0 14.395.5103 991.4 16.379.5056 990.3 17.271.2135 989.2 19.463.5089 987.5 21.853.6049 986.5 22.846.8002 984.0 26.239.1780 980.7 31.83

29、3.9068 982.0 31.638.5958 986.2 22.552.2848 989.3 17.764.9067 991.5 17.178.9213 992.3 19.394.2122 993.9 16.5109.541 994.6 16.0114.229 995.2 16.4119.211 995.3 15.5129.583 996.2 13.1140.157 997.1 13.1150.400 997.6 15.3160.646 998.0 12.9176.022 998.8 12.3191.403 999.4 11.1附录3 台风“龙王”路径相关数据编号 经度 纬度 气压 风速1

30、 143.5 19.6 996 187 143.1 20 992 2013 142.6 20.3 990 2319 141.7 20.9 990 2325 141.1 21.5 970 3531 140.1 21.7 960 4037 139.2 22 955 4043 138.5 22.1 950 4549 137.4 22.3 950 4555 136.8 22.4 940 5061 136.1 22.4 935 5067 135.6 22.4 935 5073 135.2 22.4 935 5579 134.7 22.4 920 6085 133.8 22.3 920 6091 133.

31、1 22.3 920 6097 132.1 22.1 935 50100 131.6 22.1 935 50103 131 22 935 50106 130.6 22.1 935 50109 129.9 22.2 935 50112 129.3 22.2 935 50115 128.5 22.3 930 55118 127.7 22.4 920 60121 127 22.5 920 60122 126.9 22.5 920 60123 126.6 22.5 920 60124 126.3 22.6 920 60125 126 22.6 920 60126 125.7 22.7 920 6012

32、7 125.4 22.8 920 60128 125.2 22.8 920 60129 124.9 22.9 920 60130 124.6 23 920 60131 124.3 23 920 60132 124 23 920 60133 123.9 23.1 920 60134 123.6 23.1 920 60135 123.5 23.2 925 60136 123.2 23.3 925 60137 123 23.3 930 60138 122.7 23.3 925 60139 122.5 23.4 930 55140 122.3 23.5 930 55141 122 23.6 930 5

33、5142 121.6 23.7 930 55143 121.2 23.8 940 50144 121 23.9 940 50145 120.6 24 940 50146 120.5 24.1 940 50147 120.2 23.8 950 45148 120 23.9 950 45 149 119.9 23.9 950 45150 119.8 23.9 950 45 151 119.6 23.9 950 45152 119.4 24 950 45 153 119.2 24 960 40154 119.1 24.1 965 38155 119 24.2 965 38156 118.9 24.3

34、 965 38157 118.8 24.4 965 38158 118.6 24.5 970 35 159 118.4 24.5 970 35161 118.1 24.4 980 30162 117.9 24.5 980 30163 117.7 24.6 988 25164 117.5 24.6 990 23 165 117.3 24.6 995 20166 117.2 24.8 995 20167 117 24.9 995 20169 116.8 25.1 1000 15附录4 问题二时间与气压关系代码：x=1:69;plot(x,shuju(:,3),'r')hold on

35、plot(x,data(101:169,4)时间和风速关系代码：x=1:69;plot(x,shuju(:,4),'r')hold onplot(x,data(101:169,5)广义回归神经网络（GRNN）代码：clc;clear allclose allnntwarn off;% 载入数据data=xlsread('F:wenti2shuju');% 载入数据并将数据分成训练和预测两类l=100;shuju=;while l<169p_train=data(1:l,1);t_train=data(1:l,2:5);p_test=data(l+1,1);

36、t_test=data(l+1,2:5);% 交叉验证desired_spread=;mse_max=10e20;desired_input=;desired_output=;result_perfp=;indices = crossvalind('Kfold',length(p_train),4);h=waitbar(0,'正在寻找最优化参数.')k=1;for i = 1:4 perfp=; disp('以下为第',num2str(i),'次交叉验证结果') test = (indices = i); train = test

37、; p_cv_train=p_train(train,:); t_cv_train=t_train(train,:); p_cv_test=p_train(test,:); t_cv_test=t_train(test,:); p_cv_train=p_cv_train' t_cv_train=t_cv_train' p_cv_test= p_cv_test' t_cv_test= t_cv_test' p_cv_train,minp,maxp,t_cv_train,mint,maxt=premnmx(p_cv_train,t_cv_train); p_cv_t

38、est=tramnmx(p_cv_test,minp,maxp); for spread=0.1:0.1:2; net=newgrnn(p_cv_train,t_cv_train,spread); waitbar(k/80,h); disp('当前spread值为', num2str(spread); test_Out=sim(net,p_cv_test); test_Out=postmnmx(test_Out,mint,maxt); error=t_cv_test-test_Out; disp('当前网络的mse为',num2str(mse(error) pe

39、rfp=perfp mse(error); if mse(error)<mse_max mse_max=mse(error); desired_spread=spread; desired_input=p_cv_train; desired_output=t_cv_train; end k=k+1; end result_perfp(i,:)=perfp;end;close(h)disp('最佳spread值为',num2str(desired_spread)disp('此时最佳输入值为')desired_input;disp('此时最佳输出值为&

40、#39;)desired_output;% 采用最佳方法建立GRNN网络 net=newgrnn(desired_input,desired_output,desired_spread);p_test=p_test'p_test=tramnmx(p_test,minp,maxp);grnn_prediction_result1=sim(net,p_test);grnn_prediction_result=postmnmx(grnn_prediction_result1,mint,maxt)'shuju=shuju;grnn_prediction_result;l=l+1;endshujuplot(shuju(:,1),shuju(:,2),'r')hold onplot(data(:,2),data(:,3),'g')附录5 问题三广义回归神经网络（GRNN）代码：indices = crossvalind('Kfold',length(p_train),4);h=waitbar(0