摘要：水庫是形狀不規(guī)則的容體，礦體和自然地形上的工程土方是形狀不規(guī)則的巖土體，彼此形態(tài)類似，容（體）積的計算方法相同，用等高線法計算時，都是按照礦體幾何學的方法，假定其為無數(shù)直線圍成的截頭錐狀體，采用梯形公式或截錐公式計算容（體）積。這種方法對計算體的幾何形狀認定不當，對計算公式的模型圖形認識有誤。應(yīng)用數(shù)值積分的方法計算庫容，可以找到識別庫容幾何形狀的方法，導出各種幾何形狀的容積計算式，查明現(xiàn)行各種計算式的模型圖形，糾正傳統(tǒng)計算方法中的錯誤。研究表明：庫容的幾何形狀有凸形、直形或凹形三類形狀，各水庫之間或一座水庫的不同高程之間都不盡相同;幾何形狀不同，容積計算公式也就各異，沒有通用的計算式;若用錯公式，則當相鄰兩計算剖面面積之差為40%時，計算的容積將有0.468%～11.111%的偏差。

關(guān)鍵詞：水庫容積;礦體;工程土方;幾何形狀;計算模型;數(shù)值積分;容（體）積計算;

中圖分類號：TV221.1 文獻標志碼：A

水庫容積、礦體和工程土方體積是水資源開發(fā)利用、礦產(chǎn)開采、工程建設(shè)必需的基本資料。水庫是形狀不規(guī)則的巨大容體，礦體和基于自然地形的工程土方是形狀不規(guī)則的巖土體，彼此形狀類似，體積的計算方法也相同。根據(jù)地形測繪資料或地質(zhì)勘探資料計算容（體）積的方法大致有3類：橫向計算的有垂直剖面法（斷面法、平行斷面法、不平行斷面法）;豎向計算的有水平剖面法（等高線法、等值線法、地形法）;把容體切割成眾多棱柱體計算的有方格法、三角網(wǎng)法以及數(shù)字高程模型法（DEM）。對于庫容或軸狀礦體、山丘土體等塊狀巖土體，多用水平剖面法;對板狀礦體或道路、渠道、河道工程土方等條狀巖土體，則用垂直剖面法。這兩種方法是最常用的方法，基本原理相同，均假定兩個剖面之間的體積是截錐體或角壔體（Prizmatoid），用梯形公式、截錐公式或角壔體公式計算容（體）積。長期以來，由于對這些不規(guī)則物體的幾何形狀及體積計算公式的計算模型缺乏正確的認識，因此在計算中陷人某些誤區(qū)，筆者在前人探討礦體計算的基礎(chǔ)上，根據(jù)對地形法計算庫容的分析研究，找到了識別庫容幾何形狀的方法，導出了各種幾何形狀的庫容計算公式，這些方法和原理同樣適用于礦體和工程土方體積的計算。

1 傳統(tǒng)計算方法及存在的問題

形狀不規(guī)則的礦體，是無法用規(guī)則體積的計算公式計算其體積的。用水平剖面法計算時必須在等高線形狀的基礎(chǔ)上將其轉(zhuǎn)換成符合其總體幾何形狀的、可計算的概化模型。但是，概化模型具有半虛擬的性質(zhì)，其四周界面是無法直接測量到的，于是鮑曼（BaymaH）就假定（1909年）：兩個水平剖面（等高線）之間的礦體是四周界面（圍巖）由無數(shù)斜率不同的直線所圍成的截頭錐狀體，從而以立體幾何的方法推導出了這種錐狀體的體積計算式，即鮑曼公式。索布列夫斯基（Cобулевский）則假定（1934年）：兩平行剖面之間的礦體是由若干不同傾斜度的平面所圍成的角壔體，并用同樣方法推出了角壔體的體積計算式，亦即索布列夫斯基公式[1-2]。這兩個公式本質(zhì)相同，實為一式。鮑曼公式計算較為復(fù)雜，故被逐步舍棄，形式簡便的梯形公式和截錐公式被引入礦體、工程土方以及庫容的計算中，并規(guī)定或建議：相鄰兩剖面等高線內(nèi)面積之差小于40%時用梯形公式，大于40%時用截錐公式;或一般庫容計算用梯形公式，“精確”或“嚴密”計算時用截錐公式[3-4]。

文獻[5-6]用數(shù)值積分的方法推導了鮑曼公式和截錐公式，并與梯形公式相比較，認為此三式都是錐狀體的計算式，但梯形公式計算結(jié)果偏大，鮑曼公式則更具普遍性，比較精確。同時，還推出了馬鞍形礦體的體積計算式。文獻[7]同意以上論述及規(guī)定，但認為礦體周圍的豎向界面大都呈凸形弧面狀，采用拋物線公式（Simpson公式）計算體積比較切合實際，若采用立方拋物線公式（辛普森3/8法則）計算則更加精確。文獻[8]認為庫容剖面的幾何形狀應(yīng)有凸形、直形及凹形等多種（相應(yīng)的體積為碗狀、截錐狀及馬鞍狀），凸形庫容適用梯形公式計算。這些文獻都對計算體形狀的假定及傳統(tǒng)的計算公式提出了質(zhì)疑。

近年來，隨著測繪技術(shù)及計算技術(shù)的發(fā)展，在庫容和土方量計算中使用了許多計算軟件，替代了手工計算，其中有些軟件是按照地形法或斷面法原理編制的，但是其建模中所用的體積計算公式仍如舊規(guī)，所以迄今為止，用水平剖面法或垂直剖面法計算庫容、礦體或工程土方時，仍然遵循百年來傳統(tǒng)的觀點與方法。而由于這些“假定”“規(guī)定”“建議”沒有弄清計算體的幾何形狀，其假設(shè)的形狀未必符合實際，同時對計算公式的模型又不夠了解，混淆了平面求積與空間求積不同的概念，因此造成了體積計算中的許多錯誤，降低了計算的精確度。用數(shù)值積分的方法計算水庫容積并分析研究，不僅可以了解各種計算公式的計算模型，而且可以找到識別庫容幾何形狀的方法，從而糾正傳統(tǒng)計算方法的錯誤，建立比較正確、合理的計算方法與步驟。

2 庫容幾何形狀與數(shù)值積分計算

2.1 根據(jù)水庫高程—面積曲線計算庫容

水庫各高程（z）的面積值（A）是庫區(qū)地形測繪必須提供的基本資料，根據(jù)這些資料就可以直接計算庫容。圖1表明，水庫的高程—面積曲線不是單一、連續(xù)的函數(shù)關(guān)系，而是不連續(xù)的函數(shù)關(guān)系。說明水庫庫容整體而言是一個不規(guī)則的體積，但在函數(shù)連續(xù)的區(qū)間卻是有一定規(guī)則的，可以用數(shù)值積分的方法計算庫容，即將整座水庫按高程—面積曲線函數(shù)的連續(xù)區(qū)間分成幾段，分別計算其庫容，然后累加得到總庫容，從而解決了不規(guī)則體積的積分計算問題。

在水庫高程—面積曲線函數(shù)的連續(xù)區(qū)間，任一微薄層的容積為dV=Adz，高程a～b的容積為

V=∫_a^bAdz （1）

式（1）中若以一個代數(shù)多項式A=a₀+a₁z+a₂z²+…+a_nz_n代替A=f（z），水庫地形圖等高線內(nèi)面積依次為A₀、A₁、A₂、…、A_n[圖2（a）]，高程間距為h，且A（0）=A₀，A（1h）=A₁，…，A（nh）=A_n，則可求得高程—面積關(guān)系方程式不同方次n的體積計算式，如當：1h，則有則有

按牛頓一柯特斯（Newton-Cotes）公式[9]還可得到以下方次的計算式：

由此可以根據(jù)水庫高程—面積曲線方程的方次選擇相應(yīng)的庫容計算式。

2.2 各種幾何形狀的體積計算

在礦體開采或工程土方的挖掘中，不僅要知道其數(shù)量，還要知道其形狀，所以習慣按形狀計算其體積。在庫容計算中了解庫容的幾何形狀也有助于水庫特性的研究及容積的正確計算。

圖2（a）為水庫地形示意圖，O為最低點在圖上的投影。若通過O點以任意θ角作一截面，則可得到此水庫的剖面。剖面可能有A、B、C三種類型，即凸形（碗狀）、直形（截錐狀）及凹形（馬鞍狀）[8]。若通過O點建立柱面坐標，則其體積為

式（7）中剖面線函數(shù)ρ（z，θ）=f（z），也可用一個代數(shù)多項式ρ（z，θ）=a₀+a₁z+z₂z²+…+a_nzⁿ或ρ²（z，θ）=a₀+a₁z+a₂z²+…+a_nzⁿ代替，將不同方次n的方程式代入式（7），即可導得3類剖面形狀的體積計算式，分述如下。

2.2.1 A型（凸形碗狀，n<1）

若n=1/4，即ρ²（z，θ）=a₀+a₁z^1/2，則1h間的體積為

若n=1/2，即則1h間的體積為[8]

2.2.2 B型（直形截錐狀，n=1）（3個剖面）間的體積為

亦即索布列夫斯基公式：

若取兩根等高線計算1h間的體積，則就是鮑曼公式[5-6]：

若式（11）中ρ₂（θ）=kρ₁（θ），k為常數(shù)，則式（11）可改寫成常用的截錐公式（棱臺公式）：

2.2.3 C型（凹形馬鞍狀，n>1）

若n=3/2，的體積為[7，9]則有

按牛頓一柯特斯（Newton-Cotes）公式[9]還可得到以下方次的計算式：

由此可以根據(jù)水庫高程—面積曲線方程的方次選擇相應(yīng)的庫容計算式。

2.2 各種幾何形狀的體積計算

在礦體開采或工程土方的挖掘中，不僅要知道其數(shù)量，還要知道其形狀，所以習慣按形狀計算其體積。在庫容計算中了解庫容的幾何形狀也有助于水庫特性的研究及容積的正確計算。

圖2（a）為水庫地形示意圖，O為最低點在圖上的投影。若通過O點以任意θ角作一截面，則可得到此水庫的剖面。剖面可能有A、B、C三種類型，即凸形（碗狀）、直形（截錐狀）及凹形（馬鞍狀）[8]。若通過O點建立柱面坐標，則其體積為

式（7）中剖面線函數(shù)ρ（z，θ）=f（z），也可用一個代數(shù)多項式ρ（z，θ）=a₀+a₁z+z₂z²+…+a_nzⁿ或ρ²（z，θ）=a₀+a₁z+a₂z²+…+a_nzⁿ代替，將不同方次n的方程式代入式（7），即可導得3類剖面形狀的體積計算式，分述如下。

2.2.1 A型（凸形碗狀，n<1）

若n=1/4，即ρ²（z，θ）=a₀+a₁z^1/2，則1h間的體積為

若n=1/2，即則1h間的體積為[8]

2.2.2 B型（直形截錐狀，n=1）（3個剖面）間的體積為

亦即索布列夫斯基公式：

若取兩根等高線計算1h間的體積，則就是鮑曼公式[5-6]：

若式（11）中ρ₂（θ）=kρ₁（θ），k為常數(shù)，則式（11）可改寫成常用的截錐公式（棱臺公式）：

2.2.3 C型（凹形馬鞍狀，n>1）

若n=3/2，的體積為[7，9]

若則2h間的體積為[6]

若剖面線為冪函數(shù)ρ（z，θ）=a₀+a₁z²，則可得到1h間容積的計算式為[8]

若n=2，的體積為

上述公式中式（3）、式（4）、式（5）與式（9）、式（10）、式（14）傳統(tǒng)上稱之為梯形公式、角壔體公式（拋物線公式、辛普森公式、索布列夫斯基公式）及立方拋物線公式（辛普森3/8法則）。式（11）與式（10）實為一式。式（13）又稱為棱臺體公式。式（11）中的A_h/2為h/2處的面積。式（12）、式（15）、式（16）中的T（1，2）、T（2，3）、T（3，1）為面積改正數(shù)[6]。式（2）～式（6）及式（8）～式（17）兩類計算式同列于表1，以供對照比較。這兩種計算庫容的方法同樣適用于垂直剖面法計算庫容或河道、渠道及道路等工程土方體積[8]，其原理相同。

2.3 庫容幾何形狀的識別

庫容的幾何形狀無法由工程測量直接測定，但可以通過間接的方法判斷識別。式（1）為根據(jù)水庫高程—面積曲線積分計算庫容，可以求得高程—面積曲線方程不同方次的庫容計算式;式（7）則以庫容剖面線形狀作為變量，可以求得剖面線方程不同方次的庫容計算式;很容易證明兩者實為一式，僅是座標形式不同而已，表一中由式（1）導得的式（2）、式（3）、式（4）、式（5）、式（6）分別與由式（7）導得的式（8）、式（9）、式（10）、式（14）、式（17）相同，這樣就找到了高程一面積曲線的形狀與庫容剖面線圖形之間的關(guān)系，只要求出高程—面積曲線A=f（z）的方程式，就可以知道庫容模型的剖面線方程P=f（z）及其幾何圖形。例如，若高程—面積曲線呈直線，方程式方次n=1，則庫容剖面線方程為n=1/2的凸形曲線，幾何形狀為碗狀（饅頭狀）的凸形體，計算公式適用式（3）、式（9）;若高程—面積曲線為2次方程，則剖面線必為斜率不等的直線，體形屬直形的角壔體或錐狀體，計算公式用式（4）、式（10）。

庫容剖面線的幾何圖形是半虛擬的概化模型，無法直接測定，但可由高程，面積曲線反映出來。而高程一面積曲線是可以根據(jù)測繪資料繪制，而且不難求得其曲線方程，加以上述各種形狀庫容計算式的導出，于是就可以走出盲目假定幾何形狀、亂用體稠十算公式的誤區(qū)。

3 傳統(tǒng)方法計算容（體）積的幾個誤區(qū)

3.1 誤區(qū)一

“兩平行剖面間礦體（庫容）的形狀都是截錐體或角壔體，采用梯形公式或截錐公式計算體積”“相鄰兩剖面面積之差小于40%時用梯形公式，大于40%時用截錐公式”“一般計算用梯形公式，精密計算用截錐公式”“鮑曼公式更具普遍性，比較精確”[1-6]。

圖1表明，水庫高程—面積曲線在不同區(qū)段具有上凸、直線或下彎等多種形狀，根據(jù)高程—面積曲線形狀與剖面線形狀的關(guān)系，說明庫容的剖面形狀也相應(yīng)有凸形、直形或凹形，即體積為碗狀、截錐狀或馬鞍狀等形狀，并非只是直形的錐狀體（角壔體）一種。水庫的大斷面圖也常表明這種情況。此外，從地質(zhì)學、地貌學、土力學及河流動力學也可得到合理的解釋。礦體、山丘等巖土體同樣如此。不同形狀的區(qū)間要用不同的體積計算式，沒有普遍適用的通用公式。

以上演算表明，截錐公式是鮑曼公式的特例，是計算正截頭錐體（四周斜率相同）的公式，不是什磨特殊的精密計算公式。與“梯形公式”計算圖形又不相同，無法比較。體積計算用哪個公式合適、正確，不在于相鄰剖面面積差的百分數(shù)是多少，而在于計算公式的數(shù)學模型與計算體的幾何形狀是否一致。因此，不問計算體實際的幾何形狀，假定其都是截錐體、角壔體，一律采用梯形公式或截錐公式計算體積，顯然是錯誤的。

鮑曼公式與截錐公式是計算2個剖面間的體積，而索布列夫斯基公式（角壔體公式）是計算3個剖面間的體積，都是直形錐狀體的計算式。除截錐公式屬特例（其計算體積略偏大）外，鮑曼公式與其他兩式實為一式，只要與計算體模型相符，都一樣精確。

3.2 誤區(qū)二

“梯形公式是剖面為梯形的截頭錐體體積的計算公式”[5-6]。

梯形公式本是計算梯形平面面積的求積公式，將其改寫成體積計算公式，認為其計算的圖形是剖面為梯形的錐狀體，沒有根據(jù)，也沒見其來源的報導。而式（9）的演算則證明，體積計算中的“梯形公式”，其計算模型是具有凸形剖面的碗狀體，剖面線方程為ρ²（z，θ）=a₀+a₁z，與“梯形”毫不相干，將此式與截錐公式及鮑曼公式歸人一類，是一個誤導，混淆了體積計算與面積計算的不同概念。文獻[8]早就指出這一情況，但沒有引起普遍注意。此式是國內(nèi)外礦業(yè)、水利、土木工程界計算礦體、庫容、工程土方等體積中廣泛使用的計算式。長期以來，也有人發(fā)覺用此式計算的結(jié)果大于截錐公式或其他錐狀體計算式，也大于其他計算方法[6，10-11]，但由于對其形狀的誤解，因此始終得不到正確的解釋。“梯形”是平面圖形，沒有體積的。式（9）與“梯形”又毫不相干，將其稱為“梯形公式”，名不符實。所以，為避免誤解誤用，在體積計算中改稱“碗狀體公式”較妥。

3.3 誤區(qū)三

對天然界面呈凸形的弧面狀礦體“采用辛普森公式計算體積比較切合實際;如果采用立方拋物線公式計算，則更加精確”[7]。

辛普森公式原是計算二次拋物線下曲邊梯形平面面積的公式，所以又稱拋物線公式。體積計算中形似辛普森公式的式（10）演算證明，它是角壔體的體積計算式，計算模型是3個水平剖面（2h）間的錐狀體，垂直剖面線為直線。把它與辛普森公式混淆，以為是計算表面為凸形，縱剖面線為拋物線的體積計算式，這是又一個誤解。實際上是其水平剖面的面積值變化呈二次拋物線，而不是其垂直剖面線呈二次拋物線。

辛普森3/8法則是三次拋物線下的平面面積插值計算式，在插值計算中，其余項誤差比辛普森公式小，所以插值精度比辛普森公式高。但在體積計算中，形似辛普森3/8法則的式（14）（有稱之為“立方拋物線公式”）實際上是剖面線方程為ρ²（z，θ）=a₀+a₁z+a₂z²+a₃z³（凹形，3/2次方）的馬鞍狀體積的計算式;而形似辛普森公式的式（10），是直形錐狀體體積的計算式，兩者模型不同，不存在計算精度的可比性，當然也就不存在誰更精確的問題。用以計算凸形弧面狀礦體體積，兩式都是錯誤的。此外，兩式都名不符實，與“辛普森”、“拋物線”亦無關(guān)，所以也宜改稱。例如前者稱“角壔體公式”，后者稱“（3/2次方）馬鞍體公式”，以免誤導。

4 計算公式的相對偏差

4.1 庫容計算式的相對偏差

傳統(tǒng)方法建議采用“梯形公式”和截錐公式計算水庫容積、工程土方或礦體體積，若實際計算體與此兩公式的適用條件（計算圖形）不符，則用錯公式就會產(chǎn)生偏差[5-6]。若用式（9）、式（10）、式（13）及式（16）分別代表凸形、直形、凹形三類體積的計算式，以V₉、V₁₀、V₁₃、V₁₆表示，相鄰等高線內(nèi)面積差為Δ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.184z²及A=75.218-1.602z+0.160z²，所以此兩段庫容形狀屬直形B型，宜用式（10）計算;高程110～120m的高程一面積曲線方程為A=4.240+5.144z，因而庫容形狀為凹形A₂型，宜用式（9）計算。表中第4列為直接積分計算的數(shù)值;第5列為用式（1）根據(jù)概化模型所屬類型相應(yīng)的計算式的計算值;表中帶（）的數(shù)字為按該級實際形狀用式（9）的計算值。計算結(jié)果表明，若以直接積分計算的數(shù)值為標準，則：①用式（1）積分計算與根據(jù)概化模型所屬類型相應(yīng)的計算式計算的結(jié)果最接近，相差僅0.287%，說明按庫容類型選擇公式計算的正確性。②B型庫容若誤用梯形公式計算，則結(jié)果偏大，當相鄰面積差ΔA>60%時偏大2.178%。ΔA=26.10%～26.61%時偏大0.488%（見末列）。V₉值也比V₁、V₁₃大，所有偏差值△v與表2的理論值相符。③式（13）為直形庫容計算的特例，結(jié)果略偏大，在本例中也得到證明，∑V₁₃比∑V₁₀大0.448%。④若減小計算間距，如將h由5m改為2m，則相鄰面積差ΔA將減小，各式計算的偏差值ΔV也將相應(yīng)減小。

5 結(jié)語

（1）庫容的幾何形狀有凸形（碗狀）、直形（截錐狀）及凹形（馬鞍狀）等多種形狀，不同水庫之間或一座水庫的不同高程之間都可能不相同。形狀不同，容積計算式也就各異，沒有通用的計算式。

（2）對水庫的高程—面積曲線進行分析研究，可以正確判斷水庫的幾何形狀，按照其幾何形狀選擇相應(yīng)的體積計算式正確計算其容積;也可根據(jù)高程一面積曲線方程選擇相應(yīng)的容積計算式。

（3）應(yīng)糾正立體求積與平面求積概念相混淆的現(xiàn)象。體積計算中的“梯形公式”是凸形碗狀體的計算公式，長期被誤以為是直形錐狀體的計算式;“辛普森公式”（拋物線公式）實際上是直形錐狀體的計算式，卻被誤以為是凸形碗狀拋物體的計算公式。一些體積計算公式的名稱不當，宜改稱。

（4）傳統(tǒng)的庫容計算方法存在嚴重的錯誤，我國大部分水庫的庫容應(yīng)重新計算，以提高水庫效益。某些高校教材、部門規(guī)章中關(guān)于庫容、礦體或工程土方計算的錯誤“規(guī)定”、“建議”應(yīng)取消，重新制定正確的計算方法與規(guī)定。

（5）水平剖面法（等高線法、地形法）有正確的數(shù)學基礎(chǔ)，并有合理的物理解釋，方法簡便，應(yīng)加強研究和完善（如對計算的代數(shù)精度和誤差分析，編制計算軟件）。水平剖面法計算體積的方法和原理同樣適用于垂直剖面法（平行斷面法），用以計算河道、渠道及道路等工程土方，或板狀礦體，但是其測量方法應(yīng)該改進，以便能正確判斷計算體的幾何形狀。

（6）識別庫容、礦體及工程土方體積的幾何形狀，按其幾何形狀計算容（體）積，不僅提高了計算的精度，有很高的經(jīng)濟價值，還可供湖泊學、礦體幾何學以及水庫水文地理等研究參考。

參考文獻：

[1]雷若夫ΠA.礦體幾何學[M].北京：地質(zhì)出版社，1957：26-37.

[2]烏沙闊夫ин.礦藏幾何學[M].北京：煤炭工業(yè)出版社，1957：442-456.

[3]波塔波夫MB.徑流調(diào)節(jié)[M].北京：高等教育出版社，1956：33-37.

[4]水利部.水庫水文泥沙觀測試行辦法[M].北京：水利出版社，1979：92-97.

[5]華羅庚，王元.關(guān)于在等高線圖上計算礦藏儲量與坡地面積的間題[J].數(shù)學學報，1961，11（1）：29-40.

[6]華羅庚，王元.數(shù)值積分及其應(yīng)用[M].科學出版社，1963：26-37.

[7]朱曉嵐，何新義，陳予恒.礦體幾何學[M].徐州：中國礦業(yè)大學出版社，1987：140-144.

[8]陶祖昶.水庫、湖泊的容積計算問題[J].人民黃河，1987，9（6）：14-17.

[9]易大義，蔣叔豪，李有法.數(shù)值方法[M].杭州：浙江科技出版社，1984：64-74.

[10]羅德仁，鄒自力，湯江龍.工程土方量比較分析[J].東華理工學院學報，2005，28（1）：59-63.

[11]米鴻燕，宰建，蔣興華.靜庫容計算方法的比較分析[J].地礦測繪，2007，23（2）：1-4.