蒙特卡罗模型

2024-10-05

蒙特卡罗模型（共7篇）

蒙特卡罗模型篇1

在服务行业、交通运输等各个行业业中,排队现象司空见惯,如:车站买票,高高速路口收费,超市收银台,银行ATM取取款,及计算机程序的处理等。这些有的属属于单服务台排队模型,有的属于多服务台台排队模型,不少文章对此进行了探讨,但但是一般均为“用随机数表进行模拟,从随随机数表中任一随机数开始顺序选择10个个随机数”,而不是用随机函数来完成随机机数的产生,完全由excel完成的蒙特卡洛模模拟更是少见,本文不仅用excel的随机函数数来完成模拟,而且会对模拟结果和公式算算出的理想结果之间进行误差比较。

一、单服务台有限队列模型

1、单服务台有限队列示意图

单服务台有限队列即只有一个服务台台的排队模型,且排队规则为损失制。所谓谓损失制,也就是当顾客到达时,若服务台台被占用且排队人数已满,则该顾客离去,,不排队等待。(图1)

2、单服务台有限队列模型基本概念

(1)等待时间Wq:从顾客进入队列时时起,至开始接受服务时为止的这段时间。。

(2)逗留时间W:一个顾客在系统中的的停留时间(包括接受服务时间)。

(3)系统中顾客的期望值L:系统中顾顾客个数的期望值(包括正在接受服务的顾顾客)。

(4)队列平均长度Lq:系统中排队等等候的顾客个数的期望值。

(5)空闲概率P0:服务台处于空闲的可可能性,这涉及到服务台的有效服务强度。。

3、蒙特卡洛模拟流程图(图2)

4、模型中各个变量的理论值

在排队模型中,到达的顾客数是一个以λ为参数的Poisson流,即顾客的相继到达的时间间隔服从λ以为参数的负指数分布,顾客的服务时间服从参数为的μ负指数分布。

二、蒙特卡洛模拟原理

蒙特卡罗(Monte Carlo)模拟,是以概率论和数理统计为指导的模拟方法。

实质是运用一连串的随机数来模拟可能出现的随机现象。模拟过程与极限状态方程的具体形式无关,且与变量的分布形式无关;收敛速度与随机变量的维数无关;模拟误差可以确定。

三、【M/M/1】:【N/∞/FCFS】模型的蒙特卡洛模拟

举例:某一个小诊所,其中患者到达服从λ=4的Poisson分布,诊断时间服从μ=6的负指数分别,且有4把椅子供患者排队等候。先用排队论公式求出系统内平均人数L、队列中平均患者数Lq、患者在系统中平均逗留时间W,在队列中平均等待时间Wq,然后用模拟方法计算以上指标,并与蒙特卡洛模拟结果比较。

运用蒙特卡罗模拟法模拟排队系统并计算排队系统的运行参数:

1、由上述条件,整理出在excel中患者到达时间间隔和医生服务时间的函数

设患者相继到达的时间间隔为T,则T的分布函数为,反解得:,因y需用随机函数LAND()代替,故在excel中时间间隔的函数式为:,同理,服务时间的函数式为:

2、在excel中模拟其它变量

(1)患者到达时间为前一患者到达时间加上二者的时间间隔;

(2)模拟服务开始时间,首先确定该患者是否离开,如果离开,不存在开始时间,如果入队等待,开始时间就是到达时间与其前面所有患者离开时刻的最大值。如图3:

(3)如果患者到达时队列未满而入队等待,患者离开时间等于患者的服务开始时间加上服务时间。如果到达时队列已满而离开,则患者到达时间即离开时间,此种情况对于我们的模拟没有意义,但是为了便于计算,我们把他的离开时间等于前面患者的离开时间。

(4)患者是否离开取决于其到达时刻的系统人数,如果队列已满则离去,否则入队等候。

(5)模拟计算患者到达时的系统人数首先把患者到达时刻与其前面的所有患者的最迟离开时间进行比较,如果前者大于后者,则此时系统人数为零,否则,计算出前面所有患者的离开时间大于该到达时间的数量减去因队满而离开的患者个数即是此时的系统人数,如图4:

(6)在未曾离去的情况下,患者的等待时间为服务开始时间减去到达时间。

(7)驻留时间等于离开时刻减去到达时间。

3、对模型进行模拟,计算排队系统的运行指标。并通过公式计算理论值进行比较

在假设患者人数为1000的情况下进行模拟,某次运行指标的数值对比如表1:

根据上表进行绘图,如下图5:

4、由于在模拟中使用了随机函数,按F9键可产生新的随机数,从而是系统的运行参数发生变化,理论值有模拟值的对比表和对比图也随之变化。

四、结语

在模型【M/M/1】:【N/∞/FCFS进行蒙特卡洛模拟的过程中,模拟的次数越多,得出的结果越准确,越接近于公式算出的理想值。尽管在精度要求高的情况下,需要非常大的模拟数量,会导致计算量较大、计算时间较长。但是随着计算机运算速度的加快,蒙特卡罗方法的运用也越来越广泛。通过上述模拟过程可见,蒙特卡罗模拟法在模拟排队过程以及求解排队系统的运行参数中,具有一定的优越性和实用性。

参考文献

[1]、蒋绍忠.《管理运筹学教程》.浙江大学出版社.2006(2).

[2]、张兰江,张建福.蒙特卡罗模拟法在排队论中的应用.交通与运输.2008.12.

[3]、张建航,李宗成,宋晓峰.单服务员排队模型及其蒙特卡洛模拟.计算机应用.2006.

蒙特卡罗模型篇2

随机潮流是当前电力系统稳态分析和规划决策的重要工具,与传统确定性潮流相比,它能够考虑各种不确定因素,从而更全面深刻地揭示系统的运行特性及薄弱环节[1,2]。随着分布式电源、电动汽车的大规模接入,储能技术的发展及电力市场化改革的稳步推进,电网将呈现出较大的随机性和高维特征。因此,加快建立适应新环境的电网分析和评估方法具有重要意义[3,4,5]。

随机潮流自20世纪70年代由Borkowska[6]提出以来,国内外学者做了大量研究并提出了多种计算方法,大致可归纳为三类:模拟法、近似法和解析法。以蒙特卡洛仿真(MCS)为代表的模拟法[7]通过简单随机采样(simple random sampling,SRS)得到输入变量的样本,再对每个采样点进行确定性潮流计算,最后统计获得各状态量的分布情况。整个过程原理简单且适用性广,在保证样本规模足够大的情况下能够获得很高的精度,不足之处是仿真次数多,耗时长,一般用于评估其他方法的优劣。近似法,如点估计法[8,9]以及解析法中的快速傅里叶变换[10]、半不变量法[11]等都是利用随机变量的数字特征并结合数值方法来求解状态量的分布,有效地减少了计算规模和时间,但当考虑输入变量相关性时,公式会非常复杂且高阶矩的误差较大。文献[12]创建了一种新的序列运算理论,简化了卷积运算,但理论架构的特殊性制约了其大规模地推广应用。文献[13]中的无迹变换可直接通过求解输出变量的均值和方差得到其概率分布,避开了复杂的非线性变换。然而由于该方法以高斯分布为变换基础,具有一定的局限性。因而近年来,一些学者又开始以MCS为基础,试图通过优化抽样或简化模型等手段来解决耗时的难题,比如文献[14]采用拉丁超立方采样技术,显著地提高了采样效率。文献[15]提出扩展准蒙特卡洛方法,弥补了拉丁超立方采样不能保证序列低差异性的缺陷。而线性MCS的引入可以避免反复的潮流迭代。MCS的计算精度除了受采样个数的限制,还与输入变量的模型密切相关。目前,大多数文献均假设各输入量服从正态分布,局限性较大。文献[16]用双峰分布和双Weibull分布分别拟合负荷及风速模型,适用于一些特殊场合,但缺乏通用性。文献[17]建立了负荷的K均值模型,用聚类的方法处理不同时段的负荷数据,模型隐含地计入了节点负荷间的相关性,但聚类过程及显著性水平检验需要大量时间。

本文提出一种基于混合高斯模型(GMM)的改进MCS随机潮流计算方法,对于有测量数据的输入变量,通过GMM直接建立其精确的概率模型,而在缺失大量数据的情况下,根据输入变量特点选取常用分布进行拟合。在MCS基础上引入均匀设计抽样(UDS)得到计及相关性的输入变量样本,并利用多重线性化潮流方程保存精度、简化计算。

1 输入变量的随机模型

负荷的多样性及分布式电源的接入使得输入变量(节点注入功率)的概率密度函数(probability density function,PDF)难以用单一的标准分布准确拟合(如图1所示)。尤其在中低压网络,节点负荷的聚集程度下降,采用正态模型会有较大偏差。

GMM理论上能够平滑任何类型的输入变量分布,利用测量数据和期望最大化(expectation maximization,EM)算法可以方便地确定参数。目前,GMM作为研究热点广泛应用于机器学习、图像分割、模式识别等领域,本文主要用于模型的建立。

1.1 GMM

GMM是由若干个高斯(正态)分布线性叠加而成,对于一维随机变量X,其PDF定义为[18]:

式中:为第i个成分的分布;ωi,μi和σi分别为混合高斯函数第i个成分的比例、均值和标准差。其中,比例参数ωi必须满足归一性条件,即。

GMM所含的高斯分布个数n直接影响模型的精确程度,n越大,模型的偏差越小,但待求参数和耗时增大。图1给出了用GMM(n=4)近似实际的负荷样本分布。

1.2 参数估计与EM算法

GMM建模的关键在于高斯成分个数n及每个成分所含参数ωi,μi和σi的确定。为简化分析,假设n为定值,一般参数估计问题可直接利用最大似然函数法,如GMM的最大对数似然函数为[19]:

式中:NL为数据点xi(i=1,2,…,NL)的个数。

要找到一组合适的参数ωk,μk和σk(k=1,2,…,n),使式(3)取最大。然而由于数据点xi所属成分类别未知,无法对式(3)直接求导解方程来求得最大值,因此这里引入EM算法[20],其过程主要分为以下两步。

步骤1(期望):估计数据由每个高斯成分生成的概率,即每个成分被选中后的期望。对于数据xi,它由第k个成分生成的概率为:

步骤2(最大化):假设式(4)求出的γ(i,k)就是正确的“数据xi由成分k生成的概率”,根据式(5)至式(7)用最大似然法得到对应的模型参数。

式中:。

重复迭代上面两步,直到似然函数式(3)的值收敛为止。为了保证取到全局最优解,可先用K均值聚类法[21]得到粗略结果作为参数初值,再用EM算法进行细致迭代。

高斯成分个数n的确定可应用赤池信息量准则(AIC)[22],它是衡量模型拟合优良性的一种标准,可表示为:

式中:p为参数个数;L为似然函数。ηAIC的值越小,表明所取的模型越好,它可以权衡模型复杂度与数据拟合准确性之间的矛盾。

综上,GMM所有参数求解的过程包括:①设高斯成分个数n取2或3;②根据随机变量X的测量数据,采用EM算法得到ωi,μi和σi(i=1,2,…,n)的值;③由式(8)求得GMM的AIC值η(n)AIC;④返回②,n=n+1;⑤重复以上过程直到|η(n+1)AIC-η(n)AIC|<ξ。对于无测量数据或测量数据存在缺失的随机变量模型的建立方法,可参见附录A。

2 具有相关性的输入变量样本的产生

2.1 UDS技术

UDS是由方开泰和王元教授于1978年提出的一种新的抽样方法[23],它主要基于空间填充的思想,考虑试验点在试验范围内均匀散布。相比SRS,其优势有:①覆盖同样大小的样本空间,UDS的采样规模更小;②UDS的稳健性更好。设有m个服从[0,1]区间均匀分布的随机变量Y1,Y2,…,Ym,N为样本数,则利用UDS产生样本矩阵Ym×N的步骤如下。

步骤1:设Λ为正整数集合的一个无穷子集,n∈Λ。生成一个正整数向量H1×m=[h1,h2,…,hm],其中h1=1,1<hj<n且对任意i≠j,hi≠hj。

步骤2:从多项分布中选取m个独立同分布样本,构成向量B1×m=[b1,b2,…,bm]。

步骤3:通过SRS生成样本矩阵Em×N=[e1,e2,…,ej,…,em]T(ej=[ej1,ej2,…,ejN]),其中,随机变量e服从[-0.5,0.5]上的均匀分布。

步骤4:样本矩阵Ym×N=[y1,y2,…,yi,…,ym]T(yi=[yi1,yi2,…,yiN])中的元素yij为:

式中:{·}表示保留小数部分;i=1,2,…,m;j=1,2,…,N。

以[0,1]上均匀分布的两个随机变量y1和y2为例,令N=5,分别用UDS和SRS产生y1和y2的样本各三次,样本点分布如图2所示。可以明显看出,UDS产生的样本点在空间散布的更均匀,当采样点数量相同时,UDS比SRS覆盖的样本空间更大,采样效率更高。

2.2 边缘变换和相关性处理

产生服从给定分布的输入变量样本是MCS随机潮流计算的基础,这里借助边缘变换:设输入变量w的累积分布函数(cumulative distribution function,CDF)记作F(w),可以证明[24],u=F(w)服从[0,1]区间上的均匀分布。对于变量y∶U(0,1),其CDF为v=(y),也是[0,1]上的均匀分布。因此,输入变量w可由式(10)的变换求得,即边缘变换。附录B图B1给出了边缘变换的图形解释,它能实现随机变量从均匀分布到任意分布的变换。

式中:Φ为标准正态分布的CDF。

输入变量w1和w2(有测量数据)之间的相关性可用相关系数表示[25]:

式中:分别为w1和w2的样本均值。

而对无测量数据的输入变量,其相关系数则作为输入参数直接给定。这样便形成一个n1×n1阶的相关系数矩阵Σ=(ρij),n1为输入变量个数。

为了使由边缘变换得到的输入变量样本满足式(11)的相关性,引入中间变量v′,v′服从高斯分布,其作用是将变量的相关性变换转移到高斯域内进行。

由UDS产生的样本矩阵Yn×N=[y1,y2,…,yi,…yn]T,其中yi~U(0,1)。先由式(12)得到对应的高斯变量vi:

由于随机变量之间的相关系数在边缘变换过程中会发生改变,需要对相关系数矩阵Σ做相应修正,对于从高斯分布到均匀分布的变换,Σ的非对角元素应修正为[26]:

修正后的相关系数矩阵记为Σ*=(ρ*ij)。通过对Σ*进行Cholesky分解可以得到具有相关性的高斯样本矩阵Vcorr[27]:

式中:V=[v1,v2,…,vi,…,vn]T(vi=[vi1,vi2,…,viN]);L为上三角阵,Σ*=LLT。

将Vcorr中的样本代入式(15)(式(12)的逆向过程),得到具有相关性的均匀分布样本:

最后再做一次边缘变换,即可产生具有相关性的输入变量wi的样本:

式中:为输入变量wi的CDF。

对于GMM拟合的输入变量,上述过程可以简化。设GMM第k个高斯成分的参数为(ωk,μk,σk),ωk可理解为第k个成分被选中的概率。在由式(14)得到具有相关性的标准高斯样本vicorr=[vi1corr,vi2corr,…,viNcorr]后,根据比例系数ω大小,随机选取每个样本点所对应的高斯成分,假设样本点vijcorr(j=1,2,…,N)对应第k个高斯成分,直接应用式(17)可得输入变量wi的样本wicorr=[wi1corr,wi2corr,…,wiNcorr],从而避免了反函数求解的困难。

实际样本的相关系数矩阵Σ′与Σ之间存在一定偏差,这是由于边缘变换是非线性的,变换过程中会产生截断误差。文献[28]给出了特定分布下系数的修正方法,通常情况下偏差都很小,可忽略不计。

3 多重线性蒙特卡洛仿真

由于潮流方程的非线性,在对每个输入变量的样本点做确定性潮流计算时,需要反复迭代,这占了MCS超过95%的耗时。文献[29]采用交流线性化模型,将潮流方程在基准运行点处进行泰勒展开并忽略2阶以上的高次项,可得:

式中:W为节点注入功率;X和Z分别为状态变量和支路潮流;X0,Z0,W0分别为基准运行点处的状态变量、支路潮流和节点注入功率;S0和T0为灵敏度矩阵,S0=J0-1,T0=G0J0-1,其中J0为雅可比矩阵,。

但当输入变量的变化范围较大时,采用单点线性化模型会引起较大的截断误差,为了在提高MCS计算速度的同时不至于损失太多精度,下面引入一种多重线性化的方法[30]。

设Ptot为系统总的有功功率,有

式中:NN为PQ节点个数;Pi为节点i的负荷有功功率;NG为PV节点个数;PGj为接入节点j的发电机输出的有功功率。

由于负荷功率和新能源出力均为随机变量,可知Ptot也是随机变量,其PDF根据式(20)得到。附录C图C1大致给出了Ptot的PDF曲线,将其等间距地划分为RT个区域,RT一般取6~8。在区域R1至RT内选取相应的基准运行点,并将潮流方程在各点处分别线性化,得到

式中:Si为第i个区域对应的灵敏度矩阵;Xi0为第i个区域的状态向量:f(·)为节点功率方程;为第i个区域的基准功率向量。

的求解步骤如下:①由第1节和第2节的建模和抽样方法产生具有相关性的输入随机变量样本S;②根据式(19),计算各样本的Ptot并判断其所对应的区域;③将属于同一区域的样本取平均值,即得到各个区域的基准功率向量。

综合以上分析,采用本文所提基于GMM的多重线性MCS方法求解电力系统随机潮流的流程如图3所示。

4 算例分析

依托IEEE 30小型节点系统对所提方法的性能进行评估,算例的结构及线路参数见文献[31]。在节点25,29处分别接入两个10 MW的小型风电场(无监测数据),假设风速均服从尺度参数为8.09、形状参数为2.17的Weibull分布,且风电机组采用恒功率因数控制方式,输出特性可表示为文献[32]中的分段函数。

节点14处还接有一容量为2 MW的太阳能光伏电站,其无功出力为0,输出有功近似满足附录A式(A2)的β分布,其中形状参数α=β取0.9,光伏系统在最强光照条件下的输出功率PTmax为1.5 MW。

表1给出了用GMM拟合的负荷模型参数,一般高斯成分个数n取2或3即可达到AIC的收敛要求。表中只列出了有功功率的参数,假定负荷无功功率与有功功率之间的功率因数角恒定,有

式中:μP和σP分别为有功功率高斯成分的均值和标准差;φ为功率因数角。

除表1中节点4,7,18,21外,其余节点的负荷功率均服从高斯分布,均值即为文献[31]的负荷数据,标准差取均值的10%。考虑节点负荷之间的相关性,可将系统大致分为两个区域,区域1包含节点16至20,相关系数ρ1=0.6,区域2(节点15、节点23至26)的相关系数ρ2=0.8,两个区域间存在弱相关性(ρ1,2=0.2)。

采用本文方法对上述算例进行仿真计算,N=1 000。为了比较负荷的GMM及高斯模型对随机潮流结果的影响,利用式(23)和式(24),将表1中的节点负荷模型替换为均值和方差与原GMM相等的高斯分布:

高斯模型和GMM两种模型下节点26电压幅值及支路19-20无功功率的PDF如图4所示。

由图4可以看出,虽然两者的波动范围基本一致,但图线之间存在一定偏差,卡方检验[33]可用于拟合优度的量化分析,卡方值越小说明函数的拟合效果越好,利用MATLAB中的chi2gof函数得到GMM的卡方值分别为128.63和103.06,小于相应的高斯分布的卡方值478.07和584.52。进一步求出越限指标,则GMM对应的电压、无功功率的越限概率分别为0.76%和1.42%,而高斯模型的结果为0.42%和2.62%,从而导致低估或高估系统运行的风险。而且建模产生的误差会代入后面的潮流计算中,在误差传递作用下放大(或缩小),具体情况视方法本身而定。

输入随机变量样本的选取会直接影响MCS的计算精度。本文2.1节介绍了UDS技术及采样特点,为了评估其有效性,假设以N=20 000的MCS所得随机潮流结果为准确值,用μa和σa分别表示输出变量的准确均值和均方差,μs和σs分别表示仿真得到的输出变量的均值和均方差,则输出变量均值和均方差的相对误差计算公式可表示为:

式中:下标r表示输出变量的类型(包括电压幅值、相角,支路有功功率、无功功率)。

图5比较了不同采样规模下UDS和SRS的平均电压相对误差曲线。图中:εμmean,U为电压均值的平均相对误差;εσmean,U为电压均方差的平均相对误差。可以看出,样本容量相同时,无论是精度还是误差稳定性,UDS均优于传统的SRS,并且采样规模越小,优势越明显。

另外,对潮流方程的简化处理也会不可避免地引入截断误差。以节点25和支路10-22为例,三种潮流模型(线性、非线性、多重线性)得到的电压幅值和支路有功的CDF如图6所示(图中“非线性”即为常规MCS)。

由图6可知,采用本文的多重线性潮流模型可以获得较高的精度,其图线与简化前非线性潮流方程的结果几乎完全吻合,而单点线性模型的曲线在两端附近的误差较大,这是由于截断误差与输入随机变量样本点和基准点的偏离程度成正相关。

引入方差和的根均值(ARMS)来衡量输出变量概率分布的计算精度[34]:

式中:ξr为ARMS指标;Cr,P,i和Cr,M,i分别为所提方法和常规MCS得到的输出变量CDF上第i个点的值。

表2列出了系统部分输出变量在不同潮流模型下的ARMS值。ARMS值越大,说明所提方法与常规MCS所得输出变量的概率分布偏差越大。可以看出,多重线性化方法所得ARMS值大多小于1%,相较于单点线性化的结果,精度上有了质的提高。

在主频为2.63GHz、运行内存为2GB的intel i3计算机上,比较所提方法和MCS的计算耗时如表3所示。

当样本规模较小(N=50)时,两种方法的耗时主要由样本的生成时间组成,由于所提方法在输入变量模型和采样技术上较MCS复杂,计算时间略微长些。而N=500,1 000时,样本规模逐渐成为MCS耗时的决定因素,随着N的增加,MCS的计算时间呈线性增长,而所提方法耗时的增加并不明显且远小于相同采样规模下的MCS。

IEEE 118节点系统测试结果详见附录D。

5 结语

针对当前随机潮流方法无法同时兼顾计算精度和耗时、建模手段过于单一等不足,本文提出一种计及输入变量相关性的基于GMM的改进MCS随机潮流方法,该方法具有以下特点:①对于已知测量数据的输入变量,不管其分布类型及复杂程度,都能用GMM精确拟合,通用性较好;②利用UDS技术提高采样效率,结合边缘变换和Cholesky分解的方法可以方便地产生任意关联的输入变量样本;③建立多重线性潮流模型,在大幅提升计算速度的同时,最大限度地减小因潮流方程线性化引起的截断误差。

算例测试结果说明,在相同精度条件下,本文方法的计算效率相较于传统MCS有了显著提高,尤其适用于实际大型输配电系统在考虑新能源接入、负荷波动等不确定因素下的运行分析、安全评估及规划调度,具有较好的工程应用价值。

PCM通信系统的蒙特卡罗仿真篇3

在现代通信系统中以PCM为代表的编码调制技术已被广泛应用于模拟信号的数字传输, 数字微波通信、光纤通信、卫星通信中均获得了极为广泛的应用。

PCM抗干扰能力强, 传输性能稳定, 远距离信号再生中继时噪声不累积。PCM通信系统能实现传输和交换一体化, 可实现数据传输与数据处理一体化, 适应信息化社会对通信的要求。

1 MATLAB中的蒙特卡罗仿真模型

1.1 蒙特卡罗方法

蒙特卡罗方法也称统计试验法, 它是一种利用统计抽样理论来近似的方法。它的基本思想是:首先建立与描述与该问题有相似性的一个概率模型或是随机过程, 使它的参数等于问题的解, 然后通过对模型或过程的观察或抽样实验来计算所求参数的统计特征, 最后给出所求解的近似值。并对解的精度作某些估计。

对于通信系统来说, 不确定性主要是由于传输过程中信道噪声的干扰, 发送信号在接收端有可能进行了错误的处理, 接收信号发生错误的变化, 而其作为性能的衡量指标就是误码率, 所以, 蒙特卡罗方法在通信系统仿真中的典型应用就是仿真出通信系统的误码率, 并通过误码率进行系统的性能分析。

概率论中的大数定律是蒙特卡罗仿真方法的主要理论基础, 随机变量的抽样分析是它最主要的手段。

1.2 二进制数字通信系统的蒙特卡罗 (Monte Carlo) 仿真模型

图1中, 首先产生一系列二进制数据序列, 产生此序列的方法是通过均匀随机发生器产生一组随机数, 范围都在 (0~1) 内, 然后判断随机数, 当落在 (0~0.5) 内, 图中二进制数据源就输出0, 否则输出1, 这样就产生出所要的二进制数据序列了。图中的高斯随机数发生器产生一组高斯分布的随机数, 用于模拟信道中的噪声。

因此, 图中检测器的输入序列就包括了噪声。然后检测器内部进行判断, 当检测器输入值小于0.5时, 检测器输出0, 否则输出1。然后, 将检测器输出信号与二进制数据源进行逐个比较, 比较结果由差错计数器统计输出, 这样就得到差错数和差错数与采样值N的比值。

2 PCM通信系统

2.1 PCM基本原理

脉冲编码调制 (PCM, Pulse Code Modulation) 简称脉码调制, 是用二进制代码取代连续模拟信号的抽样值从而进行信息传输的编码技术, 即将时间和幅值都连续的模拟信号转变成时间和幅值都离散的二进制信号。从模拟信号变成数字信号经过三个步骤, 即抽样, 量化, 编码。

对模拟信号进行等时间间隔抽样, 变成时间上离散的PAM信号, 此时PAM仍为模拟信号, 再将PAM信号无限的抽样值量化成有限个可能值, 使之成为时间和取值上都离散的多进制数字信号, 最后对量化后的多进制信号进行二进制编码, 最终成为PCM信号。

为了改善小信号的量化性能, 本文采用压扩非均匀量化。ITU有两种建议方式, 它们分别是A律和μ律方式, 我国和欧洲国家采用A律方式, 但是A律压缩实现较复杂, 通常采用13折线法编码来近似A律。

2.2 PCM通信系统的仿真模型

本文PCM通信系统的仿真采用蒙特卡罗模型, 仿真框图如图2所示:

3 PCM通信系统的仿真

3.1 模块的仿真图及结果分析

首先得到需要进行抽样的模拟信号如图3所示。

模拟信号a=sin!t"

再对此模拟信号进行抽样, 抽样间隔为2π/32, 则可知在一个信号周期内, 对此模拟信号进行了32次抽样, 抽样信号图如图4所示。

得到抽样值之后应先进行量化再进行编码然后成为二进制数据流, 本文运用MATLAB语言进行仿真时将量化和编码结合在一起。取抽样信号的前十个抽样值:

抽样信号a1=sin t!1"

将该十个抽样值进行13折线量化编码后的八位二进制数分别为 (假设2048为总的量化单位) :

进行上述编码后, 将编码的二进制数据流分成八个一组, 通过译码后恢复出原来的抽样值。用MATLAB语言进行仿真的译码仿真结果如图5所示。

由图5可知, 经译码之后恢复的抽样信号图与原模拟信号的抽样信号图的形状大致一样, 只是幅值上存在一些差值, 这和理论上得到的结论是一致的。即在进行PCM调制后存在量化误差 (称为量化噪声) , 而且这种量化误差是无法避免的, 只能尽量减小。

3.2 PCM系统性能分析

由于传输过程中有量化噪声, 在传输过程中有系统噪声, 在进行系统分析时, 将散布在整个PCM系统中的噪声等效在信道中, 噪声类型为高斯随机噪声。故在用MATLAB仿真分析系统性能时, 用均值为0的加性高斯白噪声近似噪声。PCM系统的性能指标主要有信噪比和误码率。故采用蒙特卡罗仿真模型进行仿真。

用MATLAB编写的不同信号强度的量化信噪比的仿真结果如图6所示:

不同幅度的模拟信号的量化信噪比仿真图

由图6可知, 小信号进行非均匀量化时, 量化噪声功率的均方根值与信号强度成比例关系, 即信号强度越大, 信号的量化信噪比也越大, 且信号强度在一定范围内的量化信噪比基本保持不变。进行非均匀量化时量化噪声对大、小信号的影响差不多, 这就改善了小信号的量化信噪比。

用MATLAB语言编写的信号经过系统传输后的误码率与量化信噪比关系的仿真结果如图7所示:

从图中两条曲线可看出, 理论值和实际仿真的结果基本一致。从图可知, 当信噪比越大时, 系统的误码率越小。所以在保证一定的信噪比条件下, 系统的误码率可以达到所要求的指标。

摘要：PCM技术是一种数字调制技术, 它将模拟信号转换成数字信号, 广泛用于电话系统, 数字微波通信系统, 光纤通信系统, 卫星通信系统中。PCM通信系统最主要的部分就是模拟信号的数字化, 它包括抽样, 量化, 编码。本文构建PCM通信系统基本模型, 运用MATLAB语言进行系统仿真, 画出仿真波形并在此基础上进行性能分析。

关键词：数字通信系统,PCM,MATLAB语言,仿真

参考文献

[1]樊昌信, 曹丽娜.通信原理[M].第六版.国防工业出版社, 2006.

[2]吴伟陵, 续大我, 庞沁华.通信原理[M].第三版.北京邮电大学出版社, 2005.

[3]高西全, 丁玉美.数字信号处理[M].第三版.西安电子科技大学出版社, 2010.

[4]John G.proakis等著, 刘树棠译.现代通信系统 (Matlab版) [M].第二版.电子工业出版社, 2006.

[5]张珽.基于蒙特卡罗方法的通信系统误码率的仿真[J].无线通信技术, 2010.

[6]冯微玮.基于MATLAB的PCM仿真[J].科技致富向导, 2012.

[7]邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析[M].清华大学出版社, 2008.

[8]Nyquist H.Certain Topics in telegraph transmission theory trans[J].AIEE, 1928.

[9]Feldman CB, bennelt WR.Band width and transmission performance BSTJ[J].1949.

蒙特卡罗模型篇4

光线追踪算法用于计算机生成真实感图像, 在电影, 电视和辅助设计领域有着广泛的应用。光线追踪算法和另一大图像生成算法光栅化相比, 其图像生成的质量和真实性均占优, 但是速度要比光栅化算法慢很多。自从该算法被提出以来, 研究人员提出了多种对此算法的改进, 比如光子映射算法和辐射度算法。这些算法和传统光线追踪算法比起来在生成图像质量相同的情况下速度快了很多。但是这些算法只针对特定的场景布置才能起到加速效果, 在一些复杂场景特别是有大量光滑和光泽表面的情况下会出现大量误差, 并且计算速度会大幅下降。

光线追踪算法的本质是求解光传输方程。通过光能守恒定律我们可以把积分空间离散的传输方程转换为场景中面积度量的积分。求解这个积分的方法是把它转换为一个无穷维的高阶积分, 并用蒙特卡罗方法采样整个积分空间, 使积分结果按概率趋近于方程解。由于渲染场景的复杂性和光照模型的多样性, 求解这个问题的最大难点在于如何在有限的时间内使蒙特卡罗方法逼近真实解。目前, 光线追踪主要还是应用于离线渲染领域。如何在保证图像质量的情况下提高算法速度就成了亟待解决的问题。

二、问题描述

2.1光传输方程。

真实感渲染中, 核心算法之一是模拟光从光源发出, 在场景物体和介质之间进行交互, 传播, 最终被吸收或者进入摄像机的过程。从物理学的角度来看, 这个过程就是电磁波 (辐射) 与实物粒子的碰撞与反射。我们将一个表面上一点p处ω0方向上辐射出射能量定义为L0 (p, ω0) , 那么根据热辐射定律, 出射能量可以表示为入射能量在表面BRDF函数和方向ω的积分, 该积分可以用下列公式表示:

其中Le (p, w0) 是表面发出的辐射;f (p, ω0, wi) 是表面的BRDF函数, 表征物体表面对辐射的反射和吸收特性。积分空间s2是场景中的所有表面。由于积分中含有即立体角, 为了能计算这个积分, 我们希望它能统一到积分空间中。有限空间的场景中, 沿着射入一个表面的光线能找到另一个场景中的点。于是我们定义函数, 它表示点p处沿着方向找到的另一个场景中的点。如果这样的点不存在, 则t函数返回一个特殊的值。将函数t带入式 (1) 可以得到

(2) 式与 (1) 式相比, 右侧的积分函数和左侧统一, 但是依然无法直接进行计算, 因为它依然含有ω。考虑到我们在上一步的变换同样适用于L (p, ω0) , 为此我们定义

其中P'为沿着点p的出射光线遇到的场景中的点。由于场景中存在遮挡关系, 我们定义几何函数

函数表示点p和P'的可见性。如果两点之间没有遮挡, 那么V的值是1, 否则是0.将上述函数引入式 (2) 中, 我们得到了光传输方程的表面形式:

注意到上式中左边和右边都有函数L。现在我们考虑一条从光源到摄像机的完整光路径, 我们把P0处的式 (3) 带入P1处的式 (3) 得到

如果光路径上还有其他点, 也可以同样按照上述方法带入, 从而把光传输方程转化成路径形式。有了式 (4) 就可以在某些简单场景中计算数值解了。但是对于大多数复杂场景, 光路径的长度可能很长甚至是无限长, 式 (4) 就变成了一个无穷维积分, 这样的积分是无法计算数值解的。对于这样的场景, 我们用2.2中的采样算法进行计算。

2.2 Metropolis-Hastings采样算法。

其中表示以概率a接受从X状态到X’状态的变化, 否则不变。由 (5) 式得

以这种采样的方式计算积分的话, 就是MCMC算法:

其中g (x) 可以是任意的函数, xi服从分布.

式 (4) 可以简化写为

使用MCMC算法, 式 (6) 可以用下式计算:

三、算法流程与实现

本文算法的实现使用了开源项目mitusba render作为输入和输出, 核心算法流程图如下:

四、实验与结果分析

实验使用光线追踪测试场景veach, 生成的图像如下:

通过生成的图像我们可以看到, 使用MCMC算法生成的图像整体具有较好的真实感, 并且能够表现出光线多次反射对场景中物体的影响。但是, 该算法的速度比较慢, 在2450M处理器上渲染上图使用了15分钟26秒。如何提高算法的性能将是一部研究的重点。

参考文献

[1]Veach.Robust Monte Carlo methods for light transport simulation[D].斯坦福:斯坦福大学, 1997.

[2]Jakob W., S.Marschner.Manifold exploration:A Markov Chain Monte Carlo technique for rendering scenes with difficult specular transport[J].ACM Trans, Graph, 2012, 31, 4, 58:1-58:13.

[3]Jaakko L., K.Tero, L.Samuli, A.Miika, D.Fr'edo, A.Timo.Gradient-Domain Metropolis light transport[J].ACM Trans.Graph, 2013, 32, 4, 95:1-95:11.

[4]Jakob W.2012 Mitsuba v0.4 EB/OL].http://mitsuba-renderer.org.

[5]Jaakko L., K.Tero, L.Samuli, A.Miika, D.Fr'edo, A.Timo.GradientDomain Metropolis light transport[J].ACM Trans.Graph, 2013, 32, 4, 95:1-95:11.

[6]Jakob, W.Light transport on path-space manifolds[D].伊萨卡:康奈尔大学.2013.

蒙特卡罗模型篇5

质量稳定性是水泥的重要质量指标之一。许多用户对水泥质量稳定性的关注程度，甚至高于对质量指标水平的关注。水泥质量稳定性以均匀性试验得到的变异系数表示，GB12573—90《水泥取样方法》和JC/T578—1995《评定水泥强度匀质性试验方法》规定，均匀性试验采用系统抽样方法，样本容量为10。生产实践表明，该抽样方法的抽样误差较大，但是难以通过生产试验的方法确定抽样误差的数值。文献[1,2]应用蒙特卡罗方法对抽样误差进行数值模拟取得满意结果。文献[3]分析了抽样方差估计值的精度与样本数目的关系，指出抽样方差估计值的标准偏差与样本容量的平方根之积近似为一常数。

本文应用蒙特卡罗方法模拟了现行均匀性试验方法的抽样误差、扩大样本容量后的抽样误差。根据试验结果提出了对现行均匀性试验抽样方法的修改建议。

1 一般原理与模拟方法

1.1 一般原理

蒙特卡罗模拟的一般原理和方法已有介绍[4]。

在一个总数为N的检验批按系统抽样方法在其中抽取n个子样的方法是：

1)按下式计算抽样间隔k:

其中INT是取整函数。

2)确定抽样起点i:

其中RAN为随机数发生函数。式(2)随机产生均匀分布于(1,k)的自然数。

3)以i作为起点，依次在i,i+k,i+2k,...,i+(n-1)k的位置上抽取子样共n个。

1.2 模拟方法

设每个检验批(即一个编号的出厂水泥)为1 000t，按系统抽样方法在其中抽取n袋水泥。则总体个数N=20 000，样本容量(即每个出厂水泥编号的抽样数量)为n，抽样间隔k=20 000/n。设总体强度标准偏差(即一个编号内各袋之间的标准偏差，均匀性试验结果是总体标准偏差的估计值)为σ，总体强度平均值为μ。

产生N个服从标准正态分布的随机数R，即R～N(0,1)。均值为μ、方差为σ2的正态分布随机变量X可通过下列变换得到：

即X～N(μ,σ2)。

产生均匀分布于(1,20 000/n)的随机数i，取Xi,Xi+k,Xi+2k，…，Xi+(n-1)k，计算子样的标准偏差S，即完成一次模拟。连续模拟30次，计算30次模拟标准偏差估计值S的标准偏差。

2 模拟结果与讨论

按μ=50MPa(该数值大小与模拟结果无关)，σ=0.5MPa、1.0MPa、1.5MPa、2.0MPa,n=10次、20次、30次、40次、50次进行全部组合模拟，以观察不同标准偏差和不同样本容量对抽样误差的影响。σ=1.0MPa,n=10次的模拟结果见表1。

不同总体标准偏差σ、样本容量n的模拟结果的统计值见表2。

取置信概率为0.95，根据表2结果计算的不同总体标准偏差σ、样本容量n下，以标准偏差表示的抽样误差见表3。

MPa

在0.95的置信概率下，以标准偏差表示的抽样误差与样本容量的关系见图1。

表3和图1表明，当样本容量为10，总体标准偏差分别为2.0MPa、1.5MPa、1.0MPa、0.5MPa时，标准偏差表示的抽样误差分别为0.813MPa、0.676MPa、0.437MPa、0.220MPa，相对误差分别为40.7%、45.1%、43.7%、39.7%。如此大的抽样误差是无法接受的。表明现行均匀性试验方法的抽样数量偏少。

样本容量增加，抽样误差随之减小。当总体标准偏差为2.0MPa时，在样本容量10～50区间,随着样本容量的增加，抽样误差减小的速率基本保持恒定。当总体标准偏差为1.5MPa、1.0MPa、0.5MPa时，在样本容量10～50区间，随着样本容量的增加，抽样误差开始时下降较快，在样本容量为30以后抽样误差减小的速率变慢。这意味着样本容量小于30时，增加样本容量可以比较明显减小抽样误差;样本容量大于30时，增加样本容量对减小抽样误差的效果不明显。

总体标准偏差减小，抽样误差随样本容量增加而减小的速率减小。这意味着在总体标准偏差较大时，增加样本容量对减小抽样误差具有更明显的效果。

总体标准偏差增加，抽样误差显著增加。图2更加清楚地显示了这种关系。

图2表明，随着总体标准偏差的增加，抽样误差近乎呈线性增加。样本容量增加，抽样误差随着总体标准偏差增加的速率减小。这意味着当样本容量较小时，随着总体标准偏差的增加，抽样误差以更快的速度增加;当样本容量较大时，随着总体标准偏差的增加，抽样误差以较慢的速度增加。对于实际操作的意义是，在总体标准偏差较大时，增加样本容量对减小抽样误差的效果更加明显。

有关文件规定以变异系数表示均匀性试验结果。现行标准GB12573—90规定均匀性试验的样本容量为10，此时不同28d抗压强度下，以变异系数Cv表示的抽样误差与总体标准偏差的关系见图3。

当总体标准偏差为1.5MPa———这是一般水泥厂均匀性试验结果的最大值，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为1.7%、1.4%、1.1%。相对于《水泥企业质量管理规程》规定均匀性试验的变异系数不大于3.0%的要求，这个误差显然是不能接受的。再次证明现行均匀性试验方法的抽样数量偏少。

在总体标准偏差为1.5MPa和1.0MPa时，以变异系数Cv表示的抽样误差与样本容量的关系分别见图4和图5。

图4显示，当总体标准偏差为1.5MPa，样本容量为30，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.9%、0.7%、0.6%。

图5显示，当总体标准偏差为1.0MPa，样本容量为20，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.8%、0.6%、0.5%。

3 均匀性试验的适宜样本容量

总体标准偏差为1.5MPa，样本容量为30时，以变异系数Cv表示的抽样误差接近1%，这个误差依然有些偏大，但考虑到抽样成本和检验成本，以及在样本容量大于30以后继续增大对减小抽样误差的效果不显著，不宜继续增加样本容量。总体标准偏差为1.0MPa，样本容量为20时，以变异系数Cv表示的抽样误差也接近1%。

一般通用水泥28d抗压强度为40～60MPa，平均值约50MPa。当总体标准偏差为1.0MPa、1.5MPa时，变异系数分别为2%、3%。变异系数为2%对应着目前水泥厂的一般水平，变异系数为3%对应着目前水泥厂的较差水平。

建议对均匀性试验的抽样数量做如下规定：在连续生产的情况下，前次均匀性试验的变异系数不大于2%时，每次试验抽样数量为20个;当非连续生产或者前次均匀性试验的变异系数大于2%时，每次试验抽样数量为30个。

4 结论

1)当总体标准偏差为1.5MPa，样本容量为10，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为1.7%、1.4%、1.1%。证明现行均匀性试验抽样方法存在很大误差。

2)当总体标准偏差为1.5MPa，样本容量为30，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.9%、0.7%、0.6%。当总体标准偏差为1.0MPa样本容量为20，对应水泥28d抗压强度分别为40MPa、50MPa、60MPa时，以变异系数Cv表示的抽样误差分别为0.8%、0.6%、0.5%。

3)根据试验结果对均匀性试验的抽样数量作出了建议。

参考文献

[1]高志,何锡文,李一峻.分层性物质的组合取样误差与份样数目之间的关系及其Monte Carlo模拟[J].分析科学学报,1999,15(5):353-357.

[2]高志,何锡文,李一峻,等.组合取样的误差理论及其Monte Carlo模拟[J].高等学校化学学报,1999,20(12):1853-1857.

[3]高志,李一峻,何锡文,等.取样方差估计值的精度与样本数目之间的关系[J].分析化学,2001,29(2):171-174.

[4]张大康.基于蒙特卡罗方法的率值稳定性定量分析(Ⅰ)——随机检验误差对生料率值稳定性的影响[J].水泥,2007,(8):1-6.

蒙特卡罗模型篇6

关键词：三维公差综合,拟蒙特卡罗,旋量参数,广义逆矩阵,Halton序列

0 引言

公差设计不仅有助于降低产品生产成本,还影响产品质量以及装配成功率[1]。公差设计所涉及的因素较多,如何合理分配公差即公差综合是一个复杂的问题。目前有关公差综合的研究主要针对线性或二维尺寸公差的优化设计,对于计算机辅助三维公差综合的研究还不多[2,3]。

在三维CAD系统中,确定零件配合特征与装配功能要求的关系并进行公差综合方法的研究,对于辅助设计人员进行公差设计、提高产品质量有着重要意义。近些年来,一些研究人员在三维公差表示的基础上,把蒙特卡罗方法用于三维公差综合。文献[4]针对产品设计阶段预定义的装配功能要求,利用Jacobian矩阵建立与装配功能要求相关的3D公差链,并以此为基础采用蒙特卡罗法进行三维公差综合仿真研究。王太勇等[5]探讨了采用蒙特卡罗法进行三维尺寸及公差设计问题的实现方法。在特征自由度三维公差表示模型的基础上,许本胜等[6]采用蒙特卡罗法进行了三维公差综合仿真研究。然而,在进行三维公差综合仿真时,以上研究侧重于阐述蒙特卡罗法进行计算机仿真的流程,在具体的应用中简单地采用伪随机数进行仿真计算,对于蒙特卡罗方法中随机数的产生机制考虑不足,导致在大样本重复运算时会出现结果不稳定、收敛速度慢的问题,不能很好地发挥蒙特卡罗方法的应用效果。

拟蒙特卡罗方法是一种采用拟随机序列的蒙特卡罗方法,相比采用伪随机数的蒙特卡罗方法,拟蒙特卡罗方法不但具有更快的收敛速度,在仿真分析中所获得的结果也更为稳定和可靠。本文针对计算机辅助三维公差综合中仿真计算收敛速度慢以及结果不稳定的问题,在获得三维公差综合仿真等式的基础上,研究基于拟蒙特卡罗法的三维公差综合仿真方法,并通过实例进行具体的分析和阐述。

1 拟蒙特卡罗方法

拟蒙特卡罗方法的基本思想为:首先建立与问题相关的概率模型或随机过程,然后通过对概率模型或随机过程的观察或抽样试验来计算所求随机参数的统计特征,最后给出所求解的近似值。与蒙特卡罗方法不同的是,在对概率模型或随机过程进行随机抽样的过程中,拟蒙特卡罗方法采用的是拟随机序列。相较伪随机序列,拟随机序列的点列分布更加均匀,在仿真模拟时可以加快收敛速度并提高计算结果的稳定性和可靠性。

图1、图2所示分别为二维伪随机序列和二维拟随机序列。

图1、图2中的二维随机序列横纵坐标分别为区间[0,1]上的随机数。图2中的拟随机序列为二维Halton序列,与图1中的伪随机序列相比,它在函数空间的分布更加均匀。随机序列对应的点列在函数域上分布越均匀,其偏差相应越小,仿真计算时的收敛速度就越快,计算准确度也就越高。拟蒙特卡罗方法的计算误差可用Koksma-Hlawka不等式来确定[7]:

$\begin{array}{l} | \frac{1}{Ν} \sum_{k = 1}^{Ν} f (x_{k}) - \int_{Ι^{d}} f (u) d u | \leq \\ V (f) D_{Ν}^{*} (x_{1}, x_{2}, \dots, x_{Ν}) \end{array}$

其中,Id表示d维空间[0,1),随机点列x1,x2,…,xN∈Id;d、N分别为序列元素的维数和个数;V(f)是函数f在空间中Hardy-Krause意义下的有界变差函数;D*N(x1,x2,…,xN)定义如下:

D*N(x1,x2,…,xN)=DN(B;

$\begin{array}{l} Ρ) = \\ \sup_{B \in Ι^{d}} | \frac{A (B; Ρ)}{Ν} - λ_{d} (B) | \end{array}$

A(B; $Ρ) = \sum_{n = 1}^{Ν} χ_{B} (x_{n})$

其中,P是空间Id中的包含点集x1,x2,…,xN;对于任意属于空间Id的子集B,λd是d维勒贝格测度;χB是B的特征函数,A(B;P)是计算满足xn∈B(1≤n≤N)的点的个数的函数。

综上可得前N阶拟蒙特卡罗序列的误差收敛阶为O(N-1(lgN)d),而蒙特卡罗序列的误差收敛阶为O(N-1/2),因而拟随机序列的收敛速度要高于伪随机序列的收敛速度[7]。

2 基于拟蒙特卡罗的三维公差综合仿真

2.1 旋量参数和旋量矩阵

传统的尺寸与极限配合的公差表示方法已逐渐向新的产品几何技术规范(geometrical product specifications,GPS)过渡,出现了一些新的公差表示方法,如特征自由度[8]、旋量参数等[9,10]。这些公差表示方法本质上有着共同之处,即在三维CAD系统中,零件可看作是由点、线、面等基本几何特征构成的,将公差看作是构成零件几何特征微观上的三维几何变动量。为了规范和方便对公差问题的分析和描述,ISO/TC213延伸了特征的概念,定义了名义特征、名义派生特征、实际特征等特征类型,如图3所示。

从图3可看出,测量得到的实际特征为拟合特征或拟合派生特征,由于存在加工误差,其与名义特征或名义派生特征之间不可避免存在偏差,因而,从测量角度看,名义特征或名义派生特征在允许偏差范围(即公差)内是变动的,变动量可通过旋量参数来表示。对于选定的参考坐标系OiXYZ,旋量参数包括:名义特征或名义派生特征以原点Oi为中心的沿X、Y和Z轴平动方向的变动ui、vi、wi和绕X、Y和Z轴转动方向的变动αi、βi和γi。三维公差为名义特征或名义派生特征在三维空间允许的变动范围,相应地可通过旋量参数所构成的旋量矩阵来表示,记为Ti=[Diθi],其中Di=[uiviwi],θi=[αiβiγi]。

需说明的是,如果名义特征或名义派生特征沿某个方向的变动不产生新的特征,则对应的旋量参数记为零。如图3中,圆柱轴线名义派生特征沿轴线方向平动和绕轴线转动不产生新的特征,因而当选定轴线为参考坐标系OiXYZ的Z轴时,其三维公差对应的旋量矩阵为Ti=[Diθi],其中Di=[uivi 0],θi=[αiβi 0]。

2.2 三维公差综合仿真等式

产品是由若干零件装配而成的,对于规定的装配功能要求,与之相关的各零件配合面形成闭合的装配特征链,如图4所示。

图4中,产品的装配功能要求可视为装配特征链中首末配合特征之间的相对变动量,记为T0=[D0θ0],其中D0=[u0v0w0],θ0=[α0β0γ0]。装配特征链中第i个配合特征的变动量Ti(i=1,2,…,n,n为装配特征链配合特征个数,且n>1)对装配功能要求的累积影响量记为T0_i=[D0_iθ0_i],D0_i=[u0_iv0_iw0_i],θ0_i=[α0_iβ0_iγ0_i]。T0、Ti间的关系可以利用机器人学中刚体的坐标转换方法建立,为了减少坐标转换的计算工作量,将参考坐标系与装配基准坐标系各坐标轴取为同向。原点Oi在装配基准坐标系中的位置记为(Xi, Yi, Zi),利用机器人学中刚体坐标变换的知识,有

TT0_i=J0_iTTi (1)

$J_{0_i} = [\begin{matrix} 1 & 0 & 0 & 0 & Ζ_{i} & - Y_{i} \\ 0 & 1 & 0 & - Ζ_{i} & 0 & X_{i} \\ 0 & 0 & 1 & Y_{i} & - X_{i} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}] (2)$

从而:

$Τ_{0}^{Τ} = \sum_{i = 1}^{n} Τ_{0_i}^{Τ} = J [Τ_{1} Τ_{2} \dots Τ_{n}]^{Τ} (3)$

系数矩阵J称为雅可比矩阵,有

J=[J0_1J0_2 … J0_i … J0_n] (4)

由于装配特征链中配合特征数目n>1,因而J为6×6n阶长方阵。式(3)中,将Ti(i=1,2,…,n)视为未知量,对于给定的T0,式(3)有解但解不唯一,为获得唯一解,取解的最小范数作为约束条件,其现实意义为:各零件配合特征理想公差视为零,获得解与理想值之间具有范数意义下的最小偏差。记x=[T1T2 … Tn]T,y=TT0,由广义逆矩阵的知识,存在唯一的矩阵J+,满足:

$\min_{J x = y} ∥ x ∥ = ∥ J^{+} y ∥$

矩阵J+称为长方阵J的广义逆矩阵,J+可通过下式进行计算[11]:

J+=JT(JJT)-1 (5)

综合式(1)～式(5)可得三维公差综合模型如下:

[T1T2 … Tn]T=J+TT0 (6)

式(6)即为三维公差综合仿真等式。

2.3 三维公差综合仿真

拟蒙特卡罗方法的关键在于拟随机序列的生成,目前已经提出了若干拟随机序列,如Faure序列、Sobol序列以及Halton序列等。本文采用Halton序列进行仿真计算。记s维Halton序列为X=(x1,x2,…,xs),通过选定某个十进制整数m和以下步骤来生成s维Halton序列。

(1)选择s个基(base)b1,b2,…,bs,这s个基可以使用前s个素数(2,3,5,7,…)。

(2)对每个基bi,将m转换成bi进制:

$m = \sum_{k = 1}^{n_{i}} a_{k} b_{i}^{k - 1} i = 1, 2, \dots, s (7)$

其中,ni为将m转换成bi进制后的位数,ak为将m转换成bi进制后第k位上的数字。转换后得到一个随机向量(m(b1),m(b2),…,m(bi),…,m(bs)),其中m(bi)为m对应的bi进制整数。

(3)将m(bi)按位反序排列(如120反序排列为021),然后在前面加小数点得到新的值并将其转换成十进制数,记为xi,有

$x_{i} = \sum_{k = 1}^{n_{i}} a_{k} b_{i}^{- k} i = 1, 2, \dots, s (8)$

式中,ak为小数点后的第k位数字。

(4)置m←m+1,重复步骤(2)、步骤(3)直至得到所需的随机向量个数。

由上述步骤生成Halton序列,进行适当数学变换可得到服从其他分布的随机序列,如正态分布转换方法为:将s维Halton序列中各随机变量xi相加,依据中心极限定理,当s充分大时,有

$\frac{1}{\sqrt{n / 12}} (\sum_{i = 1}^{s} x_{i} - \frac{n}{2}) ~ Ν (0, 1) (9)$

通过式(9)获取标准正态分布随机变量,且当s=12时可达到较好的精度,从而要得到服从正态分布的随机变量x～N(μ,σ2),作如下变换即可:

$x = (\sum_{i = 1}^{12} x_{i} - 6) σ + μ (10)$

综合上述三维公差综合模型并利用Halton序列生成特定分布随机数的方法,基于拟蒙特卡罗的模拟仿真思路为:首先,根据特定分布函数,通过Halton序列获取装配功能要求T0各旋量参数的随机数;然后,由式(6)计算得到装配特征链中各配合特征旋量参数的值,重复计算便得到各配合特征旋量参数的样本;最后,对所得样本进行统计分析,得到样本均值、样本方差等统计量。图5所示为基于拟蒙特卡罗的三维公差综合仿真流程。

3 实例分析

图6a所示齿轮泵齿轮的精度等级为5级,齿轮分度齿圈径向跳动公差为Fr=0.016mm,综合考虑补偿热变形、保证正常润滑等因素,要求啮合轮齿侧向间隙Fc最小值为0.03mm,最大值为0.1mm。

3.1 装配功能要求与配合特征旋量参数

轮齿啮合的接触形式为线接触(图6b),装配功能要求(functional requirement of assembly,FRA)即侧向间隙可看作是两轮齿名义啮合线A1、C1空间位置发生了相对变动。装配链中首末配合特征旋量矩阵T0=[0 v0w0α0β0γ0],各旋量参数与Fr、Fc的关系为

${\begin{cases} v_{0} = F_{c} 0.03 m m \leq v_{0} \leq 0.1 m m \\ w_{0} = F_{r} / 2 - 0.016 / 2 m m \leq w_{0} \leq 0.016 / 2 m m \\ α_{0} = F_{c} / 2 r 0.03 / r \leq α_{0} \leq 0.1 / r \\ β_{0} = F_{r} / 2 b - 0.008 / b \leq β_{0} \leq 0.008 / b \\ γ_{0} = F_{r} / b - 0.03 / b \leq γ_{0} \leq 0.1 / b \end{cases} (11)$

式中,r为齿轮分度圆半径;b为齿轮宽度。

图7为与装配功能要求相关的各零件配合特征及三维公差对应的旋量参数。

3.2 三维公差综合仿真

图7中,{Oi}(i=1,2,…,6)为装配特征链中各配合特征参考坐标系,{Oi}与装配基准坐标系{O0}各坐标轴方向相同,各原点坐标为:O1=O6=(0,0,0),O2=(-b,0,r),O3=(-b,0,r),O4=(-b,0,-r),O5=(-b,0,-r)(b=19mm,r=12.5mm)。利用式(2)、式(4)得到系数矩阵J为

$J = [\begin{matrix} 0 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 0 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 & 1 & 0 & 0 & 0 & r & 0 \\ 0 & 1 & 0 & - r & 0 & 0 & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & - r & 0 & 0 & 0 & 1 & 0 & 0 & 0 & - b & 0 & 1 & 0 & 0 & 0 & - b \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & b & 0 & 0 & 0 & 1 & 0 & b & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}]$

由式(6),利用MATLAB软件计算得到齿轮泵各配合特征三维公差域内旋量参数与侧向间隙各变动量之间的关系。以分度圆特征C1[0 v6w6α6β6γ6]为例,通过计算得到以下等式:

$\begin{array}{l} v_{6} = \frac{r γ_{0} + v_{0}}{2 (11 + r^{2})} \\ w_{6} = \frac{b β_{0} + w_{0}}{2 (11 + b^{2})} \end{array}$

$\begin{array}{l} α_{6} = \frac{α_{0}}{2} \\ β_{6} = \frac{(11 + 2 b^{2}) β_{0} - 2 b w_{0}}{4 (11 + b^{2})} \\ γ_{6} = \frac{(11 + 2 r^{2}) γ_{0} + 2 r v_{0}}{4 (11 + r^{2})} \end{array}} (12)$

在进行三维公差综合仿真时,首先需对装配功能要求T0进行模拟采样。在批量生产的情况下,尤其当装配特征链中配合特征数目n较大时,可认为T0服从正态分布。这里假定Fc、Fr服从正态分布。对于Fc～N(μ,σ2),根据设计要求,取侧向间隙均值μ=(0.03mm+0.1mm)/2=0.065mm;方差σ用于描述实际结果与设计目标相偏离程度,可由6σ设计原则来确定[12],取σ=0.0117mm。同样可得到Fr～N(0,0.00272)。由式(10),利用12维Halton序列分别获取Fc、Fr的随机数,进而通过式(11)即可得到T0中各旋量参数v0、w0、α0、β0及γ0的随机数。

以旋量参数γ6为例,分别采用蒙特卡罗和拟蒙特卡罗方法,利用获得的旋量参数v0、w0、α0、β0及γ0的样本通过式(12)进行仿真计算,在不同容量样本和相同容量样本情况下进行两次实验,得到的仿真统计结果见表1和表2。

表1中,γ6(p)为采用蒙特卡罗方法计算的样本均值,γ6(q)为采用拟蒙特卡罗方法计算的样本均值。由表1可看出,在不同容量样本情况下进行多次仿真计算,采用拟蒙特卡罗方法得到的标准差S较小,说明采用拟蒙特卡罗方法获得的结果波动较小。另外,由表1也可看出,当增大样本容量时拟蒙特卡罗方法的收敛速度更快,尤其当样本容量较大时,计算结果更为稳定。

由表2可看出,在相同容量样本情况下进行多次仿真计算,采用拟蒙特卡罗方法得到的结果的波动较小,表明采用拟蒙特卡罗方法得到的结果更为稳定。

由上述分析可知,拟蒙特卡罗方法更为高效,获得的结果也更稳定和可靠,因而能更好地运用于三维公差综合模拟仿真。取样本容量M=5000,通过拟蒙特卡罗法进行公差综合仿真计算,以分度圆特征C1为例,其结果见表3。

表3中, $\bar{X}$ 一栏为几何特征C1各旋量参数的三维公差分配结果,S为样本标准差,Xmax、Xmin分别为最大值、最小值。各旋量参数最大值、最小值之差即为C1在参考坐标系{O1}中各方向的变动范围,差值越小,表明对相应方向上C1的几何变动量要求越严格;标准差对应样本分布集中程度,其值越小,样本分布越集中,相应几何变动量越难以加工保证,在设计阶段应予以重视。

4 结语

本文在利用旋量参数和旋量矩阵建立三维公差综合仿真等式的基础上,详细讨论了用拟蒙特卡罗法进行三维公差综合仿真的方法。通过将所阐述的方法应用于齿轮泵装配体的三维公差综合,得到泵体、齿轮轴及轮齿配合特征旋量参数统计量。由于三维公差设计问题的复杂性,如何将仿真计算得到的配合特征各旋量转化成具体的公差大小还有待进一步的研究。通过对实例结果的分析可以看出,相比蒙特卡罗方法,采用拟蒙特卡罗方法进行三维公差综合仿真时的收敛速度更快,获得的结果更为稳定和可靠,可以简单高效地应用于三维公差综合的仿真计算,从而验证了本文所阐述方法的有效性和实用性。

参考文献

[1]Hong Y S,Chang T C.A Comprehensive Review ofTolerancing Research[J].International Journal ofProduction Research,2002,40(11):2425-2459.

[2]王恒,宁汝新.计算机辅助公差设计与分析的研究现状及展望[J].航空制造技术,2006(3):73-75,102.Wang Heng,Ning Ruxin.Present Status and Prospect ofComputer Aided Tolerance Design and Analysis[J].Aero-nautical Manufacturing Technology,2006(3):73-75,102.

[3]孙晓燕,黄克正,张勇.计算机辅助三维公差设计的研究现状[J].工具技术,2007(41):15-19.Sun Xiaoyan,Huang Kezheng,Zhang Yong.PresentState and Review of Computer Aided 3D ToleranceDesign[J].Tool Engineering,2007(41):15-19.

[4]LaperrieáL,Kabore T.Monte Carlo Simulation of Toler-ance Synthesis Equations[J].International Journal ofProduction Research,2001,39(11):2395-2406.

[5]王太勇,熊越东,路世忠,等.蒙特卡罗仿真法在尺寸及公差设计中的应用[J].农业机械学报,2005,36(5):101-104.Wang Taiyong,Xiong Yuedong,Lu Shizhong,et al.Application of Monte Carlo Method to Dimensionand Tolerance Design[J].Transactions of the Chi-nese Society for Agricultural Machinery,2005,36(5):101-104.

[6]许本胜,王灿,景晖,等.基于特征自由度变动的公差综合建模与仿真[J].中国机械工程,2011,22(24):2981-2985.Xu Bensheng,Wang Can,Jing Hui,et al.ToleranceSynthesis Modeling and Simulation Based on DOF ofGeometric Variations of Features[J].China Me-chanical Engineering,2011,22(24):2981-2985.

[7]雷桂媛.关于蒙特卡洛及拟蒙特卡洛方法的若干研究[D].杭州:浙江大学,2003.

[8]黄美发,钟艳如.CAD系统中并行公差的建模方法[J].中国机械工程,2004,15(18):1623-1625.Huang Meifa,Zhong Yanru.A New Model for Con-current Tolerancing in CAD System[J].China Me-chanical Engineering,2004,15(18):1623-1625.

[9]胡洁,熊光楞,吴昭同.基于变动几何约束网络的公差设计研究[J].机械工程学报,2003,39(5):20-26.Hu Jie,Xiong Guangleng.Study on Tolerance De-sign Based on Variational Geometric ConstrainsNetwork[J].Chinese Journal of Mechanical Engi-neering,2003,39(5):20-26.

[10]张为民,陈灿,李鹏忠,等.基于雅可比旋量法的实际工况公差建模[J].计算机集成制造系统,2011,17(1):77-83.Zhang Weimin,Chen Can,Li Pengzhong,et al.Toler-ance Modeling in Actual Working Condition Based onJacobian-Torsor Theory[J].Computer IntegratedManufacturing Systems,2011,17(1):77-83.

[11]姚俊,张玉春.工程矩阵方法[M].2版.北京:国防工业出版社,2007.

蒙特卡罗模型篇7

随着国民经济的增长,道路交通事业也得到了迅速发展,机动车保有量不断增加,但同时城市交通拥堵,交通事故频发,交通环境恶化等现象也逐渐凸显出来。为了解决地面交通迅速发展所引发的各种问题,提出了智能交通系统概念,并日益受到来日各方的重视。其中车辆跟踪作为智能交通系统的重要组成部分,在交通监控、停车场调度、事故检测、自动导航等方面表现出了广阔的应用前景和研究价值[1]。近年来,针对车辆跟踪问题,研究人员进行了积极的探索,提出了一些相对有效的解决方案。比如文献[2]利用LDA对目标和背景进行区分,通过模型匹配方式实现车辆跟踪。由于建模和匹配的复杂性、运行环境的快速变化及模型更新速度的限制使得跟踪的效果和实时性不够理想。文献[3]中,利用前后景差分方法实现简单背景下的车辆目标跟踪,虽然方法简单,但由于每一帧都需要进行全局运算,且在已知环境下进行,因此在跟踪速度和方法的扩展方面性能较差。文献[4]利用卡尔曼滤波结合目标的区域特征,通过局部的线性预测,实现车辆简单运动轨迹下的位置估计。尽管该方法克服了全局运算的缺陷,但由于卡尔曼滤波的线性估计特性,却也相应制约了该方法在非线性运动状态下的适应能力。文献[5]中,利用均值移位(Mean Shift)算法有效克服了目标的非线性运动跟踪,但在目标跟踪时,该算法对目标位置的定位仅考虑了均值漂移向量提供的方向,而没有考虑目标实际的宏观运动。当目标的颜色分布和背景相似或有所干扰时,Mean Shift算法可能出现跟踪错误,最终丢失目标。为更好地解决非线性、非高斯的目标跟踪问题,近年来发展起来的序列蒙特卡罗(Sequential Monte Carlo)方法成为了研究热点,该方法避免了数学或模型近似的缺陷,利用加权样本实现对目标状态的后验估计,因而将其用于复杂的非线性、非高斯的目标跟踪中是有效的、可行的。

1 序列蒙特卡罗滤波技术

1.1 基本原理[6]

序列蒙特卡罗滤波(SMCF)是通过非参数化的蒙特卡罗模拟来实现贝叶斯估计,用样本而不是以函数形式对概率密度进行描述。其基本思想是利用一系列随机样本的加权和表示所需的后验概率密度,得到状态的最优估计值,其中样本即为目标状态描述的一种可能。当样本点数增至无穷大时,蒙特卡罗模拟特性与后验概率密度函数相等价,从而使滤波精度逼近最优估计。其在数学范畴内的理解:令{Xi0:k,wki}为所求后验概率密度p(X0:k|y1:k)的一个随机观测量,其中,{Xi0:k,i=0,...,Ns}为一支撑点集。与点集相应的权值为{wki,i=0,...,Ns}。Xi0:k={xj,j=0,...,Ns}表示0至k时刻的状态集,归一化权值Σiwki=1。于是k时刻的后验概率密度可以近似地表示为:

该结果可由强大数定律理论来保证近似的精确度。

这种方法的突出优点是不再受线性和高斯分布的限制,原则上适用于可用状态空间模型表示的任何非线性系统。另外,由于所采集的样本具有独立同分布特性,因此该方法也相应具备了并行处理的能力。与其它形式的递推贝叶斯滤波器相比较[7],还具有方法灵活、易于执行和近似精度高等特点。

1.2 SMCF技术执行构架

根据上述对SMCF技术基本原理的论述和分析,可将其通用的滤波步骤归纳如下:

(1)初始采样:初始时刻(k=0),根据先验分布p(x0),建立初始状态样本集{XKi}i=0,...,N~p(xk)。

(2)重要性采样:在时刻k(k>0),遵从期望分布q(xk|xjk-1,zk)采集样本XKi,使XKi~q(xk|xjk-1,zk)。

(3)计算权值并归一化:利用公式(2)~(4)进行各样本的权值计算,并进行归一化处理。

(4)重采样:根据公式求出Neff,如果Neff≥Nth,转到步骤2;否则执行重采样策略,得到新的样本集。

(5)令k=k+1,返回步骤2。

序列蒙特卡罗滤波过程如图1所示。

2 基于SMCF的车辆跟踪算法

视频车辆跟踪过程等价于在连续的视频帧间创建基于位置、形状、纹理、色彩等有关特征的对应匹配问题。本文中采用基于SMCF的色彩匹配车辆跟踪方法,该方法首先利用运动预测和SMC技术对下一帧目标可能出现的区域进行位置和尺度的抽样;然后计算各分布样本与参考目标的相似匹配程度;最后通过样本加权估计目标的后验状态来获得车辆目标的位置。

2.1 目标模型

颜色特征是目标跟踪的基本特征,其对目标的部分遮挡、旋转及尺度变化等情况都具有较好的稳定性。考虑到HSV颜色空间比RGB颜色空间更符合人的视频特性,当下正有效且广泛地应用于图像处理和计算机视觉领域中,因此本文采用在HSV空间下的颜色特征统计作为目标状态描述,颜色量化等级为Bin=8×8×8。为增强跟踪过程的鲁棒性,文中还融合了目标颜色分布的空间信息,利用HSV空间直方图构建目标模型qs,进行样本与参考目标的相似度量。目标模型qs的数学描述如下:

其中,l表示目标区域中心;lm表示第m个像素的位置;表示目标区域的大小;M表示目标区域像素点总数;k(·)表示核函数(一般选择高斯核函数);δ(·)为Kronecker Delta函数;C为归一化常数;u为特征分量。

2.2 系统动态模型

系统动态模型表述目标的状态随时间的转移过程。在SMC滤波中,目标的状态可表示为样本,系统动态模型描述各样本在帧间的传播过程。

鉴于车辆目标的刚性特点,本文采用包含运动车辆的最小方框图来描述状态变量X。X可表示为:

其中,x和y分别表示包含车辆目标的最小方框图的中心位置;h和w为该矩形方框高和宽的一半;vx和vy描述的是中心位置的运动。因此SMC滤波在某一时刻t的加权样本集可表示为,N表示样本总数;Xt(i)表示第i个样本;wt(i)表示第i个样本相应的权值。

系统动态模型表示为:

其中,Wt-1是一个高斯随机噪声矩阵;A是一个状态转移矩阵,可将其表示为:

其中,△t表示视频序列采样间隔。

2.3 系统观测模型

观测模型用来测量物体状态向量和当前观测值的符合程度。观测值表示每个样本所代表的目标的可能状态和真实状态之间的相似程度,接近目标真实状态的样本被赋予较大的权值;反之,则赋予较小的权值。

本文利用采样样本与期望目标之间的色彩分布和空间距离定义观测模型为:

其中,qy和qx分别表示样本和期望的色彩特征分布;d为衡量两者色彩分布相似度的距离;η为特征重要性参数。

基于SMCF车辆跟踪算法流程如图2所示。

3 仿真结果与实验分析

为测试本文算法的实际跟踪效果,对几组赛车视频在MATLAB7.0环境下进行了跟踪实验。仿真中计算机的配置为PIV2.4G,1G内存,采样数为80,视频帧率为25fps,图像格式为24位真彩色。

实验1:利用本文提出的车辆跟踪算法对赛道上快速行驶的车辆进行跟踪,如图3所示。算法采用帧间预测实现了目标区域的快速定位,有效地减小了跟踪处理的运算量,提高了实时跟踪效果。

实验2:分别用本文算法和均值移位算法对赛道上短时遮挡的车辆进行跟踪,如图4和图5所示。本文算法在对遮挡车辆跟踪时,因采用了运动预测和加权采样的方式,最大限度地减少了失跟和误跟现象的发生。而均值移位算法因只能通过局部极值点进行跟踪,当车辆快速行驶时,跟踪窗口不能及时更新,且在发生遮挡时跟踪窗口仍把参考目标为最优匹配的区域当作跟踪目标区域,因此在遮挡情况下跟踪窗口脱离目标,导致跟踪失败。

4 结束语

本文提出了一种基于帧间预测和特征匹配的序列蒙特卡罗滤波车辆跟踪算法,该方法通过运动区域预测减小了运算开销,并且采用样本匹配加权后验估计实现了对非线性车辆的位置及轨迹的快速跟踪。实验证明在SMCF的框架下,算法能够对无规则快速运行的车辆进行实时有效的跟踪,同时对背景复杂及部分遮挡或短时目标丢失等情况均具有较好的适应性和鲁棒性。

参考文献

[1]赵飞,张金荣.智能交通系统在城市道路安全中的应用研究[J].交通标准化,2008(1):45-48.[1]赵飞,张金荣.智能交通系统在城市道路安全中的应用研究[J].交通标准化,2008(1):45-48.

[2]NGUYEN H,SMEULDERS A.Robust tracking using foregrou-nd-background texture discrimination[J].Computer and Vision,2007,69:277-293.[2]NGUYEN H,SMEULDERS A.Robust tracking using foregrou-nd-background texture discrimination[J].Computer and Vision,2007,69:277-293.

[3]MAGEE D.Tracking multiple vehicles using foreground,back-ground and motion models[J].Image and Vision Computing,2-004,22:143-155.[3]MAGEE D.Tracking multiple vehicles using foreground,back-ground and motion models[J].Image and Vision Computing,2-004,22:143-155.

[4]SALARPOUR A,FATHI M,DEZFOULIAN M H.Vehicle tra-cking using kalman filter and features[J].Signal and Image P-rocessing,2011,2:1-8.[4]SALARPOUR A,FATHI M,DEZFOULIAN M H.Vehicle tra-cking using kalman filter and features[J].Signal and Image P-rocessing,2011,2:1-8.

[5]XU D,WANG Y,AN J.Applying a new spatial color histog-ram in mean-shift based tracking algorithm[C]//New ZealandConference on Image and Vision Computing,2005.[5]XU D,WANG Y,AN J.Applying a new spatial color histog-ram in mean-shift based tracking algorithm[C]//New ZealandConference on Image and Vision Computing,2005.

[6]杜云明,颜兵兵.加权采样的序列蒙特卡罗滤波技术[J].火力与指挥控制,2012,37(3):46-50.[6]杜云明,颜兵兵.加权采样的序列蒙特卡罗滤波技术[J].火力与指挥控制,2012,37(3):46-50.

【蒙特卡罗模型】推荐阅读：

蒙特卡罗模拟08-16

蒙特卡洛06-09

皮埃蒙特07-15

《卡罗纳》教学设计07-14

卡罗琳的祝福09-30