库容（矿体、工程土方）的几何形状与容（体）积计算
2018-09-10陶祖昶
摘要：水库是形状不规则的容体，矿体和自然地形上的工程土方是形状不规则的岩土体，彼此形态类似，容（体）积的计算方法相同，用等高线法计算时，都是按照矿体几何学的方法，假定其为无数直线围成的截头锥状体，采用梯形公式或截锥公式计算容（体）积。这种方法对计算体的几何形状认定不当，对计算公式的模型图形认识有误。应用数值积分的方法计算库容，可以找到识别库容几何形状的方法，导出各种几何形状的容积计算式，查明现行各种计算式的模型图形，纠正传统计算方法中的错误。研究表明：库容的几何形状有凸形、直形或凹形三类形状，各水库之间或一座水库的不同高程之间都不尽相同;几何形状不同，容积计算公式也就各异，没有通用的计算式;若用错公式，则当相邻两计算剖面面积之差为40%时，计算的容积将有0.468%～11.111%的偏差。
关键词：水库容积;矿体;工程土方;几何形状;计算模型;数值积分;容（体）积计算;
中图分类号：TV221.1 文献标志码：A
水库容积、矿体和工程土方体积是水资源开发利用、矿产开采、工程建设必需的基本资料。水库是形状不规则的巨大容体，矿体和基于自然地形的工程土方是形状不规则的岩土体，彼此形状类似，体积的计算方法也相同。根据地形测绘资料或地质勘探资料计算容（体）积的方法大致有3类：横向计算的有垂直剖面法（断面法、平行断面法、不平行断面法）;竖向计算的有水平剖面法（等高线法、等值线法、地形法）;把容体切割成众多棱柱体计算的有方格法、三角网法以及数字高程模型法（DEM）。对于库容或轴状矿体、山丘土体等块状岩土体，多用水平剖面法;对板状矿体或道路、渠道、河道工程土方等条状岩土体，则用垂直剖面法。这两种方法是最常用的方法，基本原理相同，均假定两个剖面之间的体积是截锥体或角壔体（Prizmatoid），用梯形公式、截锥公式或角壔体公式计算容（体）积。长期以来，由于对这些不规则物体的几何形状及体积计算公式的计算模型缺乏正确的认识，因此在计算中陷人某些误区，笔者在前人探讨矿体计算的基础上，根据对地形法计算库容的分析研究，找到了识别库容几何形状的方法，导出了各种几何形状的库容计算公式，这些方法和原理同样适用于矿体和工程土方体积的计算。
1 传统计算方法及存在的问题
形状不规则的矿体，是无法用规则体积的计算公式计算其体积的。用水平剖面法计算时必须在等高线形状的基础上将其转换成符合其总体几何形状的、可计算的概化模型。但是，概化模型具有半虚拟的性质，其四周界面是无法直接测量到的，于是鲍曼（BaymaH）就假定（1909年）：两个水平剖面（等高线）之间的矿体是四周界面（围岩）由无数斜率不同的直线所围成的截头锥状体，从而以立体几何的方法推导出了这种锥状体的体积计算式，即鲍曼公式。索布列夫斯基（Cобулевский）则假定（1934年）：两平行剖面之间的矿体是由若干不同倾斜度的平面所围成的角壔体，并用同样方法推出了角壔体的体积计算式，亦即索布列夫斯基公式[1-2]。这两个公式本质相同，实为一式。鲍曼公式计算较为复杂，故被逐步舍弃，形式简便的梯形公式和截锥公式被引入矿体、工程土方以及库容的计算中，并规定或建议：相邻两剖面等高线内面积之差小于40%时用梯形公式，大于40%时用截锥公式;或一般库容计算用梯形公式，“精确”或“严密”计算时用截锥公式[3-4]。
文献[5-6]用数值积分的方法推导了鲍曼公式和截锥公式，并与梯形公式相比较，认为此三式都是锥状体的计算式，但梯形公式计算结果偏大，鲍曼公式则更具普遍性，比较精确。同时，还推出了马鞍形矿体的体积计算式。文献[7]同意以上论述及规定，但认为矿体周围的竖向界面大都呈凸形弧面状，采用抛物线公式（Simpson公式）计算体积比较切合实际，若采用立方抛物线公式（辛普森3/8法则）计算则更加精确。文献[8]认为库容剖面的几何形状应有凸形、直形及凹形等多种（相应的体积为碗状、截锥状及马鞍状），凸形库容适用梯形公式计算。这些文献都对计算体形状的假定及传统的计算公式提出了质疑。
近年来，随着测绘技术及计算技术的发展，在库容和土方量计算中使用了许多计算软件，替代了手工计算，其中有些软件是按照地形法或断面法原理编制的，但是其建模中所用的体积计算公式仍如旧规，所以迄今为止，用水平剖面法或垂直剖面法计算库容、矿体或工程土方时，仍然遵循百年来传统的观点与方法。而由于这些“假定”“规定”“建议”没有弄清计算体的几何形状，其假设的形状未必符合实际，同时对计算公式的模型又不够了解，混淆了平面求积与空间求积不同的概念，因此造成了体积计算中的许多错误，降低了计算的精确度。用数值积分的方法计算水库容积并分析研究，不仅可以了解各种计算公式的计算模型，而且可以找到识别库容几何形状的方法，从而纠正传统计算方法的错误，建立比较正确、合理的计算方法与步骤。
2 库容几何形状与数值积分计算
2.1 根据水库高程—面积曲线计算库容
水庫各高程（z）的面积值（A）是库区地形测绘必须提供的基本资料，根据这些资料就可以直接计算库容。图1表明，水库的高程—面积曲线不是单一、连续的函数关系，而是不连续的函数关系。说明水库库容整体而言是一个不规则的体积，但在函数连续的区间却是有一定规则的，可以用数值积分的方法计算库容，即将整座水库按高程—面积曲线函数的连续区间分成几段，分别计算其库容，然后累加得到总库容，从而解决了不规则体积的积分计算问题。
在水库高程—面积曲线函数的连续区间，任一微薄层的容积为dV=Adz，高程a～b的容积为
V=∫abAdz （1）
式（1）中若以一个代数多项式A=a0+a1z+a2z2+…+anzn代替A=f（z），水库地形图等高线内面积依次为A0、A1、A2、…、An[图2（a）]，高程间距为h，且A（0）=A0，A（1h）=A1，…，A（nh）=An，则可求得高程—面积关系方程式不同方次n的体积计算式，如当：1h，则有则有
按牛顿一柯特斯（Newton-Cotes）公式[9]还可得到以下方次的计算式：
由此可以根据水库高程—面积曲线方程的方次选择相应的库容计算式。
2.2 各种几何形状的体积计算
在矿体开采或工程土方的挖掘中，不仅要知道其数量，还要知道其形状，所以习惯按形状计算其体积。在库容计算中了解库容的几何形状也有助于水库特性的研究及容积的正确计算。
图2（a）为水库地形示意图，O为最低点在图上的投影。若通过O点以任意θ角作一截面，则可得到此水库的剖面。剖面可能有A、B、C三种类型，即凸形（碗状）、直形（截锥状）及凹形（马鞍状）[8]。若通过O点建立柱面坐标，则其体积为
式（7）中剖面线函数ρ（z，θ）=f（z），也可用一个代数多项式ρ（z，θ）=a0+a1z+z2z2+…+anzn或ρ2（z，θ）=a0+a1z+a2z2+…+anzn代替，将不同方次n的方程式代入式（7），即可导得3类剖面形状的体积计算式，分述如下。
2.2.1 A型（凸形碗状，n<1）
若n=1/4，即ρ2（z，θ）=a0+a1z1/2，则1h间的体积为
若n=1/2，即则1h间的体积为[8]
2.2.2 B型（直形截锥状，n=1）（3个剖面）间的体积为
亦即索布列夫斯基公式：
若取两根等高线计算1h间的体积，则就是鲍曼公式[5-6]：
若式（11）中ρ2（θ）=kρ1（θ），k为常数，则式（11）可改写成常用的截锥公式（棱台公式）：
2.2.3 C型（凹形马鞍状，n>1）
若n=3/2，的体积为[7，9]则有
按牛顿一柯特斯（Newton-Cotes）公式[9]还可得到以下方次的计算式：
由此可以根据水库高程—面积曲线方程的方次选择相应的库容计算式。
2.2 各种几何形状的体积计算
在矿体开采或工程土方的挖掘中，不仅要知道其数量，还要知道其形状，所以习惯按形状计算其体积。在库容计算中了解库容的几何形状也有助于水库特性的研究及容积的正确计算。
图2（a）為水库地形示意图，O为最低点在图上的投影。若通过O点以任意θ角作一截面，则可得到此水库的剖面。剖面可能有A、B、C三种类型，即凸形（碗状）、直形（截锥状）及凹形（马鞍状）[8]。若通过O点建立柱面坐标，则其体积为
式（7）中剖面线函数ρ（z，θ）=f（z），也可用一个代数多项式ρ（z，θ）=a0+a1z+z2z2+…+anzn或ρ2（z，θ）=a0+a1z+a2z2+…+anzn代替，将不同方次n的方程式代入式（7），即可导得3类剖面形状的体积计算式，分述如下。
2.2.1 A型（凸形碗状，n<1）
若n=1/4，即ρ2（z，θ）=a0+a1z1/2，则1h间的体积为
若n=1/2，即则1h间的体积为[8]
2.2.2 B型（直形截锥状，n=1）（3个剖面）间的体积为
亦即索布列夫斯基公式：
若取两根等高线计算1h间的体积，则就是鲍曼公式[5-6]：
若式（11）中ρ2（θ）=kρ1（θ），k为常数，则式（11）可改写成常用的截锥公式（棱台公式）：
2.2.3 C型（凹形马鞍状，n>1）
若n=3/2，的体积为[7，9]
若则2h间的体积为[6]
若剖面线为幂函数ρ（z，θ）=a0+a1z2，则可得到1h间容积的计算式为[8]
若n=2，的体积为
上述公式中式（3）、式（4）、式（5）与式（9）、式（10）、式（14）传统上称之为梯形公式、角壔体公式（抛物线公式、辛普森公式、索布列夫斯基公式）及立方抛物线公式（辛普森3/8法则）。式（11）与式（10）实为一式。式（13）又称为棱台体公式。式（11）中的Ah/2为h/2处的面积。式（12）、式（15）、式（16）中的T（1，2）、T（2，3）、T（3，1）为面积改正数[6]。式（2）～式（6）及式（8）～式（17）两类计算式同列于表1，以供对照比较。这两种计算库容的方法同样适用于垂直剖面法计算库容或河道、渠道及道路等工程土方体积[8]，其原理相同。
2.3 库容几何形状的识别
库容的几何形状无法由工程测量直接测定，但可以通过间接的方法判断识别。式（1）为根据水库高程—面积曲线积分计算库容，可以求得高程—面积曲线方程不同方次的库容计算式;式（7）则以库容剖面线形状作为变量，可以求得剖面线方程不同方次的库容计算式;很容易证明两者实为一式，仅是座标形式不同而已，表一中由式（1）导得的式（2）、式（3）、式（4）、式（5）、式（6）分别与由式（7）导得的式（8）、式（9）、式（10）、式（14）、式（17）相同，这样就找到了高程一面积曲线的形状与库容剖面线图形之间的关系，只要求出高程—面积曲线A=f（z）的方程式，就可以知道库容模型的剖面线方程P=f（z）及其几何图形。例如，若高程—面积曲线呈直线，方程式方次n=1，则库容剖面线方程为n=1/2的凸形曲线，几何形状为碗状（馒头状）的凸形体，计算公式适用式（3）、式（9）;若高程—面积曲线为2次方程，则剖面线必为斜率不等的直线，体形属直形的角壔体或锥状体，计算公式用式（4）、式（10）。
库容剖面线的几何图形是半虚拟的概化模型，无法直接测定，但可由高程，面积曲线反映出来。而高程一面积曲线是可以根据测绘资料绘制，而且不难求得其曲线方程，加以上述各种形状库容计算式的导出，于是就可以走出盲目假定几何形状、乱用体稠十算公式的误区。
3 传统方法计算容（体）积的几个误区
3.1 误区一
“两平行剖面间矿体（库容）的形状都是截锥体或角壔体，采用梯形公式或截锥公式计算体积”“相邻两剖面面积之差小于40%时用梯形公式，大于40%时用截锥公式”“一般计算用梯形公式，精密计算用截锥公式”“鲍曼公式更具普遍性，比较精确”[1-6]。
图1表明，水库高程—面积曲线在不同区段具有上凸、直线或下弯等多种形状，根据高程—面积曲线形状与剖面线形状的关系，说明库容的剖面形状也相应有凸形、直形或凹形，即体积为碗状、截锥状或马鞍状等形状，并非只是直形的锥状体（角壔体）一种。水库的大断面图也常表明这种情况。此外，从地质学、地貌学、土力学及河流动力学也可得到合理的解释。矿体、山丘等岩土体同样如此。不同形状的区间要用不同的体积计算式，没有普遍适用的通用公式。
以上演算表明，截锥公式是鲍曼公式的特例，是计算正截头锥体（四周斜率相同）的公式，不是什磨特殊的精密计算公式。与“梯形公式”计算图形又不相同，无法比较。体积计算用哪个公式合适、正确，不在于相邻剖面面积差的百分数是多少，而在于计算公式的数学模型与计算体的几何形状是否一致。因此，不问计算体实际的几何形状，假定其都是截锥体、角壔体，一律采用梯形公式或截锥公式计算体积，显然是错误的。
鲍曼公式与截锥公式是计算2个剖面间的体积，而索布列夫斯基公式（角壔体公式）是计算3个剖面间的体积，都是直形锥状体的计算式。除截锥公式属特例（其计算体积略偏大）外，鲍曼公式与其他两式实为一式，只要与计算体模型相符，都一样精确。
3.2 误区二
“梯形公式是剖面为梯形的截头锥体体积的计算公式”[5-6]。
梯形公式本是计算梯形平面面积的求积公式，将其改写成体积计算公式，认为其计算的图形是剖面为梯形的锥状体，没有根据，也没见其来源的报导。而式（9）的演算则证明，体积计算中的“梯形公式”，其计算模型是具有凸形剖面的碗状体，剖面线方程为ρ2（z，θ）=a0+a1z，与“梯形”毫不相干，将此式与截锥公式及鲍曼公式归人一类，是一个误导，混淆了体积计算与面积计算的不同概念。文献[8]早就指出这一情况，但没有引起普遍注意。此式是国内外矿业、水利、土木工程界计算矿体、库容、工程土方等体积中广泛使用的计算式。长期以来，也有人发觉用此式计算的结果大于截锥公式或其他锥状体计算式，也大于其他计算方法[6，10-11]，但由于对其形状的误解，因此始终得不到正确的解释。“梯形”是平面图形，没有体积的。式（9）与“梯形”又毫不相干，将其称为“梯形公式”，名不符實。所以，为避免误解误用，在体积计算中改称“碗状体公式”较妥。
3.3 误区三
对天然界面呈凸形的弧面状矿体“采用辛普森公式计算体积比较切合实际;如果采用立方抛物线公式计算，则更加精确”[7]。
辛普森公式原是计算二次抛物线下曲边梯形平面面积的公式，所以又称抛物线公式。体积计算中形似辛普森公式的式（10）演算证明，它是角壔体的体积计算式，计算模型是3个水平剖面（2h）间的锥状体，垂直剖面线为直线。把它与辛普森公式混淆，以为是计算表面为凸形，纵剖面线为抛物线的体积计算式，这是又一个误解。实际上是其水平剖面的面积值变化呈二次抛物线，而不是其垂直剖面线呈二次抛物线。
辛普森3/8法则是三次抛物线下的平面面积插值计算式，在插值计算中，其余项误差比辛普森公式小，所以插值精度比辛普森公式高。但在体积计算中，形似辛普森3/8法则的式（14）（有称之为“立方抛物线公式”）实际上是剖面线方程为ρ2（z，θ）=a0+a1z+a2z2+a3z3（凹形，3/2次方）的马鞍状体积的计算式;而形似辛普森公式的式（10），是直形锥状体体积的计算式，两者模型不同，不存在计算精度的可比性，当然也就不存在谁更精确的问题。用以计算凸形弧面状矿体体积，两式都是错误的。此外，两式都名不符实，与“辛普森”、“抛物线”亦无关，所以也宜改称。例如前者称“角壔体公式”，后者称“（3/2次方）马鞍体公式”，以免误导。
4 计算公式的相对偏差
4.1 库容计算式的相对偏差
传统方法建议采用“梯形公式”和截锥公式计算水库容积、工程土方或矿体体积，若实际计算体与此两公式的适用条件（计算图形）不符，则用错公式就会产生偏差[5-6]。若用式（9）、式（10）、式（13）及式（16）分别代表凸形、直形、凹形三类体积的计算式，以V9、V10、V13、V16表示，相邻等高线内面积差为ΔA，按照柯西不等式代换计算，则不同ΔA值时其体积的理论偏差△V见表2。由表2可知，直形体积若误用“梯形公式”计算体积，则当相邻面积差ΔA达到40%时，体积就可能偏大0.468%～5.882%;凹形体如果误用“梯形公式”计算，则将偏大6.302%～11.111%，即使换用截锥公式计算，亦将偏大5.807%～10.593%。
4.2 计算实例
表3为某水库的计算实例。由水库高程—面积曲线求得高程100～110m及120～130m的高程—面积曲线方程式分别为A=21.003+1.624z+0.184z2及A=75.218-1.602z+0.160z2，所以此两段库容形状属直形B型，宜用式（10）计算;高程110～120m的高程一面积曲线方程为A=4.240+5.144z，因而库容形状为凹形A2型，宜用式（9）计算。表中第4列为直接积分计算的数值;第5列为用式（1）根据概化模型所属类型相应的计算式的计算值;表中带（）的数字为按该级实际形状用式（9）的计算值。计算结果表明，若以直接积分计算的数值为标准，则：①用式（1）积分计算与根据概化模型所属类型相应的计算式计算的结果最接近，相差仅0.287%，说明按库容类型选择公式计算的正确性。②B型库容若误用梯形公式计算，则结果偏大，当相邻面积差ΔA>60%时偏大2.178%。ΔA=26.10%～26.61%时偏大0.488%（见末列）。V9值也比V1、V13大，所有偏差值△v与表2的理论值相符。③式（13）为直形库容计算的特例，结果略偏大，在本例中也得到证明，∑V13比∑V10大0.448%。④若减小计算间距，如将h由5m改为2m，则相邻面积差ΔA将减小，各式计算的偏差值ΔV也将相应减小。
5 结语
（1）库容的几何形状有凸形（碗状）、直形（截锥状）及凹形（马鞍状）等多种形状，不同水库之间或一座水库的不同高程之间都可能不相同。形状不同，容积计算式也就各异，没有通用的计算式。
（2）对水库的高程—面积曲线进行分析研究，可以正确判断水库的几何形状，按照其几何形状选择相应的体积计算式正确计算其容积;也可根据高程一面积曲线方程选择相应的容积计算式。
（3）应纠正立体求积与平面求积概念相混淆的现象。体积计算中的“梯形公式”是凸形碗状体的计算公式，长期被误以为是直形锥状体的计算式;“辛普森公式”（抛物线公式）实际上是直形锥状体的计算式，却被误以为是凸形碗状抛物体的计算公式。一些体积计算公式的名称不当，宜改称。
（4）传统的库容计算方法存在严重的错误，我国大部分水库的库容应重新计算，以提高水库效益。某些高校教材、部门规章中关于库容、矿体或工程土方计算的错误“规定”、“建议”应取消，重新制定正确的计算方法与规定。
（5）水平剖面法（等高线法、地形法）有正确的数学基础，并有合理的物理解释，方法简便，应加强研究和完善（如对计算的代数精度和误差分析，编制计算软件）。水平剖面法计算体积的方法和原理同样适用于垂直剖面法（平行斷面法），用以计算河道、渠道及道路等工程土方，或板状矿体，但是其测量方法应该改进，以便能正确判断计算体的几何形状。
（6）识别库容、矿体及工程土方体积的几何形状，按其几何形状计算容（体）积，不仅提高了计算的精度，有很高的经济价值，还可供湖泊学、矿体几何学以及水库水文地理等研究参考。
参考文献：
[1]雷若夫ΠA.矿体几何学[M].北京：地质出版社，1957：26-37.
[2]乌沙阔夫ин.矿藏几何学[M].北京：煤炭工业出版社，1957：442-456.
[3]波塔波夫MB.径流调节[M].北京：高等教育出版社，1956：33-37.
[4]水利部.水库水文泥沙观测试行办法[M].北京：水利出版社，1979：92-97.
[5]华罗庚，王元.关于在等高线图上计算矿藏储量与坡地面积的间题[J].数学学报，1961，11（1）：29-40.
[6]华罗庚，王元.数值积分及其应用[M].科学出版社，1963：26-37.
[7]朱晓岚，何新义，陈予恒.矿体几何学[M].徐州：中国矿业大学出版社，1987：140-144.
[8]陶祖昶.水库、湖泊的容积计算问题[J].人民黄河，1987，9（6）：14-17.
[9]易大义，蒋叔豪，李有法.数值方法[M].杭州：浙江科技出版社，1984：64-74.
[10]罗德仁，邹自力，汤江龙.工程土方量比较分析[J].东华理工学院学报，2005，28（1）：59-63.
[11]米鸿燕，宰建，蒋兴华.静库容计算方法的比较分析[J].地矿测绘，2007，23（2）：1-4.
