2010储油罐的变位识别与储容罐的定位标定全国一等奖

储油罐的变位识别与储容罐的定位标定
摘要
为了解决储油罐纵向倾斜和横向偏转对罐容表读数的影响，本文分别介绍了不同变位条件下罐内储油量与油位高度及纵向倾斜和横向偏移参数αβ、之间的模型，然后根据所得理论值与实验所给数据做比较，进行误差分析；根据误 差分析结果对罐容表进行校核，从而达到准确计量目的。

针对问题一，对于纵向变位后的储油罐，我们根据实际情况，针对不同油位高度采用微元法及几何法建立了罐内储油量与纵向倾斜角α及h 之间的模型，通过将理论值与附表中实验值作容积曲线及残差图来做比较、分析误差，根据误差分析结果，代入待标定的高度序列，给出了罐容表重新标定模型
2387582.965124.06i i i i i i V V e e h h =-=-+-，其中
为了进一步证明了理论模型的准确性，最后利用Matlab 通过空间坐标变换
的方法模拟储油罐中油面的上升过程，得到了容积曲线从而证明了理论模型的准确性，并对变位参数进行了灵敏度分析。

针对问题二，我们利用微元法和几何法在第一问所建模型的基础上，将油罐分为圆柱体和两个球冠端面，把圆柱体分成五段积分，球冠端面分成两个部分，结合纵向倾斜和横向偏移参数αβ、，得到储油量与αβ、及h 的一般关系。

为了利用所建模型确定变位参数，我们采用差分最小二乘法及网格搜索算
法，得到了 2.0072, 1.903185a b ==o o ，使
21
11
([(,,)(,,)](,))n
t
i t i r i i i V h
V h V h h αβαβ++=--∆∑最小，
并以该变位参数计算了以10cm 为间隔的罐容表（见后文），最后根据附表2中实际出油值、显示高度与计算出的理论值进行比较，对模型的正确性与方法的可靠性进行了检验。

关键词：分段切割 微元法积分 网格搜索 计算机模拟 差分最小二乘拟合
一、问题重述
储油罐使用一段时间后，不可避免的会发生纵向倾斜和横向偏转等问题，从而导致罐容表发生改变，按照有关规定，需重新对罐容表进行标定，本论文要求我们解决以下问题：
（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为α=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。

请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。

（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β）之间的一般关系。

请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为
10cm的罐容表标定值。

进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。

二、模型的假设
（1）假设伸入油罐内管子的体积及储油罐的管壁所占的体积是忽略不计的；（2）假设油浮子不会因为长期浸泡变形而对读数产生影响；
（3）仅考虑系统误差，排除随机误差的影响；
，不能过大
（4）假设纵向倾斜和横向偏转的参数αβ
三、符号及变量说明
V：储油罐内油的体积(单位L)
h：通过油浮子测得的油面高度(即油位高度，单位dm)
α：纵向倾斜时的倾角(单位度)
β：横向偏转时的偏角(单位度)
r：与冠状体相对应的球的半径(1.625m)
r：在冠状体上截面圆的半径
h
下文中涉及到的其它符号均在下文中标注
四、问题的分析
针对问题一，对于无变位的储油罐，我们直接利用微元法得到储油量与油位 高度的函数关系，根据所得结果与附表1中实验值进行比较，分析出了实验结果与理论结果存在误差，并给出了误差是随油位高度呈线性增长。

附表1所给数据只有油量与油位高度，显然寻求出不同条件下油量与油位高度的函数关系是我们必须首先考虑的问题，结合题图4，即实验用的小椭圆油罐，采用微元积分方法很容易导出油量与油位高度的关系式，然后，运用得到函数解析式()V f h =计算出不同高度的储油量，运用Matlab 软件在同一图内分别作出理论计算值与实际观测值的图像，通过将理论值与附表中实验值作容积曲线及残差图来做比较、分析误差,即可分别求得无变位与变位状态下实际值与理论值的偏差，根据计算得到罐容表校正后的标定值，此外用计算机模拟的方法验证了模型的正确性，并对参数进行了灵敏度分析。

针对问题二，为了找出罐内储油量与油位高度及变位参数α、β之间的一般关系，按照储油罐纵向倾斜变位及横向偏转两种方式，分别建立了罐内储油量与油位高度和变位参数之间的函数关系式。

对储油罐纵向倾斜，我们采用微元法得到了1V 与h 、α 之间的关系式：1(,)V f h α=；对储油罐横向偏转，我们用几何法得到h 和β之间的关系式为：()h f β=，从而就建立了总体积2V 与h 、α、β之间的关系：2(,,)V f h αβ=；采用非线性差分最小二乘拟合法来确定α、β的值； 再利用实际检测数据来检验模型的正确性和可靠性，将模型计算的结果同检测数据进行误差分析。

五、模型的建立与求解
5.1 问题一解答过程：
为寻求油量与油位高度的关系，我们通过微元法逐一求解不同条件下油量与油位高度的函数关系式，如下为函数关系式的建立过程： 5.1.1当截面为题图4椭圆面 (1)油罐未发生变位时：
图1 椭圆截面图
用平行于水平面的平面截得图形为以x 、L 为边的矩形，以椭圆中心为原点，建立如图1所示平面直角坐标系，沿y
轴方向取一微元dy ，得：
()2dV S y dy
L x dy =⋅=⋅⋅⋅
由椭圆方程2222x y a b +=1，得：x = ()2
21212h b V a a h L L h b b arcsin b b b b π--=⎡⎤⎛⎫⋅⋅=⋅-⋅+⋅-+ ⎪⎢⎥
⎝⎭⎣⎦
⎰其中：a ----椭圆长半轴的长度，值为0.89m ；
b ----椭圆短半轴的长度，值为0.6m ；
L----油罐椭圆柱体部分长度，值为2.45m
(2)油罐纵向倾斜时：
图2 区域划分图
由图2可知：根据不同的油量高度，我们将油罐划分为六个区域分别加以讨论，显然区域⑥不可能出现(否则油会溢满)，故不需考虑。

下面计算在其它区域的情况，以区域③用切片法计算为例：
图3
如上图，沿着x轴正方向投影，根据投影位置的不同，又可将柱体分为两部分来计算
图4 图5
由图4：当投影轴在中心轴的下半侧时，
PQ AB x tan
α==⋅，'O P b x tan α=-⋅
将()'y O P b x tan α=-=--⋅带入椭圆标准方程22
221x y a b
+=，得：
0)
x a x =>'·
·''O MQN O P b x tan ab arccos ab arccos O Q b
S α
-⋅==扇形()'1'2O MN O P MN a b x tan S α∆=⋅⋅=⋅-⋅⋅则图4中阴影部分面积：
()1''O MN
O MQN S
S S b x tan a b arccos a b x tan b αα∆=--⋅=⋅-⋅-⋅⋅阴扇形由图5：当投影在中心轴下半侧时，
PR x tan α=⋅
，2PQ b x tan α=-⋅，'O P b PQ x tan b α=-=⋅-
'''O NQM O P O Q
x tan b
ab arccos
ab arccos
b
S α⋅-∴=⋅=⋅扇形
()'1'2MO N
O P MN a x tan b S α∆=⋅⋅=⋅⋅-()x tan b ab ab arccos
a x tan
b b S απα⋅-=-⋅+⋅⋅-阴2
综上：可求出此时油量与油位的关系为：
1
tan tan tan tan b h
L h
b L V dx dx S S ααα
α
+-=
+
⎰
⎰
阴1阴2
同理可求出在其他区域的函数解析式：
① 当1h L tan α≤<⋅0时，
01tan 0
h
L V dx S α+=
⎰
阴
②当10L tan h b L tan αα⋅≤<-⋅时，
01tan 1tan h
L h
L V dx S αα
+-=
⎰
阴
② 当01b L tan h b L tan αα-⋅≤<+⋅时，
01tan 2tan tan 1
b
h
L tan a h b L a a
V dx dx S S α+-=
+
⎰
⎰
阴阴
④当102b L tan h b L tan αα+⋅≤<-⋅时，
012tan tan h
L a h
L a
S V
dx +-=⎰
阴
⑤当022b L tan h b α-⋅≤<时，
2tan 2tan 1
02b
h
L V ab dx
S tan b h L ααπα-⎛⎫
+
⎪⎝
⎭-=⋅-⎰
阴
其中：0L —油浮子与左侧面的距离，值为0.4m ；
1L —油浮子与右侧面的距离，值为2.05m ；
由上面式子，当标定高度知道后，就可根据()V f h =唯一确定此时对应的油
量，根据计算值就可与实际标定值作曲线拟合，得到理论值与实际值的残差及相应的误差，由误差大小一方面可以检验模型建立的合理性，如果误差超过允许的范围，说明模型不合理，需要改进；同时又可以通过尽量减小误差来进行变位校核。

5.1.2 误差分析
分析附表1，无变位进油实验结束时间和出油实验开始时间非常接近，因此可以看作是连续实验，可将进油的累加油量终值看作出出油的初始值，对这两个过程分别作高度-油量曲线，发现两曲线基本重合，可知出油过程相当于进油的逆过程，因此误差分析只考虑了进油情况。

下图为无变位时理论值与实际值比较图象：
如图可见、在无罐体变位时，通过积分得到的容积与实际储油量的误差值和油位高度h呈线性关系，，我们将每次实际注入油量与容积增值做比较得到下表：(前28次)
由表可见，每次进油，实际占用罐体体积都是稍大于50L的,并且对于一定注入油量(50L)偏大值稳定于1.75L，通过分析及查找相关文献资料，我们认为产生这种误差的原因与多种因素有关，并且与储油罐内部结构有关，如测量油位装置和进出油管、罐内固体残渣、油污、焊接凸结构所占用的体积都会使得单位油占体积偏大。

在原因不明的情况下，直接用曲线或多项式弥补误差欠妥，这里因对罐内结构缺乏了解，因此不做进一步分析。

罐体倾斜时，使用上述积分模型，对不同高度h分段处理积分，就可以得到由模型计算出的理论值，和实际值作对比曲线如下：
从上图线可见理论值略高于真实值，列出其中一段对象(如流水号233-240)
由表可知：理论值总是比实际值大，显然，以原来方式标定油量是不可取的，所以罐体刻度必须重新标记，而且误差的大小在允许的范围内，也验证了模型建立的合理性。

另外地，进一步分析误差项的特性，用Matlab画出误差项的残差图, 不同于无变位下误差与浮标高度所呈的线性关系，这里误差值是以曲线的形式出现,因此利用Matlab对误差项进行多项式拟合，得到的曲线图如下：
可见二次或三次拟合曲线都能较好的拟合误差值的趋势，因此在制作罐容表时，可以将用二次拟合得到的值作为误差的补偿项，加入原式,这里取二次拟合曲线得到参数如下:
2387582.965124.06i i i e h h =-+-
其中e 为误差项，改正后的标定值应该减去对应的i e ,i h 为浮标高度。

5.1.3 变位后容量表的标定
由于产生误差的原因多样，为了使油量测量不至于偏差太大，此处用二次拟合模拟误差项，以油量与液面高度的函数关系计算出理论值减去模拟误差项作为标定值：
i i i V V e =-
如下为选取部分刻度从计量基准点算起，以1cm 为间隔对应的容量标定值,
5.1.4 计算机模拟
为了检验积分方法计算出的油量是否准确，函数关系是否合理，我们运用计算机模拟实验来验证其结果的正确性，计算机模拟过程如下： (1)当油罐无变位时
如图，先建立三维坐标系XYZ
通过计算机模拟的方法，对于每一个高度h ,生成一个平行于面xOy 的面
ABCD ，同时随机在面ABCD 产生点集V ，满足如下条件：
1）点集均匀分布在面ABCD 上。

2）面ABCD 的边界和面积是已知的。

3）面ABCD 足够大，能完整的截取罐体。

记空间中在罐体内的点的点集为E ，则由概率论知识，可以得到截面''''A B C D 的面积如下：
''''{|,}(1){|}A B C D ABCD
num x x V x E S num x x V S ∈∈=∈
其中{|}num x x V ∈表示点集V 中点的数量，计算油面''''A B C D 大小，把每一小段油位增值h ∆与面''''A B C D 的面积的乘积近似的看作储油量的体积元:
''''A B C D V h S ∆=∆⋅
通过累计起某一高度h 以下所有体积元V ∆（将h ∆取为以适当小值），近似地得到高度h 与罐容积V 关系：0
h
V V =∆∑。

(2)油罐发生变位时
当储油罐发生纵向或横向变位，为了保证1.1中油量增加的方向为重力的反
方向，即方法的通用性，我们对罐体作了纵向和横向变位的空间三维坐标变换，使得坐标系中z 轴方向始终指向油上升方向，即与油面垂直，方法如下：
A
三维坐标变换中的系数
T A A A Z Y X ),,(为某点在空间直角坐标系A 的坐标； T B B B Z Y X ),,(为该点在空间直角坐标系B 的坐标；
由于坐标系变换后坐标尺度不变，而且变换前后坐标原点重合，有：
()B A B A B A X X Y R Y Z Z ω⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦
根据图中旋转角度与油罐变位方向可得0x y z w w w α
β===
且油罐变位角度一般均为小角度，可近似的得到坐标旋转矩阵()R w ：
1
0()()()()0
11Z Y X R R R R βωωωωαβ
α
-⎡⎤
⎢⎥=⋅⋅=⎢⎥⎢⎥-⎣⎦
自此就可以得到两个坐标系间转化关系：B A X R X =⨯。

(3)、模拟过程
我们通过MATLAB 编程模拟油面高度h 逐渐上升的过程，程序流程如下：
1） 初始化当前油面某点分别在原坐标系和旋转坐标系下位置
1_a x i s r e a l 和1_axis tran ，在旋转系生成均匀分布且足够大的点集
V ，并计算面积，记为1S 。

2）
1__[0,0,];i i axis tran axis tran h -=+∆并计算对应的原坐标系位置_()_i i axis real R w axis tran =⨯,检验在原坐标系下，_i
axis real 是否超出罐体边界，若是，则到步骤4），否则到步骤3);
3）
1;i i =+并且通过逻辑运算统计点集{|,}num x x V x E ∈∈数目，并带入计算式(1)，计算出当前体积元i V ∆，返回步骤1）。

4）
通过累加运算，求不同油面下罐容积0
h
V V =
∆∑。

(4)、模拟结果和灵敏度分析
利用Matlab 画图，将积分模型结果，计算机模拟数值，真实实验值，画在一个坐标系上，可以得到：
可见计算机模拟值曲线处于真实值和积分模型值曲线之间，说明通过计算机模拟得到的解是可信的。

借助计算机模拟，我们就可以画出不同a值下的模拟值曲线，从而进行灵敏度分析：
可见a值越大，高度-容积曲线是越平矮的，同时总容积是减少的，因为油面会更早接触的浮标口而溢出。

5.2 问题二求解：
问题二求解同样需要先计算出油量与液面高度的函数关系，求圆柱体部分的容积与问题一中椭圆截面相似，只需将a 、b 都替换为R 即可，于是得到圆形截面柱体部分体积计算如下： 5.2.1 柱体部分体积计算
令 ()
12
R xtan S R arccos R R x tan R αα-=⋅--⋅
(
)222
x tan R S R R arccos R x tan R R απα⋅-=-⋅+⋅⋅-其中：1S 表示投影区域在中性轴下方时，弓形阴影的面积，
2S 表示投影区域超过中性轴时，弓形阴影部分面积
（1）当20h tan l α≤<⋅时，
1tan 10
h
l V S dx α+=⎰
（2）当21tan h R tan l l αα⋅≤<-⋅时，
12tan 1tan h l h
l V S dx αα
+-=⎰
（3）当12R tan h R tan l l αα-⋅≤<+⋅时，
1122tan tan tan tan R h l h
R l V S dx S dx ααα
α
+-=+⎰
⎰
（4）当212R tan h R tan l l αα+⋅≤<-⋅时，
1
22
tan tan h
l h
l V S dx αα
+-=⎰
（5）当122R tan h R l α-⋅≤<时，
222tan 1an 2
t R
h l R h V R l S dx tana α
απ-⎛⎫2-=⋅-+ ⎪⎝
⎭⎰
其中：R =1.5m （圆柱截面的半径）；
12m l =（油浮子离圆柱较近端点的距离）
； 2l =6m （油浮子离圆柱较远端点的距离）
5.2.2计算冠状部分的容积
（1）当20h tan l α≤<⋅时，0V =右
图1
如上图，左侧部分体积可化分成上、下两个部分来计算，由于液面在中心轴上方和下方时，下面部分体积计算完全一样，于是我们不妨只考虑液面在上面情况，如图所示：
D
F
G
E r-H
1
h
r r
H
图2
图3
由图2可知： 12h h l tan α=+⋅，1OD h R =- 在ODE ∆中，由勾股定理得：2
2
2
h OD r r +=
h r =∴由图3得：
2
()DFG
h h
r H
r arccos r H S S S r ∆-=-=⋅--⋅
阴扇2
01()h h h r H
arccos
r H
dz V r r ⎛
-∴=⋅--⋅ ⎝
⎰左下
计算上面部分体积时，由于α非常小，可以认为11O N O M =，那么图中阴影部分面积可近似为扇形面积，21O N S α=⋅⋅1
2
扇，只需求出ON 即可求出上面部分体积，如图1：1()O N r H =--，
在1
OOA ∆中，1O A =11
O A
O A
V dx S -=∴⋅⎰左上扇
则 V V V =+左左上左下
112
02()O A
O A
h h
h r H
arccos
r H dz dx S r r -⎛
-=⋅--⋅ ⎝
⋅+⎰⎰
扇
（2）当212tan h R l tan l αα⋅≤<-⋅时，此时左右两侧冠状体都有油量，需计算两侧冠状的体积，左侧情况与上面完全相同，这里不在赘述，计算右侧体积时与左侧类似，同样划分为上下两部分，运用切片法求得：
(
)0
2
2h h h r H arccos r H dz r r V ⎛
-=⋅-- ⎝⎰1右 11
2O A
O A V dx S -=⋅⎰右扇
其中 22h h l tan α=-⋅，222
2()h r R h r +-=，211
2
O N S α=
⋅⋅扇
1O A =
1()O N r H =--
则
V V V =-右1右2右
(
)110
2
2h O A O A h h r H arccos r H dz r r dx S -⎛
-=⋅--⋅- ⎝⋅⎰
⎰扇 （3）当222R l tan h R α-⋅≤<时，左侧冠状体完全充满液体，
()
2π33
H r H V ⋅⋅-∴=左
右侧冠状体积算法与上面计算相同。

5.2.3 横向偏转时测量高度与实际高度关系
图4 图5
（1）当液面在中心线下方时，如图4有：
OC R AC R h =-=-,()OD OC cos R h cos ββ∴=⋅=-⋅
0()h R h cos R β∴+-⋅=,即0()h R R h cos β∴=--⋅
（2）当液面在中心线上方时，如图5有：
OD AD R h R =-=-()OC OD cos h R cos ββ∴=⋅=-⋅
0()h R cos R h β∴-⋅+=，即：0()h R R h cos β∴=--⋅
由此可见；在两种不同情况下，0h 与h 的关系表达式是一致的，都为：
0()h R R h cos β=--⋅，根据测量值h ，就可求出液面对应竖直方向液面高度
0h ，从而求出相应的油容积。

5.2.4 确定变位参数,αβ
分析原附表2中所给进、出油量数据，通过将显示储油量作差值得到“显示的出油量”，我们发现，显示的出油量与实际的出油量数据不相符合，给出其中前10组数据如下：
可见当前显示储油量不准确，不能真实的反应油量的存贮和变化，而且通过
上表可见实际出油值是总是大于显示出油量，因此该误差绝不是由系统误差造成的，应该是储油罐发生纵向或横向变位导致罐容表发生变化，应该计算新的罐容表。

根据前面已求出的变位储油罐容积公式，我们采取最小二乘拟合的方式确定变位角a 和b 。

但当前罐容表显示储油量信息不准确，因而它不能作为确定a 和
b 的拟合对象，即我们无法通过对“浮标高度－显示存储油量”关系获得正确的
a 和
b ，所以我们转而拟合“浮标高度－实际出油量”的关系如下：
2111
,min {[(,,)(,,)](,)}n
t i t i r i i i V h V h V h h αβ
αβαβ++=--∆∑
其中(,,)t i V h αβ表示理论变位油罐模型在参数,a b 和高度i h 下的值；
1(,)r i i V h h +∆表示从显示高度1i h +降到高度i h 时，实际的出油量，通过读附表2
获得。

这里我们采用网格搜索的方式寻找变位参数α和β，即通过不断缩小搜索范
围和步长逼近最优解，步骤如下：
1）记ft 为搜索步长值，初始值为1，限制,(0,8)αβ∈且1;i = 2）在当前步长下搜索目标函数最小值，记最优值为[]i
i αβ，当最优值变化量
小于一定值eps ，转到步骤3）,否则1;i i =+并减小步长值,此处2
ft
ft =。

重新限制搜索范围至[,]i i αβ周围，转向2）； 3）输出当前[,]i i αβ，结束迭代
将上述获得理论变位油罐模型带入，通过Matlab 编程求解上述拟合模型，我们得到的a
和b 为： 2.0072 1.903185βα
=。

，=
即该地下储油罐纵向倾斜2.0072度，横向偏转1.903185度。

给出罐体变位后油位高度间隔为10cm 的罐容表标定值如下：
5.2.5 模型可靠性检验
为了验证以上模型以及罐容表的正确、可靠性，将求得的,a b 带入，求出理论出油量1(,,)(,,)t i t i V h V h αβαβ+-序列，并与实际储油量1(,)r i i V h h +∆序列比较分析得到：
平均误差很小，而最大误差却有4.63%，逐个对应项分析发现，这种情况属于极少数个例，大部分项误差非常小，在603个流水号中，误差大于1%的有80个，误差大于2%的只有20个，误差大于3%一共只有9个（包含下表最大项）
由上表可见，通过我们的模型计算出的显示油量下降量和实际出油量，每次几乎都是相等的，误差量相比原表大大减少，累加起来，得到的罐容表也一定会是准确可信的。

六、模型的评价与改进
模型优点：
（1）本论文在正确理解了题意的基础上，建立了科学、合理的数学模型，为油罐的正确标定提供了理论依据。

（2）模型的理论性很强，大部分公式都是严格运用高等数学中积分法积出。

模型缺点：
（1）模型算出的是理论值，考虑到实际问题，可能存在系统误差、偶然误差等，我们没有考虑，如果在此基础上进行校核，就会使得计算结果更符合实际情况。

（2）问题二中，考虑到倾斜角 很小这个条件，我们在计算体积时采用分段方
式计算，其中用切片法求体积时将切片（不规则）部分面积近似用扇形代替，考虑到这样的等效对总体积影响很小，认为是合理的。

参考文献
[1]管冀年,赵海，卧式储油罐罐内油品体积标定的实用方法，计量与测试技术，第31卷第3期:21-22页，2004。

[2]田铁军，倾斜立式罐部分容积的计算，现代计量测试，第7卷第4期：39-44页，1999
[3]姜启源，谢金星，叶俊，数学模型（第三版），北京：高等教育出版社，2003
[4]韩中庚，数学建模方法及其应用，北京：高等教育出版社，2005
[5]何国伟，误差分析方法（第一版），北京：国防工业出版社，1978。