在matlab中模拟平稳随机过程71907

收稿日期: 2000203209 作者简介: 彭 兢(1972 - ) ,男 ,新疆阜康人 ,博士生 ,100083 ,北京。在 MATLAB中模拟平稳随机过程彭 兢 金长江(北京航空航天大学 飞行器设计与应用力学系) 摘 要: 叙述了平稳随机过程理论的基本定义 ,以及运用实谱分解方法模拟平稳随机过程的原理。 全面分析了在微型计算机上使用 MAT LAB 软件模拟平稳随机过程的注意事项 ,包括:信号源的选取 ,仿真算法和相应参数的设置 ,由于功率谱密度的不同定义而需要进行相应的调整等 3 个方面。 结合大气紊流和航空母舰在海上的随机运动的 2 个实例 ,给出了在 MAT LAB 中模拟类似过程的方法 ,并对结果进行了检验。 分析表明在 MAT LAB 中能够较好地模拟平稳随机过程。关 键 词: 平稳过程; 功率谱密度; 紊流; 舰船运动中图分类号: O 211。 61文献标识码: A 文 章 编 号: 100125965(2001)0520585204 随着微型计算机的发展 ,计算机模拟技术日益完善。MAT LAB 语言和软件包于 1980 年由美国学者 Cleve Moler 等人推出后[1],受到控制界研究者的普遍重视 ,风靡了全世界 ,尤其在教学和科研方面已经得到了广泛应用。

可以使用multisim 交互式地搭建电路原理图，对模拟、数字、模拟/数字混合电路进仿真。5、应用moldflow注塑仿真软件，对微非球面模具进行注塑过程的模拟仿真，优化了工艺参数，降低了注射过程的试模成本、提高了注塑效率。通过共性关键技术平台构建，铸造过程模拟仿真，各种加工过程模拟仿真，结构件的cad/cae技术，均质调控铝熔体纯净化和提高铸锭冶金质量。

2 线性系统中的平稳随机过程线性定常系统如图 1。x( t)G(s)y( t)图 1 线性系统 文献[2]证明 :如输入 x( t) 是平稳随机过程 ,输出 y ( t)也是平稳随机过程。 并且数学期望满足E[ y( t) ] = E[ x( t) ]∫式(5)中 , h(λ ) (λ ≥ 0) 为 G( s) 的脉冲响应函数。同时二者的功率谱密度有如下的对应关系 :Sy(ω) =| G(iω) |2Sx(ω)式(6)中 , G(iω)为传递函数 G( s)的频率响应。如将一个白噪声信号 x ( t) 通过线性传递函数 G(s) ,得到有色噪声信号 y( t) ,如图 1。 根据式(6) ,因为白色噪声的功率谱为常数 ,令它为 1 ,则 :Sy(ω) =| G(iω) |2= G3(iω) G(iω)式(7)中 , G3(iω) 为 G(iω) 的复共轭。 G( s) 被称∞0h(λ ) dλ(5)(6)(7)2001 年第27卷 第5期10 月北 京 航 空 航 天 大 学 学 报Journal of Beijing University of Aeronautics and AstronauticsOctober 2001Vol。

27 No15为成形滤波器的传递函数形式。 工程上常应用此方法来模拟平稳随机过程[3]。2 几个重要的注意事项2。 1 信号源考虑在 MAT LAB 中用模拟一个图 2 的系统。图 2 连续输入随机线性系统假定图 2 中 a 为常数 ,白噪声输入的功率谱密度为 1 ,均值为 0 ,则根据式(5) 和式(6) 可知 ,输出信号为平稳随机过程 ,均值为 0 ,功率谱密度为Sy(iω) =| G(iω) |2Swn(ω) =1a2ω2+ 1(8)式(8) 中 , Swn(ω) 为白噪声信号的功率谱密度平稳随机过程 matlab。 根据式(4) ,输出信号 y 的方差应为12 π∫ 对图 2 进行仿真 ,取 a = 1 ,采用 Euler 定步长算法 ,步长0。 01 s ,仿真时间 0～200 s,随机数种子17 ,得 到 输 出 信 号 y 的 均 值 和 方 差 分 别 为- 0。 0060 ,0。 0048。若取步长0。 001 s,仿真时间、 随机数种子和仿真算法 不 变 , 输 出 y 的 均 值 和 方 差 分 别 为0。

0013 , 0。 0005。而输出 y 的理论均值为 0 ,方差为 0。 5 ,可见输出 y( t)的模拟结果的均值基本符合理论预测 ,但方差却不正确。进一步的理论分析证实了这一点。设对于如图 2 的系统 ,在一个计算步长内输入信号为一常数 xk,直接对系统离散化 ,可得系统的模型为yk+1 = e-Δ t/ ayk + (1 - e-Δ t/ a) xk式(10)中 ,Δ t 为计算步长。 因此E[ y2(1 - e-Δ t/ a)2E[ E2又因为 E[ y2是相互独立的 ,所以 E[ xkyk] = 0。 此时(1 - e-Δ t/ a)2(1 - e- 2 Δ t/ a)若Δ t/ a → 0 ,将(12)式进行幂级数展开 ,得σ2y = Ry(0) =∞- ∞1a2ω2+ 1dω =12 a(9)(10)k+1] = e- 2 Δ t/ a+ 2e-Δ t/ aE[ xkyk] +k](11)k +1] = E[ y2k] =σ2y,同时由于 xk和 ykσ2y ==1 - e-Δ t/ a1 + e-Δ t/ a(12)σ2y =Δ t/ a + o(- Δ t/ a)22 !(- Δ t/ a)22 !2 - Δ t/ a + o(13)当Δ t → 0 时 ,σ2说明输出信号的方差与计算步长Δ t 有关。

以上分析说明 :对于随机输入连续系统 ,使用MAT LAB 的 SIMULINK工具箱进行模拟时 ,一般的线性系统离散化方法不再适用 ,必须进行修正。并且从式(13) 中得出这样的近似结论 :仿真随机输入连续线性系统时 ,对普通随机输入信号进行加权 ,乘以1/Δ t ,使输出信号的方差趋于理论值平稳随机过程 matlab。针对这一问题 ,MAT LAB 从 4。 2 版本以后 ,在SIMULINK的工具箱中提供了一个新的白噪声源:有限 带 宽 的 白 噪 声 BWN (Band2Limited WhiteNoise) 。 它有 2 个参数:噪声功率 (Noise Power) 和采样时间(Sample Time) 。 有限带宽的白噪声实际上是将离散白噪声先加权乘以阶保持后得到的连续噪声信号。例如 ,对于一个功率谱密度为 1 ,均值为 0 ,方差为 1 的白噪声序列 ,若用于采样时间为0。 01 s的系统 ,仿真时间0～1000 s ,随机数种子为 17 ,则实际输出白噪声信号的噪声功率约为 100 ,均值为- 0。 0113 ,方差为99。 2 ;若用于采样时间为0。 001 s的系统 ,仿真时间0～200 s,随机数种子仍为 17 ,则实际输出的白噪声信号的噪声功率约为1000 ,均值为0。

04 ,方差为995。 5。2。 2 仿真算法MAT LAB 中常用的数值积分算法 ,如 Runge2Kutta 法 ,Adams 法 ,Euler 法等 ,用于连续系统的计算结果精度很好。 但对于随机输入系统 ,如果数值模拟采用的时间步长小于生成白噪声的时间步长 ,或者采用变步长数值积分法 ,例如 Runge2Kutta法模拟时 ,在每个计算步长内要对输入函数多次求值 ,这都会导致输入信号变化 ,使得模拟结果与理论值不同。 因此 ,在MAT LAB 中应采用定步长法仿真 ,如 Euler 法 ,并且采样时间不得小于白噪声模块的采样时间。2。 3 功率谱的定义文献[3]中指出 ,实际使用功率谱密度的定义不统一 :如傅里叶变换及其逆变换中的系数规定有时不统一 ,单边功率谱和双边功率谱的不同等。例如文献[3]中有12 π∫则均方根、 自相关函数与功率谱密度之间满足x = Rx(0) =∫MAT LAB 中采用了式(1)～(4)的定义,因此在模拟y→Δ t2 a。1/Δ t ,再经过零Sx(ω) =∞- ∞Rx(τ )e- iτωdτ(14)σ2∞- ∞Sx(ω) dω(15)685北 京 航 空 航 天 大 学 学 报 2001 年如功率谱定义如式(14)的随机信号时,还必须乘以2 π ;如果采用的是单边功率谱定义,则乘以π 。

3 算 例3。 1 舰尾紊流中的自由大气紊流分量舰尾紊流的自由大气紊流分量是与飞机相对位置无关的分量 ,可白噪声进行滤波产生。 在 x ,y , z 3 个方向的分量分别用 u , v , w 表示。 其空间频谱形式如下式(16)所示[4]。5。 661 + (30。 48Ω)226。 59[1 + (121。 92Ω)2][1 + (304。 8Ω)2][1 + (40。 64Ω)2]21 + (30。 48Ω)2Su(Ω) =Sv(Ω) =Sw(Ω) =(16)式(16) 中的单位为(m/ s)2/ (rad/ m) ;Ω 为空间频率 ,单位 rad/ m。 且σi =∫谱 ,因此输出还必须乘以 π 。根据式(16)和式(6) ,可使用如图 3 的模块图进行模拟。∞0Si(Ω) dΩ ,是单边功率图 3 模拟自由大气紊流的 SIMULINK模块图图 3 中的 3 个增益均为 π 。 模拟条件为 :有限带宽白噪声的随机种子依次为 29713 , 78912 ,29321,噪声功率为 1 ,采样时间均为1 s ,仿真算法为 ode1 ( Euler) 法 , 步 长 为 1 s , 起 止 时 间 为0～32767 s。

模拟结果为 :输出信号的均值为 Eu=0。 0174 , Ev= - 0。 001 , Ew= 0。 0126;均方根为σu=0。 545 ,σv= 0。 520 ,σw= 0。 320。 而理论值为 :均值Eu= Ev= Ew= 0; 均方根σu= 0。 540 ,σv= 0。 515 ,σw=0。 323。为进一步考查输出是否满足要求 ,还应判断其功率谱密度是否一致。 由于在 MAT LAB 中使用PSD 函数进行谱估计时 ,功率谱密度也是按照式(2)来定义的 ,因此得到的结果还必须除以因子π ,然后才能和本例中的单边功率谱的理论值进行比较。 得到其功率谱密度后 ,与理论值进行比较 ,结果如图 4 ,实线为模拟结果 ,虚线为理论值。图 4 自由大气紊流功率谱密度模拟结果和理论值比较从模拟结果的统计特性比较和功率谱密度曲线的比较来看 ,模拟结果和理论值一致。3。 2 航空母舰的海上运动模拟航空母舰在海上的运动对于研究舰载飞机在航空母舰上的起飞着舰有十分重要的意义。航空母舰的随机运动是海浪随机运动的作用的结果 ,在研究舰载飞机起飞着舰时航空母舰的运动时 ,通常采用平稳随机过程的理论和方法 ,用白噪声输入一个线性的传递函数来模拟[5]。

文献[5]中给出了模拟 ESSEX级航母甲板运动的模型。 航空母舰的运动状态和模型如下: ES2SEX 级 , 航 速 14。 8 km/ h , 航 向 左 舷 60° , 浪 高5。 18m,海况主要为浪涌。 纵摇运动 (和垂荡运动h 对应的传递函数分别为0。 00583 s2[ξ1,ω1][ξ2,ω2]G θ =(17)Gh =0。 353568s( s + 0。 04)[ξ3,ω3][ξ4,ω4]式(17) 、 式(18) 中 , [ξ ,ω ] = s2+ 2ξωs + ω2;ξ1 =0。 2;ω1= 0。 55;ξ2= 0。 3;ω2= 0。 64;ξ3= 0。 2;ω3=0。 4;ξ4= 0。 2;ω4= 0。 55。 式 (17) 中 , G θ 的单位是rad ,其分子和文献[5]中的模型的分子之间相差一个因子π/ 180。因而 ,模拟使用的模块图如图 5 所示。

由于所采用的功率谱定义如式(1) 、 (4) ,模拟时不必再乘以比例因子 π 。(18)785第 5 期 彭 兢等 :在 MAT LAB 中模拟平稳随机过程图 5 模拟航空母舰海上运动的 SIMULINK模块图模拟条件 :有限带宽白噪声的随机种子为56789,噪声功率为 1 ,采样时间为 0。 1 s ,模拟方式为ode1(Euler)法定步长仿真 ,步长也为0。 1 s ,起止时间为0～1023 s ,输出采样时间为1 s。 模拟结果 :均值 E θ = 0。 002°, Eh= - 0。 001 m;均方根为σθ =1。08°,σh=1。79 m。 文献[5]中给出的理论值为:均值 E θ =0°, Eh=0m;均方根σθ =1。0°,σh=1。68m。图 6 为模拟结果的功率谱密度 Sθ 及 Sh和理a 纵摇频谱b 垂荡频谱图 6 航母海上运动功率谱密度模拟结果和理论值比较论值的比较。 由统计特性和图 6 中功率谱密度曲线的比较可以看出 ,模拟结果与理论值一致。

4 结 论以上 2 个具体实例说明:按照本文所述的方法 ,只要合理选择参数 ,并严格遵守功率谱密度的定义 ,用 MAT LAB 软件包中的 SIMULINK工具箱可以很好地模拟平稳随机过程 ,例如大气紊流和舰船运动等。 只要选取的随机数种子合适 ,模拟结果的数目足够多 ,模拟结果的统计特性功率谱密度曲线都与理论值一致。模拟中还要注意 2 点:一是选取信号源的采样时间和仿真步长时还要考虑传递的带宽 ,采样时间与系统的带宽匹配 ,模拟结果才合理;另一个是以上所描述的具体过程均为 MAT LAB 5。 1 中实现 ,而 MAT LAB 中的命令和工具箱随版本不同而变化 ,具体使用中需要仔细验证。参 考 文 献[1] 薛定宇。 控制系统的计算机辅助设计 — — — MAT LAB 语言及应用[M]。 北京:清华大学出版社 ,1996。[2] 汪荣鑫。 随机过程[M]。 西安:西安交通大学出版社 ,1987。 43～103。[3] 肖业伦 ,金长江。 大气扰动中的飞行原理[M]。 北京:国防工业出版社 ,1993。166～173。[4] 飞行力学杂志社。

军用规范 — — — 有人驾驶飞机的飞行品质(MIL2F28785C)的背景资料和使用指南[S]。 陕西富平:飞行力学杂志社 ,1985。[5] Durand T S, Teper GL。 An analysis of terminal flight path controlin carrier landing [R]。 AD 606040 , 1964。 72～77。Simulation of Stationary Random Process with MATLABPENGJing JIN Chang2jiang(Beijing University of Aeronautics and Astronautics , Dept。 of Flight Vehicle Design and Applied Mechanics)Abstract : The basic definitions of stationary random process and the principles of using spectrum decomposingto simulate the stationary random process were given。

The three important factors to be considered when simulatingstationary randomprocess with MAT LAB including selection of signal source , simulating algorithm and correspondingparameters , and regulation caused by different definition of the power spectral density。 The procedure of simulatingstationary random processes with MAT LAB was illustrated by two examples: air turbulence and aircraft carrier mo2tion at sea。 The tested simulation results indicate that stationary random process can be simulated very well withMAT LAB。Key words : stationary processes; power spectral density; turbulence ; ship motion885北 京 航 空 航 天 大 学 学 报 2001 年