时间：2022-09-15 07:05:20
【关键词】水文频率;计算;新型
目前，在各类水利水电工程、交通建设以及海港工程项目中，都需要指定设计标准的水文设计值，而这样的设计值往往通过对水文频率的计算获得。水文设计值结果对工程的设计、投资和工期计算都有着重要意义，所以对于设计值必须确保经过多方面分析、检查。对于水文设计值及相关的统计参数，并不是将资料进行简单的统计计算，比如某一项水利工程能够提供50年的资料，那么就需要考虑到是否有百年一遇的可能性，甚至更长。接下来就根据实际的水文频率计算理论展开探讨。
1 常规的水文频率计算方法
在以往的水文频率计算中较为普遍的计算方法为矩法与线性矩法[1]。这两种方法的具体计算过程有着很多的相似之处，都是在频率分布模型和经验频率制定后，对估计参数进行计算的大方向，其中矩法是依据多年依赖的常规矩，以变数X的r幂数来进行计算，而线性矩法则是通过近期所提出的线性矩，采取变数X的一次幂和一定的概率作为权重概率矩，然后再一次进行线性组合。
这两种方法如果出现了X估算有误差的情况出现，经过高次幂后会使得误差值更大，所以通常不会单独地使用一种方法，而是将两种方法进行结合使用，这样将误差控制在一个尽量小的范围内。在这两种方法中，对于平均值的计算结果是相同的，所以常常以计算值为主，不再作变动，但因为资料、方法等方面还存在着一些不确定因素，所以计算所得的参数值往往只能作为估计的初值，具有一定的参考价值却达不到具体应用的水准，这样对于后续设计以及施工工作的展开造成了不小的难度，这些数值往往需要通过实践和空间上的综合平衡分析后才能够得到实际应用，这样的做法即会造成计算的不准确，在耗时上也会更多。
2 水文频率具体计算方法和评选标准
上文分析了现在普遍使用的水文频率计算方法的实际情况，下面就水文频率的具体计算进行深入分析。
首先是对水文系列的统计参数进行估算，然后利用所选取的频率分布模型进行频率曲线线型计算，计算其中各种频率的设计值[2]。
目前我国的设计人员在各类水文频率的计算方法中，都是从假定水文特征量入手，假定水文特征量X服从某一水利工程的已知分部条件，拟合出X的理论频率曲线F（x），之后根据估算的参数和频率曲线，推求出可以满足设计标准的设计值xp，这样的算法在同一样板上，因为估计方法的不同，也会使计算结果不同。
其具体的计算公式如下：
其中为设计值的理论值。
这样的计算方法在实际中的应用往往会出现总体的期望破坏概率的情况，期望概率大大超过设计频率时有发生，这样就无法保证水文设计值的准确，说明了我国目前的水文频率计算方法还不能使设计对象达到安全标准。
3 水文频率新型结算方法探究
下面笔者对利用MATLAB软件对水文频率进行计算的方法开展分析，就MATLAB软件的概念，以及整个的计算流程做出深入探讨。
3.1 MATLAB概念
MATLAB是上个世纪八十年代由美国MathWorks公司推出的高效运算软件，在经过数十年的发展后，MATLAB软件具有科学计算、试图交互系统、动态系统仿真等一些先进的计算能力。其主要的优势体现在操作简单、函数运算丰富、图形符号功能全面、科学智能选择运算方法以及多种语言支持[3]。
3.2 MATLAB的参数计算
将水文频率计算引入MATLAB软件计算中，一般采用的是P-Ш型分布，水文变量X的三个统计参数为、、，另一组的参数用、、表示，通过这两组的参数来推算出水文设计值。
在MATLAB计算中，平均值的函数表示为：
P-Ш型分布属于伽马分布，数以对于P-Ш型分布的两个统计参数和在MATLAB计算中需要达到P-Ш型分布统计参数的目的，P-Ш型分布的概率函数公式为：
在该公式中说包含的和两个参数在MATLAB中的古估算函数为：
和的极大似然估计值分别与向量对应。
最终，、以及的参数值计算由以下公式得出：
为的极大似然估计值，为的极大似然估计值。
3.3 MATLAB的理论频率计算
在得出了、和得数值后可以进一步的金粟庵理论概率，其具体结算公式如下：
其中MATLAB的P-Ш型分布函数定义为：
利用上述的公式将变量的频率，在P-Ш型分布函数的逆函数中表示为：
在将上述的所有参数估计与理论频率计算带入MATLAB软件进行计算编程后，为了达到水文频率的自动化计算目标还需要进行频率计算的图解配线，这需要通过人机交互的方式来进行展开。
3.4 MATLAB的图解配线
首先，在计算机图解的配线上需要考虑到一个核心问题，即海森概率格纸的绘制[4]。下面就对海森概率格纸的绘制进行进一步的分析。
在以往的水文计算中，正态累计的分布函数表达式为：
将此表达式带入到MATLAB中可以将此变形为：
海森概率格纸在=0.01%时，横坐标设定为0;在=99.9%时，横坐标设定为1;=50%时纵坐标为0，这样就得到了概率的图解配线。
4 MATLAB的实际应用和其优点
4.1 MATLAB的计算步骤
上文对MATLAB的水文概率计算做出了具体的分析与计算机编程，下面就MATLAB的实际应用步骤进行分析。主要的计算方式分为6个步骤来完成：
步骤一：将基本的数据资料输入计算机，包括了水文的变量值x，系列长度N。
步骤二：根据上述的海森图解配线来开展海森概率格纸的绘制。
步骤三：计算出基岩频率，在概率格纸上将对应的坐标点标出。
步骤四：通过上述步骤，将、以及的参数值计算出来
步骤五：计算出理论概率，并计算相应的坐标点，绘制理论频率曲线
步骤六：进行屏幕配线，观察见过点与理论频率的拟合层度，从而确定参数的优劣情况。
4.2 MATLAB的优点
利用MATLAB软件进行编程计算使整个水文频率计算过程简捷、高效、准确，给后续工作的开展提供了良好的基础;对MATLAB软件的学习过程也比较简单，很容易上手，对于那些习惯了传统的水文频率计算方法的人来说，并不会带来太大的困难。
5 结束语
现代水利工程的建设离不开水文频率的计算，在计算方法上如今也出现了多种样式，本文通过对目前较为普遍使用的矩法与线性矩法进行分析，总结这两种方法的不足之处，同时提出了利用MATLAB软件的新型水文频率计算方法，对于今后的水文频率计算有参考作用，将水文频率计算与现代计算机结合是未来的发展趋势，如何使水文频率计算能够更加准确、可靠，从而让水文参数值完美地符合工程设计的需要，是每一个从业人员都要进行思考的，也为将来更加简捷高效进行水文频率计算工作提供新的思路。
参考文献：
[1]金光炎.水文频率分布模型的异同性与参数估计问题[J].水科学进展，2010（4）.
关键词：液压；凿岩机；设计研究；建议
中图分类号：TD421.2 文献标识码：C
Abstract： The paper analyzes various design and research literature on hydraulic rock drills in China， enumerates major arguments in some related papers， points out the imperative problems， and puts forward academic and policy-related advice on the R&D of the hydraulic rock drills industry in China.
Keywords： hydraulic； rock drill； design & research； opinion
1 我国液压凿岩机设计研究历史回顾
我国液压凿岩机的研制起步非常早，距今已近50年。我国液压凿岩机技术研究也很早，二十世纪七十年代，我国长沙矿冶研究院等就翻译了大量国外液压凿岩机专利和论文，七十年代后期，我国作者撰写了大量液压凿岩机的技术论文，直到现在，关于液压凿岩机的论文与专利仍然不断出现在各类期刊与网络上。
我国液压凿岩机设计研究方面的论文，有博士论文，硕士论文，学士论文，有期刊论文，也有网络论文，在整个液压凿岩机技术论文中，占据了很大一部分。
估计有20个高等学校（其中包括3所师范大学）和近10个研究院所的作者发表过液压凿岩机技术论文，其中大部分是关于液压凿岩机设计研究方面的论文。专利文献更是与液压凿岩机结构设计有关。可见，液压凿岩机设计研究在我国并非冷门。
我国液压凿岩机的技术论文与设计研究论文大量产出。与此形成鲜明对比，我国液压凿岩机产品发展很缓慢，独立开发的产品，很少能在市场站住脚。
2 我国液压凿岩机文献的分类
2.1 研究型论文
2.1.1 设计计算型论文
此类论文数量多。
设计计算型论文包括了仿真计算，优化设计，活塞间隙优化设计，蓄能器设计优化，换向阀设计优化，活塞缓冲设计，钎尾缓冲设计，冲击频率调节方法等等。
博士、硕士、学士的论文几乎全部属于此类。高校与研究院所作者的论文大都属于此类。此类论文数量多。
2.1.2 实验研究型论文
实验研究型论文包括实验方法介绍，实验设备和仪器介绍，实验项目介绍，测试结果与测试曲线分析，误差分析等等。
此类论文中介绍实验原理与实验方法的较多。专门介绍液压凿岩机测试结果与曲线，并且进行分析的论文较少，大多散见于设计计算论文中，作为计算正确性的证明。
高校与研究院所作者有些论文属于此类，此类论文数量不多。
2.1.3 材料热处理研究型论文
材料热处理研究型论文包括活塞、钎尾的选材与热处理，铜套的选材等等。关于活塞、钎尾的选材与热处理的论文较多，关于其他零件的论文数量少。
2.1.4 加工工艺型论文
加工工艺型论文包括凿岩机主要零件的机加工的设备，工装夹具，量具与测量方法，工序与工艺步骤，怎样保证尺寸精度与位置精度，等等，大概是技术保密的原因，此类论文数量极少。
2.2 综述型论文
2.2.1 产品品种型号综述型论文
产品品种型号综述型论文包括当时国内外液压凿岩机产品的厂家，系列，型号，主要性能参数等等。
此类论文数量较多，但大多是介绍阿特拉斯、山特维克公司的液压凿岩机产品，介绍国内产品的较少。
2.2.2 产品历史发展综述型论文
产品历史发展综述型论文包括介绍国内外液压凿岩机历史发展，产品型号推出的年代，市场情况等等， 此类论文数量不多，但大多是介绍蒙特贝塔、阿特拉斯、山特维克公司的液压凿岩机产品历史，介绍国内产品历史的较少。
2.2.3 产品技术综述型论文
产品技术综述型论文包括液压凿岩机结构类型，技术进展，等等。此类论文数量少而且内容重复较多。
2.2.4 产品市场调研分析指导型论文
产品市场调研分析指导型论文包括当时国内使用的液压凿岩机产品的主要型号，市场保有量，年销售量与销售额，备件销售量与销售额，市场的细分，将来市场的趋势预测，市场对液压凿岩机的性能与技术要求。
此类论文很重要，对于国家的产业政策和企业的产品规划指导意义很大。
此类论文数量极少，几乎没有，倒是在网络上各种市场调查公司的广告不少，只提供提纲，具体内容要付费才能提供。
由于我国缺乏市场液压凿岩机产品销售数据统计，要写出高可信度的此类论文，必然要花费巨大的人力财力，只靠个人或一个单位的力量难以完成，最好有行业组织与行政力量的参与。
2.3 使用维修型论文
2.3.1 产品使用方法型论文
产品操作方法型论文包括凿岩机冲击压力，旋转压力，推进压力的优化匹配，钻头直径与旋转速度的匹配，操作方法等等。
此类论文数量不多，有北京科技大学高澜庆教授等的论文“液压凿岩机主要工作参数对凿岩速度影响的试验研究”，有广东省水利水电第二工程局梁明华论文“液压凿岩机旋转速度与凿孔直径的关系” [1]，是凿岩机使用实践总结，指导意义更大。
2.3.1 产品维修与故障分析型论文
产品维修与故障分析型论文多数为水电工程局、铁路隧道工程局的技术人员所写。此类论文来源于生产实践，言之有物，参考价值很大，数量不多，列举如下：
煤科总院北京建井所黄园月、李耀武、郭孝先“液压凿岩机的故障分析与防治”；
广东省水利水电第二工程局何雄彬“HDl35A和COPl238ME型液压凿岩机工作原理及常见故障处理” [2]；
广东水电二局股份有限公司李同明“阿特拉斯ROCD7钻机使用中易出现的问题与改进” [3]；
李强“Atlas 1838型凿岩机六种常见故障的排除”[4]；
张兆钦“COP1238凿岩机技术特点及使用维护” [5]；
田华军“HL500型液压凿岩机的日常维护保养” [6]。
2.4 产品介绍型论文
2.4.1 产品性能介绍型论文
产品性能介绍型论文多为介绍阿特拉斯、山特维克公司产品的论文，也有少量介绍国内海卓公司产品，乐清采矿机械厂的液压凿岩机产品的论文。
2.4.2 市场应用报道型论文
国外作者关于液压凿岩机市场应用报道型论文很多，国内此类论文数量很少，有少量关于YYT26支腿式液压凿岩机工程案例，使用成本效益分析的报道。
2.5 专利文献
最近的十几年关于凿岩机械的专利有六十多件，大都是液压凿岩机械的结构设计方面的专利。尚未见到有重大影响，或产生效益的专利。
3 我国液压凿岩机技术研究的不足
3.1 凿岩机文献的五多五少现象
通过对我国液压凿岩机文献的分析，我们发现有五多五少现象。
（1） 设计论文多，实验与工艺论文少；
（2） 研究论文多，使用维修论文少；
（3） 性能介绍型论文多，市场应用报道型论文少；
（4） 历史发展综述论文多，市场调研分析指导型论文少；
（5） 专利多，影响小。
3.2 缺少学术交流讨论
我国大概有20多年没有开过全国性液压凿岩机技术方面的交流会或研讨会。
我们众多的论文作者，都是在自说自话，没有交流，没有讨论，更没有争论。论文提出的论点，基本没有人跟进，进行理论与实践两个方面的证明或证伪，或补充，只是不加评价的引用。
我怀疑，很少有人认真阅读这众多的论文，如果不是博士生，硕士生，大学生因为做毕业论文的需要，阅读的人就更少了。我国关于液压凿岩机论文的影响力实在有限。
3.3 设计理论研究与产品生产脱节
我国大概可称得上液压凿岩机设计研究论文数量大国，但是我国自行开发的液压凿岩机产品多数从市场上消失，现在还能在市场上站住脚，有一定销售量的产品，多数为1980年代技术引进的产品。这种情况至少说明我们的设计研究和产品生产是脱节的，或者说明我们的设计理论还不成熟，还要不断改进。
4 我国液压凿岩机设计理论的主要论点
4.1 关于凿岩机冲击机构设计理论
4.1.1 三段法理论
认为活塞往复运动一个周期是由活塞回程加速、回程减速、冲程加速三个阶段组成的。4.1.2 设计变量理论
将一个无量纲数称为设计变量，这个无量纲数可以是：
（1） 冲程时间与活塞运动周期之比；
（2） 活塞回程加速度与冲程加速度之比；
（3） 活塞行程与可能最大行程之比；
（4） 活塞后腔受压面积与前腔受压面积之比，很明显，第四个无量纲数最直观，最易检验。
4.1.3 优化设计理论
根据以下5个设计目标，求得一个最佳设计变量。
（1） 蓄能器容积变化的最大值最小；
（2） 蓄能器隔膜震动次数最小；
（3） 能量利用效率最高；
（4） 冲程时，最大瞬时流量最小；
（5） 回程时，最大瞬时流量最小。
4.2 关于换向阀设计理论
4.2.1 换向阀中位正开口理论
换向阀在中位时，换向阀的高压窗口与低压窗口都有微小开启，叫正开口。认为正开口有益于冲击机构性能改善。主张应该采用正开口。
4.2.2 换向阀最优行程理论
换向阀行程的选取受到阀口通流面积，换向时间与阀耗油量三重制约，存在一个优化点。
4.3 关于蓄能器设计理论
4.3.1 高压蓄能器与活塞运动的最佳耦合理论
蓄能器的蓄油与排油量不仅与蓄能器的容积有关，而且与系统的工作频率有关，如果蓄能器固有频率选择不合适，即使加大蓄能器容积，也不能增大蓄能器的蓄油与排油量。
在蓄能器结构参数固定，凿岩机进油压力确定的情况下，可以改变蓄能器充气压力来改变谐振频率，使得系统处于最佳工作状态。
这个理论是极有理论研究价值与学术价值的，是一个动态的理论，是北京科技大学首先进行研究的，可惜无人跟进，这方面的论文太少。
4.3.2 蓄能器一次振动理论
在活塞运动一个周期内，蓄能器隔膜只有一次振动。这个理论已经包含在第4.1.2小节的优化设计理论中了，这是一个静态的理论。
4.4 关于信号孔位置
4.4.1 活塞回程换向信号孔位置的计算
这个理论是由北京科技大学与中南大学提出的，这是极其重要的计算，缺了这个计算公式，无法进行凿岩机的图纸设计。
4.4.2 活塞冲程换向信号孔位置的确定
这个理论是由北京科技大学与中南大学提出的，这也是极其重要的计算。既不能换向太早，使得撞击钎尾时，活塞已经被减速，也不能换向太晚，造成活塞二次打击钎尾。
4.5 钎尾反弹缓冲动力计算
研究这个理论的有浙江大学张新等人，广东工业大学机电工程学院刘智等人，中南大学机电工程学院赵宏强等人，杨国平等人，难能可贵的是，张新做了实验研究，有实验曲线图。
4.6 其他理论
（1）活塞密封间隙与长度优化理论；
（2） 活塞空打缓冲结构计算；
（3） 冲击频率无级调节自动换挡理论。
以上3种理论均有，就不一一论述了。
5 我国液压凿岩机设计研究亟待解决的问题
5.1 换向阀的结构尺寸问题
（1） 换向阀中位负开口的优缺点分析；
（2） 最优负开口量的计算与实验验证；
（3） 换向阀最优行程理论与实验验证。
5.2 活塞与缸体结构尺寸问题
（1） 活塞前后腔面积的确定及优化；
（2） 活塞回程换向信号孔位置的确定；
（3） 最优冲程换向信号孔位置计算与实验验证。
5.3 蓄能器与活塞运动的最佳耦合理论
研究蓄能器隔膜的动态频率响应，蓄能器最佳充气压力的计算。
5.4 钎尾反弹缓冲系统动力计算与实验验证
不但研究缓冲机构本身，进行静力学计算，更要研究钎尾反弹缓冲系统的动态响应，将缓冲机构，钎具与缓冲液压油路作为一个系统，研究缓冲活塞的频率响应，缓冲液压系统是否能有效吸收钎尾反弹的能量，又是否能迅速将钎头重新抵紧岩石。
5.5 凿岩机反打系统动力计算
既要研究反打机构本身，进行静力学计算。更要将反打机构与反打液压油路作为一个系统，研究反打活塞的位移，运动速度与冲击能量。
5.6 零件的气蚀与腐蚀问题
（1）活塞前后腔气蚀问题的计算与解决方法；
（2）凿岩机壳体联接平面点蚀的原因与预防。
6 我国液压凿岩机设计研究的学术建议
6.1 反求法
对国外高端液压凿岩机结构尺寸的分析，用反求法研究凿岩机冲击机构的设计计算公式与方法。反求法基本属于归纳法，反求法需要统计多台高端液压凿岩机的结构尺寸，统计数据越多，代数计算公式越接近实际，设计指导作用越大。这方面的工作似乎还没有人去做。
6.2 验证法
将我们国内的设计理论，如用最佳设计变量理论，换向阀最优行程理论，蓄能器设计理论，分析国外高端液压凿岩机的结构尺寸、充气压力等参数，检验是否符合我们的理论。如果大致符合，则验证了我们理论的正确性，如果相差很大，则要寻找原因。
6.3 代数计算公式研究
6.3.1 液压凿岩机设计研究方法的分类
理论分析与公式计算，实验研究，计算机数值仿真是凿岩机设计研究的三驾马车。
三段法，最佳设计变量等属于理论分析计算。
理论分析与公式计算是凿岩机设计研究的基础，是实验研究与数值计算的的基础与指导。
6.3.2 液压凿岩机代数设计计算公式的建立
公式计算必须将凿岩机的物理（实际）模型简化为力学模型，再进一步简化为数学模型，再简化为代数公式。简化必须是合理的，不能与力学模型有大的矛盾和冲突。
在不断的简化中，必然有失真，这时就要用经验系数去校正。这方面我们还有许多工作要做，对每一结构类型的凿岩机，都可以研究出一套基本通用的计算公式。
6.3.3 液压凿岩机代数设计计算公式的输入输出
代数公式计算输入的是冲击能，冲击频率，冲击末速度，进油压力， 输出的是活塞质量，活塞前后受压面积，活塞行程，信号孔位置，阀的结构尺寸，阀行程，开口量，凿岩机进油流量等等。
6.4 强化实验研究
6.4.1 液压凿岩机实验的分类
在理论分析指导下的实验研究是必不可少的，我们的力学模型，数学模型是否正确，简化是否合理，都需要实验验证。数值仿真计算就更需要实验验证了。
这里所说的实验主要是凿岩机内部机理实验，也可是型式实验，而不是出厂实验。
6.4.2 液压凿岩机实验的规范
实验的各种条件，包括样机，液压系统，测量方法，仪器仪表，都应该是明确的，实验的结果应该是真实的，可重复的，重复实验的误差应该进行分析。实验的时间、地点与参加人应该注明。
6.5 数值计算研究
6.5.1 数值计算的定义
数值计算就是虚拟样机技术，因此数学模型要尽可能逼近力学模型与物理模型，数值计算需要输入的数据很多，需要详细的凿岩机图纸数据，否则不能称之为虚拟样机。
数值计算又叫动态仿真，我国在数值计算方面的论文太多了。有基于AMESim的仿真，基于Simulink的动态仿真，基于MATLAB的计算机仿真，准匀加速度法仿真计算，键合图方法仿真，等等。
6.5.2 液压凿岩机数值计算输入输出的基本要求
冲击机构数值计算中，进油流量必须是输入值，而压力是输出值。
数值计算输出结果不能仅仅是冲击能，冲击频率等，必须能输出活塞运动速度、位移，换向阀的速度、位移，活塞前后腔压力的曲线。并且能够描述出液压油的空化与气蚀现象。
数值计算的结果曲线应该用实验曲线验证，未经实验验证的仿真计算是不可信的。
数值计算的结果应该是真实的，不能弄虚作假，其他人也可以重复运行程序。
6.5.3 液压凿岩机数值计算程序的基本要求
数值计算的程序应该模块化，数值计算需要花费大量的精力与时间，不是一个人短时间能够完成的，因此数值计算的程序应该模块化，可以由一组人员分工，进行编写。
数值计算程序应该界面友好，参数的改变，应在界面中进行，而不应在程序中进行。
数值计算的模型与程序应该是持续改进的，和所有的计算机软件一样，应该有版本号。液压凿岩机数值计算程序应该通用工具化，不要搞成专用工具。不能只有作者自己用，换一个人就不能用。
6.5.3 液压凿岩机数值计算程序的难点
虚拟样机技术是一项浩繁的工程，但是并不是遥不可及。内燃机的燃烧过程，牵涉到化学，燃烧学，热学，力学，都可以做到数值仿真，并且在发动机设计中起到重要作用。冲击机构的数值仿真也一定能做到，但是这要求有一个精干的团队和一个好的实验条件。
数值计算程序都是针对某一个特定的图纸的，因此通用性较差。不要指望适合于单腔回油的凿岩机，也能适合非单腔回油的凿岩机。也不要指望适合于芯阀结构的凿岩机，也能适合套阀结构的凿岩机。这也是数值计算应用的一个难点。我们至少应该做到，对某一相同结构类型的液压凿岩机，具有通用性。
7 我国液压凿岩机设计研究的政策建议
7.1 以企业为主体，产学研结合
过去几十年，我国液压凿岩机设计研究是以高校和研究所为主体的，和企业产品开发联系较少。事实证明，我国液压凿岩机设计研究进步不大，只有以有实力的大企业为主体，厂学研结合，学校与研究所自己不搞产品生产，只为社会和企业提供知识与技术，各自发挥自己的优势，才是我国液压凿岩机技术与产品的发展道路。
7.2 以行业协会为主体，组织技术交流
与凿岩机的生产、使用有关的行业有凿岩机械气动工具、煤炭、矿山、铁路、水电、冶金等等，在1980年代，矿山机械行业组织过液压凿岩机技术交流会， 1990年代初期， 煤炭建井行业组织过液压凿岩机技术交流会。
8 结语
（1）对我国液压凿岩机研究文献进行了分类，5大类，13个小类；
（2）分析了我国液压凿岩机技术研究的不足，3个方面的不足；
（3）列举了我国液压凿岩机设计理论的主要观点，6大类，13个观点；
（4）列举了我国液压凿岩机设计中亟待解决的5大类共11个问题；
（5）提出了我国液压凿岩机设计研究的5条学术建议；
（6）提出了我国液压凿岩机设计研究的2条政策建议。
参考文献：
[1] 梁明华.液压凿岩机旋转速度与凿岩直径的关系[J].工程机械，2001，32（6）.
[2]何雄彬.HDl35A和COPl238ME型液压凿岩机工作原理及常见故障处理[J].广东水利水电，2001（8）.
[3] 李同明.阿特拉斯ROCD7钻机使用中易出现的问题与改进[J].四川水利，2004（3）.
[4] 李强.Atlas 1838型凿岩机六种常见故障的排除[J].工程机械，2005，36（3）.
Abstract： Numerical calculation and measurement of prestressed tendon is the key point of construction quality control of prestressed structure， and tension numerical calculation is the foundation and the most important. Whether the tension numerical calculation is accurate or not plays a decisive role in the quality and service life of the prestressed component. Tensioning numerical calculation involves the content， more complicated and difficult steps. Taking the tension numerical calculation of negative prestressed tendons （low relaxation strand） of prestressed reinforced concrete box beam on A2 highway as an example， this paper introduces the numerical calculation， calculation method， calculation steps， main points and points of attention of the tension moment of negative bending moment prestressed reinforcement of box beam roof in bridge engineering.
P键词：箱梁；顶板；负弯矩；预应力筋；张拉；数值
Key words： box beam；roof；negative bending moment；prestressing tendons；tension；numerical
中图分类号：U445.4 文献标识码：A 文章编号：1006-4311（2017）14-0089-03
0 引言
预应力混凝土构件广泛应用在各种建筑工程中。预应力筋的张拉是预应力混凝土构件施工中的核心技术，对预应力混凝土构件的质量起着决定性的作用。预应力筋的布设位置、形式多样，施工难度各不相同，但施工精度要求高、注意事项多，尤其是负弯矩预应力筋的施工更有许多不同的特别要注意之处[1]。先简支后连续箱梁被广泛应用到大中型桥梁中，它是一种介于简支结构和传统连续结构之间的桥梁上部结构形式[2]，但桥梁结构中负弯矩张拉不被大家重视，影响了桥梁的安全和使用寿命[3]。预应力筋的张拉数值计算是预应力筋张拉质量控制的关键和重点。准确计算负弯矩预应力筋理论伸长量是控制张拉施工质量的基础[4]。本文以某高速公路A2标预应力钢筋混凝土箱梁顶板负弯矩预应力筋（低松弛钢绞线）的张拉数值计算为例，介绍了桥梁工程中箱梁顶板负弯矩预应力筋的张拉数值计算内容、计算方法、计算步骤、要点和注意事项。
1 工程概况
某高速公路A2标段主线全长7.000公里（起讫桩号为K210+000～K217+000），共有大桥949.6米/5座，大桥上部结构选用标准跨径为25米、30米的先简支后连续装配式预应力混凝土连续箱梁。
箱梁顶板负弯矩预应力筋设计采用低松弛高强度钢绞线，采用BM15-5、BM15-4扁锚体系。预应力筋孔道采用金属波纹扁管。负弯矩束布设及材料技术参数详见表1。
2 计算准备工作
计算准备工作是为负弯矩束相关数值计算收集必须的数据和依据。主要包括对预应力筋（钢绞线）取样检测、张拉设备检定、测定影响预应力筋计算长度的结构尺寸等。
2.1 预应力筋取样检测
预应力筋取样检测主要是通过对购入施工现场的钢绞线按规定方法和频率制取试样。送有资质单位检测，对材料的质量进行验证、复核，对比产品出厂合格证的数值；同时为预应力筋张拉相关数据计算提供钢绞线准确的必须数据（如面积AP、弹性模量Ep等）。
2.2 张拉设备标定校验
负弯矩预应力筋张拉采用电动高压油泵（含压力表）、穿心式单根千斤顶各2台。压力表表面最大读数为张拉力的1.5～2.0倍，标定精度应不低于1.0级[5]。
在张拉前张拉机具应在有相应资质的检验计量单位进行仪器仪表的标定，为张拉提供依据[6]。千斤顶与压力表必须配套校验，以确定张拉力与压力表之间的关系曲线。根据大量的试验数据分析研究，张拉力与压力表读数之间具有一元线性方程关系，可以用一元线性回归方程表示如下：
Y=a×X+b （1）
式中：Y：压力表读数（MPa）；X：张拉控制力（kN）；a为系数（无单位）；b：修正值（MPa）。a、b的值因压力表、千斤顶及组合不同而不同。
通过式（1）可以计算负弯矩预应力筋张拉表见值，进行分级测量、计算预应力筋伸长值。
2.3 测定影响预应力筋计算长度的结构尺寸
张拉预应力筋必须安装锚具和千斤顶，锚具、千斤顶占用一定的预应力筋长度，其中千斤顶工具夹片到锚垫板之间的预应力筋也被张拉、伸长，对预应力筋的伸长值计算和测量直接产生影响，故必须测量这部分预应力筋的长度，也即工作锚具的厚度B、千斤顶顶压器外端与工具夹片内端之间距离D（如图1）的和。
3 张拉力计算
3.1 预应力筋张拉控制应力（бcon）
根据设计规定，钢绞线张拉锚下控制应力为бcon =0.75fpk=0.75×1860=1395MPa。由于施工工艺、获取质量控制数据等的需要，张拉必须分行程进行，бcon是唯一的依据，如张拉需要行程为10%бcon、20%бcon、100%бcon，相应的锚下应力为139.5MPa、279.0MPa、1395MPa。据此，可进行预应力筋张拉端的张拉力P、预应力筋张拉理论伸长量植的计算。
3.2 预应力筋张拉端的张拉力（P）
预应力筋张拉端的张拉力P（N）可由下式计算求出：
P=бcon×Aρ×n×c （2）
式中：бcon：预应力筋张拉控制应力（Mpa）；Aρ：预应力筋的截面面积（mm2）；n：同时张拉预应力筋的根数，此处n=1.0；c：超张拉系数，不超张拉时为1.0。
预应力筋张拉端的张拉力P不是固定数值，张拉力P随每批钢绞线实际送检测得或钢绞线生产厂家提供的截面面积Aρ值不同而不同。如某批钢绞线检测得Aρ=140mm2则计算P得：P=1395×140×1×1.0=195300（N）。
各个阶段的张拉控制力依据设计给定锚下张拉控制应力带入式（2）可分别求得（见表2）。
3.3 张拉压力表读数计算
张拉压力表读数是进行张拉施工控制的依据。压力表读数用各阶段张拉控制力数值带入式（1）求得。必须注意每个压力表读数都必须由张拉力带入式（1）求得，不得以某次计算得到的压力表读数再按张拉行程按比例计算其他行程压力表读数。因为关系方程中的相关系数和修正值是关键所在。尽管两者计算结果相差有时非常小。这是经常被忽视的一个问题。现以一次实际标定得关系方程来计算压力表读数并作比较以更直观地说明此问题，计算及比较结果详见表3所示。
4 预应力筋计算长度（L）
预应力筋计算长度（L）既不是设计图纸中包含工作长度的束长，也不是锚固端之间的束长。是张拉受力时参与伸长的分段预应力筋长度之和，对于两端张拉的预应力筋来说等于两张拉端千斤顶工具夹片内端之间预应力筋的长度。包括锚固端之间的束长（T）、工作锚具厚度（B）、千斤顶顶压器外端与工具夹片内端之间距离（D），也即预应力筋计算长度L=T+2B+2D。此处预应力筋为直线且两端张拉，可以以预应力筋的中点作为计算截面，即以（T+2B+2D）/2长度值带入相关公式计算每端理论伸长值，然后合计求和取得整束预应力筋的理论伸长值。
5 理论伸长值（ΔL）计算
5.1 计算依据
按文献[5]，预应力筋的理论伸长值ΔL按下式分段计算叠加。
5.2 伸长值计算
由上述式（2）、（3）及设计图纸相关数据，分别计算得30m、25m箱梁顶板负弯矩束张拉理论伸长值。现列出30米箱梁顶板负弯矩束张拉理论伸长值如表5所示。
5.3 伸长值控制范围
根据文献[5]规定张拉采取“双控制”，预应力筋采用应力控制方法张拉，以伸长值进行校核，若理论伸长值和实际伸长值之间的差超过±6%时应立即停止张拉，找到原因并解决后，方可继续张拉[7]。预应力筋伸长值的控制范围为：
6 张拉实际伸长值计算
根据文献[5]张拉程序为0初应力（量伸长量初读数L1）бcon持荷5分钟（量伸长量终读数L3）锚固。实际施工时，在施加应力前，钢绞线处于松曲状态，此时不能作为钢绞线伸长值的起点，而钢绞线恰好被作用至绷紧状态的应力较小且不易掌握。为了准确测算预应力筋实际发生的所有伸长值，文献[5]提出钢绞线在初应力作用下的非弹性伸长值可以用公式求出或采用相邻阶段推算法获得。根据工程实践检验，大都采用相邻阶段推算法。此法具有简便易行、数据准确可靠等优点。具体的操作方法如下：
按010%бcon20%бcon100%бcon三个阶段，共量出三个伸长值，即：
第一阶段010%бcon，量伸长量初读数L1；
第二阶段10%con20%бcon，量伸长值为L2；
第三阶段20%con100%бcon，量伸长量终读数为L3。
由第一、二阶段可以推算出钢绞线非弹性伸长值ΔL''=L2-L1。张拉实际伸长值为：
ΔLs=（L3-L1）+ΔL''=（L3-L1）+（L2-L1）=L3+L2-2L1。
实际伸长值理论上还应考虑锚环的压缩量，由于此值很小而一般予以忽略不计[8]。
7 结论
影响预应力筋张拉数值和施工质量的因素很多，如孔道平顺、位置正确等，负弯矩管道精确定位可有效保障后续施工顺利进行和工程质量[9]。
张拉力、压力表读数、钢绞线理论伸长值等的计算必须认真、准确、严谨，要有专人计算、另外有人复核。张拉有专人记录、测量伸长值，现场实测值与理论值进行比较，对张拉应力进行校核[7]。张拉应力及压力表读数计算错误可能导至断丝、伸长率超出控制范围等事故。检查张拉应力及压力表读数是否计算错误直至重新送检钢绞线。
必须在预应力筋的张拉设备标定及原材料检验、张拉理论伸长值计算与伸长值范围确定、以及现场张拉控制与记录等各个环节加强全面质量意识，认真做好自检、互检、验收。务必做到计算测量精确、操作方法正确和施工工艺控制严格。
箱梁顶板负弯矩预应力筋的张拉数值计算内容、计算方法、计算步骤、要点和注意事项等同样适用于其他形式梁板顶板负弯矩预应力筋、以及正弯矩预应力筋，只是预应力筋数量、张拉设备型号不同相关参数须做相应调整。
参考文献：
[1]刘小力，叶松，熊明祥，等.箱梁顶板负弯矩预应力筋扁形施工[J].安徽建筑，2010，17（3）：56-57.
[2]田耀娟.调整张拉槽口确保箱梁负弯矩钢束准确定位施工[J].山西建筑，2014（28）：176-177.
[3]春星.箱梁负弯矩预应力筋的张拉施工控制[J].中国公路，2013（10）：124-125.
[4]刘瑞定.桥梁负弯矩预应力筋理论伸长量及油表读数计算的研讨[J].四川水泥，2015（4）：213.
[5]中交第一公路工程局有限公司.JTG/TF50-2011公路桥涵施工技术规范[S].北京：人民交通出版社，2011.
[6]张月明.探讨桥梁预应力T梁的施工技术与质量控制[J].黑龙江交通科技，2014，4（242）：86，88.
[7]熊明祥.龙潭湾大桥施工项目管理[D].西安：长安大学，2014.
关键词：计算思维；原码；运算器；计算机组成原理；计算机教育
中图分类号：TP311 文献标识码：A 文章编号：1009-3044（2015）28-0077-02
Again Discuss Design idea of the original code on the Computational Thinking
SUN Li
（Chenxian College of Southeast University， Nanjing 210088， China）
Abstract： In accordance with the Computational Thinking process，Design idea of the original code is discussed.Why appear the original code in the computer component principle.The original code can solve what arithmetic problems. Operational characteristics of the original code is related and analyzed from the aspects of addition， subtraction， multiplication and division.The final，Summary the of design idea of the original code.on the Computational Thinking process.
Key words： computational thinking； The original code； ALU； computer component principle； computer education
在《基于计算思维的计算机组成原理课程研究与教学》项目进行中，要围绕计算思维培养的思路选择适用的计算思维教学案例，在运算器设计这一章，我们选取原码作为计算思维教学案例之一。
计算思维是我们在遇到问题时，考虑用计算工具在问题可解的情况下，思考如何运用计算语言描述解决问题的过程。这个思考过程的结果是可在计算工具运行并由计算语言表述的有先后顺序的序列。
计算机组成原理课程全书讲述的核心问题是计算机是怎么设计的，第二章为了解决运算器设计问题，出现了原码。真值转换为机器数表示是运算器设计的第一步，真值转换为机器数不是编码完成即可，而是为用转换后的机器数编码在运算器完成各类运算且能简化运算器设计，真值的定义为带“+”“-”符号的数，机器数的定义是符号数字化的数（特别指二进制数据），那么真值转换为机器数表示，只要把符号数字化、其余位不变化就能最简单地完成机器数编码，我们把这种编码就叫原码，所以在机器数编码设计中原码的优势的是非常明显的，原码应该是最重要的一种编码，那原码把握住这个机会了吗？如果原码能满足运算要求，可能别的编码就不会再有机会出现。
下面，以最简单的核心运算：加减乘除为例，讨论一下原码运算性能。（注：本文数据如未作特殊说明，皆为二进制数据）。
1 原码加减法
加法运算是一个重要的基本运算器设计单元，后继很多运算器设计都与其有关。
原码加减法运算可参考的原型只有使用几千年的加减法竖式，但这种加法的特点是符号和数值分开计算，不能符号和数值一起参加运算。
要想进行加法要先讨论参加运算的两个数值的符号，如果符号相同做加法，如果异号做减法，那就要完成加法器设计必先完成减法器设计。
再看减法，要想进行减法要先讨论参加运算的两个数值的符号，如果符号相同做减法，如果异号做加法，那就要完成减法器设计必先完成加法器设计。
综上分析，加法器和减法器的设计，二者出现了互为条件的死锁，原码加减法运算器设计思路到此止步。
2 原码乘除法
原码乘除法的设计原型是：笔―纸乘除法方法，然后加法器器件为基础，增加一些硬件从而实现乘除法运算器设计。因此原码进行乘除法运算要相对简便些，这是因为乘积的符号可以根据被乘数与乘数的符号来确定，其规则是：两数相同，乘积符号为正；两数异号，乘积符号为负，除法类同。因此，原码乘除法实际上是两个无符号数相乘，乘积结果再加上符号位。原码乘除法的设计思路可行，此处我们只以原码乘法为例说明之，因篇幅原因除法略。
2.1原码乘法原型
在笔―纸乘法方法中，乘积由部分积相加得到。部分积的个数由乘数B的位数决定，当乘积B的个位数值是1时，部分积为被乘数左移后的值，当乘积B的个位数值是0时，部分积为0。
以上是人工乘法过程，为了便于计算机实现，部分积相加得出乘积的过程可以改为：每得出一个部分积就进行一次加法，并将被乘数左移改为部分积右移，这样只要用和被乘数、乘数位数一样的n位加法器就能完成2n位加法运算了。
2.2原码一位乘法
由于原码乘法被转化为一系列的加法和右移操作，n位相乘需要n次加法周期，原码一位乘法具体操作是：每次通过乘数B的末位值来判断是否需要加法，当乘积B的末位数值是1时，执行加被乘数的操作，当乘积B的末位数值是0时，不执行加被乘数的操作；然后部分积直接右移一位（高位用进位位填充），移出的一位放入乘数B中的最高位，同时B也右移一位，把乘数移出（判断过的）一位，即刚刚比较过的B的末位移出丢弃，然后对新的B的末位进行判断，重复以上过程直到处理完乘数B的所有位。
原码一位乘法中，每次判断加或不加之后，部分积都要进行右移操作。为了存放部分积，部分积寄存器P初值得为0，结果值即乘积高n位存放在部分积寄存器中，低n位存放在乘数寄存器B中。
2.3原码二位乘法
2.3.1原码二位乘法思想
为了提高乘法的运算速度，可以采用二位同时乘或多位同时乘的方法，二位乘法就是每次处理同时考虑乘数二位，根据它们的不同组合（00、01、10、11）一步求出两位的部分积，只需增加少量的逻辑线路，就可以将乘法的速度提高一倍，从而提高乘法的速度。
2.3.2对于原码两位乘位算法的几个思考
1）部分积与被乘数采用三个符号位，这样表示的原因是：当加法处理+X或-X操作时，可能会影响到中间结果符号位的最低一位，这样最低一位符号位就不能代表结果真正的符号位，当加法处理+2X操作时，可能会影响到中间结果符号位的中间一位，这样中间一位符号位就不能代表结果真正的符号位，而无论做怎样的操作都不能影响到中间结果符号位的最高位，所以中间结果符号位的最高位能正确标识结果的符号位，我们在每次处理二位后都要进行移位操作，当同步右移二位时，符号位移出的空位是用中间结果符号位填充的，当部分积与被乘数采用三个符号位时，就能保证每次移位操作的简单与正确。
2）在实际运算中CX操作以+[CX]补完成，这是因为此时运算为乘数与被乘数的数值部分的运算，而正数的原码与补码相同，所以可以看成是两个正数原码的运算，也可以看成是两个正数补码的运算，结果依然为正数，所以可以把运算中CX操作以+[CX]补完成；
3）原码二位乘法算法中数值位N的奇偶问题的讨论
有关奇偶问题的处理，不同书中方法不一，这个奇偶处理算法在设计运算器时还是有点复杂的，也是教学过程中同学容易出错的地方。个人觉得如果奇偶处理放在运算进行之前，会简化运算器设计的，就是当发现数值位N为奇数时，我们在不影响数据实际值的情况下直接增加一位变成偶数位参加运算就可以了，最后只进行T是否是1的判断和处理就结束了。
3 结论
原码在运算器设计中，最基本的加法运算器设计不能实现，原码乘除法运算器的设计以加法器器件为基础，原码乘除法运算器性能尚可。原码加减运算器留下的设计空白是必须填补的，这给了设计新的机器数编码的机会，这也是原码进行运算器设计的缺陷。最后总结原码设计的计算思维过程见图一。
参考文献：
[1] 陈国良.计算思维导论[M].北京：高等教育出版社，2012.
关键词：微分方程 数值解法 双语教学 有限差分法
中图分类号：G642 文献标识码：A 文章编号：1673-9795（2014）04（b）-0060-02
微分方程数值解法就主要研究如何通过离散算法将连续形式的微分方程转化为有限维问题，如代数方程组，进而来求解其近似解[1]。它以逼近论、数值代数等学科为基础，探讨有效的微分方程数值解法。主要包括求解区域网格划分、离散方程的建立、方程性能分析、近似解收敛性分析等环节。探索微分方程数值解法是有积极而重要的科学意义的，这是因为：（1）在实际应用中，我们只关心方程在某个范围内对应于某些特定的自变量的解的取值或近似值；（2）绝大多数情况下，无法找到方程的解析解，即使解析解存在也不一定能表示为显式解。微分方程数值解法在计算物理、化学、流体力学航空航天等很多工程领域具有广泛的应用。目前已发展成为一门计算技术学科，其核心理论内容也成为高校计算数学和应用数学等专业的核心基础专业课程之一[2]。
1 双语教学的必要性
现代社会的高素质专业人才不仅要具备扎实的专业知识，还须具备流利地应用英语进行沟通和交流的能力。双语教学是教育部积极倡导的一种课堂教学模式，在2001年公布的《关于加强高等学校本科教学工作提高教学质量的若干意见》中指出要“积极推动使用英语等外语进行教学”[3]，主要是在课堂教学过程中采用母语和以英文为代表的多种语言教学。其目的就是为了跟上经济全球化的步伐和迎接科技革命的挑战。对高新技术领域中的诸如信息技术、生物技术、金融、法律等专业，力争三年内，外语教学课程达到所开课程的5%～10%[3]。2005年，在教育部颁布的《关于进一步加强高等学校本科教学工作的若干意见》中进一步要求高校要“以大学英语教学改革为突破口，提高大学生的国际交流与合作能力”，进一步明确了要“提高双语教学课程的质量并扩大双语教学的课堂数量”[4]。可见，国家教育部门对高校采用双语教学给予了相当的重视和期望。
微分方程数值解法既有数学上严密的逻辑性、独特的理论结构体系，又在各种工程计算中有着重要的应用，因此是联系纯数学理论和工程应用的桥梁和纽带。另一方面，很多数值计算软件开发平台和帮助文件都是用英文开发的，而数值微分各种理论算法又可以直接用伪代码表示，如何对数学专业英语很娴熟，那么应用这些数值计算软件就得心应手，亦可以熟练与国际同行交流。再者，该课程一般在高年级开设，通过大学两年的英语教学积累，大部分同学已经达到了大学英语四级水平，可以较容易的阅读数学专业文献。同时，高年级的同学对数学基础理论知识，如数学分析、高等代数、数值分析、常微分方程、偏微分方程等有了较好的掌握，继续接受方程的数值解的概念和理论是顺理成章的事情。因此，无论是实际工程需要还是学生自身素质，对微分方程数值解进行双语教学都是可行的、必须的。本文拟结合重庆理工大学信息与计算科学专业课程的设置，对微分方程数值解法的双语教学模式进行探讨，以寻求适合我校数学专业课程的双语教学模式。
2 课堂教学模式探讨和上机实验
课堂理论教学是学习《微分方程数值解法》的主要方式，务必引起足够重视。大学教育离不开课堂教学，而课堂教学离不开讲授。理论是科学的基础，理论是创新的基石，只有掌握了理论结果和相关概念，才能进一步有所创新。在教材选取上面，我们选取了李荣华、刘播等主编的《微分方程数值解法》[1]作为主要参考教材。选用Arieh Iserles主编的《A first course in the numerical analysis of differential equations》为主要辅助教材（网站下载）[5]。该英文版教材作者英文功底深厚，相应的概念、定理、定义表达简洁容易理解，阅读该教材有种阅读英文科技小品的感受，对提高学生的英文水平非常有帮助。授课采用计算机多媒体辅助教学。首先让学生阅读中文教材以熟悉所学概念定理等内容，同时对所学的算法知识、理论知识也有一定的了解。上课PPT采用全英文书写，采用中文授课。当用英文表示所学概念时，老师给出其相应的中文含义，由于学生先期对该概念有了一定的预习，那么接受英文概念则不是太困难。只要教师及时对这个英文专业词汇进行解释，学习过程中则不会存在太大的困难。
英文概念词汇有助于学生获悉如何用英语表达我们常见的数学概念和定义定理等内容。同时也有助于学生进一步理解数学概念内涵和激发学生学习英语的热情。例如，第一章中对于常微分方程的向量场的概念，如果采用英文Vector field则更容易理解。对于Euler 法的重要基础地位，英文教材描述颇有味道：In a deep and profound sense， all the fancy multi-step and Runge-Kutta schemes are nothing but a generalization of the basic paradigm （yn+1=yn+hf（tn，yn），n=0，1，…）[5]。这句话既强调了Euler迭代公式的基础地位，进一步说明多步法（multi-step）和龙格-库塔法（Runge-Kutta）的新奇性和实用性。虽然Runge-Kutta法是Euler法的推广，但是其理论推导在短时间内不容易弄清楚，主要困难在于需要学生了解数值积分的代数精度概念、误差收敛阶，多元函数的Taylor展开，即如何灵活应用未知函数y（t）的各阶导数与右端函数f（t，y）的偏导数之间关系来对参数ki进行Taylor级数展开。
在实践教学方面，教育部对高校本科教学工作的若干意见中重点强调了要进一步加强实践教学，注重学生创新精神和实践能力的培养，切实提高大学生的实践能力，切实加强实验等实践教学环节[3～4]。所以，微分方程数值解法的计算式实验环节也需引起足够重视。通过计算机编程，有助于学生更好的理清各种算法的运算步骤，深入理解算法内涵，对掌握微分方程数值解法的学习方法能起到重要的作用。
3 存在的问题和总结
在教学伊始，学生的学习积极性并不高涨。主要是因为同学们接受新鲜事物有一个过程，心底里认为使用英语教学没有必要，课前预习不充分，不愿意花精力去记忆消化英文概念和理解英文句法。为达到较理想的教学效果，还需要学生在思想上高度重视。国外原版英文教材价格太贵，并且教材内容比我们教学大纲要多，我们必须有针对性地选择重点章节讲解，并不能面面俱到。受师资水平和学生英文水平限制，我们目前上课还无法使用英语口语教学。一是授课教师没有在国外高校进行过改门课程的讲授，口语不纯正；二是学生的专业数学概念词汇少和听力理解。这就要求在平时教学过程中，师生都要有目的的加强练习，及时发现问题并提出可行的解决方案并不断积累经验。
参考文献
[1] 黄振侃.数值计算-微分方程数值解[M].北京工业大学出版社，2006.
[2] 李荣华，刘播.微分方程数值解法[M]. 高等教育出版社，2009.
[3] 教育部.关于加强高等学校本科教学工作提高教学质量的若干意见[Z].2001.
【关键词】数值分析；教学内容；分析设计
中图分类号：G64 文献标识码：A 文章编号：1006-0278（2014）07-215-01
《信息与计算科学专业教学规范（试行稿）（见《大学数学》第19卷1期2003）》中规定了信息与计算科学专业主要的培养目标是培养具有良好的数学基础和数学思维能力，掌握信息与计算科学的基本理论、方法和技能，受到系统的科学研究训练，能解决信息科学与工程技术实际问题的高级专门人才。数值分析作为已经成为物理学、力学、计算机应用、航空航天专业学生的学习课程，提高数值分析课程的教学质量，已成为当前高校教育改革的焦点。
一、教学过程中的教学方法探索
（一）突显专业特色
由于各专业在专业内涵培养目标、培养要求、专业发展方向、课程体系等多个方面都存在本质的区别，就必须对数值分析这门课程的课程性质、教学内容、教学目标、教学要求、教学课时数等方面有所考虑.因此，在深化信息与计算科学专业数值分析课程教学改革的同时，不能照搬照套其他专业数值分析教学模式，而应从更高的角度、更宽的视野、更高的标准来探索信息与计算科学专业数值分析课程教学改革，使之自成体系、独具一格。
（二）深化课程教学改革
深化数值分析课教学改革，须改变重视理论教学、忽视实验环节的现状、做到理论与实验并重，实现齐头并进、协调发展，以使学生加深理解数值分析理论课的教学内容，更好地掌握其中的数学原理和方法，更有利于提高学生算法设计能力。同时加强信息与计算科学专业数值分析实验课教学改革，构建有专业特色的实验课教学模式，是促进实现这一培养要求的重要途径之一。
（三）消除教学条件限制和教学模式的影响
数值分析这门课程对教学条件和教学模式、教学水平有很高的要求，所以受教学条件限制和教师教学模式的影响很大，使该课程的教学出现重理论轻实践的弊端。学生课上听完老师讲解的算法推导、误差分析等理论后，课下很少进行实际计算，学生学完该课程后只记得过程有些枯燥乏味的复杂公式，而不会学以致用、解决实际问题，使教学没有达到预期目标。在教学中，可充分利用投影进行程序使用后演示。
讲解中设置一些有层次的问题，引导学生积极思考，加入课堂学习中来，课后注重辅导，及时解决学生中存在的问题。
二、数值分析课程的学习方法
数值分析是一门理论与实践相结合应用很强的学科，学习《数值分析》的方法是与其他学科不同，故方法得有所改变，才会有好的学习效果。
（一）确定目标
在学习数值分析这门课程的时候，首先要明确学习目的。在介绍课程时要明确要求学生为自己定下一个切实可行的目标，并为这个目标的实现制定短期、中期、长期的学习计划，并实时进行对照，找到其中的优点或者是计划中存在的问题，在制定计划的过程中做到心中有数，量力而行。该课程是学习MATLAB软件的最佳平台，也是数学建模必备的知识之一。
（二）了解内容
上课之前，应引导学生粗略地了解一下课本中出现的内容。目前，在数学教学中流行的所谓“任务驱动”学习方法，就是指先有结果，再研究实施策略的学习方法.在任务驱动教学中，打破了常规教学方法中由浅入深的基本顺序，每一章节的知识点都是通过几个有代表性的案例来学习的，甚至包括认识程序。让学生先体会到效果，从而增加学习兴趣。用这种方法来学习数值分析，尤其是一些视窗界面的计算程序，往往可以达到事半功倍的效果。
（三）正确利用书籍
与数值分析这门课程有关的参考资料书有很多，在学习之前，由于数值分析这门课程偏重于应用，故而，在学习数值分析的理论知识后，一定要在理论知识的指导下，进行必要的实践和练习，使得学习的知识能有的放矢，并知道具体的知识点运用的具体方向和社会运用的价值，要实现这些，就必须正确利用有关的书籍，不仅丰富学生的知识库，扩宽学生的视野和知识面，而且对学生以后的就业和理论知识的深造有很大益处。
（四）注重实践
数值分析在模仿别人练习的基础上有想法，并尽力去寻求解决办法。数值分析中的误差是衡量我们的计算是否有效的标准.如果在误差允许的范围内，则算法是有效的，否则就是一个无效的问题求解。大多数人认为数值分析仅学习理论知识，很难和程序设计结合起来。通过上面的论述，大家已经能够初步地认识到这个学科是与程序设计紧密联系才能够体现它的重要性。
参考文献：
[1]华东师范大学数学系.数学分析（上册）第三版[M].北京：高等教育出版社，2009：24-28.
关键词：并行计算；OpenMP ；OpenCL；GPU；MPI；自然对数
中图分类号：TP312 文献标识码：A 文章编号：1009-3044（2013）14-3415-05
1 概述
e作为数学常数，是自然对数函数的底数。它就像圆周率π和虚数单位i，是数学中最重要的常数之一，也是第一个被获证为超越数的非故意构造的数。自然对数底e是数学分析中使用非常广泛的无理数之一，它在金融、数学等领域[1-2]有着重要的应用。
对于自然对数底的计算，尚晓明[3]介绍了三种在教学中的常见的计算方法，而张新仁，徐化忠[4]介绍了计算机串行计算自然对数底的方法和实现。虽然自然对数底的值已经内置在很多的计算软件中，但是对于自然对数底的并行计算方法缺少相关的介绍。对自然对数底进行并行结构分析，不仅可以了解这一重要的数学常数，还可以分析不同的并行计算模式的优缺点。
2 e值计算算法分析
对于自然对数底e的计算公式有很多，在文中只讨论一种计算公式（1）。
2.1 串行算法分析
2.2 并行算法分析
考虑到式（1）的结构特点，它是由有限个多项式构成，在计算上有一定的独立性，结合多CPU 的架构，在这里讨论两种计算e值的并行算法。
2.2.1 多项式分段并行算法
循环并行化是使用OpemMP来并行化程序的重要部分，它是并行区域编程的一个特例。使用parallel for编译指导语句能将for循环中的工作分配到一个线程组中，而每个线程组中的每一个线程将完成循环中的一部分。对于算法3.2.2，还可以利用openMP 的编译指导语parallel for 进行简化。这样可以简化程序，方便串行程序到并行程序的转换。
3.3 MPI实现
3.4 OpenCL实现
随着GPGPU在超级计算环境中的日益普及，它在应用领域表现出优秀的性能功效比和性价比优势。OpenCL作为GPGPU的并行语言的一种，OpenCL [6]是一个为异构平台编写程序的框架，此异构平台可由CPU，GPU 或其他类型的处理器组成。OpenCL基于kerels函数运行，执行时也是拷贝复制执行模式。参考上面三种语言的并行实现计算自然对数底e，可以很快的得到所对应Kernel函数，在kernels中要注意合理利用线程ID号来分配多项式的起始位置。
4 分析与讨论
上文中，我们分析了自然对数底e 的并行计算结构和具体语言实现。为了分析四种并行计算语言的优缺点，我们在特定的设备（AMD 2*1）上进行数值试验。由于自然对数底E 的计算精度不是关注的重点，在这里不考虑计算机的精度截断（即计算机精度不够，导致无用计算）。
4.1 两种并行算法比较
5 结论
本文通过自然对数底e值并行计算的四种语言实现，分析了四种并行语言的实现特点和计算优势。由实验结果可以看出，CPU 和GPU 的并行计算还是有一定差异的，需要根据各自的特点进行协调配合，可达到更好的计算效率。
参考文献：
[1] 陈仁政.不可思议的e[M].北京：科学出版社，2005：1-88.
[2] BRIAN J M.e：数中之大师[J].余敏安，译.数学译林.2007（3）：213-226.
[3] 尚晓明.超越数π与e的计算方法及其应用[J].焦作大学学报，2012（1）：73-74.
[4] 张新仁，徐化忠.自然对数的近似计算[J].山东电大学报，2012（3）：64.
关键词：数学模型；缓坡方程；实用性
中图分类号：TQ018文献标识码： A 文章编号：
近岸水域对于当今人类活动具有特别重要的意义。外海波浪由外海传入近岸浅水地区时，受水深、地形、底摩擦、障碍物、水流等因素的影响，会发生变形、折射、绕射、反射和破碎等各种波浪变形现象。波浪研究方法主要包括理论分析、模型实验、现场观测资料研究以及数值模拟的方法等。许多理论如Stokes波理论、浅水波浪理论、波浪折射、辐射和绕射理论以及高阶Boussinesq方程理论得到的空前发展。随着计算技术的迅速发展，利用数学模拟探讨波浪传播变形的规律已成为主要研究手段。
缓坡方程基于缓坡条件的假设，是势波理论三维Laplace方程的一种简化近似形式，它将三维问题转化为二维问题，能够模拟出波浪的折射、绕射和反射联合作用等现象。Eckart(1952)第一个提出了在浅水中传播的缓坡方程。Berkhoff(1972)[1]在缓坡假设下，根据势波理论用小参数展开的方法推导出适合不同水深的反映波浪折射绕射联合作用的方程，即著名的缓坡方程。许多学者针对缓坡方程的缺陷，提出各种不同的类型的改进和推广。其中Dingemans(1997) [2]对海底作用项进行修正。
本文在以往参与实际工程项目以及前人做的各项工作及得出的各项成果的基础上，将进行以下几个方面的研究工作：
1.以Dingemans(1997)推导的方程为基础建立缓坡方程，建立综合考虑海底陡坡影响、底摩擦、波浪破碎、非线性色散关系的综合缓坡方程求解模型。
2.使用ADI 法得到无条件稳定的差分方法，在处理边界条件时，引入虚拟边界层。
3.对比不同试验地形上的物理模型试验结果和数值计算结果，验证模型应用性。
4.将模型应用于实际工程中，为设计提供参数。
1.理论模式
1.1控制方程
依照线性波理论，在推导缓坡方程过程中一个常用的假设就是波速势在垂直方向上的为，沿着坐标（x, y）的变化通过h很弱的表现出来，根据该假设，可得到如下方程：
（1）
Dingemans(1997)对海底作用项进行修正，为简单起见，这里直接给出Dingemans改进后的方程简洁表达式。方程的具体表达式和推导可以参考相关文献【2】。
（2）
式中，为势函数，为相速，为群速，k为波数，线性色散方程为，h为水深。
1.2缓坡方程的扩展
以上给出的方程在方程的实践应用中往往还需要考虑一些影响因素, 如海底陡坡影响、底摩擦、波浪破碎、非线性色散关系等。
1.2.1抛物化缓坡方程
在速度势中引入慢变时间变量，对单色波有：
（3）
式中：为波陡，或，对于陡变地形有。
代入方程（2）中，可得新的方程为：
（4）
1.2.2考虑能量扩散
在波浪传播过程中由于风能、底摩擦以及波浪破碎的影响，通常在缓坡方程中添加能量耗散项，即在时间关联型方程中加入项（Pan，2000[3]）。
根据Kirby(1984)势函数结构，分解时间关联型缓坡方程，从分解方程实部可以得到：
（5）
（6）
（7）
式中：波能，，为波幅，、为波浪破碎项和底摩擦项。
(1)底摩擦
（8）
式中：为波数，为水深函数， 为波能，为角频率，为底摩擦系数，它的取值通常为0.01-0.02之间.
(2)波浪破碎影响
（9）
（10）
式中：取1.0，为调整波高。
对于波浪破碎指标，本文根据工程实际需要，根据规范选取合适的计算方式。
1.2.3色散关系
在计算波浪相关数值中，以水深和波浪周期来确定波长是需要考虑的问题之一。波长一般通过色散关系来求解，线性波色散关系一般表示如下：
（11）
式中：为角频率，为波数，为水深函数。
在实际计算过程中，式（12）的精度较差，又因为需要迭代求解，本文为了简便原因，决定采用Hunt公式进行计算。
1.3拓展后的方程
通过对方程的扩展，得到如下形式的方程：
（12）
2. 方程离散和求解
方程的数值离散采用ADI方法，该格式能无条件稳定。将时间步长分两步，分别对x，y方向进行迭代求解。ADI方法虽然无条件稳定，但收敛的速度太慢，为了加快收敛速度，引入松弛因子，每次迭代完毕执行如下：
（13）
式中。
为了更好的处理边界条件，需要用到Kirby的三维近似方法，但它有很大的自身限制性。这样，就在各边界处虚设一边界层，运用二维的方法对边界进行处理。
3. 模型验证
为了验证本文模型的可用性, 采用了Berkhoff经典试验地形进行了数值计算, 并将数值计算结果和实验数据进行对比。
Berkhoff地形是斜坡和椭圆型暗礁相结合的地形，Berkhoff等（1982）对该地形进行了波浪传播实验，并取得了8个断面的实测资料，该试验作为复杂地形上的经典试验之一，常被用于波浪传播数学模型的验证。该物理模型水深地形图见图3所示，计算范围20 x 25，斜坡梯度为1:50，斜坡梯度方向与波浪入射边界方向夹角为20º，入射波周期=1.0，入射波高=0.0464。本文数值试验参数为：计算区域为20 x 25，入射波周期=1.0，入射波高=0.0464，时间步长取0.05，空间网格步长取0.1，上边界为入射边界，下边界为消波边界，左右边界为全反射边界，消波系数为0.3。
图3 Berkhoff物理模型水深地形图
图4断面1、3、5、7波高比较（注：“”为试验值，“—”为缓坡方程计算值）
由以上的相对波高等值线图4可以看出，数值模拟结果和物理试验吻合良好。在平底部分，对波高的影响很小，随着斜坡的变化，以及椭圆地形的存在，发生了波浪的折射和绕射。
4. 模型应用
广东某港，为工程需要而考虑建防波堤，以保护内港区，由于我们只关心港池内的波况, 为节约计算时间, 本文只取如图5中3600m1800m 的矩形区域作为计算域。因东南向波浪为主波向, 现仅计算该波向50年一遇重现期下的极端高、设计高水位。经过推算，口门等深线-6处的波高H=2.8m，周期T=7.8s，现取时间步长=0.5s，计算网格步长v x = v y = 5m。计算所得各波面如图6、图7所示。由图可见,入射波浪由东南向向前传播, 遇防波堤后部分反射,在防波堤趾处,入射波浪绕射进入港池,受防波堤掩护作用, 港池内波高较小。
图5 港区数值计算区域地形以及防波堤位置图
图6SE向50年一遇极端高水位等比波高等值线图 图7SE向50年一遇设计高水位等比波高等值线图
关键词：插值算法；改进的IDW；IDW；Kriging；Spline；ArcGis Engine
中图分类号：P208 文献标识码：A 文章编号：1007-9599 （2013） 02-0000-02
目前在GIS在海洋环境监测中的应用越来越广泛，然而实际监测过程中要得到覆盖全部海域的各种环境要素的监测数据是不容易现实的，通常情况下有关监测部门选取有限的具有代表性的位置进行监测。如果要对整个海域的环境状况进行评估，那么就需要空间插值（又称空间估计）。
1 空间插值算法的概述
1.1 空间插值的依据
空间插值是根据已知空间或分区空间的数据推求任意空间数据的方法，其理论依据“地理学第一定律”：空间位置上越靠近的点具有相似特征值的可能性越大，反之则越小。也就是说在推求未知点的特征值时，距离该点越近的已知点对该点的影响越大，即权重系数越大，反之则越小。
1.2 空间插值算法的类型
空间插值算法需要的数据可以是离散数，也可以是分区数据。其目的在于将这些非连续的数据经过插值后形成连续的数据曲面，便于对各个空间的特征值进行比较，从而得出对实际生活、生产有价值的分析结果。空间插值算法可分为内插法和外推法，所谓内插法是指通过已知空间数据推求同一空间的其他点的特征值的计算方法；所谓外推方法是指通过已知空间数据求其他区域的点的特征值的方法。该系统建立在ArcGis engine平台里，其中主要集成了反距离加权插值算法、克里金插值算法、样条函数插值算法、趋势面插值算法、差变函数插值算法。本文将对上述前2种算法进行研究，并对改进反距离加权插值算法。
2 算法研究
2.1 反距离加权插值算法
2.2 克里金插值算法
3 结论
通过研究IDW插值、Kriging插值、Spline插值、改进的IDW的原理和实现，影响插值效果的因素很多，比如插值方法、已知特征值的分布、采样数目。
IDW算法受已知区域内样本点数目及其权重的共同影响，由于其操作简单，容易实现，目前仍然被广泛使用。由于该方法过于依赖权重系数，所以其受非均匀分布的特征值影响较大，易产生偏差。Kriging算法有坚实的统计学和概率论的基础，可以得到每个估计值的可靠度，计算出每个估计值的误差大小，能够推求出更加准确的特征值。其不足在于该方法计算量大、速度慢、耗费资源量大。Spline算法比较适用于有大量密集点内插和可用于逐渐趋于表面化的情况，但是其误差难以估计。改进的IDW算法相比于普通IDW，前者的插值结果在均方误差、偏差均值、标准差、绝对偏差均值、绝对偏差最大值和绝对偏差最大值等方面的数值均小于采用后者所得的插值结果，也就是说引入变异函数的改进IDW算法的插值效果优于普通IDW的插值效果。通过研究发现Kriging插值算法优于IDW插值算法，而IDW插值算法优于Spline算法。然而对于对空间插值而言， 没有绝对最佳的方法， 只有根据实测样本数据进行充分的分析，选择最佳的方法。
参考文献：
[1]林振山，袁林旺，吴得安.地理建模[M].气象出版社，2003.
关键字：重力坝；双滑面；稳定性；可靠度；抗力角
中图分类号：TV642 文献标识码：A
重力坝的抗滑稳定性一直是重力坝设计的关键问题，潘家铮提出被动抗力法、主动抗力法、和等安全系数法来计算其稳定性。在重力坝实际工程中，影响重力坝稳定的各项因素存在差异性和不确定性，所以采用定值的安全系数法存在一定的局限性。因此，许多学者寻求用可靠度理论来进行重力坝的可靠度评价 。随着双滑面重力坝的可靠度方法研究加深，对于双滑面的可靠度计算，通常使用静力法求解，此时的抗力角取值就具有非常重要的影响力。在安全系数的研究中，部分学者都对抗力角进行了研究 。但是，在可靠度方法中，对于抗力角的取值以及影响的研究还不多，本文通过对比安全系数与考虑可靠度理论的等K′法对抗力角的敏感性，分析并得到抗力角的取值建议。
1 考虑可靠度的安全系数法（等K’法）
通常的可靠度分析往往是通过坝基岩体的破坏准则或者静力平衡来建立坝基抗滑稳定极限状态方程，然后用JC法进行可靠度计算，但是这种方法往往需要大量复杂计算，而且得到可靠度不能和现今的安全系数指标对应起来，因此一种考虑可靠度的安全系数法（等K’法）被常晓林等提出来。
水利水电工程结构可靠度设计统一标准（GB50199-94）规定，采用各参数与分项系数相结合的方法，水工结构承载能力极限状态基本组合的设计式可统一写为：
对于双滑面重力坝，如图1，其等K′法的表达方法如公式（3）。
由公式1可以看出，等K′法需要对中间变量求解，其可由得到，实际计算中通常假设抗力角，然后求出中间变量进行计算，实际上不同抗力角取值对中间抗力有产生影响，从而对等K′法的结果产生影响。
2 抗力角对其影响以及取值分析
目前抗力角选取分为以下几种：①取；②取为岩体内摩擦角；③取与主滑方向平行。第一种的取值，比较保守。后两种取值未经多个实际工程来验证，规范上也未推荐。所以，抗力角的取值没有定论。
在以上的双滑面重力坝等K′法分析中，要求解就需要对抗力角进行假设，然后求出抗力Q，从而求解。本文认为，通过对比等K′法与通用的等安全系数法，先研究两种方法对抗力角的敏感性，再进一步分析其取值变化规律，得到抗力角的建议取值。
3 实例分析
某双软弱夹层重力坝，工程安全等级为一级，简化后的抗滑稳定计算剖面见（图2），均值0. 4，变异系数0. 3；均值0. 07，变异系数0. 4；均值1. 0，变异系数0. 26；均值0. 9，变异系数0. 42；排水孔中心线处扬压力系数0. 25，变异系数0. 30；坝体混凝土容重2.4；基岩容重2.7。
采用常规的等安全系数法计算安全系数K，以及采用等K′法计算K′，他们对抗力方向角的敏感性结果如下（表1），（图3）。
由计算结果可以分析得到：等K法安全系数K以及等K′法的K′，随着抗力角的正大逐渐增大，对抗力角都比较敏感，当时，两者最小，工程上偏于安全，随着抗力角的增大，两者逐渐增加。但是两者间的差距随着抗力角增大减少，表示两者差异性对抗力角敏感性降低，两种方法差异性降低。另外，两种方法对比后显示，在150m的坝高下，K′数值小于K，表示前者工程上用于设计偏于安全。对于本例题，对于抗力角的取值，取0°显然偏于安全，在实际工程中往往过于安全，通过曲线可以看出，在抗力角取值25°之后，曲线斜率增大，表示抗力角敏感度加强，当抗力角取值在10°～25°之间时，取值可以是较为安全的，可以认为是较为合理的。
4 结论
本文基于可靠度理论对双滑面重力坝的抗滑稳定性分析，指出抗力角取值对其不确定性，并结合等K′法与通用的等安全系数法对某重力坝算例分析，指出不同方法对抗力角的敏感度，以及得到抗力角的合理取值，得到以下结论：
（1）以等K′法计算坝基稳定性，能够得到唯一的值，该值结合了坝体参数的不均一性，与可靠度理论结合，是可行的。
（2）在抗力角取值时，0°是过于安全，随着抗力角取值增大，得到的K′以及K值也随着增大，但它们对抗力角的敏感性减低。
（3）在实际工程应用中，抗力角的取值不尽相同，要根据抗力角的敏感性分析，对不同工程条件，不同计算方法选择不同的抗力角取值。
参考文献
[1]潘家铮.重力坝设计[M].北京：水利电力出版社，1987.PAN Jiazheng. Design of dames [M]. Beijing： Water Conservancy and Electric Power Press， 1987. （china）.
[2]中华人民共和国国家标准编写组.GB 50199-94水利水电工程结构可靠度设计统一标准[S].北京：中国计划出版社，1994.
关键词：汉土什公式；正逆问题；条件数；公式性态；误差传递；近似关系；变化曲线
中图分类号：P641；TV 211.12 文献标志码：A 文章编号：16721683（2015）05095904
Performance analysis of unsteady flow equation with first kind of leakage recharge
FANG Bei1，LIU Yuanhui1，GUO Jianqing2
（1.College of Science，Chang′an University，Xi′an 710064，China；
2.School of Environmental Science & Engineering，Chang′an University，Xi′an 710051，China）
Abstract：When the original data has random errors，the Hantush equation for the first kind of leakage aquifer system can transfer the errors when it is used in the calculations of the direct and inverse problems.The conditional numbers to characterize the performance of Hantush equation have been derived，and the approximate relationship between the errors from the calculated results and errors from the original data is determined for the direct and inverse problems.On the basis，the transfer effects of different u and r/B values on the random errors from the original data are obtained，the variation curves of conditional numbers for characterization of the transfer effects are plotted，and the scopes of u and r/B values are determined when the direct and inverse problems have “bad performance”.
Key words：Hantush equation；direct and inverse problems；conditional numbers；equation performance；error transfer；approximate relationship；variation curve
在第一类越流系统中，通常采用汉土什公式分析抽水试验数据，确定越流含水层的水文地质参数[1]，其形式简单，故在实际计算中得到了广泛应用。把在抽水过程中，采用已知含水层水文地质参数计算含水层的水位降深值的问题称为正计算问题；而采用抽水过程中观测到的水位降深值，反求含水层水文地质参数的问题称为逆计算问题。在正逆计算问题中总会有不确定因素出现，导致代入公式的参数存在一定的误差，从而由汉土什公式计算出的结果也具有不确定性[2]。掌握汉土什公式的性态及公式本身对误差的传递作用[3]对解析解的正逆计算具有重要意义。Jiu等曾采用灵敏度分析的方法对含水层参数的误差进行了研究[4]，McElwee等曾分析了地下水模型的灵敏性[56]，这些分析都采用定量化的方法讨论了计算问题时公式的性态，但却得出了定性的结论。郭建青等用误差分析的方法讨论了泰斯公式的性态[78]及一维河流水质方程解析解的性态[910]，把“条件数”作为衡量公式性态的标准，用以判断公式本身对原始数据的随机误差的敏感程度[11]。这种分析方法的原理简单，易于理解和计算。本文将利用文献[7]的思路计算汉土什公式的条件数，并对汉土什公式的性态进行初步分析，最后达到定量的描述该公式对数据误差的传递作用的目的，进而依据结论可以尽量避免“病态”条件下的计算，以免计算结果“失真”。
1 基本原理
1.1 基本公式
函数随机误差的计算[12]可表述如下。
设f（x）是关于n个随机变量x1，x2…xn的函数，在每个变量xi的均值附近作二阶泰勒展开，则有
f（x）=f（x）+∑ni（x）（xi）（xi-xi）+12∑ni∑nj2f（x）xixj（xi-xi）（xj-xj）
（1）
由于xi是随机变量，故x-xi也是随机变量[13]，从而f（x）是随机函数。对式（1）两端分别取数学期望，则有：
f（x）=f（x）+12∑ni∑nj2f（x） xixjσxixj
（2）
式中：σxixj为随机变量xi和xj的协方差。如果xi和xj相互独立时，就有σxixj=0，从而可得
f（x）=f（x）
（3）
用式（1）减去式（2）得
f（x）-f（x）=∑nif （x） xi（xi-xi）+12∑ni∑nj2f（x） xixj（xi -xi）（xj-xj）-12∑ni∑nj2f（x） xixjσxixj
（4）
式（4）两端分别自乘并取数学期望后得：
σ2f≈∑ni∑njf（x）xif（x）xjσxixj
（5）
式（3）和式（5）是以下分析的基础。
1.2 汉土什公式的均值形式
在第一类越流系统中，通常用于求解水文地质参数的汉土什公式[14]如下：
s=Q4πtW（u，B）
（6）
其中，
W（u，B）=∫∞uexp（-y-r24yB2）ydy
（7）
u=r2μ4Tt
（8）
式中：s为水位降深值（m）；Q为抽水流量（ m3/d）；T为含水层导水系数（m2/d）；μ为承压含水层的弹性释水系数，无量纲；1/B为含水层的越流补给因子，无量纲；r为径向距离（m）；t为时间（d）。
利用汉土什公式进行正逆问题计算时，代入公式的原始数据总会存在着随机误差。在正计算问题中，把T、μ、1/B看作是随机变量，则由式（6）计算得到的水位降深值s也是随机变量；在逆计算问题中，把水位降深s的观测值看作是随机变量，则由式（6）反求出的T、μ、1/B也是随机变量。
2 分析与讨论
在式（3）中，令
f（x）=f（T，μ，1/B），x1=T，x2=μ，x3=1/B则有
s=f（T，μ，1/B）
（9）
表明当随机变量T、μ、1/B相互独立时，汉土式公式的均值形式与确定形式是一致的，则在此不需进行过多的讨论。
在式（5）中，令
f（x）=f（T，μ，1/B），x1=T，x2=μ，x3=1/B
则有
σ2s=φ 2 Tσ2T+φ 2μσ2μ+φ21/Bσ21/B
（10）
其中，
φ T=sT
（11）
φ μ=sμ
（12）
φ 1/B=s1/B
（13）
式中：σ2s，σ2T，σ2μ，σ21/B分别为s、T、μ、1/B的方差。
式（10）反映了随机函数的方差与随机变量方差间的关系。为了分析与讨论方便，引入如下无量纲变量：Cs=σss，CT=σTT，Cμ=σμμ，C1/B=σ1/B1/B
则式（10）可以写成如下形式：
C2s=ψ 2TC2T+ψ 2μC2μ+ψ21/BC21/B
（14）
再由各式之间的关系可以导出：
ψ TTsφ T=-exp（-u-r24uB2）W（u，B）-1
（15）
ψ μμsφ μ=-exp（-u-r24uB2）W（u，B）
（16）
ψ1/B1/Bsφ1/B=-r22B2∫∞u1ydy
（17）
2.1 正问题的分析与讨论
式（14）可以理解为汉土什公式正问题计算时，第一类越流含水层系统中的参数T、μ、1/B的相对误差的均方差分别为CT、Cμ、C1/B，由公式（6）计算出的水位降深值s的相对误差的均方差为Cs。在数值分析原理[15]中称φ T、φ μ、φ1/B的绝对值为条件数，其值大于10时，公式本身对代入数据的误差起“放大”作用，如果条件数很大时，表明即使代入数据的误差很小，也会导致计算结果的误差较大甚至“失真”，此时称该计算问题为“病态”问题，其他情况下都可以看作是“良态”的；当条件数小于10时，公式本身对代入数据的误差起“缩小”作用。
图1（a）为利用式（15）绘制的相应于不同的u和r/B值的条件数|ψ T|变化曲线。可以看出，当r/B一定时，|ψ T|随着u值的增大先减小后增大，当u足够小时|ψ T|趋近于1.0，即公式本身对T的误差的“放大”作用较弱，此时的正计算问题可以看作是“良态”的；当u
图1（b）为利用式（16）绘制的相应于不同的u和r/B值的条件数变|ψ μ|化曲线。可以看出，当r/B一定，u足够小时，|ψ μ|趋近于0，即此时在正计算问题中，公式（6）本身对μ所具有的随机误差几乎没有影响；|ψ μ|的极小值随着u的减小而沿着r/B增大的方向移动，即当u值较小且r/B值较大时，公式本身对μ所具有的随机误差的影响最小；随着u的增大而增大，当u>1.0时，|ψ μ|的值多大于1.0，在这种情况下，利用公式（6）计算出的结果具有很大的误差，此时正计算问题属于“病态”问题。
图1（c）为利用式（17）绘制的相应于不同的u和r/B，条件数变化|ψ1/B|曲线。可以看出，当r/B一定时，|ψ1/B|随着u的增大而减小；当u一定时，|ψ1/B|随着r/B的减小而减小，即|ψ1/B|的极小值随着u的增大而沿着r/B减小的方向移动。当r/B
图1 正问题条件数变化曲线
Fig.1 Variation curves of conditional numbers of direct problems
2.2 逆问题的分析与讨论
公式（6）的逆计算问题的条件数，可以通过对式（14）作适当的变形得到：
σ2T=1σ 2T（σ2s-φ 2μσ2μ-φ 21/Bσ 21/B）
（18）
σ2μ=1σ2μ（σ 2s-φ 2Tσ2T-φ 21/Bσ21/B）
（19）
σ21/B=1σ21/B（σ2s-φ 2Tσ2T-φ 2μσ2μ）
（20）
由上式可知，在利用汉土什公式确定越流含水层水文地质参数时，当其中两个参数已知时，其所具有的误差对水位降深的观测误差起到了消减作用，故为保守和讨论方便起见，仅把水位降深的观测误差看作是影响计算结果误差的独立随机变量。以下分别对σμ=0，σ1/B=0；σT=0，σ1/B=0和σT=0，σμ=0三种情况进行讨论。
（1）σμ=0，σ1/B=0的情况。
在这种情况下，式（14）可以写为
σT=1φ Tσs
（21）
其无量纲形式为
CT=1ψ TCs
（22）
式（22）可以理解为逆计算问题中，当水位降深观测值的相对误差的均方差为Cs时，由公式（6）反求出的含水层参数T的相对误差的均方差为CT，则相应的条件数为|1/ψ T|。
图2（a）为利用式（15）的倒数形式绘制的相应于不同的u和r/B的条件数|1/ψ T|变化曲线。可以看出，|1/ψ T|随u的减小而趋近于1.0，此时公式（6）对原始数据s的随机误差的“放大”作用较弱；|1/ψ T|随r/B的增大而减小，当r/B的值较大时，|1/ψ T|值较小，此时原始数据s的随机误差对计算结果误差的影响最小；当0.5时|1/ψT|随着u的增大而减小，故在0.5附近，|1/ψ T|取得极大值。在u=0.5附近|1/ψ T|的值相当大，即利用公式（6）反求T时，即使原始数据s具有很小的随机误差，也会导致计算结果具有很大的误差甚至“失真”，此时的计算问题属于“病态”问题。
（2）σT=0，σ1/B=0的情况。
在这种情况下，式（14）可以写为
σμ=1φ μσs
（23）
其无量纲形式为
Cμ=1ψ μCs
24）
上式为逆计算问题中，水位降深观测值的相对误差的均方差Cs与利用式（6）反求出的含水层参数μ的相对误差的均方差Cμ间的近似关系，则相应的条件数为|1/ψ μ|
图2（b）为利用式（16）的倒数形式绘制的相应于不同的u和r/B的条件数|1/ψ μ|变化曲线。可以看出，当u值较大时，|1/ψ μ|值较小且大都小于10，此时利用式（6）反求μ时，公式本身对原始数据s的误差起到了“缩小”作用；当u值较小或r/B较大时，|1/ψ μ|值多大于10，即表明在原始数据s具有很小的随机误差时，利用式（6）反求的结果误差较大甚至“失真”，此时的计算问题属于“病态”问题。
（3）σT=0，σμ=0的情况。
在这种情况下，式（14）可以写为
σ1/B=1φ1/Bσs
（25）
其无量纲形式为
C1/B=1φ1/BCs
（26）
上式表示逆计算问题中，水位降深观测值的相对误差的均方差Cs与利用公式（6）反求出的含水层参数1/B的相对误差的均方差C1/B间的近似关系，则其相应的条件数为|1/ψ1/B|。
图2（c）为利用式（17）的倒数形式绘制的相应于不同的u和r/B的条件数|1/ψ1/B|变化曲线。可以看出，|1/ψ1/B|随u的增大而增大，随r/B的增大而减小。当r/B>05时，|1/ψ1/B|值始终小于10，此时式（6）对原始数据s的误差起到“缩小”作用；|1/ψ1/B|的极小值随着r/B的增大而沿着u减小的方向移动，此时公式本身对原始数据s的误差传递几乎没有影响；当r/B较小时，|1/ψ1/B|值趋近于无穷，此时即使原始数据s具有较小的随机误差，利用式（6）反求出的1/B的误差也可能较大甚至“失真”，即此时的计算问题属于“病态”问题。
图2 逆问题条件数变化曲线
Fig.2 Variation curves of conditional numbers of inverse problems
3 结论
通过以上分析和讨论，可以得出如下结论：利用汉土什公式计算时，若原始数据具有随机误差，则利用公式（6）计算得到的结果也会具有一定的误差，即汉土什公式对代入数据误差具有传递作用；在正计算问题中，当u值较大时，T、μ的误差对正计算问题的影响是“病态”的，而当r/B较大时，1/B的误差对正计算问题的影响是“病态”的；在逆计算问题中，当u=0.5及附近时，计算μ和T的问题属于“病态”的，而当r/B较小时，计算1/B的问题属于“病态”的。其它条件下的计算问题均属于“良态”的。在利用汉土什公式进行正逆问题计算时，应尽量避免“病态”问题下的计算，以免计算结果产生较大的误差甚至“失真”，以便保证计算结果的精度。
参考文献（References）：
[1] 朱学愚，钱孝星.地下水水文学[M].北京：中国环境科学出版社，2005.（ZHU Xueyu，QIAN Xiaoxing.Groundwater hydrology[M].Beijing：Chinese Environment Science Press，2005.（in Chinese））
[2] 项彦勇.地下水力学概论[M].北京：科学出版社，2011.（XIANG Yanyong.Fundamental theory of groundwater mechanics[M].Beijing：Science Press，2011.（in Chinese））
[3] 张诚坚，高健，何南忠.计算方法[M].北京：高等教育出版社，1999.（ZHANG Chengjian，GAO Jian，HE pututational methods[M].Beijing：Higher Education Press，1999.（in Chinese））
[4] Jiao J J.Sensitivity features of aquifer parameters and their implications on parameter estimation[D].Ph.D.Thesis，Dept.of Civil Engingineering，Univ.of Birmingham.1993.
[5] McElwee，C D.and M A Yukler.Sensitivity of groundwater models with respect to variations in transmissivity and storage[D].Water Resour.1978，14 （2）：451459.
[6] McElwee，C D.Sensitivity of groundwater models[D].In：Advances in Transport Phenomena in Porous Media.Edby J.Bear and M.Y.Corapcioglu.NATO Adv.Study Inst.Ser.E.v.1987，128：751817.
[7] 郭建青，李云峰，王洪胜.泰斯公式性态分析与误差估计方法[J].煤田地质与勘探，1999，27（4）：3842.（GUO Jianqing，LI Yunfeng， WANG Hongsheng.The performance analysis of theis equation and error estimation method[J].Coal Geology and Explorati on，1999，27（4）：3842.（in Chinese））
[8] 郭建青.泰斯公式对数据误差传递作用的初步分析[J].工程勘察，1988（1）：4448.（GUO Jianqing.The analysis transfer of the heis equation in datas' error[J].Engineering Investigation，1988（1）：4448.（in Chinese））
[9] 郭建青，李云峰，钱会.一维河流水质方程解析解性态的初步分析[J].水文，1999，8（6）：3033.（GUO Jianqing，LI Yunfeng，QIAN Hui.Analysis on the peerformance of onedimension stream water quality equation's analytical solution[J].Hydrology，1999，8（6）：3033.（in Chinese））
[10] Knopman，D S and C I Voss.Behaviour of sensitivities in the onedimensional advectiondispersion equation：implications for parameter estimation and sampling design[D].Water Resour.1973，23（2）.
[11] 薛禹群，谢春红.地下水数值模拟[M].北京：科学出版社，2007.（XUE Yuqun，XIE Chunhong.Numerical simulation for groundwater[M].Beijing：Science Press，2007.（in Chinese））
[12] 沙定国.误差分析与测量不确定度的评定[M].北京：中国计量出版社，2003.（SHA Dingguo.Evaluation of uncert ainty in measurement and analysis of error[M].Beijing：China Metrology Publishing House，2003.（in Chinese））
[13] 田铮，秦超英.随机过程与应用[M].北京：科学出版社，2009.（TIAN Zheng，QIN Chaoying.Stochastic process and application[M].Beijing：Science Press，2009.（in Chinese））
论文关键词：住房公积金，见角进元，excel，缴交基数，函数
一、住房公积金
住房公积金是单位及其在职职工缴存的长期住房储蓄金，实际上是一种住房保障制度函数，是住房分配货币化的一种形式论文范文。其具有强制性、互助性、保障性，单位为职工缴存的住房公积金是职工工资的组成部分，单位为职工缴存住房公积金是单位的义务，享受住房公积金政策是职工的合法权利。职工个人及单位（国家机关、国有企业、城镇集体企业、外商投资、城镇私营企业及其城镇企业、事业单位、民办非企业单位、社会团体）缴纳的住房公积金函数，实行专户存储，归职工个人所有。
现在缴纳住房公积金已经成为了势在必行，但住房公积金的缴交基数却也是需要经常调整的，由于住房公积金缴交基数的计算有一个特殊的要求函数，就是计算时要求“见角进元”，用手工操作的话，如果单位人数过多，工作量将相当大函数，很有谱可能出现计算差错论文范文。于是经过摸索，发现在excel中可以运用函数简单的解决此问题。
二、准备工作
1、首先，上一年度的工资表复制并转帖到同一个excel文档，并命名为住房公积金调整表函数，主要字段有职工号、姓名、身份证号码、调整前工资、调整后工资、调整后缴交基础（个人8%、单位20%、合计和备注）等字段。制作好“住房公积金缴交基数调整表”如下图1：
2、计算时要用到的两个公式Roundup 和 Rounddown
⑴Roundup函数是将数字按照指定的位数进行无条件进位，即远离零值，向上舍入数字论文范文。语法：ROUNDUP(number,num_digits)，其中Number 为需要向上舍入的任意实数函数，Num_digits 四舍五入后的数字的位数。num_digits 大于 0，则向上舍入到指定的小数位。 如果 num_digits 等于 0，则向上舍入到最接近的整数论文范文。 如果num_digits 小于 0，则在小数。例如：=ROUNDUP(86.9,0) 将 86.9 向上舍入函数，小数位为 0 ，值为87 ；=ROUNDUP(4.16159,3) 将 4.16159 向上舍入，保留三位小数 ，值为4.162。
⑵Rounddown函数是将数字按照指定的位数进行无条件舍去函数，即靠近零值，向下（绝对值减小的方向）舍入数字论文范文。语法：ROUNDDOWN(number,num_digits)，其中Number 为需要向下舍入的任意实数。Num_digits 四舍五入后的数字的位数。 如果 num_digits大于 0，则向下舍入到指定的小数位论文范文。 如果num_digits 等于 0函数，则向下舍入到最接近的整数。如果num_digits 小于 0，则在小数点左侧向下进行舍。例如：=ROUNDDOWN(86.9,0) 将 86.9 向下舍入，小数位为 0 ，值为86 ；=ROUNDDOWN(4.16159,3) 将 4.16159 向下舍入函数，保留三位小数 ，值为4.161论文范文。
三、嵌入函数，开始计算
基本的计算思路是用rounddown函数无条件舍去“调整后工资基数”与“比率”的乘积得分为数值，使其乘数保留至角位函数，然后再用roundup函数把此值无条件进位至元。
具体操作第一步，是在已经建好的调整表中，在F7中输入“=Roundup(Roundup(E7*8%,1),0)”，按回车键得到邓仁斌的个人调整后缴交数为440函数， 然后再选中F7，把鼠标放在F7的右下角，待鼠标图标变成十字架就往下拖，也即往下填充就可以函数，至此，“见角进元”的问题就解决了。如下图 2：
第二步，同理在G7中输入“=ROUNDUP( ROUNDDOWN(E7*20%,1),0)”,在按回车键，得到公司帮邓仁斌缴纳的住房公积金是1100函数，然后再由G7往下填充，计算如下图3：
第三步，在H7中输入“=SUM(F7+G7)”，按回车键即可得到邓仁斌有自己和公司总共缴交住房公积金为1540函数，然后再由H7向下填充即行，到此整个计算过程和步骤也结束了，最终计算结果如下图4：
四、注意的问题
1、在制作“住房公积金缴交基数调整表”时，会用到Vlookup函数提取职工月工资，利用Isna和If函数计算月平均工资，本文笔者并未对此作详细讲解，但在用这些函数时应要注意函数语法和用法问题。
2、在计算住房公积金缴交基数时，有一个要求“见角进元”，而解决这个问题也要用到两个函数Rounddown和Roundup两个函数，同样也要注意到他们的语法和用法。
3、如果只用一个函数不能解决一个问题，在函数嵌套中，我们可以把几个函数的功能结合起来解决用单个函数无法解决的问题，是的实际工作中遇到的许多问题都能迎刃而解，大大提高工作中的效率。就如本文中应注意的，如果只用rounddown或roundup函数不能解决“见角进元”这个问题，那么只能把二者结合起来才可以得到正确的结果。
