数值积分范文
时间：2023-04-02 19:02:20
导语：如何才能写好一篇数值积分，这就需要搜集整理更多的资料和文献，欢迎阅读由公务员之家整理的十篇范文，供你借鉴。
篇1
关键词：高斯数值积分方法；积分变量；积分区间；非圆弧拱；多拱梁法
1问题的提出
高斯数值积分方法是一种节点很少、精确度很高的方法，它的特点是节点不等距，计算精度很高，一般利用正交多项式的有关关系式来确定其节点位置和系数.当节点为n时，其代数准确度可达2n-1次.如节点数为3,则求积公式对于任意5(=2×3-1)次多项式都是准确的，这样的精度完全可应用于拱坝程序中的拱段及梁段的计算.笔者曾在圆弧拱多拱梁法程序中，广泛采用3节点高斯数值积分，取得成功[1].但将这一方法推广应用于非圆弧拱多拱梁法程序，却有一定的难度.
由计算数学可知[2]，一般定积分式与高斯积分式的变换形式为：
式中:l=(b-a)/2，为积分限变换系数；n为高斯节点数；ξi为高斯积分节点坐标；gi为与ξi对应的高斯积分系数.
如所周知，拱圈形、载常数计算公式为：
(1)
(2)
式中:S为弧长.
当采用高斯积分公式且为等截面圆拱时，上式相应改为：
(3)
(4)
式中:φ为与弧长S对应的中心角；r为中心半径；E为坝体弹模；Ii,MLi为与高斯节点i对应的截面惯性矩及静定力矩.从上式可见，由于等截面圆拱r为定值，存在简单的dS=rdφ,S=rφ的关系，所以积分变量由S变为φ，积分区间由弧段0～S变为中心角0～φ.在等截面圆拱计算中，许多参数(如坐标及静定力系等)可直接由中心角用显式求得.积分变量由S改为φ后，可大大简化计算.
然而，对于非圆弧拱，问题要复杂得多.因为非圆弧拱的曲率半径和曲率中心处处都在变化，往往不能用简单的显式来表示某一拱段中心角与弧长的关系.并且某些曲线计算弧长也很麻烦，甚至不易用显式求得.如何采用适宜的积分变量和积分区间，是迫切需要解决的问题.
2不同曲线积分变量区间的选取
2.1基本资料5种非圆弧拱示意见图1.几种非圆弧拱曲线方程及有关公式见表1.
图15种非圆弧拱平面示意
表1几种非圆弧拱曲线公式
类型抛物线对数螺旋线椭圆双曲线
曲线方程y=x2/2Rρ=ρ0eaψ
a=cosβ
曲率半径r=(R2+x2)3/2/R2r=Reaψ
r=a2b2x
[(a-y)2/a4+x2/b4]3/2r=a2b2x
[(a+y)2/a4+x2/b4]3/2
弧长S=1/2{xA/cosψA+RLn[(sinψA+1)/cosψA]}S=R/a(eaψA-1)
dSReaψdψ
备注R为拱冠曲率半径β为切线角,R为拱冠曲率半径a为长轴之半,b为短轴之半a为实轴之半,b为虚轴之半
2.2弧的微分公式表1所列弧的微分(dS)算式，除了对数螺旋线稍简单外，其余3种曲线的算式都比较复杂.若从曲线的一般性质来看，当曲线方程y=f(x)时，则弧的微分为：
式中:y′=tgφ，为函数y=f(x)在点x的导数，φ为过点(x,y)的切线与x轴的交角，不难看出，此角与该点的中心角相等.将y′=tgφ代入上式：
所以dS=dx/cosψ(5)
将式(1)作为抛物线、椭圆、双曲线3种曲线弧的微分一般公式，比表1中所列dS算式要简捷得多.式(1)也适用于其他一阶导数存在的任何曲线.
2.3积分变量区间的选取根据各类曲线的性质，选取3种积分变量区间，分述于后.
(1)对于抛物线、椭圆、双曲线3种曲线，采用式(1)所列弧的微分一般公式.此时积分变量为x，积分区间为与S相应的x变化区间.高斯积分时各项均应乘以1/cosφi,φi为与节点i对应的中心角.如A1算式改为:
(6)
(2)对于对数螺旋线，采用表1所列算式：dS=Reaφdφ.此时积分变量为φ，积分区间为与S相应的φ的变化区间，各高斯积分项i均应乘以Reaφi，如A1算式为：
(7)
上式除了高斯积分项乘以eaφi外，其余与圆弧拱相似.
图2五心拱平面示意
(3)对于五心拱，其中弧段为圆弧，算法与等截面圆拱相同，当其边弧段为变截面时，则中心拱弧线为非圆弧曲线，且不便用显式表达.仔细考察该段曲线，发现它与以R3=(RM+RD)/2为半径的圆弧很相近(RM,RD分别为边弧外、内半径)，此圆弧中心O3在边弧起始截面外弧中心O1与内弧中心O2联线中点处，如图2所示.
边弧段计算时，积分变量为φ，积分区间为与边弧S相应的φ的变化区间φ3，如A1算式为：
(8)
上式形式上与圆弧拱一样.
3算例
以上述3种积分变量区间的选取方式，计算各类曲线的半拱弧长，举例于下.
设采用3节点高斯数值积分方法，节点坐标及高斯积分系数列于表2.如以中心角φ为积分变量，以Δφ为积分区间，则与节点i相对应的φi=Δφ(1+ξi)/2.又如以水平坐标x为积分变量，以Δx为积分区间，则与节点i相对应的xi=Δx(1+ξi)/2.
表2节点坐标及高斯积分系数
节点号i节点坐标ξi高斯积分系数gi
1-0.7745966910.555555582
200.888888896
30.745966910.555555582
3.1求抛物线、椭圆、双曲线拱半拱弧长例1：设抛物线拱拱冠曲率半径R=140m，拱端中心角φA=45.32°，拱端坐标xA=141.5726m，求半拱弧长S，可以采用两种方法.一种是按表1所列S的公式直接计算，这是理论积分后的公式，是精确的.另一种是采用高斯数值积分方法，以x为积分变量，xA为积分区间.
(1)按理论公式
S=1/2{xA/cosψA+RLn[(sinψ\-A+1)/coxψA]}
算得S=162.921371m.
(2)按高斯数值积分方法，应有
算得gi/cosφi=2.301610646,S=(141.5726×2.301610646)/2=162.922485m.
该数值积分值与理论计算值162.921371m相比，仅相差0.001114m，相对误差仅为6.8×10-6.
例2：设椭圆拱拱冠曲率半径R=164.51m，拱端中心角φA=51°，坐标xA=141.5843m，长轴之半236.9m，短轴之半197.42m，求半拱弧长.例3：设双曲线拱拱冠曲率半径R=152.71m，拱端中心角φA=41°，坐标xA=141.58421m，实轴之半954.41m，虚轴之半381.77m，求半拱弧长.
上两例因两种曲线无理论积分公式直接用显式计算弧长S，只能与表1中dS算式的高斯数值积分值相比.用dS算式直接数值积分时，积分变量、积分区间与前述方法一样，仍为x及xA，但每一高斯积分项不用除以cosφi，而是乘以dx前的算式等，可见后一算法较繁.两例成果S及比较见表3.
表3椭圆拱、双曲线拱计算成果比较
类别例2椭圆拱例3双曲线拱
原dS算式成果164.143661158.610046
dS=dx/cosψ成果164.143737158.610031
两种算法差值/m0.0000760.000015
相对误差4.63×10-79.457×10-8
由表3可见两种算法成果非常接近.
3.2求对数螺旋线拱半拱弧长例4：设对数螺旋线拱拱冠曲率半径R=154m，拱端中心角φA=46°，切线角β=52°，a=ctgβ=0.781285626,以中心角φ为积分变量，积分区间为拱端中心角φA，求半拱弧长S.
(1)按理论公式计算
S=R/a(eaψA-1)=171.9726625m
(2)按高斯数值积分计算
两种算法成果差值为0.0000061m，相对误差仅为3.547×10-8.
3.3用近似方法求五心拱边例5：设五心拱半拱边弧夹角10°，外半径RM=290m，内半径RD=193.514m，边弧起点拱厚6.466m，拱端厚度8.466m，求半拱边弧长S.
(1)用较精确的计算公式S=φ3×(RA+2×R3)/3.式中3=(RM+RD)/2;φ3为以O3为近似中心的边弧夹角，以弧度计；RA为边弧拱端点与O3联线的长度，见图2.算得RA=241.56728m，φ3=0.206893551弧度，R3=241.757m，从而算出S=50.00488034m.
(2)用近似计算公式S＝φ3×R3=50.01796429m.
两种算法所得的差值为0.013m，相对误差为2.6×10-4.
由以上成果可知，以R3为半径，以φ3为中心角所得的圆弧与边弧的中心弧很近似，因此在高斯数值积分计算时，可以采用较简单的积分变量φ及积分区间φ3.
4结语
综上可知，高斯数值积分应用于非圆弧拱时，需针对非圆弧拱曲线的性质和特点，选取适宜的积分变量与相应的积分区间，这样常可收到事半功倍之效.如对于抛物线拱、椭圆拱、双曲线拱，采用dS=dx/cosφ的通式，既避免了各种曲线的繁复计算，也便于程序规格化.对于对数螺旋线拱，则利用对数函数微分、积分都简单以及极坐标方程弧的微分的特点，积分变量选取φ而不选取x，既大大简化了计算，也很容易利用圆弧拱的算法稍加变换.对于五心拱，在控制误差足够小的前提下，边弧线采用近似圆弧，大大简化了计算.
将3节点高斯数值积分方法推广应用于各种非圆弧拱坝多拱梁法程序，对于提高拱坝程序的计算速度和精度，有很大价值；对于将来推广应用高斯数值积分于各类复杂结构的分析、计算，也有一定的启发和借鉴作用.
参考文献
篇2
关键词：曲面面积 数值计算 数值微分 积分
中图分类号：O172 文献标识码：A 文章编号：1672-3791（2014）03（c）-0238-02
二重积分的数值计算方法有很多，但是在实际应用中，曲面面积的很重要，而曲面面积计算的数值方法却不多，目前还没有找到一种高效、精确的计算其表面积方法。文献[1][2]模型的建立是基于多网格化下小区域内曲面积近似等于平面面积，因此计算结果存在一定误差，且计算精度不易分析。为了减小误差，提高精度，我们建立利用积分中值定理和数值微分公式，建立一个新的计算表面积的数值计算公式―― “四点”插值算法。
1 单元构造和数值计算公式
已知曲面函数为，则考虑曲面在矩形区域内的表面积。对区域进行分割，首先考虑如图1网格单元区域：
利用积分中值定理[3]
则
若，充分小，则由偏导数的连续性有：
，
于是
由三点数值微分公式[4]
，
于是
2 误差估计
其中
由三点数值微分公式[4]
，其中
由二阶泰勒展开：
，其中
于是
其中
同理
所以
3 复化公式
计算矩形区域内函数的表面积，在格网化区域计算表面积。首先对区域进行划分，把目标区域划分成个方格，则有：
，
取如图2的方格，则在每个方格上应用表面积的近似计算公式，只需计算4个信息点。
4 算例分析
例：曲面函数在矩形区域内的表面积。
其表面积计算的精确值为：
在相同的分割网格下：
“四点”插值算法节点数：
三角形法需要的节点数：
数值计算结果如表1。
5 结语
通过实验的matlab仿真，可知基于本文的方法求解曲面面积的算法误差和传统的“三角形法”误差虽然都是，但本文方法的误差是“三角形法”的，计算时间是“三角形法”的二十分之一。由此可以看出本算法需要信息点少，精度较好，运算速度快，具有较大的实用价值。
参考文献
[1] 陈吉龙，武伟，刘洪斌.DEM在林地表面积计算中的应用研究[J].西南农业学报，2008，21（5）.
[2] 魏东，张秀程.基于递归算法的三维地形面积计算方法研究[J].沈阳：沈阳工业大学信息科学与工程学院，2007（3）.
篇3
关键词 地层热阻 热响应测试地埋管换热器
中图分类号：TK172 文献标识码： A
Thermal resistance calculation on vertical ground heat exchanger
Li Jinghui1Wang Jiankui2 Lu Lin2 Fang Xugen1
1 Zhejiang Construction Division Building Energy Technology Co., Ltd.
2 Zhejiang Academy of Building Research & Design.ltd
Abstract: Using the data of grock-soil thermal response test, andaccordingto , this paper presents the numericalcalculationofthe vertical ground heat exchanger's thermal resistance calculation formula,andalso analyzeshow heat exchanger resistance affect the performance of heat exchanger. The calculation and analysis have a certain reference valueinguidingthe design and construction of ground heat exchanger.
Key word: ground heat resistance，heat responsetest，ground heat exchanger
0 引言
地层热阻是决定土壤源热泵合理设计和科学应用的核心因素，是影响地下埋管换热器传热性能及土壤源热泵系统节能与经济性的重要原因。换热器的热阻工程实际计算及热阻的影响对于设计人员设计及方案的优化有着重要的影响。
1 地埋管换热器的热阻计算
地埋管地源热泵系统地埋管换热器的设计计算主要有采用单位井深换热量、专业软件和依据《地源热泵系统工程技术规范》[1]（以下简称《规范》）作为目前指导地源热泵工程设计和施工最重要的规范，在其中第4.3.5条明确指出“竖直地埋管换热器的设计也可按本规范附录B的方法进行计算。
但是《规范》附录B中对地层热阻及短期连续脉冲负荷引起的附加热阻的计算公式的表述存在一定的歧义，也不完整，在实际应用中会产生较大的误差甚至错误[2]。本文以供热工况下的换热器的热阻为例各热阻值对进行计算分析。
1.1岩土热响应测试数据
表1 岩土热响应测试基础数据表
项目 测试孔 项目 测试孔
钻孔深度（m） 100 钻孔直径（mm） 135
埋管形式 单U型 埋管材质 PE管
埋管内径（mm） 26 埋管外径（mm） 32
钻孔回填材料 原浆 细沙 主要地质结构 粉质粘土
岩土热响应测试测试报告结果：埋管区域的平均综合导热系数为1.667 W/ m·℃，平均容积比热为2.034×106 J/m3·℃。岩土体初始温度20.1℃。
1.2值计算
公式1
式中：
——传热介质与U形管内壁的对流换热热阻，（m·℃/W）;
——U形管的内径（m）；
——传热介质与U形管内壁的对流换热系数[W/ m2·℃]。
换热器温差＜20℃，根据热物性测试报告定性温度按6℃，查表得=56.3×10-2W/ m·℃，ν=1.553×10-6m2/s ，=11.60 =1547×10-6N·s /m2
d：管内径=26mm管内媒质流速按0.9 m/s计算；=0.9×0.026/1.553×10-6=15067＞104；=0.023×15067×11.6=4020；K==4020×0.563÷0.026=87045 W/ m2·℃；=1÷(3.14×0.026×87045)=1.41×10-4 m·℃/W
1.3 值计算
公式2
公式3
式中：
——U形管的管壁热阻，（m·℃/W）;
——U形管导热系数[w/(m·℃)]；
——U形管的外径(m)；
——U形管的当量直径(m)；对单u形管，n=2；对双U形管，n=4。
=0.42~0.45 w/(m·℃)；=0.032 m； ——0.0452m
=0.0539 m·℃/W
1.4计算
公式4
式中：
——钻孔灌浆回填材料的热阻(m·℃/w)；
——灌浆材料导热系数[w/(m·℃)]；
——钻孔的直径(m)。
=1.47w/(m·℃)， =0.135m，= 0.0809m·℃/w
1.5计算
对应于单个钻孔：
公式5
公式6
式中：
——地层热阻(m·℃/w)；
——指数积分公式；
——岩土体的平均导热系数[w/(m·℃)]；
——岩土体的热扩散率m2/s；
——钻孔的半径（m）；
——运行时间（S）；
根据文献[2-4]对进行修正简化为下式：
公式7
= 公式8
公式9
=1.667w/(m·℃)，=/cρ，cρ=2.034×106J/m3℃，则=0.8196×10-6m2/s，钻孔的半径=0.0675（m）,运行时间=68×24×3600=5.8752×106（S）；z=0.0006
==0.3712 m·℃/w
1.6计算
公式10
式中：
——短期连续脉冲负荷引起的附加热阻(m·℃/W)；
——短期脉冲负荷连续运行的时间（s）。
根据文献[2-4]修正简化为下式：
公式11
=10×3600s
=0.1297 m·℃/w
2 热阻影响分析
根据《规范》中公式B.0.2-4中供热运行份额Fn取0.2。根据以上计算得各热阻值及比例见下表，竖直埋管换热器分项热阻计算分析。
表2 竖直埋管换热器分项热阻计算分析表
项目 ×Fn ×(1-Fn)
热阻值
（m·℃/w） 0.001 0.0539 0.0742 0.1038 0.1184
份额（%） 0.04 15.39 21.18 29.60 33.79
根据以上数据可知，单孔U型管回填材料的热阻、地层热阻和短期连续脉冲负荷引起的附加热阻为构成换热器热阻的主要部分，U型管内壁的对流换热热阻很小，对换热器热阻影响几乎可以忽略。
3 结论
3.1对竖直地埋管换热器各热阻值进行了详细的计算，结合相关文献对地层热阻及短期连续脉冲负荷引起的附加热阻进行了简化计算。
3.2 计算结果表明，单孔U型管回填材料的热阻、地层热阻和短期连续脉冲负荷引起的附加热阻为构成换热器热阻的主要部分，对于优化换热器设计方案具有一定的参考价值。
参考文献
[1] 中国建筑科学研究院.GB 50366-2005 地源热泵系统工程技术规范 [S].北京：中国建筑工业出版社,2009版
[2] 雷建平,於仲义. 关于《地源热泵系统工程技术规范》地层热阻计算式的商榷及应用分析[J], 暖通空调,2009,39(6):27-30
[3] 余延顺. 土壤耦合热泵系统地下埋管换热器传热模型的研究[J].暖通空调,2005,35(1):26-31
篇4
在计算机上运用数值分析解决实际问题的过程为：实际问题→数学模型→数值计算方法→程序设计→上级计算求出结果[1]。数值实验的设计应当充分体现这一过程，同时也应当充分体现数学建模思想. 而目前地方高校在数值分析课程的教学与实验中，普遍存在重理论、轻实践、纯数学十足的问题，部分即使避免这些问题，做到利用计算机进行可视化教学和算法编程实践，但学生也缺乏分析问题、解决问题的能力，因为学生的算法与程序多是百度而来，缺乏思考. 因此，数值分析课程的数值实验因当因时而异，因专业而异，紧跟时展，做到充分吸引学生的注意力，激发学生的兴趣，调动学生积极性，让学生积极主动的去想办法解决问题，而非被动的接受。这样的数值实验才能有效培养学生应用能力。
数值分析的内容广泛，包含插值拟合、数值积分与数值微分、数值代数、微分方程数值解法、非线性方程与方程组的数值解法，最优化等等。下面以插值拟合、微分方程数值解法、非线性方程数值解法、最优化的相关实际问题为例，研究数值实验的设计。
2 结束语
通过以上案例发现：如果知道解决实际问题的数值分析方法，利用MATLAB能很方便的求出实际问题的结果。但是对于多数学生而言，由实际问题无法得到数学模型，更无法知道相应的数值分析方法。因此在学生有一定数值分析基础后，才能引入这些实际问题的数值实验；同时选取的数值实验必须为学生深入研究预留了空间，因为在讲解相关的知识背景，详细分析问题，建立模型后，对学生进行的是分层次指导解决这些问题（即或利用MATLAB函数命令简单编程计算，或利用MATLAB的专用工具箱计算，或设计算法流程利用MATLAB编程计算）；这样每一个学生都能参与到实验中来，每一个学生在实验中都有收获；最后撰写实验报告，阐明实验的目的、要求、过程、收获等.如此，通过这些数值实验，学生既自己动手解决实际问题，培养了其综合应用能力，又让其充分体会到数值分析的魅力，为大学生数学建模竞赛和后续专业学习打下坚实基础.当然，随着技术的进步，不同专业的数值分析应用实验需要教师不停的探索和研究。
参考文献：
[1] 徐翠薇，孙绳武.计算方法引论[M]. 3版.北京：高等教育出版社，2007.
[2] 吕国英.算法设计与分析[M]. 2版.北京：清华大学出版社，2009.
[3] 张德丰等. MATLAB数值分析[M]. 2版.北京：机械工业出版社，2012.
[4] 杜廷松.关于《数值分析》课程教学改革研究的综述和思考[J]. 大学数学，2007，23（2）：8-15.
[5] 张涛，. 数值实验在《数值分析》教学中的重要性分析[J]. 长江大学学报：自然科学版，2009，6（1）：371-372.
篇5
关键词：质子交换膜燃料电池；数值模拟；两相流；水传递
中图分类号：U473.4文献标文献标志码：A文献标DOI：10.3969/j.issn.2095-1469.2012.06.03
质子交换膜燃料电池（Proton Exchange Membrane Fuel Cell, PEMFC）可以看作是一个正极不断生成水的发电装置。如果正极生成的水不能有效排出，就会造成正极水淹。一方面造成催化层氧气供应不足，此时原本在高电位下进行的氧气得电子产生还原反应，就会被低电位下进行的H+得电子还原反应取代，电池端电压由正转负，形成所谓“负差效应”；另一方面催化层活性区域减小，导致催化剂利用率降低。目前大部分工作放在系统设计和工程方面，这些策略常不可避免地会带来严重的寄生功率损耗。对小功率便携式PEMFC而言，由于所容许的空间和能量有限，不可能采用主动型排水装置。对于大功率的燃料电池系统，虽然允许采用增加气体流速，强制排除水分措施[1-5]，但是，这些方法只能排除双极板流场中的水分，无法有效解决夹藏在多孔电极孔隙中的水，而PEMFC生成的第1滴水恰恰在气体多孔电极的催化层中。
多孔电极内液态水量过多，液态水会淹没电极减少电极有效面积，同时堵塞传质通道严重阻碍氧气传输，使阴极氧气供应不足，浓差极化增大，从而导致电池性能大幅度下降[6-9]。正极水淹造成催化层氧气供应不足的核心问题是：非极性的O2分子在极性水分子中的溶解度太低，在PEMFC工作温度（60～80?℃）下，由于水的蒸汽压增大，O2分压减小，O2在水中的溶解度还会进一步减小。我们基于“相似相溶原理”，在电极微孔中预先加入对O2有很高溶解度的非极性硅油（DMS），代替孔隙作为气体不可侵占的传输通道，确保氧气的供应不会因为水的堵塞而中断，而未被硅油占据的孔隙则作为反应生成水的固定排除通道，从而实现了气体传输与水排出通道的有序分离[10]。图1为传统气体电极（Conventional Gas Pt/C Electrode, CGE）和抗水淹电极（Anti-Flooding Pt/C Electrode, AFE）示意图，图2展示了CGE催化层内氧气的通道为空气和液态水，AFE催化层内氧气的通道为空气、DMS和液态水。AFE的优势主要体现在以下几方面：（1）氧气在DMS内的溶解度是水中的10倍，DMS提供了输送氧气的专用通道。（2）DMS主要填充在孔径在20?～70?nm之间的孔隙中，因而解决了最难克服的水淹问题，该尺度内的传质是传统的对流的方法难以解决的。（3）气体多孔电极催化层微孔而非双极板流场的水淹问题才是影响燃料电池性能的首要问题，而AFE直接解决了发生在多孔电极内的水淹问题。
本文在文献[11]和[12]模型基础上，针对抗水淹电极内DMS在多孔电极内的作用，发展一个三维两相流的数学模型，用于研究抗水淹电极单电池，然后讨论多孔电极内液态水饱和度分布和电流密度，分析抗水淹电极内DMS对燃料电池输出特性的影响。
1 几何模型与数学方程
图3中试验与数值计算都是基于图3（a）所示的蛇形流道。图3（b）的单电池由阳极极板、阳极流道、阳极扩散层、阳极催化层、质子交换膜、阴极催化层、阴极扩散层、阴极流道和阴极极板依次组合而成，图3（c）为单电池的计算域，电极活性面积为20?mm×21?mm。三维数值计算域是根据流场板的流道外形来建立流道的实体几何域作为流场计算的流体域。由于流场板为固体导电区域，只有电传导和热传导的计算，不涉及多组分的流体传输和各种源项的计算，因此流场板的网格尺寸相对大些，而催化层是电化学反应发生的区域，涉及多种源项的计算和多组分流体的传输，所以对催化层网格进行了局部加密。
本文的数值模拟工作基于一个稳态、多组分、非等温和三维的两相流传递模型，模型考虑水的相变，模型简化了液态水的生成和传输过程，将液态水饱和度分布描述为水蒸汽过饱和度的简单函数，这样就可以通过气体扩散传质方程得到的水蒸汽分压力直接计算出液态水饱和度分布。液态水饱和度分布的结果将被用来修正反应气体扩散系数和电化学反应速率方程，这样模型就能够模拟液态水生成堵塞传质通道和淹没电极有效面积的影响。
流场边界条件：阴阳极入口设为质量入口边界，先设定进口边界上总的质量流量，然后再指定不同组分的摩尔分数；阴、阳极出口设为压力出口，设定相对压力为0?Pa。流场板和流道的交界面以及流场板和扩散层的交界面设为壁面边界条件；流道和扩散层的交界面、扩散层和催化层的交界面以及催化层和质子交换膜的交界面均设为内部连续性边界，其余边界为壁面边界条件。
温度场边界条件：模型外部边界面及流道入口和出口均设定为燃料电池的工作温度；内部交界面设为连续性边界条件。
电势场边界条件：阴极流场板电位边界设为电池输出电压，阳极流场板电位边界设为0?V，内部各交界面设为连续性边界条件，其余边界设为绝缘边界条件。
数学模型应用了以下假设：（1）气体混合物为理想气体。（2）电池内流体的流动为层流。（3）多孔介质为各向同性，即在同一层中具有相同的特征参数，如渗透率、组分扩散系数和黏性系数等。（4）流道
足够长，即流动在出口处达到充分发展状况。（5）催化层内催化剂和DMS含量比例恒定且均匀。
质量守恒方程：
，
x、y和z 3个方向的动量守恒方程：
，
组分守恒方程：
，
温度控制方程：
，
式中：为孔隙率；ρ为密度；为速度矢量；Sm是由于电化学反应发生而导致的混合气体质量源项；ui为x、y和z 3个方向的速度；为动力学粘度；p为流体压力；Su为动量源项；Ci为组分体积浓度； Dm,i为组分扩散系数；Si是由于电化学反应发生而导致的各组分质量源项；Cp为等压比热容；T为热力学温度；keff为有效导热系数；ST为能量源项。各方程具体源项和参数计算方法参见文献[11]和[12]。
以下是针对抗水淹电极建立的局部有效扩散系数修正方程。
传统多孔电极里组分气体的有效扩散系数由孔隙率和水的饱和度s共同决定，水的饱和度定义为多孔电极内局部液态水的体积与孔隙体积的比值。
，
，
式中：Dm,i为组分气体的扩散系数；VWL为局部液态水的体积；Vpore为孔隙体积。
有效介质理论[13]用于估计有效扩散系数，对于重叠的球状颗粒
．
多孔电极液态水的饱和度，通过减小扩散面积和增加扩散路径的弯曲阻碍了气体的扩散
g(s)=(1-s)x ．
如果孔隙的结构不随着水饱和程度改变，根据有效介质理论，x=1.5。Mezedur等[14]针对他们模型中数值计算得到x=0.71。Martys[15]用数值方法研究了多孔介质中液体以球状形式存在的影响问题，通过网格计算得到液态水的分布，结果表明x=2。但对于接触角为90°的情况，x远远大于2。因此多数理论分析中x一般大于2[16]。所以x的值是与具体的多孔介质的结构有关，是不同燃料电池多孔电极的物性参数，本文通过数值仿真结合试验结果推算出x值为2。
针对抗水淹多孔电极，多孔电极里组分气体的有效扩散系数由孔隙率、水的饱和度S和DMS在多孔电极里的体积分数Soil共同决定。假设DMS均匀分布在催化层内，在式（8）基础上发展为式（9）。
，
，
，
式中：为氧气在DMS内的扩散系数；koil为扩散系数的比值；Voil为局部DMS的体积。
氧气在DMS内的扩散系数由Wilke-Chang公式（12）[17]计算
，
式中：为溶质A在溶剂B里的扩散系数；MB和分别为溶剂B的摩尔质量和粘度mPa・s；T为绝对温度 K；VbA是溶质A在沸点时的摩尔体积 （O2: 25.6?cm3/ mol）；是与溶剂B相关的因子，DMS是1。在温度333?K时，氧气在DMS里的扩散系数是1.0×10-5?cm2 /s。
氧气在水蒸汽里的扩散系数由式（13）计算
，
式中：Dij为双组元气体扩散系数。当 p0=1.01×105?Pa，T0=352.3?K时，DO2-H2O=0.352?cm2/s；p0=1.82×105?pa，T0=333?K时，氧气在水蒸气里的扩散系数DO2-H2O=
0.182?cm2/s。则由式（10）可得koil=5.5×10-5，所以在式（9）中koil Soil可以被忽略掉。
根据试验测量的DMS填充的面密度值（如2.5?mg/cm2）乘以填充的面积4.2?cm2计算出DMS的填充总质量，然后将填充总质量除以DMS的密度就可以得到DMS的体积Voil，使用式（11）就可以计算得到DMS在多孔电极里的体积分数Soil。模拟计算DMS填充的面密度值2.5?mg/cm2对应的Soil值为0.333。
将发展后的气体有效扩散系数修正式（9）绘制成曲线图，如图4所示。图中展示了电极分别是抗水淹电极（其中DMS填充的面密度值2.5?mg/cm2）和传统电极（其中DMS填充的面密度值0?mg/cm2）时液态水饱和度对组分气体扩散系数的影响，单个曲线体现的是随着液态水在多孔介质里的增加，因此阻碍了气体的扩散，所以有效扩散系数减少。曲线间的不同体现的是，多孔电极填充了DMS后，由于DMS占据了部分气体通道同样也会阻碍气体的扩散，所以气体的有效扩散系数整体降低了。
电极反应速率：电极反应速率采用Butler-Volmer动力学方程。
，
，
式中：Ra和Rc分别为阳极和阴极电极反应速率，下标a、c分别代表阳极、阴极；Aa和Ac为催化层比表面积；i0,a和i0,c为交换电流密度；和分别为氢气和氧气的当地摩尔浓度；和分别为氢气和氧气的参考摩尔浓度；和为浓度指数，对于阳极=0.5，对于阴极=1；和为电化学反应能量传递系数。
在阴极催化层内反应生成的液态水不仅占据了传质通道而且减少了催化剂的有效表面积，传统多孔电极通过式（16）修正，考虑液态水对催化剂有效表面积的影响，如图5所示。
k(s)=1-s．
为了描述液态水和DMS对电化学反应速率的影响，水的饱和度S和DMS在多孔电极里的体积分数Soil共同修正阴极的电化学反应方程，在式（16）基础上发展为式（17）。
k(s)=[Soil+(1-Soil)(1-s)]y ，
式中：电化学反应速率修正函数指数y为3.5。
边界条件和数值计算方法：在流道入口处解动量平衡方程时采用质量流量边界条件，入口质量流量、组分浓度和温度已知；在流道出口处解动量平衡方程时采用压力出口边界条件，在该面上的所有变量（压力除外），梯度都为0；在流道进出口之间的固体壁面（除燃料电池各部分交界面）采用组分流量为0，无滑移壁面边界条件，并假定极板外壁面温度恒定为电池工作温度，通过给定两极板处的电势求电流。模型仿真求解使用FLUENT12.0，应用求解压力耦合方程组的半隐式方法SIMPLE(Semi-implicit Method for Pressure-linked Equations)[18] 对所有的守恒方程进行迭代求解，直到连续两次迭代的各个物理量相对误差为10-5。
2 抗水淹电极模型验证
为了通过数值模拟计算来研究抗水淹电极PEMFC，使用上述考虑了DMS影响得到的式（9）和式（17）来分别模拟传统电极PEMFC和抗水淹电极PEMFC。模拟计算中传统电极PEMFC的Soil值为0，抗水淹电极PEMFC中的DMS含量为2.5?mg/cm2。
本文根据单电池试验的实际情况，设置了仿真数学模型的物理参数和运行工况列于表1，其它参数主要从文献[11]得到。
图6对比了PEMFC试验值和模拟计算值。图中显示了不论是传统电极还是抗水淹电极的电池极化特性曲线，模拟计算和试验值都吻合得较好，模拟计算和试验结果一样都显示出了活化极化区、欧姆极化区和浓差极化区；还显示出传统电极PEMFC和抗水淹电极PEMFC性能的差异，试验中单电流密度值大于1.32 A/cm2后阴极为抗水淹电极的单电池性能要好于传统电极单电池性能，模拟计算结果中当电流密度值大于1.3 A/cm2后阴极为抗水淹电极的单电池性能要好于传统电极单电池性能；在输出电压为0.35 V时，试验中抗水淹电极较传统电极PEMFC的电流密度值提高了0.87 A/cm2，模拟计算中抗水淹电极较传统电极PEMFC的电流密度值提高了0.74 A/cm2。
图7和图8展示的是当输出电压为0.35 V时阴极催化层中间截面的电化学反应速率云图。由图7可见传统电极电化学反应速率在图的右边低左边高，也就是说产生的电荷分布右端低左端高，整体分布不均匀；由图8可见抗水淹电极电化学反应速率在局部上成高低的交错变化，但是在整体上产生的电荷分布较均匀，没有出现图 7中产生的电荷电池的一端低而另一端高的现象。这说明阴极电极填充了DMS后，能够使电荷分布均匀，从而改善电池输出性能。
3 抗水淹电极性能预测
为了进一步研究DMS填充量对电池输出特性的影响，应用上述的抗水淹电极模型对不同的DMS含量的燃料电池极化特性进行模拟计算。模拟计算中设置了7种DMS含量，从0?mg/cm2(即电极不含DMS的传统电极)到6?mg/cm2。将不同的DMS面密度值转换成DMS的在孔隙里体积分数后，分别代入到式（9）和式（17），得到不同DMS含量时所对应的扩散系数修正方程和电化学反应速率修正方程，见表2。
图9展示的是将不同DMS含量的扩散系修正函数的曲线图，由图可以看出随着DMS的含量由0?mg/cm2增大到6?mg/cm2的过程中，有效扩散系数的值出现整体下降，这反映的是DMS对气体扩散的阻碍；在每种DMS含量下，随着液态水饱和度的增大，气体有效扩散系数不断下降，这反映的是液态水对气体扩散的阻碍，所以计算模型同时反映了DMS和液态水共同对气体扩散的阻碍效果。
图10展示的是不同DMS含量的电化学反应速率修正方程的曲线图。由图可以看出随着DMS的含量由0?mg/cm2增大到6?mg/cm2的过程中，催化剂的有效面积整体提高了，这反映的是DMS起到增大催化剂有效面积的作用；在每种DMS含量下，随着液态水饱和度的增大，催化剂有效面积不断减少，这反映的是液态水使催化剂有效面积减少，所以计算模型同时反映出了DMS和液态水共同对催化剂有效面积的影响。
综合图9和图10的分析可知，增大DMS含量，一方面增大了催化剂的有效面积，另一方面也增大了对气体的阻碍作用，前者有利于电池输出特性的提高，后者不利于电池输出特性的提高。因此，有必要对DMS的含量进行选择，采用表2中不同DMS含量时的扩散系数修正方程和电化学反应速率修正方程，对质子交换膜燃料电池进行了模拟计算。
图11展示的是不同DMS含量时PEMFC的极化特性。对比DMS含量为0?mg/cm2、1?mg/cm2、2?mg/cm2、3 mg/cm2和4?mg/cm2各极化特性曲线发现，随着DMS含量的增大，在浓差极化段的输出电压不断增大，即随着DMS含量的增大PEMFC的输出特性得到了提高。但是当DMS含量为5?mg/cm2而输出电压为0.45?V时，极化特性曲线出现了剧烈的下降，并且电流密度值有倒退的趋势，这说明增大DMS含量对燃料电池输出特性的负面影响开始占主导地位，即由于气体扩散能力的降低，使多孔电极更早地出现了传质极限，因此输出电压急剧下降，当DMS含量为6?mg/cm2时输出电压为0.5?V时就出现了电压的急剧下降，可见当DMS含量超过4?mg/cm2后在高电流输出时电池性能开始下降。为了更清楚地对比不同DMS含量所对应的电流密度值，取输出电压为0.35?V时不同DMS含量所对应的电流密度值绘制成图12，图中显示开始电流密度值随着DMS含量的增加而增加，当DMS含量超过4?mg/cm2后电流密度值开始下降。这说明抗水淹电极内DMS的含量不是越大越好，DMS的含量对于提高PEMFC性能存在着最优值，在本次试验电池所对应的仿真模拟计算中DMS含量的最优值为4?mg/cm2。
文中采用DMS作为电极多孔介质中氧气额外的固定传输通道，是在电极微孔中预先加入进去的，它不参与化学反应，不会被消耗，因此其生命周期与PEM同步，具有耐久性。
4 结论
篇6
关键词: 图像去噪；全变分模型；数字TV滤波器；扩散系数
中图分类号:TP391.41
文献标识码:A文章编号:1672-8513(2010)03-0211-05
A New Numerical Realization Algorithm Based on the Total Variational Model
WANG Yiyan
(Department of Physics and Engineering Technology，Sichuan University of Arts and Science，Dazhou 635000，China)
Abstract:
Based on the classical total variational restoration model，this research proposes a new numerical realization algorithm of digital TV filter for image denoising from the viewpoint of local diffusion. The proposed algorithm can select the diffusion coefficients adaptively according to the local gradient of the image，so as to reduce the diffusion in edge regions while enhancing the diffusion in flat regions. It is used as an anisotropic filter for the noise image. So the proposed algorithm can not only suppress noises very effectively but also preserve image edge details very well. Experimental results show that this algorithm has better performance and higher computation efficiency than the conventional digital TV filter.
Key words: image denoising; total variational model; digital TV filter; diffusion coefficient
图像在获取或传输过程中不可避免地会遭受噪声的干扰，从而导致图像质量下降，对图像的后续处理产生严重的影响.为了消除这种影响，就必须在后续处理前，对图像进行去噪预处理.因此，图像去噪一直是图像处理领域的研究热点.其难点在于：去噪的同时要求保留图像的边缘等细节信息.目前，基于偏微分方程的变分模型方法是图像去噪的一种流行方法，其高质量的处理效果已引起了国内外广大学者的重点关注［1-6］.Rudin等［1-2］提出的Total Variation(TV),模型是其中的典型代表.后来，Chan等［4］在文献［1］的基础上提出了数字TV滤波器，并证明该滤波器是基于TV模型的最小能量泛函在半点格式下［7］的等价数值求解格式.由于数字TV滤波器是建立在TV模型的基础上，因此它具有良好的边缘保持特性.同时，数字TV滤波器的设计特点还使之具有特殊的信息反馈机制.国内肖亮等［5］也对数字TV滤波器的设计及性能进行了深入的分析与研究，验证了数字TV滤波器的上述性能.然而，传统数字TV滤波器本质上是一种倒数梯度加权滤波器［5］，对图像去噪后的像素值是通过原来噪声图像中对应点的周围像素的梯度值平均加权得到，因此，它是一种线性滤波器，在去噪时容易导致图像模糊.本文根据全变分最小的能量泛函，从局部扩散的角度出发，设计了一种新的应用于图像去噪的数值算法.该算法可以根据图像的局部梯度信息自适应地选取扩散系数，使边缘处扩散较弱，在平坦区域扩散较强.实现了对噪声图像的各向异性滤波处理，在去除噪声的同时能够较好地保护图像的细节信息.实验结果表明，该算法去噪后的图像比传统数字TV滤波器具有更好的客观评价指标和主观视觉效果，同时运算效率也更高.
1 全变分模型及其数字滤波器
1.1 全变分模型
令u为原始清晰的图像，u0为被噪声污染的图像，即u0=u+n.式中n是均值为0，方差为σ2的高斯噪声.即：E(n)=0，E2(n)=σ2.（1）
全变分模型是由Rudin等［1-2］提出的，并且现在是图像恢复中最成功的方法之一.Rudin等认为，有噪声图像的的全变分比无噪声图像的全变分明显大，最小化全变分可以达到去除噪声的目的.全变分模型是一种变换，可用下式表示：
TV(u)=∫ΩΔudxdy.（2）
其中Δ u=uxuy，Δu=ux2+uy2是图像的梯度场，Ω为图像的定义域，(x，y)∈Ω为像素点的二维坐标.
并假设所有可观测的图像具有有界变分，有界变分空间定义为：
艳：基于全变分模型的新型数值实现算法
BV(Ω)={u:u∈L1(Ω)，并且TV(u)
BV(Ω)是在BV范数uBV=uL1+TV(u)下的Banach空间.
由（1）式可知，TV范数最小化应满足以下2个条件：
∫Ωudxdy=∫Ωu0dxdy，1Ω∫Ω(u－u0)2dxdy=σ2.（4）
其中，Ω为图像区域Ω的面积.对于常数c，有TV［u+c］=TV［u］.因此，第1个条件通常总是成立的［4］.于是，从u0中恢复u的问题就演变为在条件1Ω∫Ω(u－u0)2dxdy=σ2下，求TV(u)=∫ΩΔudxdy的最小化问题.通过引入参数λ，可以定义全变分最小的能量泛函：
minuJ(u)=∫ΩΔudxdy+λ2∫Ωu－u02dxdy.（5）
上式中第1项为图像的正则化项，第2项为数据保真项.λ为Lagrange乘子，是平衡正则化项和保真项之间的权重参数.
由于目标函数J(u)是u的凸函数，极值存在的一阶充分必要条件为：Δ J(u)=0.然而由于当Δu=0，一阶条件不能满足.为了克服这种退化问题，一般对图像正则化项引入小扰动进行修正：
Δuε=ux2+uy2+ε2.（6）
只要ε保持足够小，便不会影响噪声滤波的性能.所以式（5）改写为：
minuJ(u)=∫ΩΔuεdxdy+λ2∫Ωu－u02dxdy.（7）
这样，就建立了图像恢复的全变分模型.
1.2 数字TV滤波器
由一阶最优条件，式（7）可归结为具有黎曼边界条件的欧拉-拉格朗日方程：
－Δ•ΔuΔuε+λ(u－u0)=0.（8）
记Ω为Ω的边界，n为边界Ω的法向量，黎曼边界条件为：
un=0，x∈Ω . （9）
即式（8）成立时J(u)最小.由梯度下降法得：
ut=Δ•ΔuΔuε+λ(u0－u).（10）
数字TV滤波器［4］实质上就是对上述极值问题的数值逼近.对于数字图像，设α是待处理的中心像素点，像素值为uα，β等是α的十字邻域，如图1所示.
图1中hαα(u)和hαβ(u)等为加权系数，分别定义为：
hαα(u)=λλ+∑β~αωαβ(u)，hαβ(u)=ωαβ(u)λ+∑β~αωαβ(u),（11）
其中：ωαβ(u)=1Δαuε+1Δβuε.（12）
显然，hαα(u)+∑β~αhαβ(u)=1.从式（11）和（12）可以看出，加权系数与中心像素和它邻域像素的正则化梯度的倒数和有关.定义中心像素α和它的十字邻域的正则化梯度分别为：
Δαuε=∑β~α(uβ－uα)2+ε2，
Δβuε=(uβ－uα)2+ε2.（13）
结合式（11）~（13），则可得到TV数字滤波器，记为Fε，λα(u)：
Fε，λα(u)=∑β~αhαβ(u)uβ+hαα(u)u0α.（14）
由式（14）可知，TV数字滤波器包含2个参数：
1） 正参数ε称为规整化参数，为一个很小的正数，目的是避免图像平坦区域梯度为0.ε取值很小，当ε≤10－4时，对滤波器影响很小.
2） 正参数λ称为匹配参数，即Lagrange乘子，决定了滤波器如何在边缘保护和平滑噪声之间取得平衡，其选取依赖于噪声水平σ.噪声越大，λ越小，平滑作用越明显，边缘损失越大；反之则λ取值越大，边缘保护越明显.
TV数字滤波器采用迭代策略对图像中的每个像素点进行处理，整个滤波过程如下：
初始化u(0)，通常取u(0)=(u(0)α1，u(0)α2，…，u(0)αN)
For k=1，2，….
For α1α2α3…αN
根据式（11）~（13）分别计算以α为中心点的Δαuε及其相邻Δβuε等4点、ωαβ(u)、hαβ(u)、hαα(u)，然后利用式（14）计算u(k)αj=Fε，λαj(u(k－1)).
特别地，由于图像在迭代过程中不断变化，因此，需要适时更新λ值［4］：
λ=1σ21Ω∑α∈Ω∑β~αωαβuβ－uαuα－u0α,（15）
其中，σ2是图像的噪声方差，需要估计，Ω是图像的面积.通常λ被初始化为1σ2.
由上述分析可知，数字TV滤波器实质上是一种倒数梯度加权滤波器［5］，对图像去噪后的像素值是通过原来噪声图像中对应点的周围像素的梯度值平均加权得到.因此，它是一种线性滤波器，在去噪时容易导致图像模糊.
2 本文算法
针对数字TV滤波器的不足，我们对全变分模型给出了一种新的结合梯度加权最速下降法和半点格式的数值迭代算法.
首先由变分法［8］可推导出式（5）对应的欧拉-拉格朗日方程为：
－Δ•ΔuΔu+λ(u－u0)=0.（16）
由梯度下降法，引入人工时间变量t，则可得如下方程：
ut =Δ• ΔuΔu + λ(u0 - u),
unΩ=0,u(x,y,0)=u0 (x,y) . （17）
令η为图像的梯度方向，ξ为与η垂直的切线方向（即边缘方向），分别表示如下：
η=ΔuΔu=(ux，uy)u2x+u2y，ξ=ΔuΔu=(－uy，ux)u2x+u2y.（18）
定义uξξ和uηη分别为u在ξ和η方向上的二阶导数，表达式如下：
uηη=u2xuxx+2uxuyuxy+u2yuyyu2x+u2y,uξξ=u2yuxx－2uxuyuxy+u2xuyyu2x+u2y . （19）
则式（17）可化为η－ξ方向上的偏微分方程(PDE)：
ut=1Δu•uξξ+λ(u0－u).（20）
由式（20）可看出，边缘方向对应的扩散系数为1Δu，它依赖于图像的局部梯度信息.在不同区域，扩散程度不一样.在边缘处，Δu较大，扩散系数较小，因此沿边缘方向的扩散较弱，从而有利于保护边缘；在平坦区域，Δu较小，扩散系数较大，因此在平坦区域的扩散能力较强，从而可以有效去除噪声，整个过程相当于对噪声图像进行各向异性滤波处理.因此，该方法能在去噪的同时保护图像的边缘细节.
采用半点格式的中心差分来离散化式（20）中的PDE，为了利用像素周围更多的信息，本文选择8邻域系统，如图2.
由一阶和二阶导数的差分格式可以得到Δu，uξξ的离散形式，进一步对式（20）利用时间步进法求解，其迭代形式为：
un+1i，j=uni，j+δt•1Δu•uξξni，j+λu0－uni，j.（22）
在实际计算过程中，设图像大小为M×N，每个像素位置点(i，j)∈1，…，M×1，…，N，采取对称边界条件为：u0，j=u1，j，uM+1，j=uM，j，ui，0=ui，1，ui，N+1=ui，N.并且为了避免上述方程中分母处的梯度为0，将梯度提升为Δuε=Δu2+ε2，其中，ε>0.Lagrange乘子λ可通过对式（16）两边同时乘以(u－u0)，然后在整个图像区域Ω上积分求得：
λ=1σ2Ω∫ΩΔ•ΔuΔu•(u－u0)dxdy.（23）
3 实验结果与分析
为了验证本文算法的有效性，我们进行了大量的仿真实验.采用平均绝对误差(MAE)和峰值信噪比(PSNR)作为去噪性能的客观评价尺度，MAE和PSNR分别定义如下：
MAE=1MN∑i，jyi，j－xi，j，（24）
PSNR=10 lg(M×N)max(xi，j)2∑i，j(yi，j－xi，j)2，（25）
式中，x是大小为M×N的原始图像，y是去噪后的图像.MAE值越小越好，而PSNR值越大越好.
以Woman灰度图像为测试对象，对其加入均值为0，不同噪声方差（依次为5，10，15，20）的高斯噪声，将本文算法与数字TV滤波器［4］进行了去噪比较.图3给出了方差为20时，Woman图像去噪后的视觉效果；表1为不同噪声条件下，2种方法去噪前后的客观评价指标对比；表2为2种算法在去噪过程中所需时间,运算在Matlab 6.5平台上进行.
从实验结果来看，本文算法去噪效果、客观评价指标和运算效率方面都比传统数字TV滤波器具有明显的优势.
4 结语
本文对全变分模型提出了一种新的数值实现方法，与传统数字TV滤波器不同，该算法从局部扩散的角度，充分利用了图像想素点的梯度信息，自适应的选取扩散系数，对噪声图像进行各向异性滤波处理.实验结果表明，本文算法在去除噪声的同时能较好的保护图像细节，与传统数字TV滤波器相比，本文算法在去噪效果、客观评价指标和运算效率等方面均具有明显的优势.
参考文献:
［1］RUDIN L，OSHER S，FATEMI E.Nonlinear total variation based noise removal algorithms［J］.Physica D，1992，60(1-4): 259-268.
［2］RUDIN L，OSHER S.Total variation based image restoration with free local constraints［C］//Proc 1st IEEE Int Conf Image Process，Austin，TX: IEEE Press，1994，1: 31-35.
［3］COBSON D，VOGEL C.Convergence of an iterative method for total variation denoising［J］.SIAM Journal of Numerical Analysis，1997，34(5): 1779-1791.
［4］CHAN T F，OSHER S， SHEN J.The digital TV filter and nonlinear denoising［J］.IEEE trans on Image Processing，2001，10(2): 231-241.
［5］肖亮，吴慧中，韦志辉，等.基于总体变差模型的数字滤波器设计及其性能研究［J］.信号处理，2003，19(3): 247-251.
［6］庄红林，王文斌，段鹏，等.整体变分最小化方法在图像处理中的应用［J］.云南民族大学学报：自然科学版，2008，17(2): 176-180.
篇7
【关键词】应用型人才；模块化；创新能力；教学改革
0 引言
应用型人才培养模式的具体内涵是随着高等教育的发展而不断发展的“应用型人才培养模式是以能力为中心，以培养技术应用型专门人才为目标的”[1]。本科应用型是本科层次教育，既有着普通本科教育的共性，又有别于普通本科的自身特点，它更加注重的是实践性、应用性和技术性。
随着现代科学技术的发展和计算机的广泛应用，科学计算已经成为平行于理论分析和科学实验的第三种科学手段。而《数值分析》这门课程是科学计算的基础和核心，逐渐成为本科生和研究生的必修课。是理论与实践结合紧密的一门科学。那么在我们讲授《数值分析》课程的过程中，将理论分析和实验教学相结合，不断地改进我们的教学方法，吸取其他数值分析教育工作者的改革措施和经验[2-4]激发学生学习的兴趣和热情， 提高学生对算法的实际应用能力。
《数值分析》是一门介绍科学计算的核心理论和基本方法的数学课程[5]，分析《数值分析》课程的特点，根据在该课程的具体教学实践中的体会，对该课程的内容进行分类整理，论述以MatLab 为教学平台，并对以数学建模为实验内容的教学方式的改革作初步探索。
1 《数值分析》课程的特点及传统教学中的不足
《数值分析》课程有不同于其他数学专业课的特点， 有其自身的理论体系。数值分析要面向计算机。数值分析是为了解决实际问题而产生的科学，每个算法除了在理论上可行外， 还要通过数值试验证明是可行的。这对学生的编程能力以及数学软件的应用能力提出了较高的要求。
基于《数值分析》课程以上的这些特点可以看出《数值分析》是一门与计算机使用密切结合的，实用性和实践性都很强的课程。但在传统的教学中，只注重讲授数值方法的原理以及误差估计等理论。另外，在教学中由于缺乏对数值实验的重视以及学生数学软件运用能力和编程能力的限制， 使学生不能把这些算法在计算机上实现并应用到解决实际问题中去。这几年我们一直在不断地对教学过程进行思考和探索，不断地总结教学经验，克服教学过程中的不足，不断地对教学方法进行改革。
2 对《数值分析》课程教学改革的几点体会
2.1 数值分析的模块化教学内容改革
“数值分析”是研究各种数学问题求解的数值计算方法及其理论的一门课程[5]。在理、工科数学类科目教学体系中，“数值分析”起着承上启下的作用，是高等理工科院校的重要基础课程，同时也具有培养学生创新思维、创新能力的特点。数值分析的教学目的不是让学生仅仅会利用己有的算法去解决某些问题，而应有理有据地、系统地传授其基本思想、基本方法和基本原理。数值分析的宗旨应是提高学生数学素养和使学生掌握实用算法并重。这就是说要使学生能够知其然，更得知其所以然，并能让学生在遇到新问题时有能力依靠数值方法去解决，注重培养学生举一反三的能力。因此，理论方法和应用都不可忽视，这就要求对教学内容、环节等重新优化设计。
经过我们的分析发现《数值分析》课程主要内容可以总结为三个方面：①离散问题的连续化， 即根据离散数据构造一个简单且易于计算的函数代替原有的复杂函数或数据， 例如函数插值、函数逼近、数据拟合等内容；②连续问题的离散化，即在一系列离散的点上求解数值解，例如非线性方程求根、数值积分、微分方程数值求解等内容；③对于方程组的求解以及矩阵的特征值计算等内容。在教学过程中，我们要从宏观上把课程内容分析给学生，让学生沿着这种思路把课本中的内容系统化， 把各个分散的内容都建立联系进行模块化教学。模块化的教学结构能够更好地适应和启发引导学生积极思维，大胆创新。
2.2 引入MatLab 软件为主的计算机平台，利用现代化的教学手段
在教学过程中，根据本课程的特点，不能一味地只做理论分析，而是要打破传统的教学体系，采用多媒体辅助教学设备，引入采用数值解法、解析解法和图形可视化相结合的方法。
在实际教学中， 我们在利用传统的教学手段的同时，引入现代化教学手段―――电子教学，将一些推导过程、实例应用利用数值计算软件，以可视化的方式来表现数值计算的过程，通过分步实现、误差数据、几何图形直观地展示给学生， 从而使学生更容易理解抽象的理论， 也能让学生学习到如何将算法在计算机上实现。我们将计算机软件MatLab 引入数值分析的各个教学内容之中，在教学中充分利用MatLab 软件为主的计算机平台。MatLab 是最为普遍和功能强大的科学计算工具，我们针对具体的数值算法，编写MatLab 数值算法程序和图形可视化程序， 进行数值数据和图形的比较。引导学生理解所学的内容，激励学生的兴趣，培养他们自己动手解决问题的能力和创新能力。
2.3 从实际问题出发，将数学建模案例和工程应用实例引入课堂教学
《数值分析》是一门应用很强的课程，以数学建模案例和工程应用实例引导学生运用数学方法解决实际问题，因此选择与授课核心内容相对应的数学建模案例和工程实例进行教学具有必要性和可能性.例如用2011年全国大学生数学建模竞赛的C题“山东省职工工资预测问题”作为实例讲解曲线拟合.职工工资增长有两个基本模型：一个是指数增长模型，一个是阻滞增长模型.按照两种不同的模型，使用非线性最小二乘法分别拟合自1978年至2010年的工资数据，预测从2015年至2035年的工资数据，并将预测数据与实际数据比较，可以让学生形象直观地用预测结果分析两个模型的优劣.又如网络流量、核辐射量的测量与计算与数据拟合理论，根据山区地形采样点数据来绘制地形图与插值理论，导弹追踪问题与微分方程数值解，核废料的妥善处理、油罐刻度设计问题与方程组的求根问题等，都是把数值分析与数学建模和工程实例相结合的典型例子，会给《数值分析》课程的教学带来蓬勃
3 结语
应用型本科院校基于模块化的数值分析教学，增强了课堂教学的直观性，优化了课堂教学与实验教学，能真正实现教与学的良性互动。基于现代化的教学手段的数值分析教学，培养学生应用数值计算方法解决实际工程问题的能力，开拓学生思维，培养具有创新能力与素质的应用型现代化建设人才。
【参考文献】
[1]庄华洁，周金其.本科应用型人才培养模式的研究与实践[J].高等教育研究，2004（6）.
[2]李大美，等.基于创新能力培养的计算方法课程改革[J].中国大学教学，2007.
[3]张光辉，任敏.MATLAB平台上《数值分析》课程教学的几点思考[J].甘肃联合大学学报（自然科学版），2012，26（5）：103-105.
篇8
关键词：建筑照明；智能控制；ZigBee技术
中图分类号：TM923 文献标识码：A
随着国民经济的快速发展，能源消耗成为国家和社会关注的热点问题，同时随之而来的节能减排工作成为研究的重要内容之一。据不完全统计在现代建筑中照明系统对能源的消耗占的比例较大[1-2]。因此，建筑电气照明节能技术的理论及应用研究成为了建筑电气综合节能研究内容之一。本文以学校智能照明控制系统设计为例，研究并设计出一套基于Zigbee技术智能照明控制系统，该系统详细的介绍了控制系统的总体设计方案、硬件电路设备选型、Zigbee无线网络的组建以及软件系统设计等内容。本文所提方案和做法对于建筑电气照明节能起到一定的参考和借鉴作用。
1 智能照明控制系统总体设计方案
针对学校的具体情况本文提出一种基于ZigBee 技术的智能照明控制系统。该系统主要由调光控制电路、ZigBee 节点模块、ZigBee 无线通信网络及远程控制上位机四个部分组成。智能照明控制系统结构图如图 1 所示。
从图中可以看出，整个控制过程如下：信息采集电路将所采集的信息经单片机处理后传递给ZigBee 收发节点模块将其转换为无线电信号发送到 ZigBee 无线网络中，经过ZigBee 收发协调器模块传输给上位机，在上位机中将获得相关的现场采集数据。
图1 智能照明控制系统结构图
2 智能照明硬件系统设计
本文的智能照明控制系统主要以单片机为核心，单片机采用AT89S52芯片，系统同时还设有参数修改、数据采集与处理以及联动控制等功能，其结构如图 2 所示。从图2可以看出，该控制电路主要由光照传感器模块、红外传感器模块、ZigBee节点模块、LED 灯管驱动电路模块、保护控制电路、输入/输出等相关部分构成。当系统检测到室内光照度后，首先将其数字化处理，然后送给单片机。单片机根据所检照度值，将所产生的脉冲宽度调制信号送到驱动电路中，从而控制流过 LED 灯管上的电流，以调整室内亮度使之与其相匹配。其硬件选型如表1所示。
图 2 智能照明控制电路结构示意图
表1 智能照明控制硬件选型表
序号 硬件名称 硬件选型
1 单片机 AT89S52
2 通信网络节点 CC2430
3 光照传感器 TSL2561 光传感器芯片
4 红外传感器 BIS0001传感器芯片
3 ZigBee无线网络的组建与智能照明控制系统软件设计
3.1ZigBee无线网络的组建
由于本文的智能照明系统设计对象主要指学校，而且在系统中主要以传感器作为主要接入对象，其数据传输量相对而言不大，采用一般的网络拓扑结构进行数据通信即可。由于学校的直径一般在几公里范围内，然而ZigBee 节点数据传输距离仅有 100m 左右，因此为了增加系统数据的传输距离，在整个网络布置时需增加ZigBee 路由器。正是由于基于以上考虑，本文所设计的无线网络拓扑结构采用簇树型网络拓扑结构。该网络结构具有成本低、实现容易、可扩充、易操作等优点。本文所设计的智能照明控制网络拓扑结构如图 3 所示。从图3中可以看出，学校网络协调器和各楼宇子网路由器是其中非常重要的两个设备，网络协调器完成网络建立、协调整个网络运行。相关各个教学楼路由器主要负责数据传输作用，目的是为保证将现场所采集的数据经过长距离传送后到达协调器。教室中的各个终端控制节点模块，主要负责现场的数据采集和初步处理以及对各个LED控制管理。
图3 智能照明控制系统网络拓扑结构
3.2智能照明控制系统软件设计
智能控制系统程序流程图如图 4 所示，从图4看出，系统首先进行初始化，之后可以通过参数的录入设定系统所需要的光照值，现场采集信息系统含有两个部分，一部分采集光照强度，另一部分通过人体红外传感器采集是否有人员的信息。上述这些信息进入到单片机中执行相应的算法程序，本文所采用的自动控制算法为比例-积分-微分（PID）算法，通过该算法的运算输出合适的控制量以减少偏差可以实现恒定光照控制，通过几次循环反馈系统进入较稳定的状态，教室的光照强度基本上稳定在预先设定的数值左右。本文所采用的调光方式为 PWM方式，由于单片机系统不带PWM模块，因此需要模拟一个PWM信号，本文利用数字电路中的定时器来模拟产生一个 PWM 信号。
图4 智能照明控制系统程序流程图
4结论
本文将ZigBee 技术应用于学校智能照明控制系统中，对智能照明控制系统的软、硬件进行设计，其中重点对系统方案、智能照明控制节点以及ZigBee无线网络的组建进行了详细地阐述，本系统可以根据室内光照调节照度，从而使学校教室内具有很好的照度。同时本系统采用ZigBee方式，减少了布线铺设，使得工程的施工更加容易，达到国家所提倡的节能减排，该系统的设计方案为其它智能照明控制系统提供了参考与借鉴。
参考文献
[1]任红.绿色照明与建筑照明节能设计[J].电气应用，2009（7）：21-25.
[2]陈岁生，卢建刚 . 基于智能检测技术的室内照明系统 [J]. 机电工程技术， 2008，37（5）：53-55.
[3]朱明硕，陈校.浅谈照明节能技术的应用[J].电气应用，2010（3）：54-58.
[4]乐钰唯.照明控制技术的发展及应用现状[J].照明工程学报，2010（3）：14-18.
篇9
[关键词] 旅游收入 增加值 旅游收入乘数
一、旅游乘数的研究回顾
1882年,经济学家巴奇霍特(Bagehot)分析了紧缩产业对经济中其他产业引起的负面影响后,不少经济学家便开始对乘数理论展开了研究。1931年,美国经济学家卡恩(Kahn)首次把乘数概念引入经济理论,他的乘数计量的是总就业量的增量和与之相关的投资品行业的初级就业量的增量之间的比。1936年,凯恩斯(Keynes)对其进行补充,形成了乘数模型。之后现代乘数理论大多沿着凯恩斯的乘数模型和里昂惕夫投入-产出模型发展。1982年马西森(Mathieson)和沃尔(Wall)提出了旅游乘数的概念,即“旅游乘数是这样一个数值,最初旅游消费和它相乘后能在一定时期内产生总收入效应”。
1998年,张凌云在《南开经济研究》上发表的《旅游业乘数效应的几个问题》是我国旅游业乘数研究的早期成果。西方学者弗莱彻(Fletcher)和斯尼(Snee)鉴定了六种旅游乘数,包括产出乘数、销售乘数、收入乘数、就业乘数、政府收入和进口乘数。左冰(2002)、石培华(2003)、吴国新(2003)在各自的论文中研究了我国的产出乘数和就业乘数。但是,更多的学者偏向研究乘数作用对区域旅游的影响,李志清(2001)对上海、李兴绪(2004)对云南、谭佩琴(2006)对澳门、王瑜(2006)对福建、张小利(2007)对西部等。乘数计算方面,我国学者田里(2002)、李天元(2002)、罗明义(2004)等研究旅游乘数的时候,在各自出版的旅游经济学书籍里面均提出了不尽相同的公式。匡林(1996)着重介绍了旅游业乘数的类型和计算,同时介绍了经济分析模型。对于旅游收入以李江帆(1999)为先,认为计算增加值代替总收入有其合理性,其后有魏小安(2000)、宋子千(2001)对此继续深入讨论。智瑞芝(2003)、邹炜(2008)在其论文中分别对黑龙江和云南旅游产业进行了增加值的计算研究。
二、旅游收入乘数的测算
1.增加值法的引入
由国民经济核算理论上看,增加值是生产单位从总产出的价值中间除去其中所含货物、服务消耗的价值之后的余额,它代表的是该生产单位会及各种生产要素在生产过程中新创造的价值(高敏雪,2006)。众所周知,国民生产总值GDP的核算就是从增加值的角度出发,衡量一个国家或地区在一定时期内新创造的价值。在宏观上,增加值算法可以避免重复计算,加总结果可以表现整个经济生产的真实成果。旅游收入是旅游经营者从旅客们在旅游过程中购买服务、货物来满足他们旅行中的游览、住宿、饮食、娱乐、购物等方面获得的货币收入。但是旅游总收入中旅游产品所消耗的中间产品,并不是旅游业所生产的,而是与旅游相关的其他部门生产的。因此,严格上来说,若直接用旅游收入计算在国民经济中所占份额的话,可能会扭曲旅游消费对国民经济的实际贡献。所以,应该从旅游总收入中扣除实物产品和中间消耗的价值,来得到旅游增加值。在实际问题中,由于统计的难度和问题,我国并没有相关单独的统计数据。旅游活动发生在乘坐交通工具、购物、用食、娱乐活动之中,数据上很难区分是属于日常生活还是旅游,这样由旅游产生的经济价值被统计作为相关的交通、商业、餐饮和其他服务业的价值中。难以独立统计,这给相关研究带来了难度,不过,运用旅游增加值剥离测算方法,可以将相关产业中属于旅游创造的价值提取出来。
按照李江帆(1999)对旅游增加值的界定,认为旅游主要涉及交通、邮电、商业、饮食业和社会服务业这五个第三产业部门。这些部门的产品只有一部分投入旅游消费,其增加值需要按一定的比例剥离出,计入旅游增加值中。这些含有游客消费的服务行业提供的增加值中旅游者消费部分所占比重称为该行业的旅游消费剥离系数。利用旅游消费剥离系数分别计算出各个行业的旅游增加值,然后将这些加总就得到该年的旅游总收入。计算步骤如图1所示。在计算的时候基于了两个假设,第一,假定行业内部的增值率相同,也就是比如在餐饮消费中不论是花在旅游上还是在其他方面带来的销售额增长都相同。第二,假定旅游消费剥离系数不随时间的变化而变化。该假设一方面由于统计局限,旅游收入中的行业构成无法得到各年的数据,另一方面也是为了化简计算。这样我们可以认为交通、邮电、商业销售、饮食业和社会服务业这五个行业的旅游消费剥离系数为24.8%、7.2%、3.0%、19.5%和54.1%。再根据统计年鉴上每年相关数据计算各个行业的增加值,乘上剥离系数,得到各个行业的旅游增加值,最后加总获得该年旅游总收入,即:∑(i行业旅游消费剥离系数×i行业增加值)=某年全国旅游增加值。具体计算结果在下表2中。
2.模型的建立和求解
根据何佳梅(2005)关于旅游收入乘数测定的方法,按照柯布-道格拉斯生产函数Y=A・Ka・Lβ来构造关于旅游收入的生产函数GDP=A・Ia,其中GDP为国内生产总值,I为增加值算法下的旅游收入,A为方程的常系数,而a则实际上表示的是旅游收入产出弹性。
由经济学的知识可以知道,旅游收入产出弹性指的就是旅游收入I的变化引起GDP产出的变化,即:
那么对旅游收入的生产函数GDP=A・Ia,根据微积分可知,a表示的就是旅游收入产出弹性,即:Ei=a
给旅游收入生产函数GDP=A・Ia两边取自然对数则可得到:
Ln(GDP)=Ln(A)+a*Ln(I)
利用最小二乘法(OLS),就可以估计出a的值,也就是旅游收入产出弹性。在估算的时候,首先,不直接取GDP的值,而将其指数化,以便保证数据的可比性。G是以1978年为基年计算的GDP指数,替代GDP进行估算。第二,在公式的右边加入时间趋势变量t,用来剔除时间趋势。最后估算公式调整为:
Ln(G)=Ln(A)+a*Ln(I)+b*t
Ln(G)=c+a*Ln(I)+b*t
利用Eviews5.0软件,运用OLS进行参数估计:
Ln(G)=6.92-0.095*Ln(I)+0.10*t
(17.61)(-1.82)(15.60)
R2=0.9966,调整R2=0.9958,F-statistic=1320.488,Prob(F-statistic)=0.000
旅游收入产出弹性为负,这与经济事实不相符合,因而怀疑方程具有多重共线性。根据Klein判别法进行分析,给出解释变量间的简单相关系数矩阵,发现其中一个简单相关系数大于R2,说明模型中存在多重共线性。因此,必须对模型进行修正。利用修正Frish法进行修正得到:
由对比可见,修正后旅游收入产出弹性都得到了比修正前更好的拟合
数据来源:相关年度的《中国统计年鉴》
3.测算结果小结
从模型测算出来乘数的结果上来看,呈现自1996年逐年下降,但是2003年突然增加,再到2004年下降后到2007年逐步上升的走势。自1998年我国出于刺激消费扩大内需的愿望,中央经济工作会议将旅游业列为我国国民经济发展新的增长点后,各省市把旅游业列入支柱产业或主导产业的名单中。在这样的宏观背景下,旅游收入增长率在随后几年超过了国内生产总值的增长率,造成乘数的逐渐变小。但是2003年我国遇上非典的非常时期,旅游业在一定程度上遭到打击,总收入相较于上年锐减了20%左右,但是国内生产总值的增长却仍然保持了一贯的高水平,这样计算出的乘数自然会突增。到2004年旅游业基本恢复到2003年以前的水平,2004年开始我国的国内生产总值继续保持高增长,这也就造成了旅游收入乘数的回升。
三、对旅游收入乘数的认识
测算乘数可以从宏观上反映某行业部门对整体经济的影响程度,体现各部门相互联系,任何部门最终需求的变动都会引起整个经济中产出、收入、就业等变动,即现代经济的特点。一方面,对于乘数较大的行业可以选择着重或优先发展,因为其可较大程度的拉动整体经济,另一方面要看清乘数作用的发挥有一定的局限性。即如旅游业是个产值贡献大、关联性强,俗称“兴一业,旺百业”的行业,其发挥联动增长作用的前提之一是社会上有一定数量可利用的闲置资源、存货、非充分就业等,这样才可以保证有扩张的后劲。对于乘数具体数量上的看待同样要带有理性的眼光,不论多精密的模型所讨论的总是基于若干假设,简化不可量化或者很难度量的条件,可能造成与实际上的偏差。如闫敏(1999)在验证旅游产业和经济发展水平之间关系的文章中,就涉及到了23各部门,这种复杂的产业交织的关系增加了问题的复杂性,要完全的从相关产业中剥离出旅游业的增加值几乎不可能。再如2003年的旅游收入乘数,数值上均比前后两年的大,但我们不可以直接得出当年旅游业发展突出的好的结论,还要看当年的具体宏观背景以及相应发生的不可抗拒的事件。
参考文献:
[1]叶秀霜:旅游经济学[M].北京:北京大学出版社,2005:152
[2]朱沁夫:旅游经济学[M].长沙:湖南大学出版社,2005:154～155
[3]师守祥:旅游业乘数研究辩证[J].旅游学刊,2007,(10):30～33
[4]匡 林:关于旅游乘数理论的几个问题[J].华侨大学学报(社会科学版),2001,(1):9～11
[5]李江帆 李美云:旅游产业与旅游增加值的算法[J].旅游学刊,1999,(5):16～19
篇10
关键词：告警数据 Hadoop Spark
1 引言
随着电信网络的不断演进，全省数据网、交换网、接入网设备单月产生告警原始日志近亿条。以上告警通过网元网管、专业综合网管、智能网管系统[1]三层收敛，监控人员每月需处理影响业务或网络质量的告警事件为20万条，但一些对网络可能造成隐患的告警信息被过滤掉。如何从海量告警数据中获取与网络性能指标、运维效率相关的有价值的数据，对于传统的关系型数据库架构而言，似乎是一个不可能完成的任务。
在一般告警量情况下，ORACLE数据处理能力基本可以满足分析需求，但当告警分析量上升到亿级，如果采用传统的数据存储和计算方式，一方面数据量过大，表的管理、维护开销过大，要做到每个字段建索引，存储浪费巨大；另一方面计算分析过程耗时过长，无法满足实时和准实时分析需求。因此必须采用新的技术架构来分析处理海量告警信息，支撑主动维护工作显得尤为必要，为此我们引入了大数据技术。
2 分析目标
（1）数据源：电信运营商网络设备告警日志数据，每天50 G。
（2）数据分析目标：完成高频翻转类（瞬断）告警分析；完成自定义网元、自定义告警等可定制告警分析；完成被过滤掉的告警分析、TOPN告警分析；核心设备和重要业务监控。
（3）分析平台硬件配置：云计算平台分配8台虚拟机，每台虚机配置CPU16核；内存32 G；硬盘2 T。
3 制定方案
进入大数据时代，行业内涌现了大量的数据挖掘技术，数据处理和分析更高效、更有价值。Google、Facebook等公司提供可行的思路是通过类似Hadoop[2]的分布式计算、MapReduce[3]、Spark[4]算法等构造而成的新型架构，挖掘有价值信息。
Hadoop是Apache基金会用JAVA语言开发的分布式框架，通过利用计算机集群对大规模数据进行分布式计算分析。Hadoop框架最重要的两个核心是HDFS和MapReduce，HDFS用于分布式存储，MapReduce则实现分布式任务计算。
一个HDFS集群包含元数据节点（NameNode）、若干数据节点（DataNode）和客户端（Client）。NameNode管理HDFS的文件系统，DataNode存储数据块文件。HDFS将一个文件划分成若干个数据块，这些数据块存储DataNode节点上。
MapReduce是Google公司提出的针对大数据的编程模型。核心思想是将计算过程分解成Map（映射）和Reduce（归约）两个过程，也就是将一个大的计算任务拆分为多个小任务，MapReduce框架化繁为简，轻松地解决了数据分布式存储的计算问题，让不熟悉并行编程的程序员也能轻松写出分布式计算程序。MapReduce最大的不足则在于Map和Reduce都是以进程为单位调度、运行、结束的，磁盘I/O开销大、效率低，无法满足实时计算需求。
Spark是由加州伯克利大学AMP实验室开发的类Hadoop MapReduce的分布式并行计算框架，主要特点是弹性分布式数据集RDD[5]，中间输出结果可以保存在内存中，节省了大量的磁盘I/O操作。Spark除拥有Hadoop MapReduce所具有的优点外，还支持多次迭代计算，特别适合流计算和图计算。
基于成本、效率、复杂性等因素，我们选择了HDFS+Spark实现对告警数据的挖掘分析。
4 分析平台设计
4.1 Hadoop集群搭建
基于CentOS-6.5系统环境搭建Hadoop集群，配置如表1所示。
4.2 Spark参数设置[6]
Spark参数设置如表2所示。
4.3 数据采集层
数据采集：由于需采集的告警设备种类繁多，故采取分布式的告警采集，数据网设备、交换网设备、接入网设备分别通过IP综合网管、天元综合网管、PON综合网管进行采集，采集周期5分钟一次。采集机先将采集到的告警日志文件，通过FTP接口上传到智能网管系统文件服务器上，再对文件进行校验，通过Sqoop推送到Hadoop集群上。
4.4 逻辑处理层
（1）建立高频翻转告警监控工作流程
先将海量告警进行初步删选，通过数量、位置和时间三个维度的分析，得出高频翻转类告警清单列表，最后由专业工程师甄别确认，对某类告警进行重点关注和监控。
（2）差异化定制方案
按组网架构细分，针对核心重要节点的所有告警均纳入实时监控方案；
按业务网络细分，针对不同业务网络设计个性化的监控方案；
按客户业务细分，针对客户数字出租电路设计个性化的监控方案。
4.5 数据分析层
Spark读取Hive[7]表的告警数据，然后在Spark引擎中进行SQL统计分析。Spark SQL模K在进行分析时，将外部告警数据源转化为DataFrame[8]，并像操作RDD或者将其注册为临时表的方式处理和分析这些数据。一旦将DataFrame注册成临时表，就可以使用类SQL的方式操作查询分析告警数据。表3是利用Spark SQL对告警工单做的一个简单分析：
5 平台实践应用
探索运维数据分析的新方法，利用大数据分析技术，分析可能影响业务/设备整体性能的设备告警，结合网络性能数据，找到网络隐患，实现主动维护的工作目标。
5.1 高频翻转类告警监控
首先制定了高频翻转类告警分析规则，将连续7天每天原始告警发生24次以上定义为高频翻转类告警，并基于大数据平台开发了相应的分析脚本，目前已实现全专业所有告警类型的分析。表4是全省高频翻转类TOP10排名。
5.2 核心设备和重要业务监控
目前以设备厂商或专家经验评定告警监控级别往往会与实际形成偏差，主要表现在以下几个方面：监控级别的差异化设定基于已知的告警类型，一旦网络重大故障上报未知的告警类型就无法在第一时间有效监控到；同一类型的故障告警出现在不同网络层面可能影响业务的程度是完全不同的；不同保障级别的客户对故障告警监控的实时性要求也是不同的。
通过大数据分析平台对差异化监控提供了灵活的定制手段，可根据告警关键字，分专业、地市、网管、机房、告警频次等维度自主定制需要的告警数据，实现日、周、月、某个时间区等统计分析。
应用案例：省NOC通过大数据分析出一条编号为CTVPN80113的中国平安大客户电路在一段时间内频繁产生线路劣化告警，但用户未申告，省NOC随即预警给政企支撑工程师，政支工程师与用户沟通后，派维护人员至现场处理，发现线路接头松动，紧急处理后告警消除、业务恢复。
5.3 被过滤告警分析
全省每天网络告警数据300万条～500万条，其中99%都会根据告警过滤规则进行过滤筛选，把过滤后的告警呈现给网络监控人员。过滤规则的准确性直接影响告警数据的质量。一般来说告警过滤规则可以从具有丰富运维经验的网络维护人员获得，但是这个过程非常繁琐，而且通过人工途径获得的告警过滤规则在不同的应用环境可能存在差异，无法满足网络维护的整体需要。采用大数据技术对被过滤的告警进行分析可以很好地完善过滤规则，让真正急迫需要处理的告警优先呈现给维护人员及时处理，真正做到先于客户发现故障。表5是动环专业被过滤的告警情况分布。
5.4 动环深放电分析
动环网管通过C接口采集蓄电池电压数据，在停电告警产生之后，电压数据首次下降到45 V，表示该局站电池出现深放电现象，通过计算这一放电过程的持续时间，记为深放电时长，该时长可以初步反映电池的放电性能。一个局站每天产生几十万条电压等动环实时数据。
在告警数据分析的基础上，实现对蓄电池电压变化数据的分析，提醒分公司关注那些深放电次数过多和放电时长过短的局站，核查蓄电池、油机配置、发电安排等，并进行整治。利用Spark SQL统计了一个月内抚州、赣州、吉安三分公司几十亿条动环数据，分析了其中深放电的情况如表6所示。
6 结论
本文利用HDFS+Spark技术，实验性地解决告警数据存储和分析等相关问题：一是通过数据分析，从海量告警数据中发现潜在的网络隐患；二是结合资源信息和不同专业的告警，最终为用户提供综合预警；三是转变网络监控思路和方式，通过数据汇聚、数据相关性分析、数据可视化展示，提高了网络监控效率；最后还扩展到对动环实时数据、信令数据进行分析。
从实际运行效果来看，HDFS和Spark完全可以取代传统的数据存储和计算方式，满足电信运营商主动运维的需求。
参考文献：
[1] 中国电信股份有限公司. 中国电信智能网管技术规范-总体分册[Z]. 2015.
[2] Tom white. Hadoop权威指南[M]. 4版. 南京： 东南大学出版社， 2015.
[3] RP Raji. MapReduce： Simplified Data Processing on Large Clusters[Z]. 2004.
[4] Spark. Apache Spark?[EB/OL]. [2016-11-27]. http：///.
[5] Matei Zaharia， Mosharaf Chowdhury， Tathagata Das， et al. Resilient Distributed Datasets： A Fault-Tolerant Abstraction for In-Memory Cluster Computing[J]. Usenix Conference on Networked Systems Design & Implementation， 2012，70（2）： 141-146.
[6] S鹏. Apache Spark源码剖析[M]. 北京： 电子工业出版社， 2015.
[7] Hive. Apache HiveTM[EB/OL]. [2016-11-27]. http：///.
[8] Holden Karau， Andy Konwinski， Patrick Wendell， et al. Learning Spark： Lightning-Fast Big Data Analysis[M]. Oreilly & Associates Inc， 2015.
[9] 员建厦. 基于动态存储策略的数据管理系统[J]. 无线电工程， 2014，44（11）： 52-54.
