Carlo方法

2024-11-04

Carlo方法（精选7篇）

Carlo方法篇1

引言

期权在金融工程中广泛用于构造各种新型金融产品, 是最活跃的衍生工具之一[2, 3, 5]。期权作为一种衍生金融衍生产品, 它依附于某种标的资产而存在, 因而它的定价决定于标的资产价格的变化。由于标的资产的价格变化是不确定的以及随机的。由此产生的期权的价格变化也定是随机的。但是一旦标的资产价确定下来, 那么期权作为它的衍生产品的价格亦将随之确定。这就是说, 若在t时刻原生资产的价格为S, 期权价格为V, 则存在函数V (S, t) 使得V=V (S, t) 。通常在期权的到期日那天, 期权的价值 (或称为期权的收益、或称为期权的价格) V (S (T) , T) 是确定的。但是期权生效日t=0那天的价值 (或称为期权价格) 是未知的。期权生效日t=0那天的价值也称为期权金, 它是期权购买者为了取得这个期权未定权益所需要付出的代价。通常, 计算目的就是要求出期权生效日t=0那天的价格, 也即V (S (0) , 0) 的值。

欧式期权定价的方法主要有三种:偏微分方程数值解法 (如有限差分法、有限元法、有限容量法等) [1, 8], 二叉树法[2, 4, 7]以及Monte-Carlo方法[9]。三种方法都有各自的优缺点。Monte-Carlo方法最主要的优点是具有一般性, 可以应用于各种金融衍生产品包括期权的定价中。在期权定价中, 有基于二叉树法的Monte-Carlo模拟, 也有基于随机过程的Monte-Carlo模拟[8, 9]。本文仅考虑后者。在国内, 利用Monte-Carlo方法来为欧式期权定价的研究工作主要是体现在标准欧式期权、欧式—亚式期权以及标准美式期权的计算上[1-10], 但是较少有作者系统讨论如何提高Monte-Carlo方法的精确性以及讨论它的优缺点。本文系统讨论各式欧式期权定价的Monte-Carlo方法的实现过程以及如何改进它们。

一、标准欧式期权的定价

在Black-Scholes期权定价模型假设条件下, 可得到计算标准欧式看涨期权的价格公式 (1) ~ (7) 。虽然有计算标准欧式看涨 (或看跌) 期权的价格公式, 但还是希望探讨标准欧式期权定价的Monte-Carlo方法。这是因为标准欧式期权是一种最简单的期权, 许多复杂的期权都是由它派生出来, 一旦掌握了标准欧式期权定价的Monte-Carlo方法, 就可以把该方法应用到更为复杂的期权价格的计算中去。下面先讨论标准欧式期权定价的Monte-Carlo方法。

设S (t) 表示在时刻时标准欧式期权标的资产的价格 (0≤t≤T) , 它是随机变化的。根据Black-Scholes期权定价模型[1, 8, 9, 10], 标的资产价格服从下列随机过程:

其中参数μ是代表标的资产的平均收益率, σ是标的资产价格变化的波动率, W (t) 是服从标准布朗运动的随机变量, d S表示标的资产价格的变化量, dt表示时间的变化量, d W (t) 表示W (t) 的变化量。如果假定μ等于无风险利率r, 则模型式 (1) 变为[1, 4, 7, 9]:

再根据风险中性理论就可得[1, 9, 10]:

这里假设欧式期权在T时刻的收益为某一确定收益函数f (S (T) ) (对于标准欧式期权, f (S (T) ) 为max (S-K;0) 或max (K-S;0) ) 。从式 (3) 可知, 期权生效日t=0那天的价值 (也即期权的价格) 为:

公式 (4) 是Monte-Carlo方法计算欧式期权价格的起点, 下面将利用Monte-Carlo方法来计算公式 (4) 中出现的数学期望。

在模型 (2) 式假设下, 有。因此有:

其中S (0) 是在t=0时刻的标的资产价格。由于W (t) 服从均值为0、方差为的T正态分布。如假设Z是服从标准正态分布的随机变量, 则可通过下式来计算式 (5) :

通过用Monte-Carlo方法分别近似式 (6) 以及式 (4) , 从而就可得到欧式期权的价格近似值。为此, 首先利用Matlab自带的程序生成服从标准正态分布的随机变量Z, 然后由Z和式 (6) 就可生成随机变量S (T) 。再通过下面的算法就可以得标准欧式看涨期权的价格的近似值:

for i=1, …, n

随机产生服从标准正态分布的随机变量zi:

end

令, 则为标准欧式看涨期权的价格的近似值。根据概率论中的大数定理, 易知:当n→+∞时。

说明1:如果把 (7) 式中的e-r Tmax (Si-K, 0) 改为e-r Tmax (K-Si, 0) , 那么可通过上面的算法可以得到标准欧式看跌期权的价格的近似值。

然而, 如果采用上面的标准Monte-Carlo方法, 计算误差可能较大。可通过下面两种方法来分别减少此方法的计算误差。

(一) 通过引入控制变量法 (简记为CV) 来减少方法的计算误差

设是通过Z (为一服从标准正态分布的随机变量) 得到的欧式期权价格, 而是通过Z* (也为一服从标准正态分布的随机变量) 得到的欧式期权价格;V*是真正的欧式期权价格, 它为一常数。如果的相关系数较大, 也即:

从而:。如果我们再引入控制变量:

由于, 考虑到式 (8) , 从而有

(二) 通过引入对称变量法 (简记为AV) 来减少方法的计算误差

设是通过Z (为一服从标准正态分布的随机变量) 得到的欧式期权价格, 而是通过-Z得到的欧式期权价格。可令:

它是通过引入对称变量得到的新估计值。它可用来近似欧式期权的价格, 并且它的方差满足:。

二、其他各种欧式期权的定价

在本节里, 主要讨论如何把第二节之中介绍的三种Monte-Carlo方法来推广到欧式—两值期权、欧式—回望期权以及欧式—亚式期权的价格计算中。

(一) 欧式—两值期权定价的Monte-Carlo法

欧式—两值期权具有两种形式:现金期权和资产期权。它们又分为看涨期权与看跌期权。现金看涨期权是指在到期日标的资产价格低于执行价格时收益为零, 当标的资产价格超过执行价格时收益为一固定数额现金Q。因而, 现金看涨期权在到期日的价格 (或收益) 为:

资产看涨期权在到期日标的资产价格低于执行价格时收益为零, 标的资产价格超过执行价格时收益为标的资产价格本身的款额, 因而资产看涨期权收益为:

由于欧式—两值期权是弱路径依赖期权, 欧式—两值期权价格的Monte-Carlo计算过程与标准欧式期权的算法过程式 (7) 基本类似。

对于现金看涨期权, 只要把式 (7) 中的改为:

对于资产看涨期权, 只要把式 (7) 中的改为:

同理, 也可为现金看跌期权价格与资产看跌期权价格的计算构造相应的Monte-Carlo计算法。

(二) 欧式—回望期权定价的Monte-Carlo法

在期权的到期日那天 (t=T) , 欧式—回望期权的收益取决于在一段特定时间内标的资产价格变化的最高价或最低价。通常有两种类型的欧式—回望期权, 它们分别是固定的回望期权和浮动的回望期权。

假设J为一段时间 (0≤t≤T) 内资产变化的最高价格, 则浮动的回望看跌期权的收益为f (S (T) ) =max (J-S (T) , 0) , 这里。相应地, 浮动的回望看涨期权的收益为f (S (T) ) =max (S (T) -J, 0) 这里。固定的回望看跌期权的收益为f (S (T) ) =max (K-J, 0) , 这里。同理, 固定的回望看涨期权的收益为f (S (T) ) =max (J-K, 0) 这里。

下面仅以浮动的回望看涨期权为例子, 采用n次离散抽样, 即得到式 (6) , 从而就可知道最小值。于是, 对浮动的欧式—回望看涨期权的价格E[e-r Tf (S (T) ) ], 可以采用以下的Monte-Carlo方法:

说明2:上面的算法可以用来计算浮动的欧式—回望看跌期权的价格, 只要将算法式 (11) 中的 (S (T) -J) +改为 (J-S (T) ) (+J=max (St1, St2, …Stn) ) 。

说明3:上面的算法也可以用来计算固定的欧式—回望看涨期权的价格, 只要将算法式 (11) 中的 (S (T) -J) +改为 (K-J) +即可。

(三) 欧式—亚式期权定价的Monte-Carlo法

欧式—亚式期权赋予期权持有者在规定的时间内以平均价格买入标的资产的权利。但是必须决定如何进行抽样, 通常有两种求抽样的方式以及有两种求平均值的方法 (见表1) 。

亚式期权又可分为浮动的执行价格亚式期权以及固定的执行价格亚式期权。当执行价格依赖于平均价格时, 这种期权就叫作浮动的执行价格亚式期权;当执行价格固定时, 就称该种期权为固定的执行价格亚式期权。各种亚式期权在交割日t=T那天的收益可以总结为:

对于固定的亚式期权:看涨期权的收益为:f (S (T) ) =max{AT-K, 0};看跌期权的收益为:f (S (T) ) =max{K-AT, 0};对于浮动的亚式期权:看涨期权的收益为:f (S (T) ) =max{S-AT, 0};看跌期权的收益为:f (S (T) ) =max{AT-S, 0}。

下面仅以固定的亚式看涨期权为例子, 也采用n次算术平均离散抽样, 这里需要先知道平均值。为计算固定的亚式看涨期权的价格, 也即, 可采用以下的Monte-Carlo方法:

说明4:上面的算法可以用来计算固定的亚式看跌期权的价格, 只要将算法式 (12) 中的改为。

说明5:上面的算法也可以用来计算几何平均抽样下的固定亚式看涨期权的价格, 只要将算法式 (12) 中的改为即可。

三、数值计算结果

在本节里, 分别记标准的Monte-Carlo方法计算得到的价格为V、通过引入控制变量而改进的Monte-Carlo方法计算得到的价格为CV、通过引入对称变量得到的改进Monte-Carlo方法计算得到的价格为AV。在编程计算中, 对看涨期权, 控制变量为标准欧式看涨期权相对应的参变量;对看跌期权, 控制变量为标准欧式看跌期权相对应的参变量。

(一) 利用Monte-Carlo方法来计算标准欧式期权的价格

假设当前标的资产价格S=100, 执行价格K=100, 短期利率r=0.1, 到期日T=1 (年) , 波动率σ=0.25。对标准看涨期权以及看跌期权, 它们价格都有准确公式。通过公式, 得到看涨期权准确价格为14.9758、看跌期权准确价格5.45954。对于看涨期权与看跌期权, 编程计算结果分别 (见表2、下页表3) 。

观察上页表2与表3不难看出:未改进的Monte-Carlo方法得到的误差较大;通过引入控制变量和引入对称变量, 都可提高Monte-Carlo方法的精确度;相比较而言, 通过引入控制变量得到的改进Monte-Carlo方法精确度较高, 通过引入对称变量得到的改进Monte-Carlo方法精确度相对较低。随着模拟次数n的增大, 各种Monte-Carlo方法精确度有相对的提高。

(二) 利用Monte-Carlo方法来计算欧式—两值期权的价格

假设当前标的资产价格S=100, 执行价格K=100, 短期利率r=0.1, 到期日T=1 (年) , 波动率σ=0.25, 现金值Q=1。对于现金看涨期权与现金看跌期权, 通过编程算得它们的价格, 计算结果分别 (见表4、表5) 。

另外, 对于资产看涨期权与资产看跌期权, 它们的价格分别 (见表6、下页表7) 。

观察表6与下页表7可以看出:随着模拟次数n增大, 标准的Monte-Carlo方法与通过引入控制变量得到的改进Monte-Carlo方法的收敛性较好。通过引入对称变量得到的改进Monte-Carlo方法的收敛性较差。当模拟次数n很大时 (如n=50 000) , 三者得到的结果非常接近。

(三) 利用Monte-Carlo方法来计算浮动的欧式—回望期权的价格

假设当前标的资产S=100, 执行价格K=120, 短期利率r=0.1, 到期日T=1 (年) , 波动率σ=0.25。对于看跌期权, 通过计算得其准确价格为14.2665;对于看涨期权, 通过计算得其准确价格为22.8089 (其公式参见[1]) 。

1. 对于看跌期权, 通过编程算得它的价格 (见表8) 。

2. 对于看涨期权, 通过编程算得它的价格 (见表9) 。

观察表8、表9不难看出:随着模拟次数m的分别增大, 三种Monte-Carlo方法的收敛性都较好。在计算看涨期权的价格时, 三种Monte-Carlo方法的误差都较小;但是在计算看跌期权的价格时, 三种Monte-Carlo方法的误差都较大 (通过增加模拟次数n, 可适当减少此误差) 。

(四) 利用Monte-Carlo方法来计算欧式—亚式期权的价格

假设当前标的资产价格S=100, 执行价格K=120, 短期利率r=0.1, 到期日T=1 (年) , 波动率σ=0.25。在下面的计算中仅仅考虑离散抽样和固定的亚式期权的计算。并且分别记算术平均抽样的亚式看涨期权为 (I) 、算术平均抽样的亚式看跌期权为 (II) 、几何平均抽样的亚式看涨期权为 (III) 以及几何平均抽样的亚式看跌期权为 (IV) 。通过编程算得它们的价格, 计算结果分别 (见表10与表11 (在计算中, 抽样次数m=5 000;n=250) ) 。

观察表10、表11不难看出:除了亚式期权 (III) 的计算结果, 三种方法的计算结果均较为接近。产生的原因可能是:在亚式期权 (III) 的计算中, 生成的随机数是伪随机数, 它们可能造成计算产生较大误差。

结论

根据Black-Scholes期权定价模型以及风险中性理论, 本文详细地讨论了各种欧式期权价格的Monte-Carlo方法。通过编程计算, 可以发现Monte-Carlo方法的主要优点是该方法具有一般性, 可推广应用于各种期权价格的计算中。缺点是精确度较低, 计算量较大。本文也发现通过减少随机变量的方差 (如引入控制变量以及对称变量) 是可以大幅提高Monte-Carlo方法的精确度。笔者相信本文研究的Monte-Carlo方法可以推广应用到其他类型的期权包括美式期权、利率期权以及债券期权等的定价中。

摘要：讨论各种欧式期权价格的Monte-Carlo方法。根据Black-Scholes期权定价模型以及风险中性理论, 首先详细地讨论如何利用Monte-Carlo方法来计算标准欧式期权价格;然后讨论如何引入控制变量以及对称变量来提高Monte-Carlo方法的精确性;最后用Monte-Carlo方法来计算标准欧式期权、欧式—两值期权、欧式—回望期权以及欧式—亚式期权的价格, 并讨论相关方法的优缺点。

关键词：Black-Scholes方程,欧式期权定价,Monte-Carlo方法

参考文献

[1]姜礼尚.期权定价的数学模型和方法[M].北京:高等教育出版社, 2003.

[2][英]Robert.Tompkins.解读期权[M].陈宋生, 崔宏, 刘锋, 译.北京:经济管理出版社, 2003.

[3]姜礼尚, 等.金融衍生产品定价的数学模型与案例分析[M].北京:高等教育出版社, 2007.

[4]马俊海.金融衍生证券定价的数值分析方法[M].杭州:浙江人民出版社, 2002.

[5]J.C.Hull.Options, futures and other derivatives, Prentice-Hall, London, 1989.

[6]P.Wilmott, S.Howison, J.Dewynne.The mathematics of financial derivatives[M].Camberidge University Press, Cambridge, 1995.

[7]P.Wilmott, J.Dewynne and J.Howison.Option Pricing:Mathematical Models and Computation[M].Oxford Financial Press, Oxford, 1993.

[8]Rudiger U.Seydel.Tools for Computational Finance (3rd Edition) [M].Springer, 2004.

[9]Paul Glasserman.Monte Carlo Methods in Financial Engineering[M].Springer, 2003.

[10]宗琮, 傅文月, 罗秋瑾, 王汉权.求解期权定价问题的一种Chebyshev—谱方法[J].福州大学学报 (自然科学版) , 2012, (2) :150-153.

Carlo方法篇2

Monte-Carlo统计迭代图像重建算法

线性衰减系数图像重建是层析γ扫描(TGS)的一个核心问题.本文从粒子输运方程出发,应用Monte-Carlo方法,提出了一种基于Monte-Carlo方法的`统计迭代图像重建算法.模拟结果表明,与一般TGS图像重建算法相比,该重建算法的重建图像误差大为减小,能够满足TGS装置的要求.

作者：张全虎隋洪志吕峰李泽作者单位：中国原子能科学研究院,放射化学研究所,北京,102413 刊名：原子能科学技术 ISTIC EI PKU英文刊名：ATOMIC ENERGY SCIENCE AND TECHNOLOGY 年，卷(期)： 37(6) 分类号：O242.2 关键词：层析γ扫描线性衰减系数 Monte-Carlo方法图像重建迭代法

Carlo方法篇3

关键词：Monte Carlo模拟；Hakanson潜在生态危害指数；敏感性分析

中图分类号： X825文献标志码： A文章编号：1002-1302（2014）06-0320-03

收稿日期：2013-09-10

基金项目：国家科技支撑计划（编号：2011BAD13B02）。

作者简介：曲春风（1981—），男，山東牟平人，硕士，实验师，研究方向为环境生态学。Tel：（0631）5688551；Email：qcf_81@163.com。

通信作者：葛长字，博士，副教授。E-mail：changzige@ouc.edu.cn。由于人类的工农业生产活动，使得大量的重金属进入到土壤中，其赋存形态因生物活动而发生改变但不能降解，致使土壤中的重金属得以积聚[1]。土壤中的重金属除污染水质外，还通过食物链对人类健康造成危害[2]。土壤中重金属的环境监测、生态危害评价等成为国内外环境、生态等领域科学工作者关注的热点[3-4]。随着“镉大米事件”的曝光[5]，我国土壤重金属污染问题更加受到重视。目前，普遍采用的土壤重金属污染的评价方法有单因子指数法、地质累计指数法、内梅罗污染指数法、污染负荷指数法和潜在生态危害指数法等[6]。Hakanson提出的潜在生态危害指数法[7]综合考虑了重金属的生态毒性及重金属区域背景值的差异，消除了区域差异影响，在国内外土壤、沉积物评价中得到广泛的应用[8-9]。该评价方法不仅受制于监测样品数量的有限性、区域污染物分布的时空不均匀性以及实测数据的误差性，而且对于重金属生态毒性的权重赋予具有一定的主观性。潜在生态危害指数法不能有效解决土壤中重金属污染评价的不确定性，即评价的风险性难以体现，而评价的风险性是决策的重要依据之一。将风险评价系统的不确定性引入到潜在生态危害指数中，将能更科学地评估土壤重金属污染的现状。本研究将Monte Carlo模拟引入到潜在生态危害指数法中，对土壤中的Mn、Cu、Pb、Zn、Cd等5种重金属的潜在生态危害及其风险性进行评估，并对各种重金属浓度的参数敏感性进行了分析。

1材料与方法

1.1基于Monte Carlo模拟的潜在生态危害评价方法

1.1.1Hakanson潜在生态危害指数法该方法由Hakanson于1980年提出，包含环境化学、生物毒理学、生态学等方面内容，常应用于沉积物及土壤的生态危害评价。基本公式为：

Cif=Ci/Cin

Cd=∑mi=1Cif

Eir=Tir·Cif

RI=∑mi=1Eir=∑mi=1Tir·Cif=∑mi=1Tir·Ci/Cin

式中：Cif为金属i污染系数，Ci为金属i的实测浓度，Cin为金属i的背景值；Cd为多种金属综合污染程度；Tir为金属i的生物毒性响应因子；Eir为金属i潜在生态风险因子；RI为多种金属潜在生态危害指数。

1.1.2Monte Carlo模拟Monte Carlo模拟因著名的摩洛哥赌城而得名，又称为随机模拟，基本思想是根据待求随机问题的变化规律，根据物理现象本身的统计规律，或者人为的构造一个合适的概率模型，依照该模型进行大量的统计试验，使它的某些统计参量正好是待求问题的解[10]。随着计算机技术的快速发展，使得快速进行Monte Carlo模拟成了可能，通过计算机模拟随机过程，进行大量模拟试验，并统计计算结果。

1.1.3基于Monte Carlo模拟的潜在生态危害指数评价利用Monte Carlo方法进行土壤中重金属潜在生态危害评估：（1）通过实测重金属浓度确定各种重金属在土壤中分布的概率密度函数；（2）通过Monte Carlo方法对概率密度函数进行随机抽样，获得重金属在土壤中的分布浓度，通过Hakanson潜在生态危害指数的方法，对土壤重金属潜在生态危害进行评估。随机模拟过程及评价结果的计算由Crystal Ball 11.1软件进行。对于不同金属的生物毒性响应因子Tir值，本研究参考经典的Hakanson指数以及国内研究成果[11-13]，设定5种重金属Tir的顺序为：Cd=30>Cu=Pb=5>Zn=Mn=1。

1.2评价区域

后湖农场位于湖北省潜江市的中部，北倚汉水、南近长江，处于江汉平原腹地；四湖主干渠和东干渠在这里交汇，并有沟通长江、汉水间的内河航运。318国道、宜黄高速公路、襄岳公路在这里立体交汇，使后湖成为江汉平原乃至整个湖北中部一个新兴的水陆交通枢纽。后湖农场包括张家窑、天新、关庙、前湖、皇装烷、流塘6个分场，拥有耕地 3 000 hm2、林地900 hm2、水面1 300 hm2；重点发展粮、棉、油、渔、猪、果六大类产品，是湖北省粮棉油高产区[14]。后湖农场作为湖北省重要的粮食产区之一，对该地区土壤中的重金属进行潜在生态危害评价，对于促进生态农业，维护人体健康都有重要的现实意义。评价区域的重金属监测数据来自于聂燕博士学位论文[14]，土壤重金属浓度背景值选取湖北省土壤元素背景值[15]。

2结果与分析

2.1潜在生态危害评价指数

5种重金属中平均值只有Mn、Pb稍低于湖北省的土壤元素背景值，Cu、Zn、Cd的平均值都超标，超标倍数为Cd为6.86>Zn为1.29>Cu为1.28；测定的重金属Mn、Cu、Zn、Pb的变异系数都不高，在土壤中分布相对比较均匀，Cd的变异系数较大，为58.63%，即土壤中的Cd分布不均匀，可能存在某些点源性的污染，统计结果见表1。表1后湖农场5种土壤重金属元素质量分数统计特征值

统计值含量MnCuZnPbCd最大值（mg/kg）539.46753.186193.35536.092.821最小值（mg/kg）480.20822.93048.55015.4900.030中值（mg/kg）506.03039.370104.33025.8201.160平均值（mg/kg）504.91039.590108.15025.8601.180CV（%）3.75020.90027.21019.87058.630概率分布模型

Weibull

（452.8，59.1，3.05）Beta

（20.64，55.17，1.88，1.54）Lognormal

（108.2，29.3，-43.7）Logistic

（25.92，2.93）Beta

（-0.27，4.36，2.79，6.11）湖北土壤背景值（mg/kg）642.030.783.626.70.172

根據Monte Carlo 模拟获得的5种重金属的概率密度分布函数，在Excel软件中，利用Crystal Ball软件将5种重金属Ci定义为假设变量，将5种重金属的生物毒性响应因子Tir定义为决策变量，将5种重金属的潜在生态风险因子Eir和多种金属潜在生态危害指数RI定义为预测值。利用Crystal Ball软件，在95%置信度的条件下，对Ci进行20 000次随机取样的Monte Carlo模拟，获得Eir和RI的预测值概率密度分布（图1）。5种重金属Eir预测值的概率分布函数与5种重金属浓度的概率分布函数类型相同（表1）中的概率分布模型，而多种金属潜在生态危害指数RI的预测值服从Beta分布。

对于重金属潜在生态危害的评价，经典的Hakanson指数评价了8种污染物[7]，本研究仅评价了5种重金属的潜在生

态危害，由于评价的重金属数量上的调整，评价指标Cif保持不变的情况下Cd需要做相应的调整。根据国内外学者风险评价指标体系的研究[16-18]，本研究采用的潜在生态危害指数评价分级标准如表2所示。表2Hakanson潜在生态危害指数评价分级标准

CfCd污染程度EirRI风险程度Cif<1Cd<8轻度污染Eir<40RI<150低潜在生态危害1≤Cif<38≤Cd<16中度污染40≤Eir<80150≤RI<300中等潜在生态危害3≤Cif<616≤Cd<32重度污染80≤Eir<160300≤RI<600较高潜在生态危害Cif≥6Cd≥32非常重污染160≤Eir<320-高潜在生态危害[4]Eir≥320RI≥600极高潜在生态危害

利用Crystal Ball软件对Eir和RI的预测结果进行统计分析，得到Mn、Cu、Zn、Pb、Cd等5种重金属的Eir平均值分别为0.79、6.47、1.29、4.85和206.2；由Crystal Ball软件计算出不同概率条件下各重金属Eir和多种金属潜在生态危害指数RI（表3），通过与潜在生态危害指数评价分级标准进行比较，得到不同概率条件下的污染程度和潜在生态危害等级（表4）。在评价区域内，Mn、Cu、Zn、Pb的污染程度很低，100%处于低潜在生态危害程度以内。Cd的潜在生态危害等级较高，处于低潜在生态危害、中等潜在生态危害、较高潜在生态危害、高潜在生态危害、极高潜在生态危害的概率分别为6.33%、880%、2403%、42.74%、18.10%，因此，对于评价区域Cd的污染应该给予足够重视。通过对评价区域多种金属潜在生态危害指数RI分析得出，该区域农田土壤重金属处于低潜在生态危害、中等潜在生态危害、较高潜在生态危害、极高潜在生态危害的概率分别为31.35%、43.92%、24.63%、0.10%，判断评价区域处于中等潜在生态危害等级。

污染等级概率（%）MnCuZnPbCd RI低生态危害1001001001006.3331.35中等生态危害00008.8043.92较高生态危害000024.0324.63高生态危害000042.74-极高生态危害000018.100.10

2.2重金属浓度的参数敏感性分析

敏感性分析就是令模型的每个参数在可能取值的变化范围内变动，预测这些参数的变动对模型输出值的影响程度，将影响程度的大小称为参数的敏感性系数，实质就是研究参数变化所引起的模型响应[19-20]。通过考察Mn、Cu、Zn、Pb、Cd 5个参数在取值范围内的变化对RI的影响程度，确定各参数的敏感性系数，通过对敏感性系数大小的分析，来判断各参数对评价结果的影响程度，重点考虑敏感性系数较大的参数。利用Crystal Ball软件对评价区域的Hakanson潜在生态危害评价的5种重金属Ci进行参数敏感性分析（图2）。评价区域对RI起主导作用的是Cd，敏感性系数达99.1%，Mn、Cu、Zn、Pb对RI影响很小，这主要是由于评价区域Mn、Cu、Zn、Pb数值较低，在评价区域都处于低生态危害程度，多金属潜在生态危害指数RI对于这4种重金属浓度的敏感性非常低。

3结论

将Monte Carlo模拟技术应用于Hakanson潜在生态危害评估中，可以很好规避由于所评价区域重金属含量的不确定性所造成的风险评估误差，对需要考虑不确定性因素的评估有一定指导意义。

从评价结果可以看出，评价区域总体处于中等潜在生态危害等级。江汉平原后湖地区作为粮棉油的中高产区，Cd处于较高的风险等级，处于较高潜在生态危害及以上的风险概率达84.87%，对于该地区土壤中的Cd污染应该引起重视。

对评价区域Hakanson潜在生态危害指数的敏感性分析表明，评价区域对多种金属潜在生态危害指数RI起主导作用的是Cd，敏感性系数达到了99.1%。

参考文献：

[1]崔德杰，張玉龙. 土壤重金属污染现状与修复技术研究进展[J]. 土壤通报，2004，35（3）：366-370.

[2]王宏镔，束文圣，蓝崇钰. 重金属污染生态学研究现状与展望[J]. 生态学报，2005，25（3）：596-605.

[3]李如忠，潘成荣，陈婧，等. 铜陵市区表土与灰尘重金属污染健康风险评估[J]. 中国环境科学，2012，32（12）：2261-2270.

[4]温雅君，孙江，高景红，等. 北京市蔬菜基地土壤重金属含量的环境质量分析与评价[J]. 中国农学通报，2013，29（14）：129-133.

[5]郭顺姬. 治理“镉大米”不容拖延[N]. 中国经济时报，2013-05-24（7）.

[6]郭笑笑，刘丛强，朱兆洲，等. 土壤重金属污染评价方法[J]. 生态学杂志，2011，30（5）：889-896.

[7]Hakanson L. An ecological risk index for aquatic pollution control：a sedimentological approach[J]. Water Research，1980，14（8）：975-1001.

[8]李如忠，潘成荣，徐晶晶，等. 基于Monte Carlo模拟的潜在生态危害指数模型及其应用[J]. 环境科学研究，2012，25（12）：1336-1343.

[9]马德毅，王菊英. 中国主要河口沉积物污染及潜在生态风险评价[J]. 中国环境科学，2003，23（5）：521-525.

[10]吉庆丰.蒙特卡罗方法及其在水力学中的应用[M]. 南京：东南大学出版社，2004：1-13.

[11]刘文新，栾兆坤，汤鸿霄.乐安江沉积物中金属污染的潜在生态风险评价[J]. 生态学报，1999，19（2）：206-211.

[12]徐争启，倪师军，庹先国，等. 潜在生态危害指数法评价中重金属毒性系数计算[J]. 环境科学与技术，2008，31（2）：112-115.

[13]何云峰，朱广伟，陈英旭，等. 运河（杭州段）沉积物中重金属的潜在生态风险研究[J]. 浙江大学学报：农业与生命科学版，2002，28（6）：669-674.

[14]聂艳.耕地质量评价的模型方法与信息系统集成及应用研究[D]. 武汉：华中农业大学，2005：143-150.

[15]国家环境保护局主持，中国环境监测总站.中国土壤元素背景值[M]. 北京：中国环境科学出版社，1990：330-483.

[16]李如忠，徐晶晶，姜艳敏，等. 铜陵市惠溪河滨岸带土壤重金属形态分布及风险评估[J]. 环境科学研究，2013，26（1）：88-96.

[17]李莲芳，曾希柏，李国学，等. 北京市温榆河沉积物的重金属污染风险评价[J]. 环境科学学报，2007，27（2）：289-297.

[18]Pekey H，Karaka

瘙塂 D，Ayberk S，et al. Ecological risk assessment using trace elements from surface sediments of Izmit Bay （Northeastern Marmara Sea）Turkey[J]. Marine Pollution Bulletin，2004，48（9/10）：946-953.

[19]蔡毅，邢岩，胡丹. 敏感性分析综述[J]. 北京师范大学学报：自然科学版，2008，44（1）：9-16.

Carlo方法篇4

天然气管道泄漏模拟是天然气管道应急决策指挥的重要依据,近年来越来越受到重视。气体的扩散受风速、风向、湍流或大气稳定度、地面粗糙度以及温度等因素的影响。传统描述流体运动的基本方程是Navier-Stokes方程和质量连续方程,以及在考虑湍流因素情况下的K-ε两方程。为了求解以上方程,一般运用数值计算的方法,包括有限单元、有限差分、简化解析解以及无网格方法,其求解计算过程涉及大量的数值计算。对于有限空间、小尺度的三维流体计算,方程求解可得出相当高精度的数值解。然而天然气管道泄漏模拟所要求的是大尺度、复杂地形、复杂气象条件下的较高精度的计算,同时还要求有较短时间的运算过程,目前很多模型很难满足。随着计算机技术以及实验物理的发展,不断完善的Monte Carlo模型开始引入大气扩散模拟中,并在天然气管道灾害应急决策指挥过程中起到重要作用。

1 扩散模型综合分析

大气扩散模型的研究最早开始于1915年,Taylor研究热在相对低温的洋流中传导,提出了著名的Taylor湍流扩散理论,后来Scrase(1930)和Best(1935)扩展了Taylor的研究,揭示了热结层以及较宽频谱波动的存在。近年来,大气扩散模型的研究发展很快,最常用的包括:

Gauss模型是实用最早的扩散模型,也是计算点源扩散浓度中最普遍的一种方法[1]。它假定在下风向气体浓度在水平与垂直方向有独立的高斯分布,属于非重气云扩散模型。Gauss烟羽模型适用于连续点源的泄漏扩散,而烟团模型适用于瞬时点源的泄漏扩散。Gauss烟团模型将泄漏源作为一系列的烟团释放到大气中。风场分布、大气稳定度在每一次的计算过程中是不变的,扩散过程中,烟团中心由轨迹确定,烟团分布服从高斯分布[2]。

重气模型用来模拟重气烟羽扩散。最常见的模型包括DEGADIS和SLAB。DEGADIS针对密度比空气重的气体在平坦地面的泄漏扩散。它所融合的Ooms喷射烟羽模型提供了气体垂直喷射模型的轨迹、浓度预测[3]。SLAB模型的泄漏源既有连续泄漏、短时泄漏,也有瞬时泄漏,还可应用于蒸气泄漏、水平和垂直泄漏[4]。

FEM3C模型是重气泄漏三维有限元扩散模型,利用计算机数值技术求解三维动态流体方程。FEM3C模型最大的特点是能计算分析复杂气象条件、复杂地形环境对气体扩散的影响。FEM3C模型采用了网格矢量化的方法加快计算速度,在计算复杂地形条件下的气体扩散时会消耗大量时钟周期。

拉格朗日模型跟随大气中粒子的扩散轨迹,建立随机行走过程。它使用一定数量的粒子来模拟物理参数的动态变化。通过运用Monte Carlo方法,建立确定的或半确定的速度变量。因此,粒子的运动即包括平均风场也包括风的湍流变化。欧拉模型主要是用固定的三维笛卡尔坐标来研究。

2 Monte-Carlo长输管道气体扩散模型

各类现有的气体扩散模型,尤其应用于长输管道事故应急决策领域内时,都存在一定的优势与不足。为了达到应急指挥决策机构过程的科学、高效、迅速等特点,综合考虑地形因素、气象因素等复杂条件影响,本文采用了基于Monte-Carlo(M-C)模型的气体扩散模拟。而气象条件中的风场预测采用了美国NOAA提出了区域大气预测模型(RAMS)。该模式是由Cotton提出的中尺度动力系统与微物理过程模式和中尺度和陆面特性模式合并而成的,动力框架是非流体静力、原始方程中尺度模式,模式的垂直坐标采用地形追随坐标,对有效利用各种探空资料提供了较好的条件。另外,针对Monte-Carlo在过去应用过程中存在的随机数“伪随机”的问题,研究基于计算机硬件信息的随机数产生。

2.1 蒙特卡罗(Monte Carlo)模型

大气扩散本身是一种随机过程,这种行为可以用大量质点随机行走的方法来模拟。如果质点对前一个时刻的运动没有记忆,而且在任何时刻各个方向的运动机会均相等,这样的质点运动轨迹成为遵循Monte-Carlo路径[5]。Monte-Carlo质点随机行走模式就是把每个污染质点当成有标志的质点,通过释放大量粒子,让它们在流场中按平均风输送,同时又用一系列随机位移来模拟湍流扩散,这样就表达了平流和湍流扩散两种作用。通过跟踪污染物例子的轨迹,可以导出轨迹的总体统计特征,从而模拟估算出污染物的空间分布和时间变化规律。

该模型能很好地反映扩散本身具有的随机特性,而且不需要附加任何特殊的假设。Monte-Carlo模式完全脱离了K理论框架,真实直观地反映了污染物在大气运动中的输送规律。既适用于平坦均匀的下垫面,又适用于复杂的地形,而且模式运行中可避免数值计算不稳定、负浓度等问题。

Monte-Carlo质点随机行走模型质点的位置是由平均风速引起的位移和湍流引起的波动共同决定的。这些质点的运动一般是在平均风速(来源于大气模型)的基础上加上一个随机变量[6]。在水平范围内,计算方程表达如下:

X(t+Δt)=Xmean(t+Δt)+U′(t+Δt)Δt,

Z(t+Δt)=Zmean(t+Δt)+W′(t+Δt)Δt (1)

式(1)中Xmean(t+Δt)为平均风速引起的位移,U′(t+Δt)、W′(t+Δt)分别为水平和垂直方向湍流风速,Δt为积分时间,一般Δ t = (Δ Z)2/ ( 8 σw2TLW),ΔZ为垂直栅格大小,TLW为拉格朗日时间尺度,σw 2为垂直速度方差。

水平方向湍流速度可由前一时刻湍流速度、自相关系数、拉格朗日时间尺度和随机数计算得到,如式(2):

U′(t+Δt)=R(Δt)U′(t)+U″(1-R(Δt)2)0.5,

R(Δt)=exp(-Δt/TLx),

U″=σuλ (2)

其中,R(Δt)为与时间步长相关的自相关系数,U′(t)前一时间湍流速度,TLu水平方向拉格朗日时间尺度,σu湍流速度的标准方差,λ为由计算机产生的(0,1)高斯正态分布随机数。

大气湍流研究表明,对流边界层中垂直方向的湍流是非均一的,并不能像水平方向那样简单表示,计算公式如式(3):

W′/σw(t+Δt)=R(Δt)(W′(t)/σw(t))+(W″/σw(t))(1-R(Δt)2)0.5+TLW(1-R(Δt))∂σw(t)/∂z,

σw(t+Δt)=σw(t)+W′(t)Δt∂σw(t)/∂z (3)

基于Monte-Carlo粒子随机行走模型,美国国家海洋与大气管理局所属的空气资源实验室和澳大利亚气象署研究中心联合开HYSPLIT(Hybrid Single-Particle Lagrangian Integrated Trajectory)软件来模拟Monte-Carlo大气扩散过程。该模式与气象文件数据/地形数据的结合、完善的模式考虑以及良好的操作界面正越来越被广泛的应用于各个领域,是目前比较先进的大气扩散模拟软件。

2.2 随机数产生

Monte-Carlo模拟的最重要的思想是通过计算机的产生特定分布的随机数来模拟现实世界中的随机现象。事实上,随机数在Monte-Carlo方法中居中心地位, Monte-Carlo法模拟的结果正确与否,关键在于随机数的质量[7]。

以往计算机最常用的产生均匀分布的随机数的算法是线性同余法,然而在大多数情况下,计算机是不能生成真正的随机数的。事实上,计算机能做的只是产生一些看上去很像完全随机的数的数字,这些数字称为伪随机数。理论上说,完全通过数学计算的方式是不可能得到理想的随机数,必须要依靠独立的外部随机事件。

HAVEGE(HArdware Volatile Entropy Gathering and Expansion)通过收集计算的相关硬件的信息熵(鼠标、键盘、硬盘、内存、网络等),这些硬件的状态信息在计算机运行的过程中是随时变动,它们成为HAVEGE的随机源。这些信息量十分巨大,至少达到100Kbits/S,而由此产生的随机数可达100Mbits/S[8]。

高斯分布随机数的产生利用Box-Muller转换,通过HAVEGE算法选取两个均匀分布的随机数v1、v2,利用Box-Muller转换式(4)可得到两个符合高斯正态分布的随机数[9]。

undefined

图1为使用高斯型随机数序列填充的结果,可见点集中分布在均值0附近, 2倍方差内集中了大部分点,3倍方差内集中了几乎所有的点,因此HAVEGE方法产生的高斯分布随机数可以达到Monte-Carlo模拟的要求。

2.3 气象条件

气象条件是气体扩散过程模拟中非常重要的因素。扩散模型一般要利用的气象条件有风速、风向、温度等。这些数据的获得途径有两种,一种是通过国家气象局和当地的气象部门通过网络实时获取。另一种方式是模拟预测。就是前文提到的区域大气预测模型(RAMS),它包含了多种物理过程及其参数化方法。如湍流混合、太阳和地球辐射、各种湿物理过程(可以描述云的形成及其与液、固态降水物质的相互作)、大气与多层土壤模式、地表植被和地表水的感热、潜热交换、地形动力作用以及积云对流参数化等。Monte-Carlo随机扩散模型利用RAMS所预测的风场(速度和风向),模拟气体扩散浓度分布。

3 M-C模型计算分析

影响气体扩散的因素很多,包括地形、气象、重力、气体性质、初始条件等,长输管道气体扩散模型计算分析主要包括在复杂条件下的扩散过程模拟,针对应急决策过程的需求特点,重点分析复杂风场和地形中管道气体的扩散机理。

3.1 风场影响

风场的主要是风速和风向,风场对天然气的扩散有非常明显的输送和稀释作用,因此,泄漏物质总是分布在泄漏点下风方向。风速对天然气扩散的影响是复杂的,不同高度的风速是不断变化的,风速的影响会加剧空气和天然气之间的传热和传质,使得天然气的扩散加剧,风速对扩散气团的迎风面和背风面的影响也不一样。无风时,扩散以泄漏点为中心,向各方向均匀扩散,质量输送以扩散为主。在有风时,主风向的平流输送作用占主导地位。一方面由于风对天然气气团的平流输送作用加剧,使得天然气气团有往下风向输送的趋势,风速越大,输送作用越显著,由于气流卷吸混合作用加强,造成下风向处的气体浓度降低。另一方面由于风速增大,引起脉动速度增大,紊流运动加剧,紊流扩散作用增大,使天然气团浓度下降,同时紊流运动的加剧也使得天然气团与周围环境的热交换变得剧烈,使扩散的过冷气体温度迅速上升,天然气团密度下降,在风的作用下更容易扩散,从而导致天然气团浓度下降。

以安阳-洛阳天然气某管道段为例,天然气的主要成分为烃类混合物,不含硫化物。泄漏事故区域初始风场分布如图2,受地形、大气湍流等因素的影响,风场(风速、风向)分布并不均匀。图4是发生泄漏事故后一小时气体浓度分布图。从开始泄漏到事后一小时,区域内风场并无多大变化,所以预测模型与ALOHA模型基本一致。泄漏事故后生7小时后,RAMS预测风场如图3。此时风向与开始泄漏初期状况相比发生了很大变化。受此影响,12小时后气体浓度分布由于受湍流、地形等因素的影响,气体浓度开始不连续,气体浓度不断稀释,在边界部分出现的浓度随机集中分布,如图5。

Monte-Carlo方法还可对粒子轨迹进行跟踪模拟。以忠武线某段含H2S天然气长输管道为例,图6模拟了事故区域内1小时浓度分布。其中所用粒子数3000,随机行走1小时。浓度分布(图6a)取决于粒子分布(图6b)。

3.2 地形条件影响

地形和下垫面的非均匀性,对气流运动和气象条件产生动力和热力的影响,尤其是当下垫面不均匀时加剧湍流运动,如山谷风、过山气流(迎坡风抬升,下坡风下沉)、海陆风、城市热岛效应等,均会改变气体的扩散条件。

由于长输管道跨越空间尺度大,一般都在上千公里的范围,所以管道经常会经过地形变化很大的区域。变化较大的地形会对气体的扩散范围、空间浓度分布、扩散时间等数据产生很大的影响,达不到应急决策的精度要求,造成决策指挥的失误。用ALOHA模型对穿越这些区域的管道泄漏进行模拟时,可能会忽视这些地形的影响,而Monte Carlo粒子随机游走模型充分结合了三维地形的特点,可满足应急决策的数据精度要求。

图7是模拟某管道泄漏事故点区域内风场的近地面分布图,受地形的影响,近地面风场会发生变化。图中可以看出,该区域内山脊地区风速比山谷要高,同时受山体坡面法向的影响,风速方向也发生了改变。Monte Carlo模型是通过计算粒子在单位时间间隔内三维空间的位移来模拟气体的扩散的,因此,受地形起伏的影响,粒子在空间运动的过程中会感受到这种地形的“强迫作用”。图8是该区域内泄漏点地面高度10m粒子的空间分布图。从图可看出,粒子在源点附近受地形影响较小,但随着扩散进行,近地面粒子受“地形强迫”作用“滞留”在与源点较短的范围内,这部分粒子占总粒子数的大部分;只有那些高空粒子受地形影响小,能扩散到较运距离,到达扩散区域的前锋。受地形影响,危险区域并不连续,主要集中在是山谷地区和泄漏点附近。

4 结论

气体扩散模拟受气象条件、环境地形等因素的影响,对于天然气管道事故应急决策而言,影响最大的是风场与地形的影响。风场对天然气的扩散有非常明显的输送和稀释作用,并且大气的湍流作用也会对气体扩散产生影响;复杂地形的“强迫作用”对扩散过程也产生重大影响。高斯扩散模型是最常用、最简单的扩散模式,但无法满足了事故灾害应急要求。Monte-Carlo模式通过大量标记粒子来描述在大气中的扩散运动,真实直观地反映了污染物在大气运动中的输送规律。该模式既体现了扩散本身具有的随机特性,而且对于复杂的地形模拟计算有很高的精度。此外,Monte-Carlo模型很容易与RAMS模型整合,预测复杂气象条件下的天然气扩散浓度分布,为应急决策提供了重要依据。

参考文献

[1]M.El-harbawi.Rapid analysis of risk assessment usingdeveloped simulation of chemical industrial accidentssoftware package.Int.J.Environ.Sci.Tech.,2008,5(1):53-64

[2]B.A.Kunkel.User's guide for the Air Force Toxic Chem-ical Dispersion Model(AFTOX).Interim report,Octo-ber 1985-December 1987.1988

[3]G.Ooms,A.Mahieu and F.Zelis.The Plume Path ofVented Gases Heavier than Air.1974

[4]D.L.Ermak.User's manual for SLAB:An atmosphericdispersion model for denser-than-air releases.LawrenceLivermore National Laboratory,1990

[5]P.A.Durbin.Stochastic differential equations and turbu-lent dispersion.NASA STI/Recon Technical Report N,1983,83:22546

[6]C.-K.S.Cheol-Hee Kim.Simulating mesoscale transportand diffusion of radioactive noble gases using the lagrang-ian particle dispersion model.Journal of EnvironmentalRadioactivity,2008,99(10):8

[7]金畅.蒙特卡洛方法中随机数发生器和随机抽样方法的研究[D].2005

[8]A.Seznec and N.Sendrier.HArdware Volatile EntropyGathering and Expansion:generating unpredictable ran-dom number at user level.2002

Carlo方法篇5

发展了一种粘弹性Monte-Carlo随机有限元法,并对固体火箭发动机药柱进行了三维随机分析.首先由Herrmann泛函出发导出了三维不可压和近似不可压粘弹性增量有限元列式.然后,同时考虑固体推进剂力学性能参数和载荷的随机性,结合Monte-Carlo技术对弹性约束的圆柱形中孔药柱进行了随机分析.为了克服直接Monte-Carlo法(DMCS)效率低的缺点,引入Latin超立方抽样技术(LHS)以提高计算效率.算例表明该方法收敛速度块、通用性强、易于工程应用.

作者：田四朋唐国金李道奎雷勇军 TIAN Si-peng TANG Guo-jin LI Dao-kui LEI Yong--jun 作者单位：国防科技大学航天与材料工程学院,长沙,410073 刊名：强度与环境 ISTIC英文刊名：STRUCTURE & ENVIRONMENT ENGINEERING 年，卷(期)：2007 34(2) 分类号：V435.13 关键词：不可压 Latin超立方抽样随机分析 Monte-Carlo模拟

★ 固体火箭发动机试验综合软件系统设计

★ 固体火箭发动机燃烧室三维流动数值计算

★ 弯桥结构动力性能的有限元计算方法

★ 六武高速公路装配式涵洞结构有限元优化分析

Carlo方法篇6

本文将ARIMA模型和Monte Carlo方法相结合, 建立一套简单而又可靠的区域煤炭消费量及其碳排放量预测方法。首先基于历年能源消费量时间序列数据建立能源消费量的ARIMA模型, 然后将通过验证后的ARIMA模型预测未来的能源消费量, 最后利用MonteeCarlo方法估算未来煤炭消费所产生的碳排放量。此外, 本文以广西煤炭消费量及其碳排放量预测为例开展了实例研究。

1 研究方法

ARIMA模型一般建模步骤包括数据序列的平稳性检验、差分处理求参数 (d) 、根据样本数据来估计模型适当的阶数p和q、模型的验证等。在本文中, ARIMA过程的处理在R语言Forecast软件包中进行[1,2]。

考虑到能源消耗预测结果以及碳排放量参数都存在不确定性, 本文采用Monte Carlo技术进行模拟。煤炭消费量的分布假设为均匀分布, 最大值和最小值分别为其预测值80%置信区间的上限和下限。单位标准煤碳量和二氧化碳排放量假设为均匀分布, 取值范围分别为0.69~0.71和2.54~2.59[3]。通过10 000次模拟, 获得碳排放量和二氧化碳排放量的平均值和80%置信区间。

2 实例研究

2.1 数据来源

本文实证研究的数据源自《广西统计年鉴2014》。自2003年以来, 广西煤炭消费总量迅速增长, 由1 770万吨标准煤增长到2013年的5 360万吨标准煤, 年均增长率达到18.3%, 在此期间年均增长率达20.3%。煤炭消费量占能源消费总量的比例比较在50%~60%之间, 截至2013年仍达55%。

2.2 ARIMA模型的验证

用1978年~2008年共31年的数据作为训练样本进行建模, 用2009年~2013年共5年的数据作为测试样本进行模型验证。根据AIC和BIC准则, 利用训练样本建立了煤炭消费量的ARIMA模型, 为ARIMA (1, 2, 0) 。然后, 用该模型预测了2009~2013年煤炭消费量 (表1) 。煤炭消费量预测值与实际值的相对误差在6.0%以内, 平均值为3.9±1.8%, 实际值都落在预测值80%置信区间的上下限以内。可见, ARIMA模型是可以满足中短期煤炭需求量预测需要的。

2.3 十三五期间广西煤炭消费量和碳排放量预测结果

时间序列模型的特点是利用新的数据动态地训练和校准。为预测十三五期间 (2016年~2020年) 的煤炭消费量, 以1978年~2013年共36年的煤炭消费数据作为建模数据, 根据AIC和BIC准则, 建立煤炭消费量的ARIMA模型, 为ARIMA (2, 2, 0) , 再利用该模型预测未来5年的煤炭消费量 (表2) 。十三五期间, 广西煤炭消费量年均增长率为6.1%, 到2020年煤炭消费量预测值的平均值为8050万吨, 80%置信区间为6735万吨~9365万吨。Monte Carlo模拟结果表明, 碳排放量和二氧化碳排放量预测值的平均值分别为5738万吨和21029万吨, 预测值的80%置信区间分别为4650万吨~6836万吨和17121万吨~24925万吨 (表2) 。对预测结果的不确定性进行认识和评估[4], 有助于提高投资决策和环境政策决策的可靠性。

3 结论

ARIMA模型和Monte Carlo方法分别是常用的时间序列分析和不确定性分析方法, 都可以在R语言等平台方便操作, 将二者相结合建立了一套简单而又可靠的区域煤炭消费量及其碳排放量预测方法。实例研究结果表明, ARIMA模型可以满足中短期煤炭需求量预测需要, Monte Carlo模拟技术提高了投资和环保决策的可靠性, 为能源消费及其碳排放量预测提供了有效手段。

参考文献

[1]吴喜之, 刘庙.应用时间序列分析[M].北京:机械工业出版社, 2014.

[2]Hyndman RJ, Khandakar Y.Automatic time series forecasting:The forecast package for R.Journal of Statistical Software, 2008, 27 (3) :1-22.

[3]涂华, 刘翠杰.标准煤二氧化碳排放的计算[J].煤质技术, 2014, 2:57-60.

Carlo方法篇7

20世纪80年代以来, 我国在水利水电工程施工中应用系统工程的方法进行施工系统分析和仿真模拟, 取得了显著成效, 计算机仿真模拟技术逐渐成为大型水电工程规划、设计、施工不可缺少的手段[1,2]。对堆石坝施工仿真的研究方面, 以往建立的仿真模型侧重于对施工全过程的静态模拟, 根据给定的施工组织设计方案寻求合理机械配置以及进行工期和工程量的计算和论证。但是, 对于大中型水利水电工程来说, 不确定因素众多, 自然条件和技术条件复杂, 使得原定的施工进度计划与实际施工过程产生较大的偏差, 以致这种静态仿真模拟失去了对实际施工的指导作用。为此, 施工过程实时仿真模拟逐渐成为了施工仿真中的新趋势。在堆石坝施工过程的实时仿真领域中, 已经有学者取得了一些初步的成果。天津大学运用自适应控制理论与方法, 建立了高堆石坝施工过程管理与控制模型;此外, 有学者研究了基于Petri网的施工过程的建模以及流程的仿真分析, 用层次化的方法建立了整体的施工过程模型, 提出了一个施工过程管理的框架;还有一些学者则针对堆石坝施工过程中的某一个方面建立了动态模型, 例如料物调配模型、交通运输模型等[3]。对施工工期的研究方面, 由于工期和成本是施工项目管理的两个重要目标, 这两个目标的控制并不是孤立的, 而是相互联系和相互制约的[4]。所以传统的工期计算一般都是利用网络计算方法, 而网络计算方法主要有:关键线路法 (CMP) 、计划评审技术 (PERT) 、决策网络技术、图示评审技术 (GERT) 等。随着计算机以及其他一些辅助工具的产生, 现在对于施工工期的计算方法又有了新的发展。四川大学涂扬举等把加速遗传算法应用于堆石坝的施工工期研究, 提出了基于加速遗传算法的堆石坝施工优化[5]。此外, 浙江大学熊鹰等也将蚁群算法运用到施工项目工期-成本优化问题中, 并且取得了很好的效果。但是, 对堆石坝有效施工工期的研究至今还没有一个确定的方法。本文提出了利用Monte-Carlo方法来模拟降雨条件下堆石坝有效施工工期, 首先分析堆石坝有效施工工期的影响因素, 然后建立基于Monte-Carlo模拟降雨条件下堆石坝有效施工工期的统计分析模型, 最后将此方法运用于工程实际中并进行归纳总结。

1 堆石坝有效施工工期的影响因素

1.1 堆石坝有效施工工期的影响因素

水利水电施工进度计划具有不确定性、随机性, 这种不确定性、随机性可归纳为主观和客观两个方面[6]。主观不确定性主要包括完成某工序所需实际施工天数估计的不确定性。客观不确定性主要是由水文、气象等因素引起的, 主要包括每月有效施工天数的不确定性和规定 (或计划) 完工日期的不确定性。对于水利水电工程施工, 一般而言, 采用轮休的办法可不考虑法定假日的停工, 但应考虑到水文和气象 (如降雨) 的影响而引起的停工。例如, 土石坝的填筑施工, 特别是黏土防渗体的施工, 降雨对其影响很大, 即每月一定标准雨日的多少, 直接影响到有效施工天数。归纳起来, 影响堆石坝施工工期的因素主要有以下几点。

(1) 技术因素。

工程规模大, 技术复杂, 无论是施工技术还是安装技术都有一定难度。例如, 堆石坝工程坝面作业, 工作面狭窄、工种多、工序多、机械设备多, 施工时须有妥善的施工组织规划, 为避免坝面施工中的干扰, 延误施工进度, 坝面压实宜采用流水作业施工;钢筋混凝土面板是面板坝的主要防渗结构, 厚度薄, 面积大, 在满足抗渗性和耐久性条件下, 要求具有一定柔性, 以适应堆石体的变形;为了避免混凝土坝不均匀沉降引起坝体开裂, 常用结构缝将坝体分段, 为了控温防裂和施工方便, 用纵缝将坝段分成若干柱状快, 在浇筑时再用临时的施工缝将柱状块分层, 以保证浇筑质量;碾压混凝土坝通常在上游面设置常态混凝土防渗层以防止内部碾压混凝土的层间渗透, 有防冻要求的坝, 下游面亦用常态混凝土, 为提高溢流面的抗冲耐磨性能, 一般也采用标号较高的抗冲耐磨常态混凝土。

(2) 进度计划安排的合理性。

由于工程初步设计阶段, 一些工作还不能深入, 工程量还不能详细计算, 故有些工作时间的估计和确定还不够准确, 各种施工方案还有待于进一步优选。

(3) 资源的充分性。

由于工程量大, 投资高, 工期长, 能否保证资源及资金的及时供应和合理使用, 也是影响工期的关键性[7]。

(4) 自然条件。

工程的地质、水文地质、气象、河流水文特征等存在许多不确定因素。例如, 在水利水电工程控制性进度计划中, 由于水文方面因素所规定的一些进度计划中的控制点, 如截流日期、拦洪日期、封堵日期等均有不同程度的不确定性, 这些控制点的不确定性对施工进度风险有直接的影响。

(5) 施工队伍及合同条件。

如何选择施工队伍和加强工程建设管理以及合同形式和条件等, 也是保证工期的重要因素[8]。

1.2 工期模型

当规定工期为Ts时, 则完工的可能性为:

$Ρ (Τ \leq Τ_{s}) = \int_{- \infty}^{Τ_{s}} f (t) d t (1)$

式中, f (t) 为工作历时为服从某一已定的概率分布 (例如β分布、均匀分布、三角形分布或其他连续分布, 也可以是离散型分布) 的随机变量。

设总工期由n个关键节点工期组成, 即T={T1, T2, …, Tn}对于规定工期Ts与其相应的完工概率为P={P1, P2, …, Pn}则完工的可能性为:

$Ρ (Τ \leq Τ_{s}) = \sum_{\begin{array}{l} Τ \leq Τ_{s} \\ i = 1 \end{array}}^{n} Ρ_{i} (2)$

不能按期完成的可能性:

$Ρ_{f} = 1 - Ρ (Τ \leq Τ_{s}) (3)$

2基于Monte-Carlo模拟降雨条件下的统计分析模型

2.1 模型假定

(1) 降雨是随机均匀分布的, 且有一定概率, 未来降雨天数与历史统计对应月的平均降雨天数应相等。

(2) 月内日降雨相互独立, 不同日的降雨互不影响, 为均匀随同分布, 某日降不同雨量的概率也相互独立。

(3) 不影响施工的降雨不考虑。

(4) 对大于半天的降雨时段作为降雨1 d考虑, 小于0.5 d的降雨时段按不降雨考虑。

2.2利用Monte-Carlo模拟仿真评估降雨条件下的有效施工天数的模型

2.2.1 搜集整理降雨资料

利用坝址地区水文站的累年各月各要素统计表, 查出一年中的每日降雨量, 然后对照《水利水电工程施工组织设计手册》中对坝体填筑停工天数标准的规定, 统计出一年中各月因雨停工的天数, 并将降雨停工天数按月制成施工停工天数柱状图, 用曲线拟合图上的数据, 拟合曲线即为每月停工天数的概率密度曲线。

2.2.2 推求正态分布N (μ, σ2) 的抽样公式

区间[a, b]范围内、期望值为μ、方差为σ2的正态分布, 记作N (μ, σ2) [a, b], 其概率密度:

$f (x) = \frac{1}{\sqrt{2 π σ}} e^{- (x - μ)^{2} / (2 σ^{2})}, - \infty ＜ x ＜ \infty (4)$

令Z= (x-μ) /σ, 则将X～N (μ, σ2) 转换成为标准正态分布变量Z～N (0, 1) 。取[0, 1]区间均匀分布随机量r1, r2, 利用二元函数变换得:

$Ζ_{1} = (- 2 l n r_{1})^{1 / 2} \cos 2 π r_{1}, Ζ_{2} = (- 2 l n r_{2})^{1 / 2} \sin 2 π r_{2} (5)$

Z1, Z2即为两个独立的标准正态分布随机变量。再由线性变换可以得到正态分布随机变量u1, u2:u1=Z1σ+μ;u2=Z2σ+μ。

$\begin{array}{l} 抽样公式即为 ∶ x_{1} = μ + σ \sqrt{- 2 l n r_{1}} \cos 2 π r_{1} (6) \\ x_{2} = μ + σ \sqrt{- 2 l n r_{2}} \sin 2 π r_{2} \end{array}$

2.2.3 用Monte-Carlo方法进行模拟仿真

采用Visual C++编程生成n个随机数, 然后对每个随机数利用抽样公式计算出每月停工天数, 将一年中12个月的每月停工天数累加起来得到一年停工天数, 由于是利用了n个随机数, 那么模拟的一年停工天数也有n组Ti (i=1, 2, …, n) 。

对一年停工天数样本计算特征值, 如均值 $\bar{Τ}$ 与方差s2。

一年停工天数样本均值:

$\bar{Τ} = \frac{1}{n} \sum_{i = 1}^{n} Τ_{i} (7)$

一年停工天数样本方差:

$s^{2} = \frac{1}{n - 1} \sum_{i = 1}^{n} (Τ_{i} - \bar{Τ})^{2} (8)$

首先, 要确定仿真运行次数n, 理论上讲仿真次数越多越好, 究竟仿真多少次才能获得一年停工天数的较好估计, 这里通过一年停工天数的样本均值 $\bar{Τ}$ 和总体均值之间的误差范围来确定。

设μ, σ2为一年停工天数的总体均值和方差, 由于T1, T2, …, Tn独立同分布, 且样本方差s2是总体方差σ2的一致估计量, 根据中心极限定理, 可知n充分大时:

$z = \frac{\bar{Τ} - μ}{\sqrt{\frac{s^{2}}{n}}} ~ Ν (0, 1) (9)$

若服从正态分布的随机变量x的概率分布函数是ϕ (x) ,

$Ρ (| \frac{\bar{Τ} - μ}{\sqrt{\frac{s^{2}}{n}}} | \leq U) = ϕ (U) - φ (- U) = β (10)$

取 $ε = U \frac{s}{\sqrt{n}}$ , 对上式进行变换, 则有:

$Ρ (| \bar{Τ} - μ | \leq ε) = β (11)$

这就是以一年停工天数的样本均值 $\bar{Τ}$ 代替总体均值μ, 在绝对误差小于ε, 在置信度为β的情况下的误差表达式。

因此, 仿真次数n的确定步骤如下:

(1) 给定置信度β, 由 $ϕ (U) = \frac{β + 1}{2}$ , 查标准正态分布表, 得U值。

(2) 给定绝对误差ε。

(3) 先做若干次仿真, 计算样本方差s2。

(4) 将ε、U、s2代入 $n = \frac{U^{2}}{ε^{2}} s^{2}$ , 便可得出满足精度要求的最小仿真运行次数。

2.2.4 绘制完工概率曲线

仿真次数满足精度后, 然后拟合出一年中停工天数的经验分布曲线和经验累计曲线, 以了解其分布规律, 具体步骤如下:

(1) 确定分组区间总范围为[Tmin, Tmax] , 其中Tmin, Tmax分别为N次仿真中一年停工天数的最小值、最大值。

(2) 根据分组区间总范围可先指定分组数然后确定组距, 或是先指定组距再确定分组数, 以便于实际分析为原则。分组数与组距满足如下关系式:

$Τ_{g} = \frac{Τ_{\max} - Τ_{\min}}{l} (12)$

式中:l为分组数;Tg为组距。

(3) 进行分组。令Ti=Tmin+Tg·i (i=1, 2, …, l) , 则分组结果为[Tmin, T1) , [T1, T2) , [T2, T3) , …, [Tl-1, Tl) 。

(4) 将仿真结果按区间统计的频数, 得一年停工天数频率直方图, 见图1。

(5) 对统计频数按区间累加, 得到一年停工天数的累计直方图, 见图2。

(6) 对一年停工天数的累计直方图进行曲线拟合, 得到一年停工天数的概率曲线, 见图3。

2.2.5 风险分析

(1) 输入指定某一年停工天数, 查询相应概率。

$p = f (t) (13)$

式中:p为停工天数的概率;t为停工天数。

(2) 给定某一概率, 查询相应一年停工天数。

$t = g (p) (14)$

3 工程实例分析

3.1 工程概况

糯扎渡水电站为云南省澜沧江中下游河段两库八级水电规划的第五级, 心墙堆石坝最大坝高261.5 m。枢纽位于云南省思茅境内。

糯扎渡水电站枢纽建筑物由心墙堆石坝、左岸开敞式溢洪道、左岸泄洪隧洞、右岸泄洪隧洞、左岸地下引水发电系统及地面500 kV开关站等建筑物组成。枢纽主要工程量为土石方明挖5 224万m3, 石方洞挖627.72万m3, 土石方填筑3 591.66万m3, 混凝土浇筑432.79万m3。

坝址地区1961～1990年各月日降雨量统计如表1。

3.2 降雨条件下的有效施工天数的模拟

(1) 按月停工天数分布。

查阅《水利水电工程施工组织设计手册》中对坝体填筑停工天数标准的规定, 并结合工程实际, 在此工程中规定:日降雨量≥10 mm时就停工。因此, 只研究表1中对≥10 mm日数的统计。

(2) Monte-Carlo模拟。

由Visual C++编程生成25个随机数, 然后利用抽样公式计算每个随机数对应的每月停工天数, 将一年中12个月的每月停工天数累加起来得到25组一年停工天数, 其计算值略。

(3) 确定仿真次数N。

设置信度β=0.96, 查标准正态分布表, 得U=2.06;给定误差ε=1 d, 先做25次仿真, 计算样本方差s2=5.026 998, 由公式 $n = \frac{U^{2}}{ε^{2}} s^{2}$ 可确定n=22次, 可看出仿真22次, 即可满足精度要求。为统计方便, 取n=25次。也就是利用上面先做的仿真结果。

(4) 绘制一年停工天数的概率曲线。

通过对Monte Carlo仿真模型的25次仿真, 得到了25个一年停工天数, 其最小值为15 d, 最大值为22 d, 均值为19 d, 均方差为2。本例中停工天数的分组组距取为1 d, 一年停工天数概率曲线如图4。这样, 若要知道某一概率所对应的一年停工天数, 只需在图5的曲线中找到相应的点即可。

4 结语

本文根据堆石坝的施工特点, 分析了堆石坝有效施工工期的影响因素并建立了工期模型。然后, 基于Monte-Carlo方法又建立了模拟降雨条件下堆石坝有效施工工期的统计分析模型, 对堆石坝有效施工工期进行了分析模拟研究。并将此模型应用于工程实践中, 为工程管理人员提供了简便的有效施工天数预测方法。需要特别指出的是, 本文是建立在月平均停工天数呈正态分布的某一特定工程基础上的, 具有一定的特殊性, 使用者可以依照本文的方法模型, 针对具体问题选择不同的概率分布类型, 从而对有效施工天数做出更为准确、更为合理的评估。

参考文献

[1]孙锡衡, 齐东海.水利水电工程施工计算机模拟与程序设计[M].北京:中国水利水电出版社, 1997.

[2]王仁超, 胡颖慧.复合坝型坝体施工单一仿真模型研究[J].中国农村水利水电, 2007, (11) :100-103.

[3]刘珊珊, 周宜红, 刘全, 等.堆石坝施工的实时动态仿真系统研究[J].系统仿真学报, 2004, 16 (11) :2 525-2 528.

[4]熊鹰, 匡亚萍.施工项目工期-成本优化问题的蚁群算法[J].浙江大学学报 (工学版) , 2007, 41 (1) :176-180.

[5]涂扬举, 马光文, 陶春华, 等.基于加速遗传算法的堆石坝施工优化[J].水力发电, 2007, 33 (1) :26-28.

[6]王卓甫, 陈登星.水利水电施工进度计划的风险分析[J].河海大学学报, 1999, 27 (4) :84-87.

[7]陈志强, 侍克斌.水利枢纽工程坝型选择影响因素的分析[J].中国农村水利水电, 2007, (8) :62-66.

【Carlo方法】推荐阅读：

写作方法(叙述方法)09-14

科学方法与创新方法07-01

安全评价方法：PsA分析方法07-06