西电雷达系统仿真2015-3

1、1 西安电子科技大学西安电子科技大学 史林史林 通信地址：西安市太白南路通信地址：西安市太白南路2号西安电子科技大学号西安电子科技大学140信箱信箱 邮编：邮编：710071 电话：电话13572155669 Email： n为什么要做系统仿真？为什么要做系统仿真？ 科学研究需要实验，进行分析、优化、评估、验证。科学研究需要实验，进行分析、优化、评估、验证。n什么是系统仿真？什么是系统仿真？ 简单地说，系统仿真就是用模型代替真实系统做实验。简单地说，系统仿真就是用模型代替真实系统做实验。n系统仿真的一般步骤系统仿真的一般步骤 n系统仿真的分类系统仿真的分类 物理仿真

2、：构造物理模型做实验，直观、形象，费用大、难修改。物理仿真：构造物理模型做实验，直观、形象，费用大、难修改。 数学仿真：构造数学模型，用计算机做实验。经济、方便、灵活。数学仿真：构造数学模型，用计算机做实验。经济、方便、灵活。 半实物仿真：数学与物理模型相结合做试验。半实物仿真：数学与物理模型相结合做试验。n数学仿真（计算机仿真）的步骤数学仿真（计算机仿真）的步骤 真实系统真实系统建立系统模型建立系统模型用模型做实验用模型做实验结果分析结果分析 建立数学模型建立数学模型设计算法设计算法编写程序编写程序运行运行结果分析结果分析 31.1离散系统数学模型及其仿真方法离散系统数学模型及其仿真方法 1

3、.1.1差分方程差分方程 描述描述n阶非时变线性离散系统输入阶非时变线性离散系统输入u(k)与输出与输出y(k)关系的关系的差分方程为：差分方程为： 给定输入序列给定输入序列u(0)、u(1)、u(k)，系统初始条件系统初始条件y(0)、y(1)、y(k-1)可计算出系统的输出相应：可计算出系统的输出相应： niniiiikubikya000)()(niniiiikubikyaky10)()()(10a1.1.2脉冲相应脉冲相应h(k) 系统的零状态相应为：系统的零状态相应为： 1.1.3系统函数系统函数H(z) 系统零状态相应的系统零状态相应的z变换为：变换为： 系统零状态相应的系统零状态相

4、应的z变换为：变换为：1.1.4状态方程状态方程 四种模型之间可相互转换四种模型之间可相互转换 离散系统的输出计算是代数运算或矩阵运算，最简单。离散系统的输出计算是代数运算或矩阵运算，最简单。niiiniiizazbzUzYkhZTzH00)()()()()()()(zUzHzY0)()()()()(iikuihkukhky)1()1()(kkkBuAxx)()(kkCxy1.2连续系统数学模型及其仿真方法连续系统数学模型及其仿真方法1.2.1连续系统数学模型连续系统数学模型n微分方程微分方程n冲击相应冲击相应n传递函数传递函数n状态方程状态方程 四种模型之间可相互转换四种模型之间可相互转换

5、连续系统的输出需求解连续系统的输出需求解n阶常系数线性微分方程，或一阶常系数阶常系数线性微分方程，或一阶常系数线性微分方程组，计算复杂，计算量大。线性微分方程组，计算复杂，计算量大。110)()(niinininiininidttudbdttydadtuhtuthty)()()()()()()()(sUsHsY)()()(tttBuAxx)()(ttCxy1.2.2连续系统数值积分法仿真连续系统数值积分法仿真 所谓数值积分仿真法就是利用数值积分计算方法求解微所谓数值积分仿真法就是利用数值积分计算方法求解微分方程和微分方程组分方程和微分方程组n一阶微分方程的数值解法一阶微分方程的数值解法 设一阶

6、微分方程为设一阶微分方程为 （1） 两边取定积分有：两边取定积分有：或或 （2） 可见，一阶微分方程的求解问题转换为定积分的计算问题。可见，一阶微分方程的求解问题转换为定积分的计算问题。0)()(,)(yaybtatytfty，dttytfdttykkkktttt11)(,)(1)(,)()(1kkttkkdttytftyty 欧拉法欧拉法 按步长按步长h取等距点取等距点 ，用矩形公式近似定积分，用矩形公式近似定积分 代入（代入（2）式有）式有 欧拉公式欧拉公式简记为简记为：误差分析：误差分析：设设 ，则称，则称 为局部截断误差。为局部截断误差。将将 在在t=tk点展开点展开Tayler级数得

7、：级数得：则局部截断误差为：则局部截断误差为： 可见，可见，h小，则局部截断误差小，但离散点数量增加，计算量将增大小，则局部截断误差小，但离散点数量增加，计算量将增大 。, 2 , 1 , 0,0kkhttkhtytftttytfdttytfkkkkttkkkk)(,)()(,)(,11)(,)()(1kkkktythftyty, 2 , 1 , 0),(1kythfyykkkk)(kktyy 11)(kkyty)()(1htytykk )(! 3)(! 2)()()()(321kkkkkktyhtyhtyhtyhyyty)()()(),()()(1111kkkkkkkkktyhtytyyth

8、fytyyty)()(! 3)(!2232hOtyhtyhkk 龙格龙格-库塔法库塔法其中：其中：局部截断误差：局部截断误差：误差估计事后估计：误差估计事后估计： 或或步长选择：步长选择：)22(643211kkkkhyykk),(1kkytfkhkyhtfkkk2,212hkyhtfkkk2,223hkyhtfkkk34,)()(511hOytykky(tk) 步长步长h计算一步计算一步 yk+1(h)，局部截断误差：局部截断误差：y(tk+1)- -yk+1(h)=ch5y(tk) 步长步长h/2计算二步计算二步 yk+1(h/2)，局部截断误差：局部截断误差：y(tk+1)- -yk+1

9、(h/2)=2c(h/2)516122)()(5511)2(11chhcytyytyhkkhkk)(1)2(1)2(11151)(hkhkhkkyyyty)(1)2(1hkhkyy)(:1)2(1khktyyYes2:hhNon一阶微分方程组的数值解法一阶微分方程组的数值解法 微分方程组：微分方程组： （3）引入向量表示：引入向量表示：则（则（3）可式表示为：）可式表示为： （4） （4）式与一阶微分方程（）式与一阶微分方程（1）式在形式上完全一样，故一阶微分方程的）式在形式上完全一样，故一阶微分方程的数值解法完全适用于一阶微分方程组的解法，只是将变量用向量替代即可数值解法完全适用于一阶微分方

10、程组的解法，只是将变量用向量替代即可 例如，欧拉公式例如，欧拉公式miytytytytytftyiimii, 2 , 1)()(,),(),(,)(0021T21)()()()(tytytytmyT020100myyyyT21mffff0yyyfy)(),()(0ttt, 2 , 1 , 0),(1kythkkkkfyyn高阶微分方程的数值解法高阶微分方程的数值解法微分方程：微分方程： （5）引入新变元：引入新变元：则则(5)式变为一阶微分方程组：式变为一阶微分方程组： 利用一阶微分方程组的解法可求解高阶微分方程。利用一阶微分方程组的解法可求解高阶微分方程。miytytytytytytftyi

11、imm, 2 , 1)()(,),(),(),(,)(00) 1() 1()( mitytyii, 2 , 1)()() 1(miytytytytytftytytytytytytyiimmmm, 2 , 1)()(,),(),(,)()()()()()()(0021213221n连续系统输出响应的计算连续系统输出响应的计算 微分方程模型：利用数值计算方法求解高阶微分方程微分方程模型：利用数值计算方法求解高阶微分方程 状态方程模型：利用数值计算方法求解一阶微分方程组状态方程模型：利用数值计算方法求解一阶微分方程组 传递函数模型：可转换为高阶微分方程，或利用拉布拉斯变换。传递函数模型：可转换为高阶

12、微分方程，或利用拉布拉斯变换。 冲击响应模型：可转换为传递函数模型，或利用数值积分。冲击响应模型：可转换为传递函数模型，或利用数值积分。n数值积分仿真法的特点数值积分仿真法的特点 优点：计算误差容易控制优点：计算误差容易控制 缺点：算法复杂，计算量大。缺点：算法复杂，计算量大。1.2.3连续系统离散相似法仿真连续系统离散相似法仿真 离散系统输出响应的计算是代数运算或矩阵运算，算法简单，离散系统输出响应的计算是代数运算或矩阵运算，算法简单，计算量小。因此，我们很自然会想到利用离散系统的输出响应来计算量小。因此，我们很自然会想到利用离散系统的输出响应来近似连续系统的输出响应，思路是将连续系统模型离

13、散化。近似连续系统的输出响应，思路是将连续系统模型离散化。n连续系统状态方程模型的离散化连续系统状态方程模型的离散化 状态方程改写为：状态方程改写为：两边左乘两边左乘 ，得，得由于由于代入状态方程有代入状态方程有)()(ttCxy)()()(tttBuAxx)()()(tttBuAxxtAe)()()(tttttBueAxxeAA)()()()()(ddttttttttttAxxeAxexexeAAAA)()(ddtttttBuexeAA两边积分两边积分于是有状态方程的解：于是有状态方程的解： （1）对对（1）式分别令式分别令t=kT、t=(k+1)T进行采样有：进行采样有： （2） （3）

14、有：有： (4)设设u( (t) )在采样间隔内为常数，即在采样间隔内为常数，即 则则(4)(4)式为式为 右端的积分与右端的积分与k k无关，令无关，令k=0k=0有：有：d)()(dd00ttdBuexeAAd)()0()(0ttttBueexexAAAd)()0()(0)(kTkTkTkTBuexexAAdTkTkTkTk)1(0)1()1()()0() 1(BuexexAA) 2 () 3 (TAed)()() 1()1()1(TkkTTkTkTTkBuexexAATkkTtkkTutu) 1( , 2 , 1 , 0),()(TkkTTkTkTkTTk)1()1()(d)() 1(u

15、BexexAATTTkTkTTk0)()(d)() 1(uBexexAA令令 ，则有，则有写成离散形式为：写成离散形式为：对输出方程采样有：对输出方程采样有：或或于是连续系统状态方程离散化后的离散状态方程为：于是连续系统状态方程离散化后的离散状态方程为：其中其中TTtdTT0)()(,)(BeeAmA)()()()() 1(kTTkTTTkuxxm)()() 1(kkkuxxm)()(kTkTCxy)()(kkCxy)()()()() 1(kkkkkCxyuxxmTtTsT1)(ILT)(AIeATTTdTdT00)()()(BBeAm例：连续系统状态方程例：连续系统状态方程 求其系统的离散化

16、状态方程求其系统的离散化状态方程解：解：)()()(tttBuAxx)()(ttCxy1001100CBANTtTsT1)(ILT)(AIeA11011000011001001ssssssAI111110111) 1(101)(1sssssssssAITtTsT1)(ILT)(AIeATTTttteeee101101TTTdTdT00)()()(BBeAmdeNNdNeeTtTTT0)(0)()()1 (0101)1(00TTTTeTNNTNeNNTNeNeNN2.1均匀分布随机数的产生均匀分布随机数的产生 统计试验中对随机数的基本要求是：统计试验中对随机数的基本要求是：n周期长周期长n统计特

17、性好统计特性好n产生速度快产生速度快 随机数产生的基本方法有：随机数产生的基本方法有：n物理法：物理法：物理产生，采样存储。如晶体管热噪声，放射粒子计数器。物理产生，采样存储。如晶体管热噪声，放射粒子计数器。 优点：真正的随机数。优点：真正的随机数。 缺点：不可重复缺点：不可重复n数学法：数学法：采用数学方法利用计算机产生。采用数学方法利用计算机产生。 优点：可重复优点：可重复 缺点：不是真正的随机数，称为伪随机数。缺点：不是真正的随机数，称为伪随机数。2.1.1平方取中法平方取中法 方法：初值方法：初值平方平方舍弃收尾各一位舍弃收尾各一位下一个数的初值下一个数的初值 公式：公式：xn+1=

18、=Int 2-kxn2 ( (mod 22k) ) yn=xn/ 22k ，yn是（是（0,1）区间均匀分布随机数。）区间均匀分布随机数。 k是字长是字长例：十进制数例：十进制数 x1=81 x12=6561 x2=56 x22=3136 x3=13 x32=0169 x4=16 x42=0256 二进制数二进制数 x1=(13)10= (1101)2 x12=(169)10=(10101001)2 x2=(10)10=(1010)2 x22=(100)10=(01100100)2 x3=(9)10=(1001)2 n缺点：若初值选取的不当，迭代后可能会成为常数，目前很少用。缺点：若初值选取的

19、不当，迭代后可能会成为常数，目前很少用。例：例： x1=(1011)2 ，迭代后为常数。，迭代后为常数。n2.1.2乘同余法乘同余法 方法：方法：xn+1=xn (mod m) yn=xn/m ，yn是（是（0,1）区间均匀分布随机数。）区间均匀分布随机数。 参数：参数：m一般取计算机的字长，如一般取计算机的字长，如16,32,6416,32,64 x0 0一般取一般取2 2b+1+1，b是正整数。是正整数。 一般取一般取2 23 3a3 3或或5 5，a是正整数。是正整数。例：乘同余，例：乘同余， =19， m=100， x0 0= =11n xnxnxn+1=xn (mod m)xn+1/

20、 m1 11 209 9 0.92 9 171 71 0.713 71 1349 49 0.492.1.3混合同余法混合同余法 方法：方法：xn+1=xn+C (mod m) yn=xn/m ，yn是（是（0,1）区间均匀分布随机数。）区间均匀分布随机数。 参数：参数： C一般取一般取4 4a+1+1，a是正整数。是正整数。 一般取奇数。一般取奇数。n从经验上看，混合同余法产生的随机数的质量比同余法好。从经验上看，混合同余法产生的随机数的质量比同余法好。2.1.4移位寄存器法移位寄存器法 周期最大为周期最大为2n, ,n是位数。是位数。例例2.1： 100001001010 0101 0010

21、 0001 100001001010 0101 0010 0001 周期周期=6=6 1 0 0 0+模模2 2加加例例2.2： 1000 1100 1110 1111 0111 1011 0101 1010 1101 0110 0011 0011 1001 0100 0010 0001 周期周期= =16 1 0 0 0+模模2 2加加2.2随机数的统计检验随机数的统计检验n什么是统计检验？为什么要做统计检验？什么是统计检验？为什么要做统计检验？ 前面介绍了几种均匀随机数的产生方法，都是采用某种数学公前面介绍了几种均匀随机数的产生方法，都是采用某种数学公式型来产生的，具有某些特定的数学规律，

22、是伪随机数。这些均匀式型来产生的，具有某些特定的数学规律，是伪随机数。这些均匀分布的伪随机数与理论上的均匀分布随机数有多大差别，是否满足分布的伪随机数与理论上的均匀分布随机数有多大差别，是否满足我们的试验要求？尚不得而知。这就需要我们设立一些能够表征其我们的试验要求？尚不得而知。这就需要我们设立一些能够表征其特性的准则，用这些准则去衡量这些随机数。如果满足这些准则，特性的准则，用这些准则去衡量这些随机数。如果满足这些准则，我们就承认它，否则就拒绝它。这个过程就称作随机数的检验，由我们就承认它，否则就拒绝它。这个过程就称作随机数的检验，由于随机数的各种参数都是用统计方法估计的，故又称作随机数的统

23、于随机数的各种参数都是用统计方法估计的，故又称作随机数的统计检验或假设检验，它依据的是最小概率原理，即在一次试验中，计检验或假设检验，它依据的是最小概率原理，即在一次试验中，概率小的事件认为实际中是不可能发生的。概率小的事件认为实际中是不可能发生的。n如何做统计检验？如何做统计检验？ 已知能表征某种分布随机数的特征参量是已知能表征某种分布随机数的特征参量是 ，根据已产，根据已产生出来的随机数样本用统计方法得到这些特征参量的估计值生出来的随机数样本用统计方法得到这些特征参量的估计值 ，将二者进行比较。如果差别不显著，就承认所产生的随机数符合，将二者进行比较。如果差别不显著，就承认所产生的随机数符

24、合要求，否则就不符合要求而拒绝之。通常用要求，否则就不符合要求而拒绝之。通常用H0表示这样的统计假设表示这样的统计假设，检验过程中是接受还是拒绝，检验过程中是接受还是拒绝H0，一般都给一个称为显著水平的临，一般都给一个称为显著水平的临界概率界概率，观测到的事件概率大于，观测到的事件概率大于，就接受假设，就接受假设H0，小于或等于，小于或等于，就拒绝假设就拒绝假设H0nppp， 21nppp21， 具体地说，我们产生出一组伪随机数具体地说，我们产生出一组伪随机数1，2，N，并假设它，并假设它的某一统计量的某一统计量x=(1，2，N)服从概率分布服从概率分布p(x)，给定一个显著水，给定一个显著水

25、平平，并令，并令=1， 称为置信度，若称为置信度，若 =P(x)，则，则x为临界值。如果为临界值。如果观测值观测值Ei E ，则认为，则认为Ei与理论值差异不显著，接受与理论值差异不显著，接受H0 ；如果；如果Ei E，则认为差异显著，拒绝则认为差异显著，拒绝H0。例如，例如， p(x)=N(0,1)，则有：，则有： p(x)x数数“好好”数数“差差”xa-xa0.80.90.950.980.990.999x1.2821.6451.9602.3262.5763.291n统计检验的步骤统计检验的步骤 （1）提出需要检验的假设，称为原假设，记为）提出需要检验的假设，称为原假设，记为H0 。 （2）

26、构造用于检验的统计量，称为检验统计量，并确定其概率分布。）构造用于检验的统计量，称为检验统计量，并确定其概率分布。 （3）给定显著水平）给定显著水平（0 1）。）。 （4）确定）确定H0的否定域，即根据检验统计量的概率分布和显著水平的否定域，即根据检验统计量的概率分布和显著水平 ， 确定使确定使H0不成立的区域。不成立的区域。 （5）根据样本观测值计算检验统计量之值。）根据样本观测值计算检验统计量之值。 （6）进行统计判断，若检验统计量之值不落入否定域，则接受）进行统计判断，若检验统计量之值不落入否定域，则接受H0，否，否 则否定则否定H0。 2.2.1频率检验频率检验 随机数的频率检验也成为

27、均匀检验，是检验随机数序列的观测频数与随机数的频率检验也成为均匀检验，是检验随机数序列的观测频数与理论频数的差异是否显著。具体地说，就是把理论频数的差异是否显著。具体地说，就是把0,1区间等分成区间等分成k个子区间个子区间，并将包含有，并将包含有N个随机数个随机数1，2，N的随机数序列，按照由小到的随机数序列，按照由小到大的顺序分成大的顺序分成k组。假设组。假设ni是第是第i组的观测频数，那么随机数属于第组的观测频数，那么随机数属于第i组的组的概率为概率为故属于第故属于第i组的理论频数为组的理论频数为构造检验统计量构造检验统计量 (1)则则z渐进服从自由度为渐进服从自由度为k1的的分布，其分布

28、函数为分布，其分布函数为kikpi, 2 , 1,1kikNNpmii, 2 , 1,kiikiiiikNnNkmmnz1212)(0,0, 0de2121)(Flim022321 zzzzkzzkkNN给定显著水平给定显著水平（一般取（一般取 =0.05或或0.01）按照下式计算出拒绝域）按照下式计算出拒绝域za。 根据随机数样本根据随机数样本1，2，N，计算检验统计量，计算检验统计量(1)式的值式的值z，将，将z 与与za进行比较，如果进行比较，如果z za，就认为差异显著，拒绝假设，就认为差异显著，拒绝假设H0；如果；如果z za，则认为差异不显著，接受，则认为差异不显著，接受H0假设。

29、假设。 例如，某长度例如，某长度N=1638的随机数序列，将其分成的随机数序列，将其分成32组，即组，即k=32。则。则mi=512，分布的自由度分布的自由度d=k1=31。当。当 =0.05时，计算出时，计算出za=44.7，根据随机数序列样本计算出的根据随机数序列样本计算出的z=17.4。由于。由于z 1000，k10，mi10。也就是说，。也就是说，在较多的间隔里，每个间隔中所落入的随机数不要太少，结果才可信。在较多的间隔里，每个间隔中所落入的随机数不要太少，结果才可信。 值得注意的是值得注意的是分布的自由度分布的自由度d30时，统计量时，统计量 渐进渐进服从服从N(0,1)分布。分布。

30、azzkkzzkde212122321122dzu2.2.2参数检验参数检验 随机数的参数检验是检验随机数特征参数的观测值与理论值的差异随机数的参数检验是检验随机数特征参数的观测值与理论值的差异是否显著，由于检验的主要参数是各阶矩，故又称为矩检验。是否显著，由于检验的主要参数是各阶矩，故又称为矩检验。 由随机数样本由随机数样本1，2，N，计算出随机变量各阶矩的观测值，计算出随机变量各阶矩的观测值一阶矩和二阶矩，即均值和方差分别为一阶矩和二阶矩，即均值和方差分别为随机变量随机变量的各阶矩和相应方差的理论值值为：的各阶矩和相应方差的理论值值为：NikikNm11NiiNm111NiiNm12211

31、1kmk22,1211kNkmKN一、二阶矩及其方差的理论值分别为：一、二阶矩及其方差的理论值分别为：构造一个统计量构造一个统计量对于对于k=1,2时（一，二阶矩），分别为：时（一，二阶矩），分别为： 根据中心极限定理，统计量根据中心极限定理，统计量zk,N渐近服从正态分布渐近服从正态分布N(0,1)。于是，给定。于是，给定显著水平显著水平，可由正态分布表查得临界值，可由正态分布表查得临界值za。如果统计量的观测值。如果统计量的观测值|zk,N|za，则拒绝假设，则拒绝假设H0；如果；如果|zk,N|50。jNiijijmNs1121111 mjNuj2.2.3概率分布检验（直方图估计）概率分

32、布检验（直方图估计） 设随机变量设随机变量x的样本序列为的样本序列为1，2，N，计算其最大和最小值，计算其最大和最小值将区间（将区间（a,b）等分成）等分成M+1个小区间个小区间ti，ti+1，其中：，其中：统计样本数据落入区间统计样本数据落入区间ti，ti+1的个数的个数ni，并计算出其频率，并计算出其频率 。则有。则有 假设假设x的概率密度函数为的概率密度函数为f(x)，则有，则有于是有于是有N,min21*1NN,max21*1110,NMMxbxabttttaNnfiiMitxtPfiii, 2 , 1 , 0,1MittxfdxxftxtPfiittiiiiii, 2 , 1 , 0

33、, )()(111iiiittfxf1)(例例2.3 混合同余法产生（混合同余法产生（0,1）区间均匀分布随机数）区间均匀分布随机数 编写编写MATLAB程序，取程序，取N=8000，产生出随机数序列。计算出样本，产生出随机数序列。计算出样本均值均值=- -0.4724，样本方差样本方差=1.5015 ，相关系数，相关系数= -0.8627。 取显著水平取显著水平=0.5，则拒绝域为，则拒绝域为1.96。经检验这组随机数样本通过了。经检验这组随机数样本通过了参数检验和相关系数检验，概率密度函数直方图估计如下图。参数检验和相关系数检验，概率密度函数直方图估计如下图。2253535115/mod)

34、1(xrxxnnnn33 统计试验法（统计试验法（Statistical Testing MethodStatistical Testing Method）又称随机抽样技术。用的最多）又称随机抽样技术。用的最多是的是的蒙特卡罗法。蒙特卡罗法。 所谓蒙特卡罗法，是一种通过对实际过程的建模、随机抽样和统计试验来求所谓蒙特卡罗法，是一种通过对实际过程的建模、随机抽样和统计试验来求解各种工程技术、数学物理、社会生活和企业管理各类不同问题的近似解的概率解各种工程技术、数学物理、社会生活和企业管理各类不同问题的近似解的概率统计方法。用蒙特卡罗法解题并不是通过真实试验来完成的，而是抓住事物运动统计方法。用蒙

35、特卡罗法解题并不是通过真实试验来完成的，而是抓住事物运动的基本特征，如统计特性，统计参数，几何特征等，利用数学方法在计算机上对的基本特征，如统计特性，统计参数，几何特征等，利用数学方法在计算机上对真实事物、过程进行仿真，即进行数字统计模拟来完成的。真实事物、过程进行仿真，即进行数字统计模拟来完成的。 蒙特卡罗法真正地用于电子信息系统领域虽然只有四十年左右的历史，但蒙特卡罗法真正地用于电子信息系统领域虽然只有四十年左右的历史，但已已经成为电子信息系统的研究与设计者的非常强有力的工具。当前，把蒙特卡罗法经成为电子信息系统的研究与设计者的非常强有力的工具。当前，把蒙特卡罗法与系统仿真结合起来，在实验

36、室完全可以复现各种复杂的雷达或电子信息系统的与系统仿真结合起来，在实验室完全可以复现各种复杂的雷达或电子信息系统的电磁环境，产生雷达或电子信息系统的各种输入信号、相干或非相干的杂波、干电磁环境，产生雷达或电子信息系统的各种输入信号、相干或非相干的杂波、干扰和噪声，对其进行信号处理、检测以产生点迹、航迹，对各类目标进行外推、扰和噪声，对其进行信号处理、检测以产生点迹、航迹，对各类目标进行外推、跟踪、数据融合，进而进行综合显示。在电子信息领域以下的几个方跟踪、数据融合，进而进行综合显示。在电子信息领域以下的几个方面蒙特卡罗法将会有广泛的应用前景：面蒙特卡罗法将会有广泛的应用前景：34应用前景：应用

37、前景： 雷达或电子信息系统电磁环境的生成；雷达或电子信息系统电磁环境的生成； 用于雷达或电子信息系统的验收与性能评估；用于雷达或电子信息系统的验收与性能评估； 用于雷达或电子信息系统的方案设计，系统性能的优化；用于雷达或电子信息系统的方案设计，系统性能的优化； 研究和选择电子信息系统抗各种有源、无源干扰手段；研究和选择电子信息系统抗各种有源、无源干扰手段； 研究或寻求各种先进的雷达或电子信息系统的体制、波形和先进的技术研究或寻求各种先进的雷达或电子信息系统的体制、波形和先进的技术； 将蒙特卡罗法和电子信息系统仿真、将蒙特卡罗法和电子信息系统仿真、EDAEDA技术结合在一起用于现代雷达或电子技术

38、结合在一起用于现代雷达或电子信息系统的自动化设计等。信息系统的自动化设计等。优点优点1 1经济性：在雷达与电子信息系统的研究和设计方面，利用统计试验法可以节经济性：在雷达与电子信息系统的研究和设计方面，利用统计试验法可以节省大量的人力、物力；省大量的人力、物力；2. 2. 有效性：利用统计试验法对雷达和电子信息系统进行研究可以大大地缩短系有效性：利用统计试验法对雷达和电子信息系统进行研究可以大大地缩短系统的研制周期；统的研制周期；3. 3. 唯一性：利用统计试验法可以完成人工难以完成或根本无法完成的工作，特唯一性：利用统计试验法可以完成人工难以完成或根本无法完成的工作，特别是对那些无法求出数学

39、解或极其复杂的问题，是目前唯一有效的方法；别是对那些无法求出数学解或极其复杂的问题，是目前唯一有效的方法；4 4蒙特卡罗法不仅能解决概率问题，也可以解决非概率问题，在大多数情况下蒙特卡罗法不仅能解决概率问题，也可以解决非概率问题，在大多数情况下都是将非概率问题转化为概率问题来进行求解；都是将非概率问题转化为概率问题来进行求解；353.1 3.1 蒲丰问题蒲丰问题 蒲丰（蒲丰（BuffonBuffon）是一位法国的科学家，他是第一个把统计试验法用来解决实）是一位法国的科学家，他是第一个把统计试验法用来解决实际问题的科学家。际问题的科学家。 他提出了用统计试验法确定园周率他提出了用统计试验法确定园

40、周率的概率统计模型。通常称其为蒲丰模型的概率统计模型。通常称其为蒲丰模型或随机投针试验模型。或随机投针试验模型。 蒲丰问题叙述如下：在一个任意的平面上，划一组平行线，线与线之间的蒲丰问题叙述如下：在一个任意的平面上，划一组平行线，线与线之间的距离为距离为2R2R，R R为一常数。然后向平面上任投一针，为一常数。然后向平面上任投一针，针的长度为针的长度为2l2l，并且有，并且有RlO RlO ，问这一针与任一直线相交的概率是多少？图问这一针与任一直线相交的概率是多少？图3 3中给出了描述这一问题的示意图中给出了描述这一问题的示意图, ,不不过图中只给出了两条平行线，图中过图中只给出了两条平行线，

41、图中M M点是针的中点，点是针的中点，x x是是M M点到最近的一条平行线点到最近的一条平行线的垂直距离，针与相交线之间的夹角为的垂直距离，针与相交线之间的夹角为。 蒲丰模型示意图蒲丰模型示意图 2RXM2L36 从投针方式和图从投针方式和图3-13-1可以看出，针的中点可以看出，针的中点M M等概率地落在长度为等概率地落在长度为R R、垂直于所、垂直于所划平行线的线段上，即随机变量划平行线的线段上，即随机变量x x在在0，1区间上是均匀分布的，其概率密度函区间上是均匀分布的，其概率密度函数为数为其它, 00,1)(RxRxf其次，夹角其次，夹角界于界于1和和1+之间的概率与增量之间的概率与增

42、量的大小成正比，即在的大小成正比，即在0，1区间上是均匀分布的，其概率密度函数为区间上是均匀分布的，其概率密度函数为0,1)( f另外，随机变量另外，随机变量x x与与是互相独立的，针与线相交与否完全由是互相独立的，针与线相交与否完全由x x和和决定，其充要决定，其充要条件为条件为Rxlx0,sin这样，这样，x x和和的联合的联合概率密度函数为概率密度函数为Rxf1),(于是，针线相交的概率为于是，针线相交的概率为dxdxppl),(0sin0dRlsin1037令令y，则有，则有RldyyRlp2)sin(10如果令如果令lR2，则得到更简单的结果，则得到更简单的结果1p可见，只要用统计试

43、验法求出概率可见，只要用统计试验法求出概率p，则，则值的估值就可求出来了，即值的估值就可求出来了，即p1上述过程就称作建模，它把蒲丰模型用一个数学模型表达出来了。上述过程就称作建模，它把蒲丰模型用一个数学模型表达出来了。 1 1 实物统计试验步骤：实物统计试验步骤： 在任一平面上按模型要求划一组平行线，线的间距为在任一平面上按模型要求划一组平行线，线的间距为2 2R;R; 向该平面任意投一根针，针的长度等于向该平面任意投一根针，针的长度等于l 2, ,投针时记录总试验次数投针时记录总试验次数N N和针和针与线相交的次数与线相交的次数m ;m ; 算事件算事件A A出现的相对频率出现的相对频率N

44、mp 38根据大数定理根据大数定理 1)(pNmp可用相对频率可用相对频率p代替概率代替概率p 最最后得到后得到的估值的估值 mNp1表表3-1 3-1 前人给出的前人给出的的估计结果的估计结果 仿真者仿真者时间时间仿真次数仿真次数 的近似值的近似值 Wolf Wolf1850185050005000 3.1596 3.1596 Smith Smith1855185532043204 3.1553 3.1553 Demorgan Demorgan18601860600600 3.137 3.137 Fox Fox1894189411201120 3.1419 3.1419 Lazzarini

45、Lazzarini1901190134083408 3.1415929 3.1415929 Reina Reina1925192525202520 3.1795 3.1795391 1 在计算机上进行统计试验的步骤：在计算机上进行统计试验的步骤： 首先，在电子计算机上产生两个相互独立的随机变量首先，在电子计算机上产生两个相互独立的随机变量y,y,的随机抽样序列的随机抽样序列, 2 , 1,Niyiiiiy,，均服从均服从0, 1区间上的均匀分布；区间上的均匀分布； 然后，分别将随机变量然后，分别将随机变量iiy,变换成变换成0, a区间和区间和0, 0, 区间上的均匀分区间上的均匀分布随机序列

46、布随机序列, 2 , 1,Nixii，即，即NiRyxiiii, 2 , 1, 在计算机上检验不等式在计算机上检验不等式ilxsin即即)sin(iiRly是否成立，如果此式成立，则认为投针试验成功，即针线相交；否则就算失败，是否成立，如果此式成立，则认为投针试验成功，即针线相交；否则就算失败，针与线没有相交；针与线没有相交； 重复步骤重复步骤N N次，得成功次数函数，即次，得成功次数函数，即iiiilxxmsin, 1),(40计算计算p p的估值的估值p和和的估值的估值),(11NiiixmNpNiiixmN1),( 最后对精度进行估计，看是否需要增加试验次数继续进行统计试验。如果满足最后

47、对精度进行估计，看是否需要增加试验次数继续进行统计试验。如果满足精度要求，则结束精度要求，则结束仿真，给出仿真结果。下表仿真，给出仿真结果。下表中给出的一组数据是在计算机上得中给出的一组数据是在计算机上得到的估计结果。从这组数据可以看出，到的估计结果。从这组数据可以看出，值的精度基本上是随着试验次数值的精度基本上是随着试验次数N N的增加的增加而增加的。而增加的。 仿真次数仿真次数 的的 估估 值值 1000 3.31125828 2000 3.16960639 3000 3.18809775 4000 3.16725106 5000 3.20718409 6000 3.21754165 70

48、00 3.01994140 8000 3.17082838 9000 3.14043478 10000 3.14087631表表 : : 值值估估计计结结果果41 综合以上步骤，给出投针试验的程序流程图，如图所示。图中综合以上步骤，给出投针试验的程序流程图，如图所示。图中N0表示满足精表示满足精度的试验次数。度的试验次数。 尽管蒲丰问题比较简单，但通过对它的分析，不难归纳出用统计试验法来尽管蒲丰问题比较简单，但通过对它的分析，不难归纳出用统计试验法来解概率问题的一般步骤：解概率问题的一般步骤： 1根据待解决问题的物理过程，首先建立描述该过程的概率模型，即对问题根据待解决问题的物理过程，首先建立

49、描述该过程的概率模型，即对问题进行建模；进行建模； 2根据要求，确定满足试验精度的试验次数；根据要求，确定满足试验精度的试验次数； 3产生产生0，1区间上的均匀分布随机数；区间上的均匀分布随机数； 4产生仿真系统所需要的随机变量，它们可能是独立的，也可能是相关的；产生仿真系统所需要的随机变量，它们可能是独立的，也可能是相关的； 5根据所选定的方法划出程序流程图，然后选定仿真语言，编出源程序，上根据所选定的方法划出程序流程图，然后选定仿真语言，编出源程序，上机进行仿真；机进行仿真； 6. 对仿真结果进行评估。对仿真结果进行评估。 这里需要说明的是这里需要说明的是, ,在进行仿真实验之前在进行仿真

50、实验之前, ,必须对计算机中所给出的必须对计算机中所给出的0,10,1区间区间上的均匀分布随机数进行统计检验。上的均匀分布随机数进行统计检验。42 图图: : 蒲丰投针试验的程序流程图蒲丰投针试验的程序流程图 启启 动动 产生产生00，11均匀分布随机数均匀分布随机数 y yi , i , i i 产生产生00， aa均匀分布随机数均匀分布随机数 x xi , i , 产生产生00，均匀分布随机数均匀分布随机数 i , i , M=m+1M=m+1 N=N+1N=N+1 计算并打印仿真结果计算并打印仿真结果 停停 机机 X =Lsini X =Lsini N N N BnSnPnPn04214

51、142 Y4111.811.824541 Y4112.324.134445 Y4412.736.844344 Y4312.449.254643 Y4312.962.164346 Y4311.473.574243 Y4212.185.684442 Y4212.698.294344 Y4312.4110.6104043 Y4010.5121.1表表 报童问题的报童问题的1010天的仿真结果天的仿真结果47N N0 01717天的利润天的利润报童使用的决策规则报童使用的决策规则1 13353355 5元元Bn=Dn-1Bn=Dn-12 23473470 0元元Bn=Bn=（Dn-1+Dn-2)/2D

52、n-1+Dn-2)/23 33583585 5元元Bn= Bn= 常数，历史上常数，历史上m m天天D D的平均值的平均值4 43663660 0元元B1=D0B1=D0B2=(D0+D1)/2B2=(D0+D1)/2Bn=(D0+D1+Bn=(D0+D1+Dn-1)/n+Dn-1)/n说明：说明：BnBn- -每天买入的报纸数量每天买入的报纸数量DnDn- -每天社会的需求量每天社会的需求量PB-PB-每份报纸的买入价每份报纸的买入价Ps-Ps-每份报纸的卖出价每份报纸的卖出价过去累计的需求量及相应的概率：过去累计的需求量及相应的概率：40 41 42 43 44 45 4640 41 42

53、 43 44 45 460.05 0.10 0.20 0.30 0.15 0.10 0.100.05 0.10 0.20 0.30 0.15 0.10 0.10表表 采用不同的决策准则时报童所获得的利润采用不同的决策准则时报童所获得的利润483.3 3.3 蒙特卡罗法在解定积分中的应用蒙特卡罗法在解定积分中的应用1 1 概率平均法概率平均法首先让我们考虑一个定积分首先让我们考虑一个定积分badxxfJ)(abxf(x)概率密度函数与积分边界关系图概率密度函数与积分边界关系图首先，根据给定的概率密度函数首先，根据给定的概率密度函数)(xf，产，产生满足该分布规律的随机变量生满足该分布规律的随机变

54、量ix，然后检，然后检验验ix是否落在是否落在区间之内，如果区间之内，如果ix处处于此区间，则认为此次试验是成功的，否则就认为是失败的。在进行于此区间，则认为此次试验是成功的，否则就认为是失败的。在进行N次独立试次独立试验之后，得出成功次数验之后，得出成功次数m，则可得到随机变量，则可得到随机变量落入区间落入区间内的概率估值，内的概率估值，imixNNmp11利用贝努力定理利用贝努力定理, , 假定事件假定事件A A出现的概率为出现的概率为p, 在在N次统计试验中事件次统计试验中事件A出现的次出现的次数为数为m,则对任意正数则对任意正数, , 有有1)(limpNmpN49显然显然, , 当当

55、N N足够大时足够大时, , 取取Nm作为上述积分的近似值是合理的，即作为上述积分的近似值是合理的，即 pp。 由此，得到统计试验步骤如下：由此，得到统计试验步骤如下： 从具有分布规律为从具有分布规律为)(xf的随机总体中抽取随机数的随机总体中抽取随机数ix； 将随机数将随机数ix与区间与区间的界的界a与与b进行比较，比较结果用一个特征指进行比较，比较结果用一个特征指标标来表示，如果满足不等式来表示，如果满足不等式bxai则则, 1否则否则0； 将比较所得到的将比较所得到的值加入一个试验成功计数器，该计数器用值加入一个试验成功计数器，该计数器用m表示，表示，我们称其为我们称其为m计数器；计数器

56、； 每次试验完毕，不管试验成功与否，均在试验次数计数器内加每次试验完毕，不管试验成功与否，均在试验次数计数器内加1 1； 在在N N次统计试验之后，计算成功计数器内的数值与试验次数次统计试验之后，计算成功计数器内的数值与试验次数N N的比值，即的比值，即为所求概率的估值为所求概率的估值Nmpp例例1：用蒙特卡罗法进行概率计算用蒙特卡罗法进行概率计算假定随机变量假定随机变量服从瑞利分布服从瑞利分布)2exp()(222xxxf0 x，试计算试计算小于任一门限小于任一门限T T的概率，这里假定，的概率，这里假定，T=1T=1，1。50 该积分的理论值该积分的理论值22021)2exp(TTedxx

57、xp用统计试验法计算的结果示于下表。用统计试验法计算的结果示于下表。 表表: : 例例1 1的计算结果的计算结果N12345678910 xi1.0420.8111.2250.5260.9272.0711.9060.6351.5861.379m0112333444P00.500.330.500.60.500.430.500.440.40N11121314151617181920 xi3.1541.8451.1331.7862.4810.4762.3911.0760.7273.000m4444455566P0.360.330.310.290.270.310.290.280.320.302020次

58、的统计试验次的统计试验: :只有只有6 6次满足次满足Txi的条件，即的条件，即30. 0p。理论值。理论值: : 393. 0p。512 2 二维随机投点法二维随机投点法下面介绍一种利用均匀分布随机数通过随机投点来计算定积分的方法。下面介绍一种利用均匀分布随机数通过随机投点来计算定积分的方法。首先，假定有一定积分首先，假定有一定积分dxxgJ)(10其被积函数满足条件其被积函数满足条件 1)(0 xg利用均匀分布随机数计算定积分原理见图利用均匀分布随机数计算定积分原理见图3-53-5。yxy=g(x)1y2y1 x2 x11 图图: : 利用均匀分布随机数计算定积分原理图利用均匀分布随机数计

59、算定积分原理图52域是由不等式域是由不等式1010yx确定的，而确定的，而又是由曲线又是由曲线)(xgy 、0X0X轴和轴和1x所限定的区域。所限定的区域。显然，该积分显然，该积分J J等于等于域的面积，因为域的面积，因为域的面积等于域的面积等于1 1。如果我们利用两个相互独立的如果我们利用两个相互独立的00， 11区间上均匀分布随机数区间上均匀分布随机数ix和和iy构成一个二维随机点构成一个二维随机点),(iiyx, 当将当将N N个随机点均个随机点均投到投到域上时域上时, 显然，这些随机点将均匀地分布在显然，这些随机点将均匀地分布在域上，那么落入域上，那么落入域中的域中的点数显然与其面积成

60、正比，而与这些点所落的位置无关。点数显然与其面积成正比，而与这些点所落的位置无关。 如果已知从一均匀分布随机总体中抽出的如果已知从一均匀分布随机总体中抽出的N N个随机点为个随机点为)( ,),(,11NNyxyx，那么它们的联合概率密度函数仍然是均匀分布的，那么它们的联合概率密度函数仍然是均匀分布的， 且有联合概率密度函数且有联合概率密度函数1),(yxf，故随机点，故随机点),(iiyx落入落入域的概率为域的概率为dxdyyxfJp),()(于是，我们得到了计算该积分面积的计算步骤于是，我们得到了计算该积分面积的计算步骤： 产生产生00，11区间上的均匀分布的随机数序列；区间上的均匀分布的