应用吸收结构法预测火山灰云漂移轨迹2

1、应用吸收结构法预测火山灰云漂移轨迹马宝君1赵谊2梁跃1李太岩1张明3 MaBaoJun1,ZhaoYi2,LiangYue1,LiTaiYan1 and ZhangMing31. 牡丹江地震台，牡丹江 1570092. 黑龙江省地震局，哈尔滨 1500903. 上海大学，上海 2004331. Mudanjiang seismic station,Mudanjiang 157009,China2. The earthquake administration of HeiLongJiang,Harbin 150090,China3. Shanghai University,Shanghai 20

2、0433,ChinaAbstract Volcanic ash clouds can poses direct threats to aviation safety,so lts necessary to forecast the trajectories of volcanic ash clouds. Depending on the date of numerical weather prediction, attracting structures and PUFF model can both do it. By using the principle of the two

3、methods, this paper using the precise trajectory atmosphere parameter data of United States from 1958 to 1997 to calculate the attracting structures of wind field.with the computational process of simplified and get rid of the dependence of numerical weather prediction data, it realize the predictio

4、n of ash cloud drift path, And comparing with attracting structures.the results show that it is feasible to track the ash cloud with the calculation method of wind field.基金项目:国家自然科学基金资助项目（41172303）；作者简介：马宝君，男，1972年生，大学本科，高级工程师，主要从事火山研究工作，E-mail：mbj885Key words wind field;attracting structures;PUFF 摘

5、 要火山灰云对航空运输安全产生直接危害，所以对火山灰云轨迹进行预测是非常必要的。吸收结构法和PUFF模型都能预测火山灰云飘移轨迹，这两种方法都依赖数值天气预报的数据，本文应用这两种方法的原理，采用美国1958年至1997年精确大气轨迹参数数据，计算风场的吸收结构，减化计算过程，摆脱对数值天气预报数据的依赖，实现火山灰云漂移路径预测，并与吸收结构方法进行了比较.结果表明，应用风场计算的方法进行火山灰云追踪是可行的。关键词风场；吸收结构；PUFF1 引言2010年4月冰岛火山喷发是对火山灰对航空安全的威胁产生了更深的认识。在上世纪90年代，九个火山灰云咨询中心（VAACs）在全世界范围内被建立，在

6、过去的十年间，许多灰云的传播和扩散模型得到发展并且被VAACs和其他中心使用。这些模型包括MEDKA1、HYSPLKT2、PUFF3、CANERM4、NAME5，结合数值天气预报中心（NWP）的数据，这些模型能够预测火山灰云在大气中的传播，为当局对灾害回应计划提供科学的依据。尽管VATD模型是非常有用的，但是这些模型在广泛的应用上存在很多限制。最大的挑战就是预测精度，主要有两个因素。第一，VATD模型使用拉格朗日粒子轨迹，这种方法的精度主要依靠NWP数据的质量，尤其是它的时间和空间分辨率。同参考轨迹相比，在传输时间上灰云粒子轨迹的平均误差超过36小时的概率能达到35%6。使用更精细网格的NWP

7、数据能够改善预测精度，但是却极其消耗时间的，时间上的延迟对于快速精确的预测都是不利的。另外一个限制广泛应用在实时监测和响应的因素是缺乏初始喷发灰云的精确描述。初始参数对于许多VATD(volcankc ash transport and dksperskon)模型的成功是非常重要的。例如最大高度，粒子大小，以及浓度分布等，然而，收集这些数据通常来讲是非常难的，尤其是在喷发后最初的几个小时内急需一些快速响应的时候。另外，这些参数值在喷发的工程中也会不断的改变。因此，利用VATD模型在最初灰云的物理特性上和风场特性上均受到限制。为了克服这些限制，Jkfeng Peng 7提出了吸收结构法，该方法致

8、力于风场的全部特性。更特殊的是，这种方法使用动态系统方法在风场中定义吸收结构。这些结构扮演着全部的吸收因子的角色，粒子会朝着它们移动。这些结构独立于粒子参数。吸收结构方法采用数值天气预报（NWP）数据，并通过天气研究和预测模型（WRF）采用大型计算机，通过定制的模拟方法得出。由于WRF软件是气象专用软件，采用NWP数据，在Unix操作系统上运行，能计算风的大小和方向，计算方法复杂且NWP 数据不易拿到。在借鉴吸收结构法PUFF模型的基础上，本文应用美国精确的大气轨道参数数据，应用龙格-库塔法计算了全球及冰岛附近的吸收结构，计算结果与彭继枫的非常接近，探索适合我国国情的火山灰云追踪方法。2 数值

9、模拟的设计方案2.1吸收结构法的基本原理7 基于烟柱初始参数的模型方法存在很大的不精确性，并且在喷发后的数小时内烟柱的参数信息很难获得， Jkfeng Peng7提出了一种基于风场数据的吸收结构方法，这种方法致力于风场的全部特性。更特殊的是，这种方法使用动态系统方法在风场中定义吸收结构，这些结构扮演着全部的吸收因子的角色，粒子会朝着它们移动，并且这些结构独立于粒子参数。 Jkfeng Peng7利用吸收结构法，选用NWP数据，用WRF软件，推算了2010年4月冰岛火山喷发后的速度场，见图2-1、图2-2、图2-3，效果很好。图2-1 风场风速大小 图2-2 风的方向图2-3 冰岛火山灰的漂移过

10、程(高程,9km）摘自Jkfeng Peng72.2 风场计算的替代方案由于对数值天气预报数据的处理非常复杂，需要使用具有大量存储空间的大型计算机进行大量的计算，并且处理过程需要一定长的时间（处理十天的数据大概需要4个多小时），所以我们在借鉴吸收结构方法的基础上，提出了大气轨道参数的计算方法，简化了处理过程。由于采用计算长时间风速平均值的方法，这种方法更具代表性。 笔者从美国国家气象中心得到了1958-1997年全球大气精确的轨道参数，该数据垂直分辩率为将大气从地球表面到同温层分成17层，水平分辩率为2.5°×2.5°1个点位数据，每月每层每点都提供1个SN向、E

11、W向风速的均值，总数据量达到了2.4G，该中心提供了1个wgrkb.exe检索软件，每执行1次只能检索1年中1个月1层中的风速，笔者应用这40年的连续观测数据得到了全球每个观测点的风资料数据，然后用龙格库塔算法对已有观测数据进行内插，最终得到全球任意点的风资料数据。2.3 龙格-库塔算法龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的，设一阶微分方程的初值问题y'=dydx=f(x,y)yx0=y0的解为y=yx。如果是等距节点，记步长h= xk+1-xk

12、根据微分中值定理y'=y'xk+h=yxk+1-y(xk)xk+1-xk=yk+1-ykh上式中xk,xk+1, 则0<<1而根据题设y'=f(x,y),y'xk+h=f(xk+h,yxk+h)故有：yk+1=yk+hf(xk+h,yxk+h)由于f(xk+h,yxk+h)是（xk,xk+1)中满足微分中值定理的点的导数y'()因此，f(xk+h,yxk+h)称为平均斜率则： yk+1=yk+h如取xk处的斜率fxk,yk作为平均斜率的近似值，则得欧拉公式，它的近似等级为Oh2。如取xk，xk+1处得斜率1, 2的平均数作为平均斜率的近似值，

13、则得改进欧拉公式，它的近似等级为Oh3。公式为：1= fxk,yk 2=fxk+1,yk+1=f(xk+1,yk+hfxk,yk) yk+1=yk+h(1+2)2由此想到，我们是否可以在xk,xk+1中多取几个点，再以它们的某种平均值作为，把精度进一步提高呢？假设在xk,xk+1中取n个点xk+1h,xk+2h,xk+nh其中01,2n1 ，斜率为1 ,2,n如果取第一个点为xk，则1=0,1=f(xk,yk)取平均斜率的计算方法为=11+22+nn则yk+1=yk+hm=1nmkm其中：m= fxk+mh,yk+hj=1m-1mjj恰当的选取公式中的常数，就可以使精度尽量

14、地高。这就是龙格库塔的一般公式，以此类推四阶龙格库塔法的一般公式为：yk+1=yk+h6(1+22+23+4)其中：1=fxk,yk2=fxk+h2,yk+h213=fxk+h2,yk+h224=f(xk+h,yk+h3)该公式的局部截断误差可达O(h5)2.3重力沉降因素在初始高度大于10Km的喷发中，半径大于100微米的较大的颗粒在数小时内将降落到地面，通常这些粒子不经历区域的或是全球的传播，因此它们的威胁是有限的。然而，更小的粒子，尤其是所谓的3级粒子8,即小于十几微米的细灰能够停留在空中数天以上，能够对航空安全产生巨大的危害9-10。在这项研究中我们主要致力于这种类型的粒子的传播研究，

15、由于这种类型粒子的直径非常小，重力沉降因素可以忽略。2.4 湍流因素Jkfeng Peng7研究表明：对于所有的湍流扩散系数值，粒子的分布集中在吸收结构上，最大的粒子浓度和吸收结构是相一致的。因为吸收结构是不依赖于湍流扩散的，所以它们能够被用于作为粒子最高浓度的指示器。湍流扩散系数的值仅仅改变粒子在吸收结构周围的分布。3 计算实例选取美国1958-1997年全球大气精确轨道参数数据，根据美国数据格式编制方法，研制了给定某一时间段某一特定地区检索各不同高度层上的风速和风向数据的软件,按月份检索东西、南北向风速资料，并按月计算19581997年风速的平均值，应用龙格库塔算法进行数据内插，然后将南北

16、向风速和东西向风速进行矢量合成计算，确定全球每点的风速和风向，以月份、经度、纬度、风速大小、风向、高程，为主要字段建立全球风场数据库。编制软件，输入火山口经度、纬度、喷发时间、飞机的飞行高度这些基本参数，用以计算火山口周边的风场，用surfer软件做出风场图，并给出风的传播方向矢量图，预测火山灰云可能移动的方向，借签puff模型3给出的公式Rk(t+dt) = Rk(t) + W(t)dt + Z(t)dt + Sk(t)dt （1）式中Rk(t)是火山灰颗粒在t时刻的位置矢量，W是当地风速，Z是一个矢量代表扰动分布Sk是重力沉降矢量，主要依赖颗粒的尺寸。由于公式（1）的主要成分风起决定性作用

17、，因此我们忽略扰动和重力沉降的影响，利用简化公式Rk(t+dt)= Rk(t) +W(t)dt计算火山灰云的位置，从而达到追踪火山灰云的目的。3.1 1958年1月全球风场计算 如图3-1，图中计算了1958年1月全球风速，上图表示风的大小，下图表示方向和大小，圆代表风头，线段长度代表风的大小，从这两张图上看，全球风场在1月份以西风为主图3-1 全球1958年1月风场（左图是风速的大小，右图为风速的大小和方向）3.2 冰岛地区1958年4月和1997年4月风场情况对比选取1958年4月与1997年4月数据，冰岛火山（西经17°37，北纬63°37）为中心计算了当时的风场情况

18、(见图3-2、图3-3)。图3-2 1958年4月冰岛附近风场（左图是风速的大小，右图为风速的大小和方向）图3-3 1997年4月冰岛附近风场（左图是风速的大小，右图为风速的大小和方向）从1997和1958年4月份风速对比情况看，虽然相隔39年，但风速的大小和方向非常相似3.3 冰岛火山灰追踪位于冰岛南部亚菲亚德拉冰盖的艾雅法拉火山，当地时间2010年4月14日凌晨1时（北京时间9时），火山开始喷发，喷发地点位于冰岛首都雷克雅未克以东125公里，岩浆融化冰盖引发洪水，附近约800名居民紧急撤离。火山喷出的火山灰还在大气层中扩散，导致冰岛、英国、德国、波兰等多国阴天。截至16日，西部欧洲多国航线

19、已经中断( Bernd Heinold （Bernd Heinold,2012）采用COSMO-MUSCAT模型模拟显示结果见图3-4。图3-4COSMO-MUSCAT模型模拟显示2010年4月16，19，21，23日冰岛火山灰云分布（摘自Bernd Heinold,2012）本文选取1958年至1997年4月份冰岛的艾雅法拉火山(西经17°37，北纬63°37)周边大气轨道参数(h=9km)，计算风速的平均值，用surfer软件绘制风场的等值线图3-5。图35冰岛周边4月份风场分布图3-5，可以看出，在冰岛南部和东南部的欧洲上空有明显的风脊存在，喷发到空中的火山灰朝着风脊

20、运动，到达欧洲上空和大西洋上空，然后逐渐向西飘散,这与媒体报导结果基本吻合。 4 讨论大多数存在的VATD模型旨在预测火山灰云在大气中的精确的发展，这种预测非常依赖输入参数的精确性，例如原始烟柱的性能，NWP风数据，湍流分布等等。与此相反，在这项研究中提出的方法不直接预测灰云在大气中的精确轨迹，而是寻找支配灰云传输的驱动力风场的性能。把风场看做一个动力系统，在风场中，这种新的方法在一般的传输过程中确定的吸收结构也能支配火山灰的传输。在对2010年4月冰岛喷发的研究中证明，火山灰云在喷发后移向吸收结构。在灰云到达这些吸收结构之后，它们跟随这些结构的变化而变化。PUFF预测中，相同的初始烟柱却有非

21、常不一样的轨迹，尤其是在长期的预测中。这就强调了灰云轨迹预测非常依赖于NWP数据的精确性，并且预测精确的灰云轨迹是非常困难的。然而来自于这两个数据集的吸收结构确实非常相似的。因此，同拉格朗日轨迹法相比，这种新的方法在速度数据上很少倾向错误，并且对于较低分辨率的NWP数据能够实现较高的精确性。因此，吸收结构代表的区域在风场中扮演吸收因子的角色，并且能够禁得起拉格朗日轨迹法所犯的错误。计算风场远比计算粒子轨迹错误少的多，因此即使是近似的、低分辨率的数据，也能够产生可靠的结果。相对于吸收结构的算法，本文采用的美国大气精确轨道参数计算的风场比较精确。在不能实时拿到风参数的前提，这不失为一个有效的方法。

22、致谢 上海大学的尹京苑老师在本文的资料收集工作中，给予了大力支持和帮助，在此表示衷心感谢。参考文献1 Piedelievre, J.P., Musson-Genon, L., Bompay, F. MEDIA-An Eulerian model of atmospheric dispersion: first validation on the Chernoby1 releaseJ. Journal of Applied Meteorology .1990,29 (12):1205-1220.2 Draxler, R.R., Hess, G.D. An overview of the Hysp

23、lkt 4 modeling system for trajectories, dispersion, and depositionJ. Australian Meteorological Magazine. 1998,(47):295-308.3 Searcy, C., Dean, K., Strknger, W. PUFF: a high-resolution volcankc ash tracking modelJ. Journal of Volcanology and Geothermal Research . 1998,(80):1-16.4 DAmours,R.,Servranck

24、x,R.,Toviessi, J.P.,Trudel,S. The operational use of the Canadian Emergency Response Model:data assimilatkon,processing,storage,and dissemination.in: OECD Nuclear Energy Agency (Ed.), Nuclear Emergency Data Management:Proceedings of an International Workshop. Organisation for Economic Cooperation an

25、d DevelopmentJ. Nuclear Energy Agency, 1998,PP:215-221.5 Jones, A., Thomson, D., Hort, M., Devenish, B.The U.K. Met offices nextgeneration atmospheric dispersion model, NAME III. In: Borrego, C., Norman, A.L. (Eds.), Air Pollution Modeling and Its Application XVIIJ. Springer, Berlin,2007,PP:. 580-589.6 Stunder, B.J.B. An Assessment of the quality of forecast trajectoriesJ. Journal of Applked Meteorology , 1996,(35):1319-1331.7 Jkfeng Peng,Rorkk Peterson. Attractkng structures in volcankc ash transport Atmospherkc EnvironmentJ.2011,dok:10.1016/j.atmosenv.2011.05.0538 Koyaguchi, T., O