2010年 全 国 大 学 生
数 学 建 模 竞 赛
储油罐的变位识别与罐容表标定
2010高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则．
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题．
我们知道，抄袭别人的成果是违反竞赛规则的， 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出．
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性．如有违反竞赛规则的行为，我们将受到严肃处理．
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全
赛区评阅编号（由赛区组委会评阅前进行编号）：
2010高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
储油罐的变位识别与罐容表标定
摘 要
加油站储油罐罐容表的精确度直接关系到加油站的经济利益，然而由于地基变形等原因，罐体的位置会发生纵向倾斜和横向偏转等变化，从而导致罐容表发生改变，影响其精度．
本文要解决的就是储油罐的变位识别与罐容表标定的问题．其中罐容表的标定，就是建立罐内油位高度与储油量的关系．对此，我们应用微积分及空间解析几何理论的相关知识，建立油罐体积函数模型V (H )．对于储油罐的变位识别问题，我们借助已建立的函数模型V (H )，用实际的油位高度确定理论储油量和变位参数值，并将理论储油量与实际测出的储油量采用最小二乘法进行拟合，然后通过拟合系数来判断模型的准确性．
对问题(1)，储油罐有无变位和纵向变位这两种情况，均要建立油罐体积积分函数模型，并运用matlab 软件求解模型，且将求解结果采用最小二乘法拟合，分析结果表明理论结果与数据模拟结果相吻合．最小二乘法拟合分析时也表明了模型求解中存在误差，从而以此为基础对模型进行修正，并得出罐容标定值表(见表一) ．
对问题(2)，同样建立建立油罐体积积分函数模型，采用离差平方和的算法并运用matlab 软件确定了变位参数α、β的值为α=20、β=4.90．以此为基础给出罐体变位后罐容表标定值表(见表三) ．
对问题一和二的模型做误差分析和修正后所得的结果显示，我们所建立的模型能很好的与实际情况相吻合，其吻合系数达到0.9996．最后我们还对模型进行了正确性验证与方法可靠性检验，并结合实际情况和应用价值对模型进行了改进与推广．
关键词 微积分；变位识别；小二乘拟合；误差；标定值
一、问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表(即罐内油位高度与储油量的对应关系) 进行实时计算，以得到罐内油位高度和储油量的变化情况．
许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化(以下称为变位) ，从而导致罐容表发生改变．按照有关规定，需要定期对罐容表进行重新标定．一种典型的储油罐其主体为圆柱体，两端为球冠体．现需要用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题．
（1）为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐(两端平头的椭圆柱体) ，分别对罐体无变位和倾斜角为α=4．10的纵向变位两种情况做了实验，实验数据如附录一所示．现需要建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm 的罐容表标定值．
（2）对于主体为圆柱体，两端为球冠体的储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数(纵向倾斜角度α和横向偏转角度β ) 之间的一般关系．利用罐体变位后在进/出油过程中的实际检测数据附录二，根据所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm 的罐容表标定值．进一步利用附录二中的实际检测数据来分析检验模型的正确性与方法的可靠性．
二、问题分析
2.1 问题一的分析
通常情况下，我们都可以通过油位计管理系统来标定罐容表，即通过测量进/出油量与罐内油位高度得到储油量的变化情况．但是储油罐在使用一段时间后，由于地基变形等原因，罐体的位置会发生纵向倾斜和横向偏转，导致罐容表测量值发生改变，这就需要我们定期对罐容表进行重新标定.
对罐容表进行重新标定前，首先我们必须判定油罐是否发生变位，即将测量的进出口油量实际值与罐体位置未发生变位时的理论值进行差分拟合，当读数误差达到一定值时，就可以判定罐体的位置是发生了变位．而对罐体无变位时罐容表的识别，可以通过对罐体已知的几何结构进行分析计算，确定罐体无变位时储油量V 与可测油位高度H 之间的函数关系V =f (H )．其次必须解决变位后罐容表如何重新标定的问题．
要解决上述问题，我们必须先建立V (H , α)的函数模型. 在建立V (H , α)的函数模型过程中, 我们参照了高等数学微积分的相关知识, 采用微元的思想得出模型.
2.2 问题二的分析
实际情况中，储油罐不单单只发生纵向倾斜，纵向和横向倾斜也应考虑，所以该问题中的情形比问题一更具有实际意义.
该问题是在问题一的基础之上增加了对横向倾角β的考虑，也就是要求我们同时考虑三个变量对储油量的影响，建立V (H , α, β).
[1]
根据事物的变化规律，针对该倾斜问题，我们发现:在两种倾斜同时发生时的结果与分步依次发生的结果是相同的，这就启发了我们可以通过分步考虑来简化模型的建立. 接着我们又考虑到，该问题中储油罐是圆柱体和球冠体这样两个特殊的对称体的组合体. 分析其几何特征可知:罐内液体不管怎么横向倾斜，其横截面均为垂直于水平面、左右对称的薄片, 也就是说横向变位对纵向变位储油量无影响.
所以为了易于模型的建立，我们假设油罐每次倾斜的完成顺序均如下图:
这样该问题中模型的建立又可以直接参照问题一中模型的建立，最后得V (H , α, β)函数模型.
模型建立的基本思路如下:
模型建立完后，得出储油量与油位高度及变位参数（纵向倾斜角度与横向倾斜角度）之间的关系模型. 我们就可以开始确定变位参数α、β的值．再将确定了变位参数α、β的值后代入模型来给出罐体变位后油位高度间隔为10cm 的罐容表标定值．最后可以再提出用实际检测数据对模型分析检验与建模方法的可靠性验证的方法．
三、模型假设及符号定义与说明
3.1模型假设
1．罐内储油不受温度压强等的影响，即储油量的体积大小只与油位高度有关；
2．油浮子为一质点，其大小可忽略不计；
3．储油罐壁的厚度很薄，可以忽略不计；
4．外界因素的改变不会影响储油罐的形状，即不会发生形变；
5．储油罐内部罐壁为理想、对称的几何图形，忽略其制造工艺带来的误差；
6．储油罐内部一系列小构建对储油量的影响忽略不计；
7．油的自身性质和蒸发损耗对储油量的影响忽略不计；
8．对油浮子与油接触时带来不可避免的仪器误差忽略不计．
3.2符号定义与说明
α:储油罐的纵向倾角；
β:储油罐的横向倾角；
l :储油罐的罐长；
V :储油量；
a :小椭圆油罐截面中椭圆的长半轴；
b :小椭圆油罐截面中椭圆的短半轴；
H :油浮子所在油面处的油位高度；
h :储油罐中任一油面位置处的油位高度；
V C (H ):测量出的储油量与油位高度关系函数；
V L (H ):实际的储油量与油位高度关系函数．
四、关于小椭圆型储油罐的模型建立与求解
为了层次清楚，我们先交待本节的结构．根据储油量的多少，以及油浮子位置的限制，对于近油位探针端下倾这种情形，分成如下几种情况进行考虑:
情形I 0≤H ≤m tan α，储油情况如图1-4所示，并建立模型Ⅱ；
情形II m tan α≤H ≤2b -n tan α，储油情况如图1-1所示，并建立模型Ⅰ； 情形III 2b -n tan α≤H ≤2b 储油情况如图1-5所示，并建立模型Ⅲ．
4.1 情形I 时，罐容与油面高度关系的模型建立
通常状况下，α在较小的情况下就会被工作人员所发现，并重新摆放储油罐，所以m tan α≤H ≤2b -n tan α时，比较常见，据此我们对此种情景在此做重点介绍．
如图1-1所示，取椭球圆柱体的中心轴为z 轴，并设该立体在过点z =0、z =l 且垂直于z 轴的两平面之间．以S 1(y ) 表示过z 且垂直与x 轴的截面面积．这时，取z 为积分变量，它的变化区间为[0, l ]；相应于[0, l ]上任一小区间[z , z +dz ]的一薄片的体积近似于底面积为S 1(y ) ，高为dz 的扁柱体的体积，即体积元素
dV 1=S 1(y ) dz ．
图1-1
以S 1(y ) dz 为被积分表达式，在闭区间[0, l ]上作定积分，便得所求立体的体积
V 1(y )=⎰S (y ) dz ． (1) 01l
接下来我们建立S 1(h ) 函数关系．
如图1-2所示，以图中椭圆柱体最左端椭球截面中心点为原点，以平行水平面和垂直水平面的方向分别为x 、y 轴建立直角坐标系，其中每片截面投影到xoy 坐标轴上的如图所示．
图1-2
图中椭圆面积公式为
x 2y 2． +
=1⇒x =22a b 现在，取纵坐标y 为积分变量，它的变化区间为[-b , y ]．相应于[-b , y ]上任一小区间[y , y +dy ]的窄条面积近似于高为dy
、底为积元素
⎛dS 1= dy ， ⎝从而得
S 1(
y )=⎰-b y ， 由图1-2有y =h -b ，代入得
S =⎰
接下来以h -b -b ． 运用MATLAB 程序，在区间[-b , h -b ]上作定为被积分表达式，积分，得所求的S 1(h ) 函数表达式为(程序见附录一)
S 1(
h )=⎰h -b
-b b h -b π=(
h -b ab arcsin +ab ． (2) a b 2
进而，我们通过建立H 与h 的函数关系，将H 引入到S 1(h ) 中，建立V 1(H )函数模
型．
选取yoz 坐标面上的截面如图1-3所示，油浮子所在处油位高度为H ，对应在y 轴上的投影点为C 点， 在z 轴上的投影长度为n ；油面上任一点的油位高度为h ，对应在y 轴上的投影为点B ，在z 轴上的投影长度为z ，即为该点的纵坐标大小；D 点为油面在y 轴上的交点．
图1-3
很明显的有线段长度关系
AB +BD =AC +CD ,
而
AB =h , BD =z tan α, AC =H , CD =n tan α,
所以
h +z tan α=H +n tan α⇒h =H +(n -z )tan α. （3） 对上述定积分公式(2)计算时，先不考虑积分限，直接对S 1(h )做不定积分，即 ⎰S (h )dh
1
=cot α
⎫h -b 1ab arcsin -πab (h -b )⎪+C ． b 2⎭
联立(1)，(2)和(3)式，因为所得结果比较复杂，为了简便起见，我们在此令k =H +n tan α-b ，j =H -m tan α-b ，所以积分后的式子为
V 1(H , α)=⎰H -m tan α-b H +n tan α-b S 1(h )dh
=ab cot α
arcsin j π⎫(m -n ) tan α⎪
b 2⎭
k ． b -ab cot α
arcsin
4.2 情形II 时，罐容与油面高度关系的模型建立
如图1-4所示，即为0≤H ≤m
tan α时的储油情形．
图1-4
此种情形下模型的建立与模型Ⅰ的建立基本相同，唯一不同的是z 轴方向上的积分上限:模型Ⅰ中的上限为罐长l ，而此处模型中油面边缘最右端与罐下壁有交界，投影到z 轴上的交点即为(0,0, H cot α+n )，所以该模型上限为H cot α+n ．即有
V 2(y )=⎰Hcot α+n 0S 2(y ) dz ． (4) 所求的S 2(h ) 函数表达式与模型Ⅰ中完全相同，即为
S 2(
h )
=⎰b h -b π=(h -b ab arcsin +ab ． (5) -b a b 2H 与h 的函数关系也为 h -b
h =H +(n -z )tan α ． (6)
同样联立(4)(5)(6)三式求解，为简便起见，我们也在此令k =H +n tan α-b ，同理得函数模型
V 2(H , α)=⎰0H +l tan α-b S 2(h )dh
=ab cot α
arcsin ⎫h -b 1-π(h -b )⎪b 2⎭0H +l tan α-b
=-ab cot α
arcsin k 1⎫-π(k +b +1)⎪． b 2⎭
4.3 情形Ⅲ时，罐容与油面高度关系的模型建立模型
如图1-5所示， 即为2b -n tan α
≤H ≤2b 时的储油情形．
图1-5
此种情形下模型的建立也与模型Ⅰ的建立基本相同，与模型Ⅰ相比，该模型相当于是模型Ⅰ中储油油体形状的立体图与一椭圆柱体的组合，所以该模型体积的求解分两部分完成，具体如下：
V 3=V 3-1+V 3-2．
其中V 3-1为椭圆柱体的体积，V 3-2为似模型Ⅰ的体积．
在求解V 3-1的函数关系式时，利用椭圆柱体的体积公式：体积=底⨯高．
求解时，为解释更加清楚，我们将yoz 坐标面上的图形截出平放如图1-6所示．
图1-6
观察图形可得如下线段关系式
AB =AC -BC , BC =EF ⨯cot α, EF =BF -BE ，
而
AC =n , BE =CG =H ．
所以最终可得椭圆柱体的高为n -(2b -H )cot α．
在求解底面时，我们直接取用椭圆的面积公式，得面积为πab ，所以有
V 3-1=πab (n -(2b -H )cot α) ． (7)
在求解V 3-2的函数关系式时，我们参照模型Ⅰ的求解过程，抓住其本质的不同之处，
仅将其积分下限换为n -(2b -H )cot α即得
V 3(y )=⎰l
n -(2b -H )S 3(y ) dz ．
所求的S 3(h ) 函数表达式也与模型Ⅰ中完全相同，即为
S 3(
h )=b =(
h -b ⎰-b a
H 与h 的函数关系也为 h -b h -b π (8) ab a r c +ab ．b 2
h =H +(n -z )tan α． (9)
同样联立(7)(8)(9)三式求解，为简便起见，我们也在此令k =H +n tan α-b ，同理得函数模型
V 3(H , α)=⎰H -m tan α-b H -b S 3(h )dh +πab (n -(2b -H )cot α)
=ab cot α
j 1⎫arcsin πm tan α⎪
b 2⎭
-ab cot α
arcsin H -b - b +πab (n -(2b -H )cot α)．
4.4情形Ⅳ时，罐容与油面高度关系的模型建立模型
如图1-7所示, 即为α=0时的情形.
图1-7
很明显，影响储油量V 的只有油面高度H ，所以我们直接建立V 与H 的函数表达式V (H ), 以下即为函数V (H )的建立：
此种情形下模型的建立与模型Ⅰ的建立基本相同, 则有
V 4(y )=⎰S 4(y ) dz ． (10) 0l
所求的S 4(h ) 函数表达式与模型Ⅰ中完全相同, 即为
b h -b π=(
h -b ab arcsin +ab ． （11) -b a b 2
由于α=0，所以液面各处H 与h 均相等，即有
h =H ． (12)
同样联立(10)(11)(12)三式求解得
bl H -b πV 4(H )=(
H -b abl arcsin +abl ． a b 2
4.5 一些补充说明
1、除了以上所建立的三种模型外，我们也考虑到了其他可能会有的情况，如图所示图1-8和图1-9，但是考虑到油浮子的测量局限性，这两种情况油浮子无法测量，所以我们在此也不做考虑．
S 2=⎰h -b
图1-8
图1-9
2、发生纵向倾斜时，可能为近油位探针端下倾，也可能为远油位探针端下倾，以上考虑的仅为近油位探针端下倾的情形．若出现远油位探针端下倾这一情形，对可测得油面高度H 的储油量函数可采用如下方法进行计算．
对于小椭圆柱体型储油罐这样的对称体，如图1-10所示，假设储油罐内有两个油浮子，分别位列储油罐内两对称的位置．并假设仅远油位探针端的油浮子可读，为H ，则另一油浮子，即近油位探针端油浮子油位高度为2a -H ，即对应的储油量函数直接套用以上模型即为：
V (H )=f (2a -H ) ．
但因远离油位探针端下倾时，微小的油位变化就会引起储油量发生很大的变化，实际工作中会很快被相关工作人员所发现，并重新放置．即远离油位探针端下倾没有太大的实际意义，所以我们在此不做进一步讨论．
图1-10
4.6 罐容与油面高度关系模型的求解
4.6.1、罐体变位后对罐容表的影响的求解
所谓罐体变位后对罐容表的影响，即考虑当的油位高度H 固定，为常数时，纵向倾角α对罐体储油体积V (H , α)的影响．可以通过理论值与实测值之间的差∆V (H , α)来判断
∆V (H , α)=V 2(H , α)-V C (H , α)．
V (H , α)是关于纵向倾角α的增函数，V (H , α)对上式拟合分析得，即α值增大时，
值增大．
4.6.2、给出α=4.10时油高间隔为1cm 所对应的一系列的罐容标定值
在使用已建立出的模型做标定之前，为确保结果的精确度，我们先采用差值拟合的方法对模型进行修正．即对油罐储油体积V (H )理论值与实际测量值的差做拟合曲线，也就是建立实测值与模型值的误差函数．
又因为采用数据拟合的方法可以反映函数曲线（面）V =f (H )反映对象整体的变化趋势，且使f (H )在某种准则下与所有实际测量值最为接近，即曲线拟合得最好．于是我们将实际测量值的数据点用matlab 拟合，在用三次拟合时，三次相前系数几乎为0，且做二次拟合时，相关系数r=0．9996，精度较高，说明拟合效果较好，故这里我们只采用中二次拟合．
现在计算误差函数．油罐储油体积无论是理论值还是实际测量值都与油位高度H 有关，所以误差函数也是H 的函数． 拟合时因进油表数值和出油表数值均为外部仪器测量，其数值较为精确，故采用进油表数值或采用出油表数值不会影响拟合效果．
用matlab [2]做二次曲线拟合[3](程序见附录二、三) 得出误差函数曲线方程为：
（1）未发生变位时，
实测值与模型值的误差函数为：
∆V (H )=V C (H )-V L (H )=1.6827⨯10-5H 2-0.15693H +17.9826．
故进行修正以后模型的函数为:
V (H )=V C (H )-∆V (H )．
用matlab 编程有修正前后拟合曲线如下图1-11所示
图1-11
（2）发生变位时
实测值与模型值的误差函数：
∆V (H )=V C (H )-V L (H )=0.00039739H 2-0.58342H +124.2537．
修正以后的模型的函数为:
V (H )=V C (H )-∆V (H )．
修正前后拟合曲线如下图1-12所示
:
V V 2现在采用相位分析法对修正后模型相似度进行检验，即用δ=cos θ=1计算得1V 2图1-12 δ=0.99，这也就证明了模型的准确性.
从上可知修正后函数模型与实际情况吻合系数较高，符合实际情况，现利用matlab 编程给出罐体变位后油位高度间隔为1cm 的罐容表标定值定标[4]表如下表一:
五、关于实际储油罐的模型建立与求解
5.1实际储油罐的模型建立
5.1.1 建立H ' (H , ) 关系式
如图2-1所示， 取出油浮子所在处的截面，并以其下端点为原心，该处切线方向为x 轴，垂直x 轴方向为y 轴，建立直角坐标系．其中平行与x 轴的直线为油面所在水平线， 所以B 点为油浮子在y 轴上的投影点，同时我们设其在y 轴正半轴上的投影高
度为H ' ．油浮子的测量高度仍然设为H ．
图2-1
由图可得如下线段关系
OB =OA +AB , AB =AC cos β, AC =CD -AD ，
而
CD =H , AD =r ，
最后可得
H ' (H , β) =(H -r )cos β+r ． (13)
5.1.2 建立V (H ' , S , α)关系式
接下来考虑纵向倾斜时，我们只需利用H ' ，结合微元积分的思想，建立函数关系式V (H ' , α, β)，最终通过(13)式将H 引入即得我们所需要的V (H , α, β)函数模型．
与问题一的思路相同，首先，我们根据储油量的多少，以及油浮子位置的限制，对于近油位探针端下倾这种情形，分成如下几种情况进行考虑:
0≤H '
我们以储油罐最下端切线方向为y 轴，以过储油罐最左端点且垂直于y 轴，并切于该点的指向上的直线为z 轴， 以y 轴与z 轴交点为原点，以过原点且垂直于y 轴和z 轴的直线为x 轴建立空间直角坐标系，如图2-2所示．同时设油面与储油罐的罐壁交点分别为点E 和点D ，观察图形得积分方向y 轴的五个区域:
0≤y ≤y E ，y E ≤y ≤1，1≤y ≤9，9≤y ≤y D ，y D ≤y ≤10．
图2-2
我们先考虑y E >1.5，即左端处一直存在一个由过点E 且垂直于y 轴的面所截出小球冠体，而对于0≤y E ≤1.5这种情况我们将在后面单独做以交代．
（1）当0≤H '
V =⎰S 1dy +⎰S 2dy +⎰0y E y E 13+H ' cot α1S 3dy ． (14)
（2）当6tan α≤H '
y E 19y D V =⎰S 1dy +⎰S 2dy +⎰S 3dy +⎰S 4dy ． (15) 0y E 19
（3）当7tan α+1.5≤H '
V =⎰S 1dy +⎰S 2dy +⎰S 3dy +⎰S 4dy +⎰S 5dy ． (16) 0y E 19y D y E 19y D 10
（4）当3-2tan α≤H '
现需要解决的就是求解D 、E 点的坐标, 来确定各模型表达式中积分上下限． 很明显点y D 和y E 是直线DE 与两个部分球面体的交点，根据初等数学的相关知识，
我们分别建立直线以及球面方程，最终联立求解得D 、E 坐标．
（一）直线方程的建立
如图2-2所示，有如下的线段关系:
OC =OA +AB +CB ，
而
OA =1, AB =2, CB =H ' cot α，
所以得
OC =3+H ' cot α．
即C 点的坐标为(0,3+H ' cot α,0)．
又已知∠ECA =α，可得其斜率为-tan α，再根据点斜式方法，可得该直线的方程为：
z =-tan αy +H ' +3tan α ． (18)
（二）球面方程的建立
如图2-3所示，设球体半径为R ．
图2-3
由图可得如下线段关系式
O ' N 2+NF 2=O ' F 2, O ' N =O ' P -NP ，
由于
O ' P =O ' F =R , NP =1，
得
(R -1)2+1.52=R 2⇒R =1.625．
则左半边球冠体球心O ' 的坐标为(0,1.625,1.5)，同理可求得右半边球冠体球心的坐标为(0,8.375,1.5)．
左半边球冠体所在球的方程：
(y -1.625)2+(z -1.5)=1.6252． (19) 2
右半边球冠体所在球的方程：
(y -8.375)
2
+(z -1.5)=1.6252． (20)
2
点E (y E , z E )既在左半边球冠体横截面圆的曲线上，同时也过直线CE ，为求E 点的坐标y E ，则可以联立方程(18)和(19)求解得：
y E =
3.25+2tan α(H ' +3tan α-1.5)
2sec α
2
同理可联立方程(18)和(20)求解得：
y D =
16.75+2tan α(H ' +3tan α-1.5)
2sec α
2
．
5.1.3建立S (H ' , α)关系式
由于弓形面积公式以及L (y )将会在以下多次用到，所以在此，我们单独将两式列出，后面就不再重复了．
设弓高为h ，则弓形的面积为：
h -r ⎫2⎛ S = π-arccos r +h -r (⎪
r ⎝⎭
如图2-2所示，有关系式
tan α=
L (y )OC -OM
=L (y )OC -y
⇒L (y )=tan α(OC -y )=(H ' cot α+3)tan α-y tan α=H ' +(3-y )tan α．
（1）当0≤y ≤y E 时，体积可看成一连串圆面该区间范围内的积分，此时记圆的面积：
22
S 1(y )=πr 12=π⎡R 2-(R -y )⎤=π⎡1.6252-(y -1.625)⎤=π(3.25y -y 2)． (21)
⎣⎦⎣⎦
（2）当y E ≤y ≤1时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h 2，
h 2=L (y )
-⎡1.5-
⎢⎣ =H '+(3-y )
tan α-⎡1.5 ． (22)
⎢⎣结合弓形面积公式得
⎛L y -1.5
S 2(y )
= π- ⎝
⎫
⎡1.6252-(y -1.625)2⎤
⎦⎣
+(
L (y )-1. -⎡⎣L y -1. ⎤⎦5
. 625
． (23)
（3）当1≤y ≤9时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h
3，
h 3=L (y )-⎡1.5
⎢⎣
=H ' +(3-y )tan α-⎡1.5，
⎢⎣得
L (y )
-1.5⎫2⎛S y =π-arccos L y -1.5 3() (24) ()⎪1.5+⎡⎣1.5⎝⎭
（4）当9≤y ≤y D 时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h
4，
h 4=L (y )-⎡1.5
=H ' +(3-y )tan α-⎡1.5
⎢⎢⎣⎣得
⎛L y -1.5S 4(
y )= π- ⎝
⎫
(25) +⎡⎣L (y )-1.5（5）当y D ≤y ≤10时，体积可看成一连串圆面在一定区间范围内的积分，此时记圆的面积：
2
S 5(y )=πr 52=π⎡R 2-(R -y )⎤=π1.6252-⎡1.625-10-y =π-67.5+16.75y -y ⎤()()． ⎣⎦⎣⎦
(26)
2
2
{}
最终联立(13)、(14)～(17)和(21)～(26)这10个公式即得(程序见附录四)
⎧y E S (y )dy +1S (y )dy +3+h tan αS (y )dy ,0≤h
3⎰y E 2⎰1⎪⎰01
⎪y E 19y D
⎪⎰0S 1(y )dy +⎰y S 2(y )dy +⎰1S 3(y )dy +⎰9S 4(y )dy ,6tan α≤h
V (H )=⎨y
19y D 10E
⎪⎰S 1(y )dy +⎰S 2(y )dy +⎰S 3(y )dy +⎰S 4(y )dy +⎰S 5(y )dy ,7tan α+1.5≤h
y E 19y D
⎪0
9y D 10⎪2
π1.52-3-h cot α+S y dy +S y dy +⎡⎤()()()⎣⎦⎰3-(3-h )cot α3⎰94⎰y D S 5(y )dy ,3-2tan α≤h
5.2 变位参数α、β值的确定
上面我们建立了罐内储油量与油位高度H 及变位参数（纵向倾斜角度α和横向偏转角度β ）之间的关系一般．现在我们需要用罐体变位后在进/出油过程中的实际检测数据确定模型中的变位参数α、β．
对于给定的油面高度H ，当α、β值不同时，理论计算出的罐内储油量不同，为与实际情况相吻合，采用如下算法来确定α、β的值：
第一步，取00
第二步，对100⨯100组α、β值，在理论模型下计算出每一组值在不同油面高度H 时罐内储油量V (H )．
第三步，计算每一组α、β值对应罐内储油量理论值与实际值的离差平方和，将对应100⨯100组离差平方和值比较取出平方差值最小时的α、β值，即:
(α, β)=min ∑(V (H )-V C )
i =1
n
2
．
另在问题A 附件2：实际采集数据表．xls 中，前一次显示油量容积值减去出油量
后与下一次显示油量容积值不相等，即显示油量容积值存在误差，故程序调用数据前，必须对出油表数据做修正，现用流水号201 205数据来说明数据修正的方法如：
积60311．43L ， 出油后理论剩油量60299．79L ，显示油量容积与出油后理论剩油量差为11．64L ，故流水号201出油后，显示油量容积修正为60448．88+11．64=60460．52L ，故任一流水号η的修正数据为：
V 修η=V 显（η+1）+V 出（η+1）．
据此得用excel [5]修正后的数据见问题A 附件2：实际采集数据表．xls 中K2---K603；
对算法用matlab 编程（程序见附录五）得α=00、β=00，此时为不发生变位的情况．故将算法修正为
第一步，取00
V C =V (H η)-V (H η+1)．
第三步，计算每一组α、β值对应罐内出油量理论值与实际值的离差平方和，将对应100⨯100组离差平方和值比较取出平方差值最小时的α、β值，即
(α, β)=min ∑(V (H )-V C )
i =1
n
2
．
同理，对算法用matlab 编程（程序见附录五）得α=20、β=4.90．
将α、β值代入模型，再用matlab 编程（附录）给出罐体变位后油位高度间隔为10cm 的罐容表标定值表如下：
5.3 正确性验证与方法可靠性检验
在对α、β值确定过程中，我们计算得到了α=00、β=00时储油量理论值数据（附录），将其与实际值做差(V (H )-V S )得差值拟合的百分差为0.23%，这就验证了模型的
正确性．据此我们提出一种正确性验证方法：
令α=00、β=00，代入储油体积函数V (H )中计算出理论的储油体积值，并与实际储油体积做离差平方和∑(V (H )-V C )，或对理论计算值与实际测量值用最小二乘法
i =1n
2
拟合，从而确定误差的大小，即模型的正确性．
六、误差分析
虽然我们建立了模型，得出了H 与V 的一一对应关系，但是模型的建立是在理想的假设基础之上的，实际上:油浮子是有一定的体积，而不是质点，如图3-1所示，其测量会有一定的误差；储油罐的厚度是存在的，所以罐内液体的长度一定比所给出的罐长要小．所以实际与理论之间必定存在一定的误差．
现在我们选择其中之一，即考虑油浮子对油位高度的测量的影响而带来的误差．根据实际情况，油位探针上的油浮子的机械外形各式各样，但我们可以将油浮子分为两类:1．储油罐无油时，油浮子与罐底接触相切（如图3-1右图）；2．储油罐无油时，油浮子与罐底有间隙，与罐壁相割（如图3-1左图）．
不考虑油的性质，无变位时，分别对上述两类进行分析，当为第一种类型时，油浮子可以检测到无油的状态，则此时的油浮子对油位高度的测量几乎无影响；当为第二种类型时，油浮子与储油罐总有一定间隙，则该间隙中的油位高度是无法测量到的，所以会对测量油位高度产生一定影响．
图3-1
七、模型的优缺点分析
本文建立的模型比较多，都是基于不同形状的储油的正截面面积在不同范围内的积分而建立起来的，有比较强的理论性及实用性，可以通过这些模型对储油罐的油位高度进行更为准确地测量，以便对储油罐内储油量进行估计，有利于油位计量管理系统的完善，其实际价值十分明显，并且对于储油罐的设计有一定指导意义．但是在建立模型时，我们忽略了一些客观因素，例如:温度、气压、罐体本身机械结构等对油浮子测量油位高度的影响，是在非常理想的状况下建立的模型，所以通过这些模型得到的理论值与实际测量值仍有一定差距．于是我们便将理论与实际的数据进行比较，求解出无变位/变位时的误差分析函数，对模型进行修正，尽量将理论值与实际测量值之间的差距减小到最小程度．
八、模型改进与推广
随着我国石油工业的发展需要，测量对于油库计量的重要性与日俱增，对测量方法和精度提出了更高的要求．在对问题的模型建立及求解过中，我们单纯地从有无变位而
造成的储油正截面不同形状的角度出发，运用微积分的思想，建立了罐内储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β ）之间一般关系的模型，是理想化的，然而结合实际情况，仍有一些影响测量油位高度的因素存在，所以我们权衡了各项影响因素后提出了:储油罐容积修正——液体静压力修正值的概念，分析油的静压力所引起的储油罐罐身产生形变而导致对计量精度的影响，并由此得出液体静压力修正值的计算公式为:∆V static pressure
πg ρR 3
=kH ，其中k =⨯10-6，这里ρ——储油罐内所盛油的
E σ
2
平均密度，g /cm 3；R ——储油罐的基本几何半径，cm ； E——储油罐罐体材料的弹性模量，N /cm 3；H ——储油罐内测量的油位高度，cm ；g ——液面在不同高度时的重力加速度，m /s 2；σ——储油罐罐体平均厚度，cm ．
在运用模型对油位高度进行标定时，除了结合误差函数外，如果再结合该液体静压力修正值∆V static pressure ，对模型进行改进，则可进一步提高测量精度，更接近真实值． 该模型研究储油罐的变位识别与罐容表标定，我们也可以将其运用于生产制造过程中原料添加量计量，潜水艇工作时蓄水量计量，洒水车盛水量计量等领域，具有较强的推广意义．
参考文献
[1]同济大学数学系．高等数学(第六版) ．北京:高等教育出版社，2007．
[2](美)Gerald Recktenwald ．数值方法和MATLAB 实现及应用．北京:机械工业出版社，2004．
[3]施浒立，赵彦．误差设计新理念与方法．北京:科学出版社，2007．
[4]黄波，王本善．中华人民共和国计量检定规程JJG266-1996．中国计量科学研究院．1996．
[5]闫志朝．Excel 函数、图表与数据分析．北京:机械工业出版社，2006． [6] 柴俊．数学实验教程(Matlab版) ．北京:科学出版社，2006． [7] 杨旭武．试验误差原理与数据处理．北京:科学出版社，2009． [8] 姜启源．数学模型(第二版) ．北京:高等教育出版社，1993．
附录
附录一 椭圆侧面面积求解 clear clc
%椭圆侧面面积求解 syms a b H y
f=2*(a/b)*sqrt(b^2-y^2)；%侧面面积积分函数 S=int(f，'y' ，'-b' ，'H-b')
附录二 变位时椭圆柱体体积求解与拟合与误差曲线拟合 clear clc
%变位时椭圆柱体体积求解与拟合与误差曲线拟合 %*************参数赋值************ a=0．89；%椭圆截面长半轴长 b=0．6；%椭圆截面短半轴长 l=2．45；%椭圆柱体长 n=0．4； m=2．05；
alpha=4．1*pi/180；%椭圆柱体纵向倾斜角度 k=H+n*tan(alpha)-b；%中间参数定义 j=H-m*tan(alpha)-b；
%*************参数赋值************ %*************模型体积求解************ if 0
V=-10^(3)*(a*b*cot(alpha)*(((1/(3*b^2))*sqrt((b^2-k．^2))+a*b*asin(k/b)+0．5*pi*(k+b+1))))； %有变位时椭圆柱体体积1 end
if m*tan(alpha)
V=-10^(3)*(cot(alpha)*a*b*(((1/(3*b^2))*sqrt((b^2-j．^2))+a*b*asin(j/b)+0．5*pi*(m-n)tan(alpha))．．．
-((1/(3*b^2))*sqrt((b^2-k．^2))-a*b*asin(k/b))))；%有变位时椭圆柱体体积2 end
if 2*b-ntan(alpha)
V=-10^(3)*(cot(alpha)*a*b*(((1/(3*b^2))*sqrt((b^2-j．^2))+a*b*asin(j/b)+0．5*pi*(m-n)tan(alpha))．．．
-((1/(3*b^2))*sqrt((b^2-(H-b)．^2))-a*b*asin((H-b)/b))))+pi*a*b*(n-(m-H)cot(alpha)))；%有变位时椭圆柱体体积3 end
%*************模型体积求解********************
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件1：实验采集数据表．xls' ，' 无变位进
油' ，'D2:D54')；
D=xlsread('F:\all\2010A\问题A 附件1：实验采集数据表．xls' ，' 无变位进油' ，'H2:H54')；
V2=polyfit(H，V-D ，2)%V-D表示体积的理论值与实际值的差 VV=V+V2；%对模型进行误差修订
%*************储油体积曲线绘制拟合************** plot(H，V ，'c*'，H ，VV ，'k-' ，H ，D ，'r--') title('有变位时储油体积拟合曲线') ； legend('理论计算储油体积曲线' ，' 修正后储油体积曲线' ，' 实际计算储油体积曲线') ； %*************储油体积曲线绘制拟合************** %*************罐容表标定值求解****************** H=0:10:3000； p=V
%*************罐容表标定值求解******************
附录三 模型检验拟合曲线绘制 clear clc
%**********模型检验拟合曲线绘制*****************
H=xlsread('C:\Documents and Settings\k01\桌面\问题A 附件2：实际采集数据表．xls' ，' 实际储油罐的采集数据' ，'E2:E603')；
D=xlsread('C:\Documents and Settings\k01\桌面\问题A 附件2：实验采集数据表．xls' ，' 实际储油罐的采集数据' ，'K2:K603')； M=polyfit(H，V ，2)
plot(H，V ，'b-' ，H ，D ，'ro') title('模型检验拟合曲线') ；
legend('修正后储油体积曲线' ，' 实际计算储油体积曲线') ； %**********模型检验拟合曲线绘制*****************
附录四 求解变位时储油罐的储油体积函数 clear clc
%求解变位时储油罐的储油体积函数 function V=volume(H)
syms alpha beta y(D) y(E) S1 S2 S3 S4 S5
%**********横向变位时，实际油位高度h 与测量油位高度H 的函数关系********** h=(H-1．5)*cos(beta)+1．5；
%**********横向变位时，实际油位高度h 与测量油位高度H 的函数关系********** %**********点D 、点E 在y 轴上的坐标值**********
y(D)=(3．25+2*tan(alpha)*(h+3*tan(alpha)-1．5) ．．．
-sqrt((3．25+2*tan(alpha)*(h+3*tan(alpha)-1．5))^2-4*sec(alpha)^2*(h+3*tan(
alpha)-1．5)^2))．．．
/(2*sec(alpha)^2)；%点D 在y 轴上的坐标y(D)
y(E)=(16．75+2*tan(alpha)*(h+3*tan(alpha)-1．5) ．．．
+sqrt((16．75+2*tan(alpha)*(h+3*tan(alpha)-1．5))^2．．．
-4*sec(alpha)^2*(h+3*tan(alpha)-1．5)^2)+8．375^2-16．25^2))．．． /(2*sec(alpha)^2)；%点E 在y 轴上的坐标y(E) %**********点D 、点E 在y 轴上的坐标值**********
%**********在不同空间范围内，储油截面图形的面积********** S1=pi*(3．25*y-y^2)；%当0
S3=(pi-acos((L(y)-1．5)/1．5))*1．5^2+(L(y)-1．5)*sqrt(3*L(y)-(L(y))^2)；%当1
S4=(pi-acos((L(y)-1．5)/sqrt(1．625^2-(y-8．375)^2)))．．．
+(L(y)-1．5)*sqrt(1．625^2-(y-8．375)^2-(L(y)-1．5)^2)；%当9
S5=pi*(3．25*y-y^2)；%当y(E)
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')．．．
+int(S3，'y' ，'3+h*tan(alpha)')%当0
while 6*tan(alpha)
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')+int(S3，'y' ，'3+h*tan(alpha)')．．．
+int(S4，'y' ，'y(E)'，'9')%当6*tan(alpha)
while 7*tan(alpha)+1．5
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')+int(S3，'y' ，'3+h*tan(alpha)')．．． +int(S4，'y' ，'y(E)'，'9')+int(S5，'y' ，'10' ，'y(E)')%当7*tan(alpha)
while 3-2*tan(alpha)
附录五 倾斜角αβ值定标 clear clc
%******************倾斜角αβ值定标*************************** for l=1:10000
for alpha=0:0．1:10
for beta=0:0．1:10，l=1:10000
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'E2:E604')；
D=xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'F2:F604')；
sum1=sum((V-D)．^2)；%求理论储油体积与实际储油体积的差方 Q(l)=sqrt(sum1)； a(l)=alpha； b(l)=beta； end end end
for l=1:10000 %比较理论储油体积与实际储油体积的和均方差大小 M=Q(l)；
if Q(l)>Q(l+1) M=Q(l+1)； alpha=a(l+1)； beta=b(l+1)； end end
alpha %输出均方差最小时α值 beta %输出均方差最小时β值
%******************倾斜角αβ值定标****************************** %******************算法修正后倾斜角αβ值定标******************** for l=1:10000
for alpha=0:0．1:10
for beta=0:0．1:10，l=1:10000
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'E2:E604')；
D=xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'D2:D604')； VC=V(k)-V(k+1)；
sum1=sum((VC-D)．^2)；%求理论储油体积与实际储油体积的差方 Q(l)=sqrt(sum1)； a(l)=alpha； b(l)=beta； end end end
for l=1:10000 %比较理论储油体积与实际储油体积的和均方差大小 M=Q(l)；
if Q(l)>Q(l+1) M=Q(l+1)； alpha=a(l+1)； beta=b(l+1)； end end
alpha %输出均方差最小时α值 beta %输出均方差最小时β值
%*****************算法修正后倾斜角αβ值定标********************
附录六 无变位进油表
28
29
附录七 倾斜变位进油表
30
附录七 问题二实际储油罐的采集数据表
31
32
33
2010年 全 国 大 学 生
数 学 建 模 竞 赛
储油罐的变位识别与罐容表标定
2010高教社杯全国大学生数学建模竞赛
承 诺 书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则．
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题．
我们知道，抄袭别人的成果是违反竞赛规则的， 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出．
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性．如有违反竞赛规则的行为，我们将受到严肃处理．
我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全
赛区评阅编号（由赛区组委会评阅前进行编号）：
2010高教社杯全国大学生数学建模竞赛
编 号 专 用 页
赛区评阅编号（由赛区组委会评阅前进行编号）：
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
储油罐的变位识别与罐容表标定
摘 要
加油站储油罐罐容表的精确度直接关系到加油站的经济利益，然而由于地基变形等原因，罐体的位置会发生纵向倾斜和横向偏转等变化，从而导致罐容表发生改变，影响其精度．
本文要解决的就是储油罐的变位识别与罐容表标定的问题．其中罐容表的标定，就是建立罐内油位高度与储油量的关系．对此，我们应用微积分及空间解析几何理论的相关知识，建立油罐体积函数模型V (H )．对于储油罐的变位识别问题，我们借助已建立的函数模型V (H )，用实际的油位高度确定理论储油量和变位参数值，并将理论储油量与实际测出的储油量采用最小二乘法进行拟合，然后通过拟合系数来判断模型的准确性．
对问题(1)，储油罐有无变位和纵向变位这两种情况，均要建立油罐体积积分函数模型，并运用matlab 软件求解模型，且将求解结果采用最小二乘法拟合，分析结果表明理论结果与数据模拟结果相吻合．最小二乘法拟合分析时也表明了模型求解中存在误差，从而以此为基础对模型进行修正，并得出罐容标定值表(见表一) ．
对问题(2)，同样建立建立油罐体积积分函数模型，采用离差平方和的算法并运用matlab 软件确定了变位参数α、β的值为α=20、β=4.90．以此为基础给出罐体变位后罐容表标定值表(见表三) ．
对问题一和二的模型做误差分析和修正后所得的结果显示，我们所建立的模型能很好的与实际情况相吻合，其吻合系数达到0.9996．最后我们还对模型进行了正确性验证与方法可靠性检验，并结合实际情况和应用价值对模型进行了改进与推广．
关键词 微积分；变位识别；小二乘拟合；误差；标定值
一、问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表(即罐内油位高度与储油量的对应关系) 进行实时计算，以得到罐内油位高度和储油量的变化情况．
许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化(以下称为变位) ，从而导致罐容表发生改变．按照有关规定，需要定期对罐容表进行重新标定．一种典型的储油罐其主体为圆柱体，两端为球冠体．现需要用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题．
（1）为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐(两端平头的椭圆柱体) ，分别对罐体无变位和倾斜角为α=4．10的纵向变位两种情况做了实验，实验数据如附录一所示．现需要建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm 的罐容表标定值．
（2）对于主体为圆柱体，两端为球冠体的储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数(纵向倾斜角度α和横向偏转角度β ) 之间的一般关系．利用罐体变位后在进/出油过程中的实际检测数据附录二，根据所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm 的罐容表标定值．进一步利用附录二中的实际检测数据来分析检验模型的正确性与方法的可靠性．
二、问题分析
2.1 问题一的分析
通常情况下，我们都可以通过油位计管理系统来标定罐容表，即通过测量进/出油量与罐内油位高度得到储油量的变化情况．但是储油罐在使用一段时间后，由于地基变形等原因，罐体的位置会发生纵向倾斜和横向偏转，导致罐容表测量值发生改变，这就需要我们定期对罐容表进行重新标定.
对罐容表进行重新标定前，首先我们必须判定油罐是否发生变位，即将测量的进出口油量实际值与罐体位置未发生变位时的理论值进行差分拟合，当读数误差达到一定值时，就可以判定罐体的位置是发生了变位．而对罐体无变位时罐容表的识别，可以通过对罐体已知的几何结构进行分析计算，确定罐体无变位时储油量V 与可测油位高度H 之间的函数关系V =f (H )．其次必须解决变位后罐容表如何重新标定的问题．
要解决上述问题，我们必须先建立V (H , α)的函数模型. 在建立V (H , α)的函数模型过程中, 我们参照了高等数学微积分的相关知识, 采用微元的思想得出模型.
2.2 问题二的分析
实际情况中，储油罐不单单只发生纵向倾斜，纵向和横向倾斜也应考虑，所以该问题中的情形比问题一更具有实际意义.
该问题是在问题一的基础之上增加了对横向倾角β的考虑，也就是要求我们同时考虑三个变量对储油量的影响，建立V (H , α, β).
[1]
根据事物的变化规律，针对该倾斜问题，我们发现:在两种倾斜同时发生时的结果与分步依次发生的结果是相同的，这就启发了我们可以通过分步考虑来简化模型的建立. 接着我们又考虑到，该问题中储油罐是圆柱体和球冠体这样两个特殊的对称体的组合体. 分析其几何特征可知:罐内液体不管怎么横向倾斜，其横截面均为垂直于水平面、左右对称的薄片, 也就是说横向变位对纵向变位储油量无影响.
所以为了易于模型的建立，我们假设油罐每次倾斜的完成顺序均如下图:
这样该问题中模型的建立又可以直接参照问题一中模型的建立，最后得V (H , α, β)函数模型.
模型建立的基本思路如下:
模型建立完后，得出储油量与油位高度及变位参数（纵向倾斜角度与横向倾斜角度）之间的关系模型. 我们就可以开始确定变位参数α、β的值．再将确定了变位参数α、β的值后代入模型来给出罐体变位后油位高度间隔为10cm 的罐容表标定值．最后可以再提出用实际检测数据对模型分析检验与建模方法的可靠性验证的方法．
三、模型假设及符号定义与说明
3.1模型假设
1．罐内储油不受温度压强等的影响，即储油量的体积大小只与油位高度有关；
2．油浮子为一质点，其大小可忽略不计；
3．储油罐壁的厚度很薄，可以忽略不计；
4．外界因素的改变不会影响储油罐的形状，即不会发生形变；
5．储油罐内部罐壁为理想、对称的几何图形，忽略其制造工艺带来的误差；
6．储油罐内部一系列小构建对储油量的影响忽略不计；
7．油的自身性质和蒸发损耗对储油量的影响忽略不计；
8．对油浮子与油接触时带来不可避免的仪器误差忽略不计．
3.2符号定义与说明
α:储油罐的纵向倾角；
β:储油罐的横向倾角；
l :储油罐的罐长；
V :储油量；
a :小椭圆油罐截面中椭圆的长半轴；
b :小椭圆油罐截面中椭圆的短半轴；
H :油浮子所在油面处的油位高度；
h :储油罐中任一油面位置处的油位高度；
V C (H ):测量出的储油量与油位高度关系函数；
V L (H ):实际的储油量与油位高度关系函数．
四、关于小椭圆型储油罐的模型建立与求解
为了层次清楚，我们先交待本节的结构．根据储油量的多少，以及油浮子位置的限制，对于近油位探针端下倾这种情形，分成如下几种情况进行考虑:
情形I 0≤H ≤m tan α，储油情况如图1-4所示，并建立模型Ⅱ；
情形II m tan α≤H ≤2b -n tan α，储油情况如图1-1所示，并建立模型Ⅰ； 情形III 2b -n tan α≤H ≤2b 储油情况如图1-5所示，并建立模型Ⅲ．
4.1 情形I 时，罐容与油面高度关系的模型建立
通常状况下，α在较小的情况下就会被工作人员所发现，并重新摆放储油罐，所以m tan α≤H ≤2b -n tan α时，比较常见，据此我们对此种情景在此做重点介绍．
如图1-1所示，取椭球圆柱体的中心轴为z 轴，并设该立体在过点z =0、z =l 且垂直于z 轴的两平面之间．以S 1(y ) 表示过z 且垂直与x 轴的截面面积．这时，取z 为积分变量，它的变化区间为[0, l ]；相应于[0, l ]上任一小区间[z , z +dz ]的一薄片的体积近似于底面积为S 1(y ) ，高为dz 的扁柱体的体积，即体积元素
dV 1=S 1(y ) dz ．
图1-1
以S 1(y ) dz 为被积分表达式，在闭区间[0, l ]上作定积分，便得所求立体的体积
V 1(y )=⎰S (y ) dz ． (1) 01l
接下来我们建立S 1(h ) 函数关系．
如图1-2所示，以图中椭圆柱体最左端椭球截面中心点为原点，以平行水平面和垂直水平面的方向分别为x 、y 轴建立直角坐标系，其中每片截面投影到xoy 坐标轴上的如图所示．
图1-2
图中椭圆面积公式为
x 2y 2． +
=1⇒x =22a b 现在，取纵坐标y 为积分变量，它的变化区间为[-b , y ]．相应于[-b , y ]上任一小区间[y , y +dy ]的窄条面积近似于高为dy
、底为积元素
⎛dS 1= dy ， ⎝从而得
S 1(
y )=⎰-b y ， 由图1-2有y =h -b ，代入得
S =⎰
接下来以h -b -b ． 运用MATLAB 程序，在区间[-b , h -b ]上作定为被积分表达式，积分，得所求的S 1(h ) 函数表达式为(程序见附录一)
S 1(
h )=⎰h -b
-b b h -b π=(
h -b ab arcsin +ab ． (2) a b 2
进而，我们通过建立H 与h 的函数关系，将H 引入到S 1(h ) 中，建立V 1(H )函数模
型．
选取yoz 坐标面上的截面如图1-3所示，油浮子所在处油位高度为H ，对应在y 轴上的投影点为C 点， 在z 轴上的投影长度为n ；油面上任一点的油位高度为h ，对应在y 轴上的投影为点B ，在z 轴上的投影长度为z ，即为该点的纵坐标大小；D 点为油面在y 轴上的交点．
图1-3
很明显的有线段长度关系
AB +BD =AC +CD ,
而
AB =h , BD =z tan α, AC =H , CD =n tan α,
所以
h +z tan α=H +n tan α⇒h =H +(n -z )tan α. （3） 对上述定积分公式(2)计算时，先不考虑积分限，直接对S 1(h )做不定积分，即 ⎰S (h )dh
1
=cot α
⎫h -b 1ab arcsin -πab (h -b )⎪+C ． b 2⎭
联立(1)，(2)和(3)式，因为所得结果比较复杂，为了简便起见，我们在此令k =H +n tan α-b ，j =H -m tan α-b ，所以积分后的式子为
V 1(H , α)=⎰H -m tan α-b H +n tan α-b S 1(h )dh
=ab cot α
arcsin j π⎫(m -n ) tan α⎪
b 2⎭
k ． b -ab cot α
arcsin
4.2 情形II 时，罐容与油面高度关系的模型建立
如图1-4所示，即为0≤H ≤m
tan α时的储油情形．
图1-4
此种情形下模型的建立与模型Ⅰ的建立基本相同，唯一不同的是z 轴方向上的积分上限:模型Ⅰ中的上限为罐长l ，而此处模型中油面边缘最右端与罐下壁有交界，投影到z 轴上的交点即为(0,0, H cot α+n )，所以该模型上限为H cot α+n ．即有
V 2(y )=⎰Hcot α+n 0S 2(y ) dz ． (4) 所求的S 2(h ) 函数表达式与模型Ⅰ中完全相同，即为
S 2(
h )
=⎰b h -b π=(h -b ab arcsin +ab ． (5) -b a b 2H 与h 的函数关系也为 h -b
h =H +(n -z )tan α ． (6)
同样联立(4)(5)(6)三式求解，为简便起见，我们也在此令k =H +n tan α-b ，同理得函数模型
V 2(H , α)=⎰0H +l tan α-b S 2(h )dh
=ab cot α
arcsin ⎫h -b 1-π(h -b )⎪b 2⎭0H +l tan α-b
=-ab cot α
arcsin k 1⎫-π(k +b +1)⎪． b 2⎭
4.3 情形Ⅲ时，罐容与油面高度关系的模型建立模型
如图1-5所示， 即为2b -n tan α
≤H ≤2b 时的储油情形．
图1-5
此种情形下模型的建立也与模型Ⅰ的建立基本相同，与模型Ⅰ相比，该模型相当于是模型Ⅰ中储油油体形状的立体图与一椭圆柱体的组合，所以该模型体积的求解分两部分完成，具体如下：
V 3=V 3-1+V 3-2．
其中V 3-1为椭圆柱体的体积，V 3-2为似模型Ⅰ的体积．
在求解V 3-1的函数关系式时，利用椭圆柱体的体积公式：体积=底⨯高．
求解时，为解释更加清楚，我们将yoz 坐标面上的图形截出平放如图1-6所示．
图1-6
观察图形可得如下线段关系式
AB =AC -BC , BC =EF ⨯cot α, EF =BF -BE ，
而
AC =n , BE =CG =H ．
所以最终可得椭圆柱体的高为n -(2b -H )cot α．
在求解底面时，我们直接取用椭圆的面积公式，得面积为πab ，所以有
V 3-1=πab (n -(2b -H )cot α) ． (7)
在求解V 3-2的函数关系式时，我们参照模型Ⅰ的求解过程，抓住其本质的不同之处，
仅将其积分下限换为n -(2b -H )cot α即得
V 3(y )=⎰l
n -(2b -H )S 3(y ) dz ．
所求的S 3(h ) 函数表达式也与模型Ⅰ中完全相同，即为
S 3(
h )=b =(
h -b ⎰-b a
H 与h 的函数关系也为 h -b h -b π (8) ab a r c +ab ．b 2
h =H +(n -z )tan α． (9)
同样联立(7)(8)(9)三式求解，为简便起见，我们也在此令k =H +n tan α-b ，同理得函数模型
V 3(H , α)=⎰H -m tan α-b H -b S 3(h )dh +πab (n -(2b -H )cot α)
=ab cot α
j 1⎫arcsin πm tan α⎪
b 2⎭
-ab cot α
arcsin H -b - b +πab (n -(2b -H )cot α)．
4.4情形Ⅳ时，罐容与油面高度关系的模型建立模型
如图1-7所示, 即为α=0时的情形.
图1-7
很明显，影响储油量V 的只有油面高度H ，所以我们直接建立V 与H 的函数表达式V (H ), 以下即为函数V (H )的建立：
此种情形下模型的建立与模型Ⅰ的建立基本相同, 则有
V 4(y )=⎰S 4(y ) dz ． (10) 0l
所求的S 4(h ) 函数表达式与模型Ⅰ中完全相同, 即为
b h -b π=(
h -b ab arcsin +ab ． （11) -b a b 2
由于α=0，所以液面各处H 与h 均相等，即有
h =H ． (12)
同样联立(10)(11)(12)三式求解得
bl H -b πV 4(H )=(
H -b abl arcsin +abl ． a b 2
4.5 一些补充说明
1、除了以上所建立的三种模型外，我们也考虑到了其他可能会有的情况，如图所示图1-8和图1-9，但是考虑到油浮子的测量局限性，这两种情况油浮子无法测量，所以我们在此也不做考虑．
S 2=⎰h -b
图1-8
图1-9
2、发生纵向倾斜时，可能为近油位探针端下倾，也可能为远油位探针端下倾，以上考虑的仅为近油位探针端下倾的情形．若出现远油位探针端下倾这一情形，对可测得油面高度H 的储油量函数可采用如下方法进行计算．
对于小椭圆柱体型储油罐这样的对称体，如图1-10所示，假设储油罐内有两个油浮子，分别位列储油罐内两对称的位置．并假设仅远油位探针端的油浮子可读，为H ，则另一油浮子，即近油位探针端油浮子油位高度为2a -H ，即对应的储油量函数直接套用以上模型即为：
V (H )=f (2a -H ) ．
但因远离油位探针端下倾时，微小的油位变化就会引起储油量发生很大的变化，实际工作中会很快被相关工作人员所发现，并重新放置．即远离油位探针端下倾没有太大的实际意义，所以我们在此不做进一步讨论．
图1-10
4.6 罐容与油面高度关系模型的求解
4.6.1、罐体变位后对罐容表的影响的求解
所谓罐体变位后对罐容表的影响，即考虑当的油位高度H 固定，为常数时，纵向倾角α对罐体储油体积V (H , α)的影响．可以通过理论值与实测值之间的差∆V (H , α)来判断
∆V (H , α)=V 2(H , α)-V C (H , α)．
V (H , α)是关于纵向倾角α的增函数，V (H , α)对上式拟合分析得，即α值增大时，
值增大．
4.6.2、给出α=4.10时油高间隔为1cm 所对应的一系列的罐容标定值
在使用已建立出的模型做标定之前，为确保结果的精确度，我们先采用差值拟合的方法对模型进行修正．即对油罐储油体积V (H )理论值与实际测量值的差做拟合曲线，也就是建立实测值与模型值的误差函数．
又因为采用数据拟合的方法可以反映函数曲线（面）V =f (H )反映对象整体的变化趋势，且使f (H )在某种准则下与所有实际测量值最为接近，即曲线拟合得最好．于是我们将实际测量值的数据点用matlab 拟合，在用三次拟合时，三次相前系数几乎为0，且做二次拟合时，相关系数r=0．9996，精度较高，说明拟合效果较好，故这里我们只采用中二次拟合．
现在计算误差函数．油罐储油体积无论是理论值还是实际测量值都与油位高度H 有关，所以误差函数也是H 的函数． 拟合时因进油表数值和出油表数值均为外部仪器测量，其数值较为精确，故采用进油表数值或采用出油表数值不会影响拟合效果．
用matlab [2]做二次曲线拟合[3](程序见附录二、三) 得出误差函数曲线方程为：
（1）未发生变位时，
实测值与模型值的误差函数为：
∆V (H )=V C (H )-V L (H )=1.6827⨯10-5H 2-0.15693H +17.9826．
故进行修正以后模型的函数为:
V (H )=V C (H )-∆V (H )．
用matlab 编程有修正前后拟合曲线如下图1-11所示
图1-11
（2）发生变位时
实测值与模型值的误差函数：
∆V (H )=V C (H )-V L (H )=0.00039739H 2-0.58342H +124.2537．
修正以后的模型的函数为:
V (H )=V C (H )-∆V (H )．
修正前后拟合曲线如下图1-12所示
:
V V 2现在采用相位分析法对修正后模型相似度进行检验，即用δ=cos θ=1计算得1V 2图1-12 δ=0.99，这也就证明了模型的准确性.
从上可知修正后函数模型与实际情况吻合系数较高，符合实际情况，现利用matlab 编程给出罐体变位后油位高度间隔为1cm 的罐容表标定值定标[4]表如下表一:
五、关于实际储油罐的模型建立与求解
5.1实际储油罐的模型建立
5.1.1 建立H ' (H , ) 关系式
如图2-1所示， 取出油浮子所在处的截面，并以其下端点为原心，该处切线方向为x 轴，垂直x 轴方向为y 轴，建立直角坐标系．其中平行与x 轴的直线为油面所在水平线， 所以B 点为油浮子在y 轴上的投影点，同时我们设其在y 轴正半轴上的投影高
度为H ' ．油浮子的测量高度仍然设为H ．
图2-1
由图可得如下线段关系
OB =OA +AB , AB =AC cos β, AC =CD -AD ，
而
CD =H , AD =r ，
最后可得
H ' (H , β) =(H -r )cos β+r ． (13)
5.1.2 建立V (H ' , S , α)关系式
接下来考虑纵向倾斜时，我们只需利用H ' ，结合微元积分的思想，建立函数关系式V (H ' , α, β)，最终通过(13)式将H 引入即得我们所需要的V (H , α, β)函数模型．
与问题一的思路相同，首先，我们根据储油量的多少，以及油浮子位置的限制，对于近油位探针端下倾这种情形，分成如下几种情况进行考虑:
0≤H '
我们以储油罐最下端切线方向为y 轴，以过储油罐最左端点且垂直于y 轴，并切于该点的指向上的直线为z 轴， 以y 轴与z 轴交点为原点，以过原点且垂直于y 轴和z 轴的直线为x 轴建立空间直角坐标系，如图2-2所示．同时设油面与储油罐的罐壁交点分别为点E 和点D ，观察图形得积分方向y 轴的五个区域:
0≤y ≤y E ，y E ≤y ≤1，1≤y ≤9，9≤y ≤y D ，y D ≤y ≤10．
图2-2
我们先考虑y E >1.5，即左端处一直存在一个由过点E 且垂直于y 轴的面所截出小球冠体，而对于0≤y E ≤1.5这种情况我们将在后面单独做以交代．
（1）当0≤H '
V =⎰S 1dy +⎰S 2dy +⎰0y E y E 13+H ' cot α1S 3dy ． (14)
（2）当6tan α≤H '
y E 19y D V =⎰S 1dy +⎰S 2dy +⎰S 3dy +⎰S 4dy ． (15) 0y E 19
（3）当7tan α+1.5≤H '
V =⎰S 1dy +⎰S 2dy +⎰S 3dy +⎰S 4dy +⎰S 5dy ． (16) 0y E 19y D y E 19y D 10
（4）当3-2tan α≤H '
现需要解决的就是求解D 、E 点的坐标, 来确定各模型表达式中积分上下限． 很明显点y D 和y E 是直线DE 与两个部分球面体的交点，根据初等数学的相关知识，
我们分别建立直线以及球面方程，最终联立求解得D 、E 坐标．
（一）直线方程的建立
如图2-2所示，有如下的线段关系:
OC =OA +AB +CB ，
而
OA =1, AB =2, CB =H ' cot α，
所以得
OC =3+H ' cot α．
即C 点的坐标为(0,3+H ' cot α,0)．
又已知∠ECA =α，可得其斜率为-tan α，再根据点斜式方法，可得该直线的方程为：
z =-tan αy +H ' +3tan α ． (18)
（二）球面方程的建立
如图2-3所示，设球体半径为R ．
图2-3
由图可得如下线段关系式
O ' N 2+NF 2=O ' F 2, O ' N =O ' P -NP ，
由于
O ' P =O ' F =R , NP =1，
得
(R -1)2+1.52=R 2⇒R =1.625．
则左半边球冠体球心O ' 的坐标为(0,1.625,1.5)，同理可求得右半边球冠体球心的坐标为(0,8.375,1.5)．
左半边球冠体所在球的方程：
(y -1.625)2+(z -1.5)=1.6252． (19) 2
右半边球冠体所在球的方程：
(y -8.375)
2
+(z -1.5)=1.6252． (20)
2
点E (y E , z E )既在左半边球冠体横截面圆的曲线上，同时也过直线CE ，为求E 点的坐标y E ，则可以联立方程(18)和(19)求解得：
y E =
3.25+2tan α(H ' +3tan α-1.5)
2sec α
2
同理可联立方程(18)和(20)求解得：
y D =
16.75+2tan α(H ' +3tan α-1.5)
2sec α
2
．
5.1.3建立S (H ' , α)关系式
由于弓形面积公式以及L (y )将会在以下多次用到，所以在此，我们单独将两式列出，后面就不再重复了．
设弓高为h ，则弓形的面积为：
h -r ⎫2⎛ S = π-arccos r +h -r (⎪
r ⎝⎭
如图2-2所示，有关系式
tan α=
L (y )OC -OM
=L (y )OC -y
⇒L (y )=tan α(OC -y )=(H ' cot α+3)tan α-y tan α=H ' +(3-y )tan α．
（1）当0≤y ≤y E 时，体积可看成一连串圆面该区间范围内的积分，此时记圆的面积：
22
S 1(y )=πr 12=π⎡R 2-(R -y )⎤=π⎡1.6252-(y -1.625)⎤=π(3.25y -y 2)． (21)
⎣⎦⎣⎦
（2）当y E ≤y ≤1时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h 2，
h 2=L (y )
-⎡1.5-
⎢⎣ =H '+(3-y )
tan α-⎡1.5 ． (22)
⎢⎣结合弓形面积公式得
⎛L y -1.5
S 2(y )
= π- ⎝
⎫
⎡1.6252-(y -1.625)2⎤
⎦⎣
+(
L (y )-1. -⎡⎣L y -1. ⎤⎦5
. 625
． (23)
（3）当1≤y ≤9时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h
3，
h 3=L (y )-⎡1.5
⎢⎣
=H ' +(3-y )tan α-⎡1.5，
⎢⎣得
L (y )
-1.5⎫2⎛S y =π-arccos L y -1.5 3() (24) ()⎪1.5+⎡⎣1.5⎝⎭
（4）当9≤y ≤y D 时，体积可看成一连串弓形面在该区间范围内的积分，设弓高为h
4，
h 4=L (y )-⎡1.5
=H ' +(3-y )tan α-⎡1.5
⎢⎢⎣⎣得
⎛L y -1.5S 4(
y )= π- ⎝
⎫
(25) +⎡⎣L (y )-1.5（5）当y D ≤y ≤10时，体积可看成一连串圆面在一定区间范围内的积分，此时记圆的面积：
2
S 5(y )=πr 52=π⎡R 2-(R -y )⎤=π1.6252-⎡1.625-10-y =π-67.5+16.75y -y ⎤()()． ⎣⎦⎣⎦
(26)
2
2
{}
最终联立(13)、(14)～(17)和(21)～(26)这10个公式即得(程序见附录四)
⎧y E S (y )dy +1S (y )dy +3+h tan αS (y )dy ,0≤h
3⎰y E 2⎰1⎪⎰01
⎪y E 19y D
⎪⎰0S 1(y )dy +⎰y S 2(y )dy +⎰1S 3(y )dy +⎰9S 4(y )dy ,6tan α≤h
V (H )=⎨y
19y D 10E
⎪⎰S 1(y )dy +⎰S 2(y )dy +⎰S 3(y )dy +⎰S 4(y )dy +⎰S 5(y )dy ,7tan α+1.5≤h
y E 19y D
⎪0
9y D 10⎪2
π1.52-3-h cot α+S y dy +S y dy +⎡⎤()()()⎣⎦⎰3-(3-h )cot α3⎰94⎰y D S 5(y )dy ,3-2tan α≤h
5.2 变位参数α、β值的确定
上面我们建立了罐内储油量与油位高度H 及变位参数（纵向倾斜角度α和横向偏转角度β ）之间的关系一般．现在我们需要用罐体变位后在进/出油过程中的实际检测数据确定模型中的变位参数α、β．
对于给定的油面高度H ，当α、β值不同时，理论计算出的罐内储油量不同，为与实际情况相吻合，采用如下算法来确定α、β的值：
第一步，取00
第二步，对100⨯100组α、β值，在理论模型下计算出每一组值在不同油面高度H 时罐内储油量V (H )．
第三步，计算每一组α、β值对应罐内储油量理论值与实际值的离差平方和，将对应100⨯100组离差平方和值比较取出平方差值最小时的α、β值，即:
(α, β)=min ∑(V (H )-V C )
i =1
n
2
．
另在问题A 附件2：实际采集数据表．xls 中，前一次显示油量容积值减去出油量
后与下一次显示油量容积值不相等，即显示油量容积值存在误差，故程序调用数据前，必须对出油表数据做修正，现用流水号201 205数据来说明数据修正的方法如：
积60311．43L ， 出油后理论剩油量60299．79L ，显示油量容积与出油后理论剩油量差为11．64L ，故流水号201出油后，显示油量容积修正为60448．88+11．64=60460．52L ，故任一流水号η的修正数据为：
V 修η=V 显（η+1）+V 出（η+1）．
据此得用excel [5]修正后的数据见问题A 附件2：实际采集数据表．xls 中K2---K603；
对算法用matlab 编程（程序见附录五）得α=00、β=00，此时为不发生变位的情况．故将算法修正为
第一步，取00
V C =V (H η)-V (H η+1)．
第三步，计算每一组α、β值对应罐内出油量理论值与实际值的离差平方和，将对应100⨯100组离差平方和值比较取出平方差值最小时的α、β值，即
(α, β)=min ∑(V (H )-V C )
i =1
n
2
．
同理，对算法用matlab 编程（程序见附录五）得α=20、β=4.90．
将α、β值代入模型，再用matlab 编程（附录）给出罐体变位后油位高度间隔为10cm 的罐容表标定值表如下：
5.3 正确性验证与方法可靠性检验
在对α、β值确定过程中，我们计算得到了α=00、β=00时储油量理论值数据（附录），将其与实际值做差(V (H )-V S )得差值拟合的百分差为0.23%，这就验证了模型的
正确性．据此我们提出一种正确性验证方法：
令α=00、β=00，代入储油体积函数V (H )中计算出理论的储油体积值，并与实际储油体积做离差平方和∑(V (H )-V C )，或对理论计算值与实际测量值用最小二乘法
i =1n
2
拟合，从而确定误差的大小，即模型的正确性．
六、误差分析
虽然我们建立了模型，得出了H 与V 的一一对应关系，但是模型的建立是在理想的假设基础之上的，实际上:油浮子是有一定的体积，而不是质点，如图3-1所示，其测量会有一定的误差；储油罐的厚度是存在的，所以罐内液体的长度一定比所给出的罐长要小．所以实际与理论之间必定存在一定的误差．
现在我们选择其中之一，即考虑油浮子对油位高度的测量的影响而带来的误差．根据实际情况，油位探针上的油浮子的机械外形各式各样，但我们可以将油浮子分为两类:1．储油罐无油时，油浮子与罐底接触相切（如图3-1右图）；2．储油罐无油时，油浮子与罐底有间隙，与罐壁相割（如图3-1左图）．
不考虑油的性质，无变位时，分别对上述两类进行分析，当为第一种类型时，油浮子可以检测到无油的状态，则此时的油浮子对油位高度的测量几乎无影响；当为第二种类型时，油浮子与储油罐总有一定间隙，则该间隙中的油位高度是无法测量到的，所以会对测量油位高度产生一定影响．
图3-1
七、模型的优缺点分析
本文建立的模型比较多，都是基于不同形状的储油的正截面面积在不同范围内的积分而建立起来的，有比较强的理论性及实用性，可以通过这些模型对储油罐的油位高度进行更为准确地测量，以便对储油罐内储油量进行估计，有利于油位计量管理系统的完善，其实际价值十分明显，并且对于储油罐的设计有一定指导意义．但是在建立模型时，我们忽略了一些客观因素，例如:温度、气压、罐体本身机械结构等对油浮子测量油位高度的影响，是在非常理想的状况下建立的模型，所以通过这些模型得到的理论值与实际测量值仍有一定差距．于是我们便将理论与实际的数据进行比较，求解出无变位/变位时的误差分析函数，对模型进行修正，尽量将理论值与实际测量值之间的差距减小到最小程度．
八、模型改进与推广
随着我国石油工业的发展需要，测量对于油库计量的重要性与日俱增，对测量方法和精度提出了更高的要求．在对问题的模型建立及求解过中，我们单纯地从有无变位而
造成的储油正截面不同形状的角度出发，运用微积分的思想，建立了罐内储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β ）之间一般关系的模型，是理想化的，然而结合实际情况，仍有一些影响测量油位高度的因素存在，所以我们权衡了各项影响因素后提出了:储油罐容积修正——液体静压力修正值的概念，分析油的静压力所引起的储油罐罐身产生形变而导致对计量精度的影响，并由此得出液体静压力修正值的计算公式为:∆V static pressure
πg ρR 3
=kH ，其中k =⨯10-6，这里ρ——储油罐内所盛油的
E σ
2
平均密度，g /cm 3；R ——储油罐的基本几何半径，cm ； E——储油罐罐体材料的弹性模量，N /cm 3；H ——储油罐内测量的油位高度，cm ；g ——液面在不同高度时的重力加速度，m /s 2；σ——储油罐罐体平均厚度，cm ．
在运用模型对油位高度进行标定时，除了结合误差函数外，如果再结合该液体静压力修正值∆V static pressure ，对模型进行改进，则可进一步提高测量精度，更接近真实值． 该模型研究储油罐的变位识别与罐容表标定，我们也可以将其运用于生产制造过程中原料添加量计量，潜水艇工作时蓄水量计量，洒水车盛水量计量等领域，具有较强的推广意义．
参考文献
[1]同济大学数学系．高等数学(第六版) ．北京:高等教育出版社，2007．
[2](美)Gerald Recktenwald ．数值方法和MATLAB 实现及应用．北京:机械工业出版社，2004．
[3]施浒立，赵彦．误差设计新理念与方法．北京:科学出版社，2007．
[4]黄波，王本善．中华人民共和国计量检定规程JJG266-1996．中国计量科学研究院．1996．
[5]闫志朝．Excel 函数、图表与数据分析．北京:机械工业出版社，2006． [6] 柴俊．数学实验教程(Matlab版) ．北京:科学出版社，2006． [7] 杨旭武．试验误差原理与数据处理．北京:科学出版社，2009． [8] 姜启源．数学模型(第二版) ．北京:高等教育出版社，1993．
附录
附录一 椭圆侧面面积求解 clear clc
%椭圆侧面面积求解 syms a b H y
f=2*(a/b)*sqrt(b^2-y^2)；%侧面面积积分函数 S=int(f，'y' ，'-b' ，'H-b')
附录二 变位时椭圆柱体体积求解与拟合与误差曲线拟合 clear clc
%变位时椭圆柱体体积求解与拟合与误差曲线拟合 %*************参数赋值************ a=0．89；%椭圆截面长半轴长 b=0．6；%椭圆截面短半轴长 l=2．45；%椭圆柱体长 n=0．4； m=2．05；
alpha=4．1*pi/180；%椭圆柱体纵向倾斜角度 k=H+n*tan(alpha)-b；%中间参数定义 j=H-m*tan(alpha)-b；
%*************参数赋值************ %*************模型体积求解************ if 0
V=-10^(3)*(a*b*cot(alpha)*(((1/(3*b^2))*sqrt((b^2-k．^2))+a*b*asin(k/b)+0．5*pi*(k+b+1))))； %有变位时椭圆柱体体积1 end
if m*tan(alpha)
V=-10^(3)*(cot(alpha)*a*b*(((1/(3*b^2))*sqrt((b^2-j．^2))+a*b*asin(j/b)+0．5*pi*(m-n)tan(alpha))．．．
-((1/(3*b^2))*sqrt((b^2-k．^2))-a*b*asin(k/b))))；%有变位时椭圆柱体体积2 end
if 2*b-ntan(alpha)
V=-10^(3)*(cot(alpha)*a*b*(((1/(3*b^2))*sqrt((b^2-j．^2))+a*b*asin(j/b)+0．5*pi*(m-n)tan(alpha))．．．
-((1/(3*b^2))*sqrt((b^2-(H-b)．^2))-a*b*asin((H-b)/b))))+pi*a*b*(n-(m-H)cot(alpha)))；%有变位时椭圆柱体体积3 end
%*************模型体积求解********************
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件1：实验采集数据表．xls' ，' 无变位进
油' ，'D2:D54')；
D=xlsread('F:\all\2010A\问题A 附件1：实验采集数据表．xls' ，' 无变位进油' ，'H2:H54')；
V2=polyfit(H，V-D ，2)%V-D表示体积的理论值与实际值的差 VV=V+V2；%对模型进行误差修订
%*************储油体积曲线绘制拟合************** plot(H，V ，'c*'，H ，VV ，'k-' ，H ，D ，'r--') title('有变位时储油体积拟合曲线') ； legend('理论计算储油体积曲线' ，' 修正后储油体积曲线' ，' 实际计算储油体积曲线') ； %*************储油体积曲线绘制拟合************** %*************罐容表标定值求解****************** H=0:10:3000； p=V
%*************罐容表标定值求解******************
附录三 模型检验拟合曲线绘制 clear clc
%**********模型检验拟合曲线绘制*****************
H=xlsread('C:\Documents and Settings\k01\桌面\问题A 附件2：实际采集数据表．xls' ，' 实际储油罐的采集数据' ，'E2:E603')；
D=xlsread('C:\Documents and Settings\k01\桌面\问题A 附件2：实验采集数据表．xls' ，' 实际储油罐的采集数据' ，'K2:K603')； M=polyfit(H，V ，2)
plot(H，V ，'b-' ，H ，D ，'ro') title('模型检验拟合曲线') ；
legend('修正后储油体积曲线' ，' 实际计算储油体积曲线') ； %**********模型检验拟合曲线绘制*****************
附录四 求解变位时储油罐的储油体积函数 clear clc
%求解变位时储油罐的储油体积函数 function V=volume(H)
syms alpha beta y(D) y(E) S1 S2 S3 S4 S5
%**********横向变位时，实际油位高度h 与测量油位高度H 的函数关系********** h=(H-1．5)*cos(beta)+1．5；
%**********横向变位时，实际油位高度h 与测量油位高度H 的函数关系********** %**********点D 、点E 在y 轴上的坐标值**********
y(D)=(3．25+2*tan(alpha)*(h+3*tan(alpha)-1．5) ．．．
-sqrt((3．25+2*tan(alpha)*(h+3*tan(alpha)-1．5))^2-4*sec(alpha)^2*(h+3*tan(
alpha)-1．5)^2))．．．
/(2*sec(alpha)^2)；%点D 在y 轴上的坐标y(D)
y(E)=(16．75+2*tan(alpha)*(h+3*tan(alpha)-1．5) ．．．
+sqrt((16．75+2*tan(alpha)*(h+3*tan(alpha)-1．5))^2．．．
-4*sec(alpha)^2*(h+3*tan(alpha)-1．5)^2)+8．375^2-16．25^2))．．． /(2*sec(alpha)^2)；%点E 在y 轴上的坐标y(E) %**********点D 、点E 在y 轴上的坐标值**********
%**********在不同空间范围内，储油截面图形的面积********** S1=pi*(3．25*y-y^2)；%当0
S3=(pi-acos((L(y)-1．5)/1．5))*1．5^2+(L(y)-1．5)*sqrt(3*L(y)-(L(y))^2)；%当1
S4=(pi-acos((L(y)-1．5)/sqrt(1．625^2-(y-8．375)^2)))．．．
+(L(y)-1．5)*sqrt(1．625^2-(y-8．375)^2-(L(y)-1．5)^2)；%当9
S5=pi*(3．25*y-y^2)；%当y(E)
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')．．．
+int(S3，'y' ，'3+h*tan(alpha)')%当0
while 6*tan(alpha)
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')+int(S3，'y' ，'3+h*tan(alpha)')．．．
+int(S4，'y' ，'y(E)'，'9')%当6*tan(alpha)
while 7*tan(alpha)+1．5
V=int(S1，'y' ，'y(D)'，'0')+int(S2，'y' ，'1' ，'y(D)')+int(S3，'y' ，'3+h*tan(alpha)')．．． +int(S4，'y' ，'y(E)'，'9')+int(S5，'y' ，'10' ，'y(E)')%当7*tan(alpha)
while 3-2*tan(alpha)
附录五 倾斜角αβ值定标 clear clc
%******************倾斜角αβ值定标*************************** for l=1:10000
for alpha=0:0．1:10
for beta=0:0．1:10，l=1:10000
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'E2:E604')；
D=xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'F2:F604')；
sum1=sum((V-D)．^2)；%求理论储油体积与实际储油体积的差方 Q(l)=sqrt(sum1)； a(l)=alpha； b(l)=beta； end end end
for l=1:10000 %比较理论储油体积与实际储油体积的和均方差大小 M=Q(l)；
if Q(l)>Q(l+1) M=Q(l+1)； alpha=a(l+1)； beta=b(l+1)； end end
alpha %输出均方差最小时α值 beta %输出均方差最小时β值
%******************倾斜角αβ值定标****************************** %******************算法修正后倾斜角αβ值定标******************** for l=1:10000
for alpha=0:0．1:10
for beta=0:0．1:10，l=1:10000
H=10^(-3)*xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'E2:E604')；
D=xlsread('F:\all\2010A\问题A 附件2：实际储油罐的采集数据．xls' ，' 实际储油罐的采集数据' ，'D2:D604')； VC=V(k)-V(k+1)；
sum1=sum((VC-D)．^2)；%求理论储油体积与实际储油体积的差方 Q(l)=sqrt(sum1)； a(l)=alpha； b(l)=beta； end end end
for l=1:10000 %比较理论储油体积与实际储油体积的和均方差大小 M=Q(l)；
if Q(l)>Q(l+1) M=Q(l+1)； alpha=a(l+1)； beta=b(l+1)； end end
alpha %输出均方差最小时α值 beta %输出均方差最小时β值
%*****************算法修正后倾斜角αβ值定标********************
附录六 无变位进油表
28
29
附录七 倾斜变位进油表
30
附录七 问题二实际储油罐的采集数据表
31
32
33