国防科技大学学报 第32卷第2期 JOURNAL OF NATIONAL UNIVERSTTY OF DEFENSE TECHNOLOGY Val.32No.22010 文章编号：1001-2486(2010)02-0068-05 基于MCMC方法的Lorenz混沌系统的参数估计· 曹小群，宋君强”，张卫民1，蔡其发2，张理论 (1.国防科技大学计算机学院，湖南长沙410073：2.中国人民解放军61741部队，北京100071) 摘要：基于贝叶斯理论，提出用马尔科夫链蒙特卡罗(MC③MC)方法来估计Lom混沌系统的未知参数。 首先导出了未知参数分布规律的后验概率密度函数：接者采用自适应Metropolis算法构造Markov链：然后截取 收敛的链序列，计算混沌系统参数的估计值。数值试验表明：该方法具有很高的估计精度，同时具有较好的抗 噪声性能。 关键词：Loem混沌系统：参数估计：马尔科夫链蒙特卡罗方法 中图分类号：0302文献标识码：A Estimating Parameters of Lorenz Chaotic System with MCMC Method CAO Xiso-qun',SONG Jun-qiang,ZHANG Wei-min',CAI Qi-f,ZHANG Li-lun' (1.Colleg f Computer.National Univ.of Defe TechologyChanga 410073,China: 2.61741 Troope of PLA.Beijing 100071,China) Abstract:Basedn Bayesian thorm,amethd is propoeed toetimate the uknown permeter of Lornz chaotic ystem uing Markov Chain Monte Carlo(MCMC)method.Firstly,the poeterior probability density functionfor unknown parmeter deduced.Sec- ondy,taking the poeterior probabilitythe invarianitition,the Adaptive Metropolis arithmisudtoorct the Makov Chain.Thiry,thevdmpeedtoccute the mathematicxpctionfthe nknow parameter.Theresuofu merical experimet show that the parameters estimated by the new method have high precisionnd the noise is filtered ffectivelyfrom obeervations at the same time. Key words:Loren chaotic ystem;parmeter estimation;Markov Chain Monte Cadlo method 自Lo©z于1963年发现第一个混沌吸引子后，混沌控制和同步大量应用于信息科学、保密通信、化 学和生命科学等领域，相应的理论和方法已经得到广泛发展1引。在实际应用中，由于混沌系统的复杂 性，某些参数不可测；或者出于保密通信的要求，系统的某些参数不可知。因此在混沌控制和同步领域， 估计混沌系统的未知参数成为一个关键问题。近几年研究人员不断提出估计混沌系统参数的新方 法刀。混沌系统的参数估计问题实质上是动力系统反问题，和正问题不同，它是通过试验或运行中的 观测数据来反求模型中的未知参数。由于测量信息不准确、不充分和系统自身的非线性等性质，因而反 问题的求解经常是不确定的，即解不一定存在；即使解存在也不唯一，在解存在唯一条件下也不稳定（解 不连续依赖于观测数据)，因此混沌动力系统参数估计问题必须采用特殊方法4-刃求解。 1MCMC方法 在贝叶斯统计方法【0-)中，认为总体分布中的未知参数向量m本身就是随机向量，在对m进行统 计推断时，首先必须对m设置一个先验分布p(m)。先验分布是在观测前所有关于未知参数信息的概 率表述。在观测发生后，使用贝叶斯公式可将未知参数的先验分布改进为后验分布。根据贝叶斯理论， 参数的后验分布、观测信息和先验分布之间有如下关系： p(mld)=p(dm)p(m)lp(d) (1) ·收稿日期：2009-10-13 基金项目：国家自然科学基金资助项目(4075064：40505023) 作者葡介：曹小群(1980一)，男，助理研究员，博士。 万方数据
国 防 科 技 大 学 学 报 筮丝鲞箜兰塑坦龇堡幽唑垦塑鲤驽翌堡婴婴堕鱼塑幽 巡：丝№：至垫!Q 文章编号：1001—2486(2010)02—0068—05 基于MCMC方法的Lorenz混沌系统的参数估计。 曹小群1，宋君强1，张卫民1，蔡其发2，张理论1 (1．国防科技大学计算机学院，湖南长沙410073；2．中国人民解放军61741部队，北京100071) 摘要：基于贝叶斯理论，提出用马尔科夫链蒙特卡罗(MCmC)方法来估计Izam伫混沌系统的未知参数。 首先导出了未知参数分布规律的后验概率密度函数；接着采用自适应Metropolis算法构造Msrkov链；然后截取 收敛的链序列，计算混沌系统参数的估计值。数值试验表明：该方法具有很高的估计精度，同时具有较好的抗 噪声性能。 关键词：L肌北混沌系统；参数估计；马尔科夫链蒙特卡罗方法 中圈分类号：0302 文献标识码：A Estimating Parameters of Lorenz Chaotic System with MCMC Method CAO x沁耐，SONG m耐，刎悯Wei．minI，CAI够舒，ZHAN(；Li—l姐I (1．C．．zdlese越compute,．National Univ．《Dd'ense Technology。Changsh 410073，China； 2．61741 Troop．0f PLA。Beijins 100071，oana) Abstln峨：Based∞Bayesian theom，a method is弘Dp嘴ed to estimate the unknown parmnetem of I正腿抛chaotic system岫唱 MaAov Chain Monte C81-10(MCMC)method．F'ustly，tIle posterior舯蝴I)r density缸吼湎for unknown patanm傩is deduced．See· ondly，峨the posteIior prob8bility锄the invariant dimlb嘶'on，the Adaptive Metropolis出onthm is used tO construct the Markov ahin8．TKImy，the伽柙nelged 8ample8 are used to calculate the m越kj倒c∞驴∞咄I dthe unknown parameters．The resul￡s 0f nu· merieal expetime出show that the parameters emmted by the new嘣kd have Iljgh弘烈施iolI and the noise is filtered effectively from damvatiom越the 8衄e time． Key wol'ds．．Lorenz ChaOliC system；parameter estimation；Matkov a柚iII Monte Carlo method 自Lorenz于1963年发现第一个混沌吸引子后，混沌控制和同步大量应用于信息科学、保密通信、化 学和生命科学等领域，相应的理论和方法已经得到广泛发展【l。3】。在实际应用中，由于混沌系统的复杂 性，某些参数不可测；或者出于保密通信的要求，系统的某些参数不可知。因此在混沌控制和同步领域， 估计混沌系统的未知参数成为一个关键问题。近几年研究人员不断提出估计混沌系统参数的新方 法【．dJ。混沌系统的参数估计问题实质上是动力系统反问题，和正问题不同，它是通过试验或运行中的 观测数据来反求模型中的未知参数。由于测量信息不准确、不充分和系统自身的非线性等性质，因而反 问题的求解经常是不确定的，即解不一定存在；即使解存在也不唯一，在解存在唯一条件下也不稳定(解 不连续依赖于观测数据)隋】，因此混沌动力系统参数估计问题必须采用特殊方法H。’求解。 1 MaMC方法 在贝叶斯统计方法[10-11 3中，认为总体分布中的未知参数向量胍本身就是随机向量，在对m进行统 计推断时，首先必须对nl设置一个先验分布p(胍)。先验分布是在观测前所有关于未知参数信息的概 率表述。在观测发生后，使用贝叶斯公式可将未知参数的先验分布改进为后验分布。根据贝叶斯理论， 参数的后验分布、观测信息和先验分布之间有如下关系： p(肼l d)=v(d m)p(m)／p(d) (1) ·收稿日期：2009．10．13 基金项目：国家自然科学基金资助项目(40775064；4啦,05023) 作者简介：曹小群(1980一)．男．助理研究员．博士。 万方数据
首小群，等：基于MCMC方法的Lorenz混沌系统的参数估计 69 其中，mT=(m1,m2,…,mw),表示长度为N的未知随机向量，d=(d,d,…,dw)是包含M个观测资 料的列向量。因为规测数据已经给出，所以p(d)是一个与m无关的常数，式(1)可写成 p(mld)op(dm)p(m) (2) 式(2)是本文利用MCMC算法估计混沌系统未知参数的出发点。其中，条件概率密度p(Im)又称似然 函数，也写作L(m),表示估计的模式参数和实际观测数据的拟合程度，值越大表明拟合程度越好， 反之越差。将先验信息和观测信息结合起来，就得到反映未知随机向量m综合信息的后验概率密度 p(m),它定义在整个解空间，表示了问题的“完全”解，即可能解集合中任意一个解出现的概率。 理论上从后验分布函数可以获得未知随机向量的任何统计特征，但在实际应用中除了极其简单的 情形，p(m1d)都无明确的数学表达式。另外，采用通常的数值积分方法（例如，Monte Carlo方法2-） 也存在极大的困难：计算量随向量维数的增多而呈指数增长。上述原因使得贝叶斯方法很难直接解决 实际问题，而近期发展的MCMC算法能有效克服该困难。 MCMC算法能对定义在高维随机向量空间M上无明确数学表达式的概率分布pP进行抽样，基本思 想是产生大量服从分布的随机向量序列{m1,m2,“,m,{,其中I为抽样数。如果向量序列满足马尔 可夫性质：向量m1的产生仅依赖于前一个向量m:,而与过去时刻i-1,i-2,“,1的状态向量m1, -2,·,m,都无关，则称此序列为马尔可夫链。马尔可夫性质的另一种表述是：若抽样算法当前访问 的是严点，则下一步访问另一点m:的概率只依赖于m,而与以前访问点无关。马尔科夫性质意味着 抽样算法完全可由转移概率矩阵P描述，矩阵元素P,表示算法在当前访问m,的条件下接着将要访问 m:的条件概率。按照构造Markov链所用转移概率矩阵的不同，MCMC方法的主要抽样算法有：Gibbs抽 样器算法W,Metropolis--Hastings算法和自适应Metropolis算法O。 本文采用自适应Metropolis算法以混沌系统未知参数的后验分布为不变极限分布来构造Markov链。 自适应Metropolis算法是Haario在2OOI年提出的一种改进MCMC抽样算法j。与传统的Metropolis- Hastings算法相比，自适应Metropolis算法不需要预先确定参数的推荐分布，而是由后验参数的协方差矩 阵来估算参数分布4。后验参数的协方差矩阵在每一次迭代后都要自适应地调整。如此，第：步参数 的推荐分布就是均值为元、协方差矩阵为C:的多元正态分布。协方差的计算公式如式(3)所示，在初 始次迭代中，协方差矩阵C:取固定值C。,之后进行自适应更新。m:是未知模式参数向量m中某个 元素在第讠次迭代的估计值。 [Cor i≤o C= (3) 5aCov(io,…,m:-1）+sel,i>io 其中，e是一个非常小的数，引进目的是确保C:不成为奇异矩阵，Cow(no,…,m:-t)表示mo,…,m:-1的 协方差阵；5是比例因子，依赖于未知随机向量空间的维数d,目的是保证接受率在一个合适的范围内； 在本文中e和s4分别取为：e=10-6,=2.41d。1为d维单位矩阵。当进行第i+1次选代时，由式 (3)可导出协方差的计算式： C1='c+流m,-(i+1)m,m+mm+e】 (4) 其中，m:-,和m:是前面i-I次和i次选代的参数均值。自适应Metropolis算法的计算流程如下： (1)设定i=0,对不同变量进行初始化； (2)随机量的生成和接受，构造Markov链： (a)利用式(3)计算协方差矩阵C: (b)产生服从正态分布的推荐参数值m~N(m,C,): (o)利用下式计算接受概率a=血{1，pdm2m} 'p(dlm.)p(m.)i (d)产生服从均匀分布的随机数4~U(0,1)； (e)若u<a,则接受m+1=m”,否则mi=m。 (3)重复上面的步骤(a)~(e),直到产生预先指定数量的样本为止。 万方数据
曹小群，等：基于MCMC方法的Lorenz混沌系统的参数估计 69 其中，肼“=(m，，m：，…，mN)，表示长度为Ⅳ的未知随机向量，d’=(d。，d2，…，如)是包含肘个观测资 料的列向量。因为观测数据已经给出，所以P(d)是一个与m无关的常数，式(1)可写成 p(m d)∞p(dl m)p(m) (2) 式(2)是本文利用MCMC算法估计混沌系统未知参数的出发点。其中，条件概率密度p(dl小)又称似然 函数，也写作￡(dl m)，表示估计的模式参数和实际观测数据的拟合程度，值越大表明拟合程度越好。 反之越差。将先验信息和观测信息结合起来，就得到反映未知随机向量廊综合信息的后验概率密度 P(flitI d)，它定义在整个解空间，表示了问题的“完全”解，即可能解集合中任意一个解出现的概率。 理论上从后验分布函数可以获得未知随机向量的任何统计特征，但在实际应用中除了极其简单的 情形，P(小ld)都无明确的数学表达式。另外，采用通常的数值积分方法(例如，Monte Carlo方法m川1) 也存在极大的困难：计算量随向量维数的增多而呈指数增长。上述原因使得贝叶斯方法很难直接解决 实际问题，而近期发展的MCMC算法能有效克服该困难。 MCMC算法能对定义在高维随机向量空间M上无明确数学表达式的概率分布P进行抽样，基本思 想是产生大量服从分布P的随机向量序列{m。，m：，…，m，}，其中，为抽样数。如果向量序列满足马尔 可夫性质：向量m川的产生仅依赖于前一个向量mi，而与过去时刻i—l，i一2，…，l的状态向量m㈠， mH，…，肼。都无关，则称此序列为马尔可夫链。马尔可夫性质的另一种表述是：若抽样算法当前访问 的是m，点，则下一步访问另一点m；的概率只依赖于m，，而与以前访问点无关。马尔科夫性质意味着 抽样算法完全可由转移概率矩阵P描述，矩阵元素船表示算法在当前访问历，的条件下接着将要访问 m；的条件概率。按照构造Ms虫ov链所用转移概率矩阵的不同，MCMC方法的主要抽样算法有：Gibbs抽 样器算法【1．】、Metropolis．Hastings算法㈣和自适应Metropolis算法【16】。 本文采用自适应Metropolis算法以混沌系统未知参数的后验分布为不变极限分布来构造Markov链。 自适应Metropolis算法是Haario在2001年提出的一种改进MCMC抽样算法【14j。与传统的MetropolisHastings算法相比，自适应Metropolis算法不需要预先确定参数的推荐分布，而是由后验参数的协方差矩 阵来估算参数分布u‘】。后验参数的协方差矩阵在每一次迭代后都要自适应地调整。如此，第i步参数 的推荐分布就是均值为厩、协方差矩阵为G的多元正态分布。协方差的计算公式如式(3)所示，在初 始i。次迭代中，协方差矩阵G取固定值co，之后进行自适应更新。畹是未知模式参数向量肌中某个 元素在第i次迭代的估计值。 G={赫。"孟小礼震 (3) 其中，￡是一个非常小的数，引进目的是确保G不成为奇异矩阵，c叫(琥。，…，晚一。)表示‰，…，畹．，的 协方差阵；妇是比例因子，依赖于未知随机向量空间的维数d，目的是保证接受率在一个合适的范围内； 在本文中e和&分别取为：￡=10～，知=2．42／d。厶为d维单位矩阵。当进行第i+1次迭代时，由式 (3)可导出协方差的计算式： G+l=掣G+孚【iWtl．1而j．1一(i+1)厩imwTi+mimTl+也】 (4) 其中，厩一，和厩i是前面i—1次和i次迭代的参数均值。自适应Metropolis算法的计算流程如下： (1)设定i=O，对不同变量进行初始化； (2)随机量的生成和接受，构造Markov链： (a)利用式(3)计算协方差矩阵G； (b)产生服从正态分布的推荐参数值m。一N(虺，G)； (c)利用下式计算接受概率口=疵{l，号舌守端)； (d)产生服从均匀分布的随机数n—u(0，1)o (e)若Ⅱ<口，则接受mi+l=m。，否则mi+l=豫。 (3)重复上面的步骤(a)一(e)，直到产生预先指定数量的样本为止。 万方数据
70 国防科技大学学报 2010年第2期 自适应Metropolis算法的最大优点是推荐分布随计算过程自动更新，不需要预先指定。同时，相比 Metropolis-Hasting算法]，由于参数不再需要分组更新，所以计算量大大诚少。 2 Lorenz混沌系统的参数估计 下面以典型的Lorenz混沌系统为例，说明利用MCMC算法对混沌系统的未知参数进行估计的有效 性。orez系统可由如下的常微分方程组表示： [dx/dt=a(y-x) dyldt=Yx-x-y (5) dz/dt=xy-bz 系统(5)的一个简单物理实现是在下方加热而在上方冷却时热对流管中的环流。其中x表示流体速 度，y和z分别为水平和垂直温差；o和y分别与流体的Prandtl数和Rayleig数成比例，b是与空间相关 的常数。当参数。=10，Y=28和b=8/3时，式(5)表示的系统呈现混沌现象。数值试验中采用四阶龙 格-库塔算法求解常微分方程，步长设置为h。观察数据的采集过程如下：先让oem系统自由演化， 在经历暂态过程之后任取混沌状态的某一点作为初值，并以此时刻作为初始时刻，然后任其演化至 时刻。这样就得到Loem系统在离散时间序列{0，h,2h,…,Th}上的标准状态(x(),y(),z()(i =0,1,2,·,T),而且这些状态量都处于系统的混沌状态。标准状态量可直接作为观测数据(x。(), y。(),(t:),也可加人随机噪声得到更真实的模拟观测资料。 本文的目标是在已知部分观测资料的条件下，采用MCMC方法估计出Lor©m混沌系统的多个未知 参数。下面利用贝叶斯公式来导出Loez混沌系统未知参数的后验概率密度函数。假设方程(5)中影 响模式状态量的3个参数(a,Y,b)都是未知的，即对三个参数进行联合反演。模式初始状态量是(x0, yo,o),观察资料为d=(x(),y(t:),()(i=1,2,…,6),表示只使用部分观察资料。根据贝叶斯 公式(1)，且不考虑分母表示的常数项，则混沌系统未知参数的后验分布p(md)=p(a,y,b|d)通过 下式进行计算： p(a,y,bld)=p(dla,Y,b)p(a,r,b) (6) 式(6)已经假设d=M(a,Y,b)+w,其中M(a,Y,b)表示Lom系统的离散正模式，w为包含观测误差的 独立分布的随机噪声，并且假定随机噪声w服从均值为零、标准偏差为a。的正态分布，即w~N(0,)。 同时假设由数值求解算子M引进的误差远小于观测误差。通过上述假设，可以得到以下形式的似然函 数： (dlm)=dlg7,b)=2xdrm-ld-Ma,ybP2a2J (7) 式中n表示观察资料数量，！·I表示欧几里得范数。 在贝叶斯推理中，首先要设定每个未知参数的先验分布。在本文中假设所有未知参数的先验分布 都满足独立的均匀分布，参数a、y和b分别满足U(σ)、U(Y)和U(b)表示的均匀分布，则总的先验分 布表示为 p(m)=U(a)U(Y)U(b) (8) 均匀分布是一种最简单的先验分布，虽然只能指定未知参数变化的上下界，但是可以缩小对参数随机抽 样的目标区域，有利于提高参数的估计精度，该结论在下面的数值试验中将得到验证。一般来说，可以 通过经验知识和历史统计确定更准确和复杂的先验分布。自适应Metropolis算法的一个优点是对于m 的任何先验分布都能够收敛于目标分布。由式(6)、(7)和(8)可得，在给定观测数据的条件下混沌系统 未知参数的后验概率密度函数为 (,74)=(-ld-M(o)U(U()U() (9) 3数值试验结果 在数值试验中，模式初始状态量(o,y0,)取为(-0.2028,2.5418,24.0873)，是orez系统混沌状 万方数据
国防科技大学学报 2010年第2期 自适应Metropolis算法的最大优点是推荐分布随计算过程自动更新，不需要预先指定。同时，相比 Metropolis．Hastjn挚算法‘嘲，由于参数不再需要分组更新，所以计算量大大减少。 2 Lorenz混沌系统的参数估计 I燃dx／dt三=孑o(y二-：x一)y ㈣ 时刻。这样就得到Lorenz系统在离散时间序列{0，h，2h，…，吼t_E的标准状态(茹(t；)，Y(t；)。：(t；))(i yo，钿)，观察资料为d=(茗。(t；)，％(ti)，％(t。))(i=1，2，…，6)，表示只使用部分观察资料。根据贝叶斯 L(d1 m)=L(dI d，7，6)=7i=：lri丽e印[一m d—M(a，y，b)0 2，(2d：)] (7) p(a，y，b d)=石￡==1：币a币[一I d一肘(a，y，6)I 2／(2do)]‘u(口)U(y)u(6) (9) 3数值试验结果 在数值试验中，模式初始状态量(聋o，yo，钿)取为(一0．2028，2．5418，24．0873)，是Lamr坛系统混沌状 万方数据
曹小群，等：基于MCMC方法的Lorm混沌系统的参数估计 71 态中一点。以该状态为初值，取积分步长h=001,设定时间积分区间为[0,8]（单位为s),然后对模式 进行积分，就可以得到系统在离散时间序列0h,1h,2h,·,800h上的标准状态变量值。试验中只取部 分状态变量值作为观测量，抽取0h,40h,80h,·,200h时刻的状态量作为n=6个观测数据。在所选状 态量上叠加高斯随机数作为观测噪声，观测噪声的均值取为0、标准偏差取为：。=0.1。未知参数σ、Y 和b的先验分布分别取如下形式： 「1/15,0<g<15 U(a)= (10) l0, otherwise f0.02,0<y<50 U(y)= 0, (11) otherwise 0.1,0<b<10 U(b)=0, (12) otherwise 以上面加人高斯噪声的状态量作为观测资料。三个先验分布概率密度取式(10)~(12)所示的表达 式，则未知参数的后验概率密度函数表示为 p(ma）=2xep-嘉2ld-M(am门以o)u)u(b） (13) 其中d=(x,y。,。)表示所选的第i个观测向量；M(m)表示对向量m=(a,Y,b)进行更新后，由 Lorenz离散模式积分到对应观测时刻得到状态量(x,y,z)的计算值。 34 31.5 33 之 30.5 32 30- 3.1 抛 29.5 3 29 29 28.5 27.5 27 27 26 26.3 100 2003004005006007008009001000 100 2003004005006007008009001000 取样状数 取样次数 图1参数b的Markov链 图2参数y的Mkow链 Fig.1 Markov chain of parameter b Fig.2 Markov chain of parameter y 确定后验概率密度函数的表达式后，就可以利用MCMC方法对未知参数进行随机抽样，预先指定 抽样次数为1O00,参数初始值取对应的先验均匀分布中的随机值。自适应Metropolis算法在进行抽样 时，每次对未知随机向量进行更新后，Lo©z模式都要在研究的时间区间上进行一次积分，计算后验概 率密度函数的值。然后，通过和上一步后验密度函数值的比较，来决定是否对未知随机向量进行更新。 三个未知参数的Makv链分别如图1、2和3所示，可以看出马尔科夫链在大约经历800多次迭代后达 到收敛。取800~1000步的链序列并采用后验均值方法计算参数估计值，分别得到。=10.001、Y= 28.021和b=2.6828。分析图中所示的试验结果可得如下结论：采用MCMC方法估计的未知参数具有 较高的准确性；在相同的取样次数情况下，参数σ具有最高的估计精度，这与指定的先验分布的准确性 有关：相对于另外两个参数，参数。的先验分布更接近真实值。采用MCMC方法估计出(a,Y,b)的值 后，代入Lo阳m混沌系统的离散模式可以得到所有状态量(x,y,z)的估计轨迹。图4、5和6分别显示了 x、y和z的估计轨迹、观测值和离散时间点上的真实值，从图中可知，估计轨迹和真实值在7.53之前非 常吻合，验证了用MCMC算法同时对多个未知参数进行估计的正确性，同时也表明了从部分较少的观 测数据成功估计出混沌系统参数的可能性。由于混沌系统的特性，即使参数估计值与真实值非常接近， 万方数据
曹小群，等：基于MCMC方法的Lomm混沌系统的参数估计 71 态中一点。以该状态为初值，取积分步长h=O．Ol，设定时间积分区间为[0，8](单位为s)，然后对模式 进行积分，就可以得到系统在离散时间序列Oh，lh，2h，…，800h上的标准状态变量值。试验中只取部 分状态变量值作为观测量，抽取Oh，40h，80h，…，200h时刻的状态量作为，l=6个观测数据。在所选状 态量上叠加高斯随机数作为观测噪声，观测噪声的均值取为O、标准偏差取为以=0．1。未知参数口、y 和b的先验分布分别取如下形式： m)-一1／．15’裂5 (10) 附，=铲’2= ㈣， u(6)：F，o：6：10 (12) 以上面加入高斯噪声的状态量作为观测资料。三个先验分布概率密度取式(10)。(12)所示的表达 式，则未知参数的后验概率密度函数表示为 p(m d)2西南expE一刍善n di一∥(m)I|2]．嘶)附)啪) (13) 其中∥=(茗：，，，：，：。i)表示所选的第i个观测向量；∥(m)表示对向量朋=(盯，y，6)进行更新后，由 Lorenz离散模式积分到对应观测时刻得到状态量(菇，Y，z)的计算值。 图1参数b的Markov链 图2参数y的MsI-kov链 Fjg．1 Msrkov chain 0f P8础m既盯6 Fig．2 M田kov cham 0f pa功蝴y 确定后验概率密度函数的表达式后，就可以利用MCMC方法对未知参数进行随机抽样，预先指定 抽样次数为1000，参数初始值取对应的先验均匀分布中的随机值。自适应Metropolis算法在进行抽样 时，每次对未知随机向量进行更新后，Lorenz模式都要在研究的时间区间上进行一次积分，计算后验概 率密度函数的值。然后，通过和上一步后验密度函数值的比较，来决定是否对未知随机向量进行更新。 三个未知参数的Markov链分别如图l、2和3所示，可以看出马尔科夫链在大约经历800多次迭代后达 到收敛。取800～1000步的链序列并采用后验均值方法计算参数估计值，分别得到d=10．001、y= 28．021和b=2．6828。分析图中所示的试验结果可得如下结论：采用McMC方法估计的未知参数具有 较高的准确性；在相同的取样次数情况下，参数仃具有最高的估计精度，这与指定的先验分布的准确性 有关：相对于另外两个参数，参数口的先验分布更接近真实值。采用MCMC方法估计出(盯，y，b)的值 后，代入Lorenz混沌系统的离散模式可以得到所有状态量(茗，y，z)的估计轨迹。图4、5和6分别显示了 菇、Y和z的估计轨迹、观测值和离散时间点上的真实值，从图中可知，估计轨迹和真实值在7．5s之前非 常吻合，验证了用MCMC算法同时对多个未知参数进行估计的正确性，同时也表明了从部分较少的观 测数据成功估计出混沌系统参数的可能性。由于混沌系统的特性，即使参数估计值与真实值非常接近， 警_诗=雹自 撼^o量要雹．‘ 万方数据
72 国防科技大学学报 2010年第2期 10.5 20 口观测值 15 估计值 真实值 95 9 8.5 65 45 1002003004005006007008009001000 20 取样次数 时间 图3参数g的Markow链 图4状态量x的比较图 Fig.3 Markov chain of parameter o Fig.4 Comparison of state variable x 口观值 45 ……估计值 d a 测值 0真实值 40 估计值 真去值月 0 30 5 0 的 馨 0 中 50 -10 15 d 0 15 10 1 2 3 4 6 2 3 4 时阀 时间 图5状态量y的比较图 图6状态量：的比较图 Fig.5 Comperison of state variable) Fig.6 Comparison of state variable 也仅能在短时间内有意义。从图4~6可知，随着系统的长时间演化，在7.58之后估计参数所表征系统 的状态量开始偏离真实状态量。 4结论 本文在贝叶斯理论的框架下，利用马尔科夫链蒙特卡罗(MCMC)算法来估计Loem混沌系统中的 未知参数。主要结论归纳如下： (1)即使观测资料中包含观测噪声，利用MCMC方法都能准确地估计出混沌系统参数；利用少量观 测资料，估计参数所表征系统的状态量和真实状态量之间的误差很小。验证了利用MCMC方法估计混 沌系统参数的有效性。 (2)在相同的随机抽样次数情况下，如果对某参数指定的先验分布准确性越高，则此参数的估计精 度也越高。 总之，MCMC方法是一种通过构造Markov链进行随机模拟的方法，通过和贝叶斯理论进行结合，构 造Markov链的过程实际上是在由先验分布界定的空间内对未知参数进行随机最优搜索的过程。本研 究的结果表明，MCMC方法和许多确定性方法刃一样能准确估计混沌系统的多个未知参数，同时具有 较好的抗噪声性能。 (下转第145页) 万方数据
国防科技大学学报 2010年第2期 图3参数口的№链 rig．3 Mazkov chain of parameter口 图4状态量髫的比较图 №．4 c0呵赫at stateⅦIiable善 图5状态量Y的比较图 图6状态量z的比较图 隐．5 c‘肋parism 0f state v越'iablc Y ng．6 c‘皿pBri8∞0f irateⅦ诎j 也仅能在短时间内有意义。从图4—6可知，随着系统的长时间演化，在7．5s之后估计参数所表征系统 的状态量开始偏离真实状态量。 4结论 本文在贝叶斯理论的框架下，利用马尔科夫链蒙特卡罗(MCMC)算法来估计lorenz混沌系统中的 未知参数。主要结论归纳如下： (1)即使观测资料中包含观测噪声，利用MCMC方法都能准确地估计出混沌系统参数；利用少量观 测资料，估计参数所表征系统的状态量和真实状态量之间的误差很小。验证了利用MCMC方法估计混 沌系统参数的有效性。 (2)在相同的随机抽样次数情况下，如果对某参数指定的先验分布准确性越高，则此参数的估计精 度也越高。 总之，MCMC方法是一种通过构造Ms／'kov链进行随机模拟的方法，通过和贝叶斯理论进行结合，构 造Markov链的过程实际上是在由先验分布界定的空间内对未知参数进行随机最优搜索的过程。本研 究的结果表明，MCMC方法和许多确定性方法C4-7]一样能准确估计混沌系统的多个未知参数，同时具有 较好的抗噪声性能。 (下转第145页) 万方数据
郭熙业，等：收发合置情况下海底混响仿真 145 4结论 文中提出了一种混合模型时域仿真方法，该方法通过设置散射体数量顾值实现了点散射模型方法 与单元散射模型方法之间的无缝切换，既克服了单元散射模型方法在散射体数量较少情况下高斯分布 假设带来的计算误差，又克服了点散射模型方法在散射体数量较多情况下计算复杂的问题。针对发射 信号为窄带信号的情况，为降低计算量，设计了先计算混响包络进而得到混响信号的方法，适用于满足 采样定理的任意采样率。最后，针对不同带宽、不同散射体密度情况仿真了混响信号，仿真结果与理论 混响在频谱特性以及统计特性方面的一致性反映了仿真结果的合理性，针对统计特性的仿真试验则验 证了混合模型方法的优越性。 参考文献： [1]Abraham D A,Lyons A P.Simulation of Non-Rsyleigh Reverberation and Clutter[J].IEEE J.Oceanic Eng.,2004,29(2):347-362 [2]徐新盛，张燕，李海森，等.海底混响仿真研究[J].声学学报，19%，23(2)：141-148. [3]蔡平，聚国龙，葛风拥，等.界面混响信号的仿真研究[J].哈尔族工程大学学报，2000,21(4)：31-35. [4]Abraha D A,Lyons A P.Novel Physical Interpretations of K-distributed Reverberation[]].IEEE J.Oceanic Eng.2002,27(4):800-813 [5]LePage K D.Statisticn of Broed-band Bottom Reverbertion Predictions in Shallowate Waveguridea[J].IEEEJ.Oceanic Eng..2004,9():330 -346. [6]杨俊龄，黄晓斌，万建伟.基于高斯序列乘积和的相千非高斯分布杂波模报新方法[】.国防科技大学学报，2006,28(5)：58-62. (上接第2页) 参考文献： [1]刘播才，梁骁明.Hcm混沌系统的广义预测控制与同步快速算法[J].物理学报，2005,54(10)：4584-4589 [2]李骊香，彭海朋，卢辉斌，等.H如m混沌系统的追踪控制与同步[J].物理学报，2001,50(4)：629-632. [3】王兴元，武相军.不确定Qm系统的参数辨识与自适应同步[J】.物理学报，2006,55(2)：65-609. [4]贺明峰，穆云明，赵立中，基于参数自适应控制的混礼问步[J】.物理学报，2000,49(5)：830-832. [5]关新平，彭海朋，李厢香，等.m混沌系统的参数牌识与控制[J】.物理学报，2001,50(1)：25-29. [6]戴栋，马西奎，李高才，等.一种基于迹传算法的混沌系统参数估计方法[J】.物理学报，2002,51(1)：24S9-2462. [7]李丽香，彭海朋，杨义先，等.基于混袍蚂蚁群算法的L混沌系统的参数估计[J刀.物理学报，2007,56(1)：51-55. [8]黄思训，伍荣生.大气科学中的数学物理问题[M].北京：气象出版社，201：411-412. [9]Gilks WR,Richardson S,Spiegehalter DJ.Markov Chain Monte Carlo in Practice [M].London:Chapman&Hall,1996. [10]张金桃，多层验前正态总体动态参数的Bas融合估计{J].国防科技大学学报，2003,25(4)：7-101. [11]幸剩被，谢红卫，张金槐.考忠验前信息可信度时的Bym估计[)J】.国防科技大学学报，003,5(4)：107-110， [12] 王同权，张若棋，沈永平，等.高能离子在材料中输运的蒙特卡罗模拟[J].国防科技大学学报，2001,23(1)：105-109 [13] 王同权，张树发，王尚武，等，电子在材料中输运的蒙特卡罗模拟[】.国防科技大学学报，2000,22(11)：81-84. [14]Smith AFM.Robert GO.Bayesian Computation Via the Gibbe Sampler and Relatod Markow Chain Monte Carlo Method[J].Joumal f Royl Statistical Society Series B,1993,55:3-23 [15]Tiemey L.Markor-chain for Exploring Posterior Distributiona[J].Armal of Satistics,994,:1701-1762. [16]Haario H,Saksman E.Tamminem J.An Admptive Metropolis Algorithm[]].Bemoulli,2001,7(2):23-242. 万方数据
郭熙业，等：收发合置情况下海底混响仿真 145 4结论 文中提出了一种混合模型时域仿真方法，该方法通过设置散射体数量阈值实现了点散射模型方法 与单元散射模型方法之间的无缝切换，既克服了单元散射模型方法在散射体数量较少情况下高斯分布 假设带来的计算误差。又克服了点散射模型方法在散射体数量较多情况下计算复杂的问题。针对发射 信号为窄带信号的情况，为降低计算量，设计了先计算混响包络进而得到混响信号的方法，适用于满足 采样定理的任意采样率。最后，针对不同带宽、不同散射体密度情况仿真了混响信号，仿真结果与理论 混响在频谱特性以及统计特性方面的一致性反映了仿真结果的合理性，针对统计特性的仿真试验则验 证了混合模型方法的优越性。 参 Ill [2】 C31 [4】 考文献： Abraham D^。Lyons AP．Simll雠Oll 0f N删eish Revelbe蒯on and C]utteor[J]．IEEE J．O嘲mlc Eng．。2004，29(2)：347—362． 徐新盛，张燕，李海森，等．海底混响仿真研究【J]．声学学报，1998。23(2)：141—148． 蔡平。粱国龙。葛凤翔，等．界面混响信号的仿真研究[J]．哈尔滨工程大学学报。2000，21(4)：31—35． ^bdm D A，Lycm AP．Novel Phygical lrII硐鹏咖删d K．dim．ilmuxi Ile嘲bcmdm[J]．IEEE J．OCe捌C E唱．，2002，27(4)：800-813． 【5】 I．ePase K D．铀ti或icB of Bmed-bend Bottom Reved删on Predi甜mm in sIIallaW-啪垤w瞅毽uid科J]．IEEE J．Oo目mic Ens．，21304，29(2)：3∞ 一346． [6】杨俊岭，黄晓斌，万建伟．基于高斯序列乘积和的相干非高斯分布杂渡模拟新方法[J]．国防科技大学学报，2006．2s(5)：58一配． (上接第72页) 参考文献： [15] [16] 刘福才。粱晓明．}悟啪混沌系统的广义预测控制与同步快速算法[J]．物理学报。加晒，54(10)：4584—4589． 李丽香，彭海朋。卢辉斌．等．}I自m混沌系统的追踪控制与同步[J]．物理学报，200l，50(4)：629—632． 王兴元，武相军．不确定Qm系统的参数辨识与自适应同步[J]．物理学报，2006．55(2)：6晒一609． 贺明峰，穆云明，赵立中．基于参数自适应控制的混沌同步[J]．物理学报。2000，49(5)：830—832． 关新平．彭海朋，李丽香．等．L肿眦混沌系统的参数辨识与控制[J】．物理学报．2001。50(I)：26—29． 戴栋，马西奎，李富才，等．一种基于遗传算法的混沌系统参数估计方法[J]．物理学报，2002，51(1)：2459—2462． 李丽香，彭海朋，杨义先，等．基于混沌蚂蚁群算法的L棚m混沌系统的参数估计[J]．物理学报，2007，56(1)：51—55． 黄思训，伍荣生．大气科学中的效学物理问题[M]．北京：气象出版社，2001：411—412． Cilk8 WR，Richardson S。sp；嘲eI}lalt贯D J．Malkov Ch抽Monte Carlo in P眦6∞[M]．h雠hI：Q嘲舶姐&Hall，1996． 张金槐．多层验前正态总体动态参数的脚嘲融合估计[J]．国防科技大学学报，2003，25(4)：99一101． 李鹏波，谢红卫，张金槐．考虑验前信息可信度时的B町∞估计[J]．国防科技大学学报，2003，25(4)：107—110． 王同权．张若棋，沈永平，等．高能离子在材料中输运的蒙特卡罗模拟[J]．国防科技大学学报，200l，23(1)：105—109． 王同权。张树发．王尚武．等．电子在材料中输运的蒙特卡罗模拟[J】．国防科技大学学报，2000，22(11)：81—84． Smith AF M．RDbelt GO．B矗廊ccq叫鲥∞Via the Gibbs Sw嘲er and Related M,lufkov al矗iII Monte Carlo眦舳[J】．Jomnal et"R删 Stati撕cal Society Series B，1993，55：3—23． Tterncy L．Markov-chah帕for Explm'吨Pom洄-ior Distributions[J]．Annals《st础击∞．1994，22：1701—1762． 胁io H，Seksmen E，Tm J．An^d叩6弛Methyls Al鲥thln[J]．Bernoulli，2001，7(2)：223—242． 1J 1J 1J 1J 1J 订幻列钉卯础刀胡明m n￡j n H ，L rL rL rl rL rL rL rl rL rL 万方数据