CN102967260A - 一种粉体物料堆的体积密度测量方法及其系统 - Google Patents

Abstract

本发明公开了一种粉体物料堆的体积密度测量方法及其系统，体积测量方法包括：将N个激光测量头并排固定于粉料车上方的架子上，当粉料车驶入架子正下方时，控制各激光测量头对粉体物料堆表面同时按照预设的相同频率进行扫描，获得各扫描点的高度数据，当粉料车的车箱驶出架子正下方时结束扫描；构建三维迪卡尔坐标系，计算各扫描点的三维坐标；利用计算机程序对扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分，得到粉料车中粉体物料堆的体积。本发明能准确真实地获取煤炭、粉料等粉体物料堆的现场体积数据，进而获取密度数据，能客观地监督运输方使用合法、安全、正规的手段进行煤炭、粉料等粉体物料的运输操作。

Description

一种粉体物料堆的体积密度测量方法及其系统

技术领域

本发明涉及图形图像技术领域，尤其涉及一种粉体物料堆的体积密度测量方法及其系统。

背景技术

为了保证粉体物料（如煤炭）运输工作的安全有序运行，防止粉体物料在运输途中发生截留、更换、以次充好等作弊现象，煤炭运输监控成为煤炭运输中一个重要的环节。对煤炭运输过程进行监控的核心是测量煤炭运输车中煤炭的体积、质量、密度等信息，然后根据这些信息判断煤炭在出发地和目的地之间是否按照正规标准进行运输。在矿山、电厂、港口及车站等有许多粉体物料堆或粉料车(如煤等)的体积和质量需要准确而快速地测量,即盘煤,多年来这一直成为技术难题。对煤炭运输过程进行监控的核心是测量煤炭运输车中煤炭的体积、质量、密度等信息。对于运输途中发生的截留问题，可以通过在出发地和目的地的地设置地磅进行很方便的质量测量，但对于更换、以次充好等作弊现象，则需要对粉体物料堆的体积和密度进行测量。

传统的盘煤方法有如下两种：一、人工盘煤。用推土机对煤堆整形成标准几何体后,人工用皮尺或手持激光测距仪进行丈量和估算，计算出煤堆的体积。但是该方法需要耗费大量的人力，测量精度低，而且只能对场地中的煤炭进行测量，无法对煤炭运输车中的煤炭进行测量。二、激光盘煤。目前大多数激光测距设备，其功能只局限于距离测量，基本不带或者只带简易的测量仪器设定程序，无法进行二次开发来满足用户需求，而根据物体进行全面的长度、分散点密度等测量的仪器，大多也不能进行计算统计，而如果进行全面的数据分析，则需要大量的人力物力与综合实施手段才能达到满意效果，且仪器相对非常庞大，无法满足煤炭行业现场需求。在其他大型项目中，会带有测量数据分析项目模块，且其项目庞大，模块众多，其使用的算法或是三重积分，或是三维重构，没有将两者相结合进行计算。

发明内容

为了准确真实地获取煤炭、粉料等粉体物料堆的现场体积密度数据，本发明提出一种粉体物料堆的体积密度测量方法及其系统。

本发明提出了一种粉体物料堆的体积测量方法，包括：

S1、撘建激光测量装置，将N个激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，将两个最远的激光测量头之间的距离调整到等于粉料车的宽度，N为大于等于2的正整数；

S2、当粉料车的车箱驶入所述架子正下方时，控制各所述激光测量头对所述粉体物料堆表面同时按照预设的相同频率进行扫描，获得各扫描点的高度数据，当所述粉料车的车箱驶出所述架子正下方时结束扫描；

S3、利用计算机程序以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴构建三维迪卡尔坐标系，计算各扫描点在所述三维迪卡尔坐标系下的三维坐标；

S4、利用计算机程序分别对每行和每列的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

进一步地，步骤S3中所述计算各扫描点在所述三维迪卡尔坐标系下的三维坐标具体包括：通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

进一步地，步骤S4中所述样条插值方法为二次样条插值。

进一步地，步骤S1中所述激光测量头个数为12。

本发明还提出了一种粉体物料堆的密度测量方法，包括，通过地磅测量出的所述粉体物料堆的质量数据，通过如上所述的粉体物料堆的体积测量方法测量所述粉体物料堆的体积数据，通过所述质量数据除以所述体积数据获取所述粉体物料堆的密度数据。

基于同一发明构思，本发明还提出了一种粉体物料堆的体积测量系统，，包括：

激光测量装置，所述装置包括N个激光测量头，所述激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，确保两个最远的激光测量头之间的距离调整到等于粉料车的宽度，其中所述N为大于等于2的正整数；

激光控制装置，用于当所述粉料车的车箱驶入所述架子正下方时，控制所述激光测量装置中各所述激光测量头对所述粉体物料堆表面同时按照预设的相同频率进行扫描，用于获得各扫描点的高度数据，当所述粉料车的车箱驶出所述架子正下方时控制所述激光测量装置结束扫描，将所述各扫描点的高度数据发送给坐标记算软件模块；

坐标记算软件模块，用于接收所述激光控制装置发送的各扫描点的高度数据，利用计算机程序分别计算各扫描点在以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴的三维迪卡尔坐标系下的三维坐标，将所有扫描点的三维坐标发送给体积计算软件模块；

体积计算软件模块，用于接收所述坐标记算软件模块发送的所有扫描点的三维坐标，利用计算机程序分别对每行和每列的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

进一步地，所述坐标记算软件模块中所述利用计算机程序分别计算各扫描点在以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴的三维迪卡尔坐标系下的三维坐标具体包括：通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

进一步地，体积计算软件模块中所述样条插值方法为二次样条插值。

进一步地，所述激光测量装置中所述激光测量头个数为12。

基于同一构思，本发明还提出了一种粉体物料堆的密度测量系统，包括：

质量测定模块，用于通过地磅测量出的所述粉体物料堆的质量数据；

体积测定模块，用于通过如权利要求6-9之一所述的粉体物料堆的体积测量系统测量所述粉体物料堆的体积数据；

密度计算模块，用于通过所述质量数据除以所述体积数据获取所述粉体物料堆的密度数据。

使用本发明所述的粉体物料堆的体积密度测量方法及其系统，能方便快捷地计算出粉体物料堆的体积和密度，特别是通过在物料运输的出发地和目的地快速测量粉体物料运输车所运输的料体物料的体积和密度，防止粉体物料在运输途中发生截留、更换、以次充好等作弊现象，保证粉体物料运输工作的安全有序运行。

附图说明

图1是本发明具体实施例一所述的粉体物料堆的体积测量方法流程图；

图2是本发明三维曲面生成示意图；

图3是本发明具体实施例二所述的粉体物料堆的体积测量系统框图；

图4到图9是本发明不同积分方式介绍的示例图。

具体实施方式

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。

实施例一

图1是本实施例所述的粉体物料堆的体积测量方法流程图，如图1所示，本实施例所述的粉体物料堆的体积测量方法包括：

S101、撘建激光测量装置；

将N个激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，将两个最远的激光测量头之间的距离调整到等于粉料车的宽度。

其中，所述激光头的数目N为大于等于2的正整数，激光头的数量可根据煤炭运输车宽度、测量环境和不同用户测量精度的不同而进行相应的调整。

以粉体物料为煤炭为例，本实施例通过将12个激光测量头设置于粉料车上方的架子上，该架子与水平面平行，各激光测量头在架子上竖直安装，确保发射的激光束垂直于水平面，图2是本发明三维曲面生成示意图，如图2所示，架子上第1个激光点与第12个激光点之间的距离与车箱宽度相等，支架与车箱宽边平行，每个激光点投射的光线与车箱底面垂直。

S102、当粉料车的车箱驶入架子正下方时，控制各激光测量头对粉体物料堆表面同时按照预设的相同频率进行扫描，获得各扫描点的高度数据，车箱驶出架子正下方时结束扫描；

例如，当有待检测的煤炭运输车进入该煤炭运输监控系统的工作区域时，12个激光头同时对煤堆自上而下进行照射扫描，并将激光在煤堆表面扫描点反馈的光信号通过模数转换器将光信号转换为数字信号，对所述数字信号进行编码和处理以获得粉体物料堆表面各扫描点的高度数据。当煤炭运输车的车尾经过激光头时，结束扫描。

S103、构建三维迪卡尔坐标系，计算各扫描点在的三维迪卡尔坐标；

利用计算机程序以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴构建三维迪卡尔坐标系，计算各扫描点在所述三维迪卡尔坐标系下的三维坐标；

通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

例如，如图2所示，各激光测量头第1次扫描的扫描点的X坐标均为0；各激光测量头第2次扫描的扫描点的X坐标均为第2次扫描距第1次扫描时车箱移动的距离；各激光测量头第3次扫描的扫描点的X坐标均为第3次扫描距第1次扫描时车箱移动的距离；依次类推，各激光测量头第最后一次扫描的扫描点的X坐标均为该次扫描距第1次扫描时车箱移动的距离，也为车箱的长度值。

第1个激光测量头所测得的所有扫描点的Y坐标均为0，第2个激光测量头所测得的所有扫描点的Y坐标均为第2个激光测量头到第1个激光测量头的距离；第3个激光测量头所测得的所有扫描点的Y坐标均为第3个激光测量头到第1个激光测量头的距离；依次类推，第12个激光测量头所测得的所有扫描点的Y坐标均第12个激光测量头到第1个激光测量头的距离，也为车箱的宽度值。

各扫描点的Z值均由步骤S102所获得各扫描点的高度数据确定，需要说明的是，测量所得的直接高度数据是指激光测量头照射到的点到该激光测量头之间的高度，该扫描点在上述三维迪卡尔坐标系下的Z坐标需要将激光测量头到车箱底部的距离减去上述直接高度数据获得。

S104、分别对对各行和各列扫描点进行样条插值形成光滑曲线，对光滑曲线组成的曲面进行二重积分。

利用计算机程序分别对每个激光测量头扫描到的所有扫描点进行样条插值形成光滑曲线，分别对每个扫描时刻所有激光测量头扫描到的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

其中，样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差，这样就避免了使用高阶多项式所出现的龙格现象，样条插值具体方法有很多，例如，包括线性样条插值、二次样条插值、三次样条插值。

本实施例以二次样条插值为例，使用Matlab函数库来实现,主要用到的Matlab中的griddata函数、mesh函数、linspace函数。例如，对10个扫描点以二次样条插值为例，所述10个扫描点坐标是

A=x1,y1,z1;x2,y2,z2;........;x10,y10,z10];

x=A(:,1);

y=A(:,2);

z=A(:,3);

N=100；//插值点数，实际使用时需确定，暂写100

x0=linspace(min(x),max(x),N);

y0=linspace(min(y),max(y),N);

[X,Y,Z]=griddata(x,y,z,x0’,y0,'v4');//插值拟合曲面。

figure;

mesh(X,Y,Z);//体积的近似值用下代码

L=sum(sum(abs(Z)));

S=abs((y0(2)-y0(1))*(x0(2)-x0(1)));

V=L*S;

所述V的值就是体积值。

依据该算法计算出粉料车内粉体物料堆的体积，进一步地，通过地磅测量出的所述粉体物料堆的质量数据，使用密度计算公式：密度=质量/体积，可求出密度。

其中，多重积分是定积分的一类，它将定积分扩展到多元函数（多变量的函数），例如求f(x，y)或者f(x，y，z)类型的多元函数的积分。

下面对各种不同的多重积分进行介绍和分析：

譬如，边长为4×6×5的长方体的体积可以通过两种方法得到：

通过函数f(x,y)=5在xy平面中的区域D，也就是长方体的底上的双重积分

∫∫_D5dx dy

或者是常函数1在长方体上的三重积分

∫∫∫1dx dy dz

令n为大于1的自然数。考虑所谓的半开n维矩形（下面简称矩形）。对于平面，n=2，而多重积分就是双重积分。

将每个区间[ai,bi)分成有限个不重叠的子区间，每个都是左闭右开。将子区间记为Ii。则，所有所有如下形式的子矩形的族

C＝I₁×I₂×···×I_n

是T的一个划分，也即，子矩形C是互不重叠的，而且它们并集为T。C中的子矩形的直径按照定义是C中最大的边长，而T的划分的直径就是划分中的子矩形的最大直径。

令f:T→R为定义在T上的函数。考虑如上定义的T的划分

T＝C₁∪C₂∪···∪C_m

其中m是正整数。如下形式的和称为黎曼和

$\underset{k=1}{\overset{m}{\mathrm{\Σ}}}f\left(P_k\right)m\left(C_k\right)$

其中，对于每个k，点Pk在Ck中，而m(Ck)是笛卡尔积为Ck的区间的边长之积。

函数f称为黎曼可积，如果如下极限存在

$S=\underset{\mathrm{\δ}\→0}{\lim }\underset{k=1}{\overset{m}{\mathrm{\Σ}}}f\left(P_k\right)m\left(C_k\right)$

其中极限取遍所有直径最大为δ的T的划分。若f黎曼可积，S称为f在T上的黎曼积分。记为

∫_Tf(x)dx

定义在任意有界n维集合上的函数的黎曼和可以通过将函数延拓到一个半开半闭矩形上来求出，其取值在原来的定义域之外为0。然后，原来的函数的积分就定义为延展的函数在矩形区域中的积分（如果存在的话）。

下文中n维黎曼积分简称多重积分。

多重积分具有很多与单变量函数的积分一样的性质（线性，可加性，单调性，等等）。而且，和单变量情况一样，可以用多重积分找出函数在给定集合上的积分。具体来讲，给定集合

和D上的可积函数f，f在定义域上的平均值为

$\stackrel{\‾}{f}=\frac{1}{m\left(D\right)}{\∫}_{D}f\left(x\right)\mathrm{dx},$

其中m（D）是D的测度。

是，积分

e＝∫∫_Tf(x，y)dx dy

是f在T上的双重积分，而若

积分

l＝∫∫∫_Tf(x，y，z)dx dy dz

是f在T上的三重积分。

注意，按常规，双重积分用两个积分号，而三重积分有三个；这只是记法上方便，也是为了通过重复积分来计算多重积分

多重积分问题的解决在多数情况下依赖于将多重积分转化为一系列单变量积分，而其中每个单变量积分都是直接可解的。包括如下几种：

1、直接检验

有时可以直接获得积分的结果，而无需任何直接计算。

2、常数

在常函数的情况中，结果很直接：只要将常函数c乘以测度就可以了。如果c=1，而且是在R2的子集中积分，则乘积就是区域面积，而在R3中，它就是区域的体积。

例如：

在D上积分f：

${\∫}_3^6{\∫}_2^{4}2\mathrm{dxdy}=\mathrm{area}\left(D\right)\·2=(2\·3)\·2=12$

3、利用可能的对称性

如果定义域存在沿着某条轴的对称性而且函数对于那个变量是奇函数，则积分为0（因为相反的两部分加起来为0）。

对于Rn中的函数，只要相关变量对于形成对称的轴是奇变量就可以了。

例一:

给定f(x,y)=2sinx-3y3+5以及T=x2+y2≤1为积分区域（半径为1的圆盘，包含边界）。

利用线性性质，积分可以分解为三部分：

∫∫_T(2sinx-3y³+5)dx dy＝∫∫_T2sin x dx dy-∫∫_T3y³dx dy+∫∫_T5dx dy

sin x和3y3都是奇函数，而且显然T对于x和y轴都是对称的；因此唯一有贡献的部分是常函数5因为其它两个都贡献0。

例二：

考虑函数f(x,y,z)=x exp(y2+z2)以及圆心在原点的半径为2的球T=x2+y2+z2≤4。该球显然是对于三条轴都对称，但是只要对于x轴积分就可以看出结果是0，因为f对于该变量是奇函数。

4、简化公式

简化公式基于简单积分区域来将多重积分转化为单变量积分的序列。它们必须从右至左计算，过程中将其它变量暂时视为常数（和偏导数的计算类似）。

5、R2中的常规区域

此种方法适用于满足下述条件的任何定义域D:

D投影到x轴或y轴任一轴，形成一个有边界的范围,以a,b代表边界值。

通过a,b两点并与

垂直的直线与D相交后的两个端点，可以用2个函数α,β定义。

x轴：

将D对x轴做垂直投影，函数是连续函数，并且D可以视为（定义在[a,b]区间上的）α(x)和β(x)之间的区域。则

${\∫\∫}_{D}f(x,y)\mathrm{dxdy}={\∫}_a^b\mathrm{dx}{\∫}_{\mathrm{\α}\left(x\right)}^{\mathrm{\β}\left(x\right)}f(x,y)\mathrm{dy}$

y轴：

将D对y轴做垂直投影，函数

是连续函数，并且D可以视为（定义在[a,b]区间上的）α(y)和β(y)之间的区域。则

${\∫\∫}_{D}f(x,y)\mathrm{dxdy}={\∫}_a^b\mathrm{dy}{\∫}_{\mathrm{\α}\left(y\right)}^{\mathrm{\β}\left(y\right)}f(x,y)\mathrm{dx}$

范例

例：可以采用简化公式的D区域。

考虑区域：D＝{(x，y)：x≥0，y≤1，y≥x²}（参看附图4）。计算

∫∫_D(x+y)dx dy

该区域可以沿x或者y轴分解。要采用公式，必须先找到限制D的两个函数和定义区间。

这个例子中，这两个函数为：

α(x)＝x²和β(x)＝1

而区间为[a，b]＝[0，1]（这里为了直观起见采用沿x轴分解）。

应用简化公式，得到：

${\∫\∫}_D(x+y)\mathrm{dxdy}={\∫}_0^1\mathrm{dx}{\∫}_{x^2}^1(x+y)\mathrm{dy}={\∫}_0^1\mathrm{dx}{[\mathrm{xy}+\frac{y^2}{2}]}_{x^2}^1$

（首先，第二个积分将x作为常数）。然后就是用积分的基本技术：

${\∫}_0^1{[\mathrm{xy}+\frac{y^2}{2}]}_{x^2}^1\mathrm{dx}={\∫}_0^1(x+\frac{1}{2}-x^3-\frac{x^4}{2})\mathrm{dx}=\·\·\·=\frac{13}{20}$

如果沿着y轴分解，可以计算

${\∫}_0^1\mathrm{dy}{\∫}_0^{\sqrt{y}}(x+y)\mathrm{dx}$

并得到同样的结果。参见图5。

R3中可(沿xy平面）分解区域的例子。

R3中的分解

这些公式可以推广到三重积分：

T是一个可以投影到xy平面的体积，它夹在α(x,y)和β(x,y)两个函数之间。那么：

${\∫\∫\∫}_{T}f(x,y,z)\mathrm{dxdydz}={\∫\∫}_D\mathrm{dxdy}{\∫}_{\mathrm{\α}(x,y)}^{\mathrm{\β}(x,y)}f(x,y,z)\mathrm{dz}$

（此定义和其它R3中的分解类似）。

变量替换

积分的极限常常不易交换（区域无法分解或者公式很复杂），这时可以采用变量替换来重写积分，令区域更加简易，从而可以用更简单的公式表达。为此，函数必须变换到新坐标系下。

例（1-a）:

函数为 $f(x,y)={(x-1)}^2+\sqrt{y};$

若采用替换x′＝x-1，y′＝y则x＝x′+1，y＝y′

可以得到新函数 $f_2(x,y)={\left(x^{\′}\right)}^2+\sqrt{y}.$

对于定义域要进行类似处理，因为原来是采用变换前的变量表达的（本例中的x和y）。

微分dx和dy要通过包含被替换的变量对于新变量的偏微分的雅可比行列式来变换。（譬如，极坐标的微分变换）。

常用的变量替换有三种（R2中一种，R3中两种）；但是，更普遍的变换可以用同样的原理来发现。

极坐标：

从笛卡尔坐标系变换到极坐标系，参见图6。

在R2中，若定义域有某种圆形对称性而函数也有某种特征，则可以采用极坐标变换（参看图中的例子），也就是说将点P(x,y)从笛卡尔坐标变换到相应的极坐标中。这使得定义域的形状改变，从而简化运算。

该变换的基本关系如下：

f(x，y)→f(ρcosφ，ρsinφ)。

例（2-a）:

函数为f(x，y)＝x+y

应用该变换得到

f(ρ，φ)＝ρcosφ+ρsinφ＝ρ(cosφ+sinφ)。

例（2-b）:

函数为f(x，y)＝x²+y²

这里有：

f(ρ，φ)＝ρ²(cos²φ+sin²φ)＝ρ²

这里使用了勾股定理（在简化操作时很有用）。

定义域的变换是根据x和y通过环厚和角度的幅度来限定ρ,φ的区间。

从笛卡尔到极坐标的区域变换，参见图7。

例（2-c）:

区域为D＝x²+y²≤4，圆周半径2；很明显，这个区域所覆盖的角度是整个圆周角，所以φ从0变化到2π，而环半径从0变化到2（内环为0的环形就是圆）。

例（2-d）:

区域为D＝{x²+y²≤9，x²+y²≥4，y≥0}，这是在正y半平面中的圆环（参看示意图）；注意φ表示平面角而ρ从2变化到3。因此变换出来的区域为矩形：

T＝{2≤ρ≤3，0≤φ≤π}.

该变换的雅可比行列式为：

$\frac{\∂(x,y)}{\∂(\mathrm{\ρ},\mathrm{\φ})}=\left|\begin{array}{cc}\cos \mathrm{\φ}& -\mathrm{\ρ}\sin \mathrm{\φ}\\ \sin \mathrm{\φ}& \mathrm{\ρ}\cos \mathrm{\φ}\end{array}\right|=\mathrm{\ρ}$

这可以通过将x=ρcos(φ),y=ρsin(φ)代入关于ρ的第一行和关于φ的第二行的偏微分中得到，所以微分dx dy变换为ρdρdφ.

一旦函数和区域的变换完成后，可以定义极坐标中的变量变换公式：

∫∫_Df(x，y)dx dy＝∫∫_Tf(ρcosφ，ρsinφ)ρdρdφ。

注意φ在[0,2π]区间中有效，而ρ测量长度，因此只能取正值。

例(2-e):

函数为f(x，y)＝x区域和例2-d相同。

从前面对D的分析，我们知道ρ的区间为[2,3]，而φ的为[0,π].函数变换为：

f(x，y)＝x→f(ρ，φ)＝ρcosφ。

最后，应用积分公式：

∫∫Dx dx dy＝∫∫Tρcosφρdρdφ。

一旦区间给定，就可以得到

${\∫}_0^{\mathrm{\π}}{\∫}_2^3{\mathrm{\ρ}}^2\cos \mathrm{\φd\ρd\φ}={\∫}_0^{\mathrm{\π}}\cos \mathrm{\φd\φ}{\left[\frac{{\mathrm{\ρ}}^3}{3}\right]}_2^3={\left[\sin \mathrm{\φ}\right]}_0^{\mathrm{\π}}(9-\frac{8}{3})=0.$

柱极坐标，参见图8：

R3中，在有圆形底面的定义域上的积分可以通过变换到柱极坐标系来完成；函数的变换用如下的关系进行：

f(x，y，z)→f(ρcosφ，ρsinφ，z)

区域的变换可以从图形中得到，因为底面的形状可能不同，而高遵循初始区域的形状。

例（3-a）:

区域为D＝{x²+y²≤9,x²+y²≥4，0≤z≤5}（也即底面为例2-d中的圆环的高度为5的"管道"）；如果采用变换，可以得到区域T＝{2≤ρ≤3,0≤φ≤π，0≤z≤5}（这是一个底面为例2-d中的矩形而高为5的长方体）。

因为z分量没有变化，dx dy dz和在极坐标中一样变化：变为ρdρdφdz。

最后，变换到柱极坐标的最后公式为：

∫∫∫_Df(x，y，z)dx dy dz＝∫∫∫_Tf(ρcosφ，ρsinφ，z)ρdρdφdz。

这个方法在柱形或者锥形区域的情况较为适用，也适用于容易分辨z区间和变换圆形底面和函数的其它情况。

例（3-b）:

函数为f(x，y，z)＝x²+y²+z而积分区域为圆柱:D＝{x²+y²≤9，-5≤z≤5}.

将D变换到柱极坐标如下：

T＝{0≤ρ≤3，0≤φ≤2π，-5≤z≤5}。

函数变为

f(ρcosφ，ρsinφ，z)＝ρ²+z

最有应用积分公式：

∫∫∫_D(x²+y²+z)dx dy dz＝∫∫∫_T(ρ²+z)ρdρdφdz;

推演一下公式，得到

${\∫}_{-5}^5\mathrm{dz}{\∫}_0^{2\mathrm{\π}}\mathrm{d\φ}{\∫}_0^3({\mathrm{\ρ}}^3+\mathrm{\ρz})\mathrm{d\ρ}=2\mathrm{\π}{\∫}_{-5}^5{[\frac{{\mathrm{\ρ}}^4}{4}+\frac{{\mathrm{\ρ}}^{2}z}{2}]}_0^3\mathrm{dz}$

$=2\mathrm{\π}{\∫}_{-5}^5(\frac{81}{4}+\frac{9}{2}z)\mathrm{dz}=\·\·\·=855\mathrm{\π}.$

球极坐标：

需要注意的是某些地区（如北美）角度标识相反。

参见图9，R3中，有些区域有球形对称性，所以将积分区域的每点用两个角度和一个距离标识较为合适。因此可以采用变换到球极坐标系；函数变换由如下关系产生：

f(x，y，z)→f(ρcosθsinφ，ρsinθsinφ，ρcosφ)

注意z轴上的点没有唯一表示，θ可以在0到2π间变化。

这个方法最为适用的区域显然是球。

例（4-a）:

区域为D＝x²+y²+z²≤16（球心在原点半径为4的球）；应用变换后得到：T＝{0≤ρ≤4，0≤φ≤π，0≤θ≤2π}。

坐标变换的雅可比行列式为：

$\frac{\∂(x,y,z)}{\∂(\mathrm{\ρ},\mathrm{\φ},\mathrm{\θ})}=\left|\begin{array}{ccc}\cos \mathrm{\φ}\sin \mathrm{\φ}& \mathrm{\ρ}\cos \mathrm{\θ}\cos \mathrm{\φ}& -\mathrm{\ρ}\sin \mathrm{\θ\φ}\sin \\ \sin \mathrm{\θ\φ}\sin & \mathrm{\ρ}\sin \mathrm{\θ\φ}\cos & \mathrm{\ρ}\cos \mathrm{\θ\φ}\sin \\ \cos \mathrm{\φ}& -\mathrm{\ρ}\sin \mathrm{\φ}& 0\end{array}\right|={\mathrm{\ρ}}^2\sin \mathrm{\φ}$

因此dx dy dz变换为ρ2sin(φ)dρdθdφ.

得到最后公式:

∫∫∫_Df(x，y，z)dx dy dz＝∫∫∫_Tf(ρsinθcosφ，ρsinθsinφ，ρcosθ)ρ²sinφdρdθdφ。

应当在积分区域为球形对称并且函数很容易通过基本三角公式简化的时候才使用这个方法。（参看例4-b）；其它情况下，可能使用柱极坐标更为合适（参看例4-c）。

∫∫∫_Tf(a，b，c)ρ²sinφdρdφ。

注意从雅可比行列式来的ρ²和sinφ因子。

注意下面例子中，φ和θ的作用反过来了。

例（4-b）:

D和例4-a相同，而f(x，y，z)＝x²+y²+z²是被积函数。

很容易变换为：

f(ρsinθcosφ，ρsinθsinφ，ρcosθ)＝ρ²，，

而从D到T的变换是已知的：

(0≤ρ≤4，0≤φ≤2π，0≤θ≤π)。

应用积分公式：

∫∫∫_D(x²+y²+z²)dx dy dz＝∫∫∫_Tρ²ρ²sinθdρdθdφ，

并展开：

${\∫\∫\∫}_T{\mathrm{\ρ}}^4\sin \mathrm{\θd\ρd\θd\φ}={\∫}_0^{\mathrm{\π}}\sin \mathrm{\θd\θ}{\∫}_0^4{\mathrm{\ρ}}^4\mathrm{d\ρ}{\∫}_0^{2\mathrm{\π}}\mathrm{d\φ}=2\mathrm{\π}{\∫}_0^{\mathrm{\π}}\sin \mathrm{\θ}{\left[\frac{{\mathrm{\ρ}}^5}{5}\right]}_0^4\mathrm{d\θ}$

$=2\mathrm{\π}{\left[\frac{{\mathrm{\ρ}}^5}{5}\right]}_0^4{[-\cos \mathrm{\θ}]}_0^{\mathrm{\π}}=4\mathrm{\π}\·\frac{1024}{5}=\frac{4096\mathrm{\π}}{5}.$

例（4-c）:

区域D是球心在原点半径为3a的球（D＝x²+y²+z²≤9a²）而f(x，y，z)＝x²+y²是被积函数。

看起来采用球极坐标变换较为合适，但是事实上，限定新区域T的变量很明显应该是：

0≤ρ≤3a，0≤φ≤2π，0≤θ≤π。

但是采用这个变换就有

f(x，y，z)＝x²+y²→ρ²sin²θcos²φ+ρ²sin²θsin²φ＝ρ²sin²θ.

应用积分公式得到：

∫∫∫_Tρ²sin²θρ²sinθdρdθ＝∫∫∫_Tρ⁴sin³θdρdθdφ

这很难求解。而如果采用柱极坐标，新的T区间为：

0≤ρ≤3a，0≤φ≤2π， $-\sqrt{9a^2-{\mathrm{\ρ}}^2}\≤a\≤\sqrt{9a^2-{\mathrm{\ρ}}^2};$

z区间可以通过将球切成两个半球并求解从D的公式来的不等式得到（然后直接变换x2+y2到ρ2）。新函数就是ρ2.采用积分公式

∫∫∫_Tρ²ρdρdφdz.

得到

${\∫}_0^{2\mathrm{\π}}\mathrm{d\φ}{\∫}_0^{3a}{\mathrm{\ρ}}^3\mathrm{d\ρ}{\∫}_{-\sqrt{9a^2-{\mathrm{\ρ}}^2}}^{\sqrt{9a^2-{\mathrm{\ρ}}^2}}\mathrm{dz}=2\mathrm{\π}{\∫}_0^{3a}2{\mathrm{\ρ}}^3\sqrt{9a^2-{\mathrm{\ρ}}^2}\mathrm{d\ρ}.$

然后应用变换

${9a}^2-{\mathrm{\ρ}}^2=t\→\mathrm{dt}=-2\mathrm{\ρd\ρ}\→\mathrm{d\ρ}=\frac{\mathrm{dt}}{-2\mathrm{\ρ}}$

(新区间变为0，3a→9a²，0)。得到

$-2\mathrm{\π}{\∫}_{{9a}^2}^0{\mathrm{\ρ}}^2\sqrt{t}\mathrm{dt}$

因为ρ²＝9a2^-t，所以

$-2\mathrm{\π}{\∫}_{{9a}^2}^0({9a}^2-t)\sqrt{t}\mathrm{dt},$

将积分限反过来，然后分配括号中的项，很容易将积分分解为可以直接积分的两部分：

$2\mathrm{\π}[{\∫}_0^{{9a}^2}{9a}^2\sqrt{t}\mathrm{dt}-{\∫}_0^{{9a}^2}t\sqrt{t}\mathrm{dt}]=2\mathrm{\π}{[{9a}^2\frac{2}{3}{t}^{\frac{3}{2}}-\frac{2}{5}{t}^{\frac{5}{2}}]}_0^{{9a}^2}$

$=2\·27\mathrm{\π}{a}^5(6-\frac{2}{5})=54\mathrm{\π}\frac{28}{5}{a}^5=\frac{1512\mathrm{\π}}{5}{a}^5.$

由于采用柱极坐标，很容易就将这个三重积分变换为简单的单变量积分。依据以上方案的基础知识，根据扫描出的原始数据，借助Matlab强大的数学计算工具，绘制出我们需要的曲面点阵图，依据前述步骤S104中所述分别对对各行和各列扫描点进行样条插值形成光滑曲线，对光滑曲线组成的曲面进行二重积分，就可以根据需要进行体积、密度、质量的计算了。

实施例二

图3是本实施例所述的粉体物料堆的体积测量系统框图，如图3所述，本实施例所述的粉体物料堆的体积测量系统包括：

激光测量装置，所述装置包括N个激光测量头，所述激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，确保两个最远的激光测量头之间的距离调整到等于粉料车的宽度，其中所述N为大于等于2的正整数；

其中，所述激光头的数量可根据煤炭运输车宽度、测量环境和不同用户测量精度的不同而进行相应的调整。

以粉体物料为煤炭为例，本实施例通过将12个激光测量头设置于粉料车上方的架子上，该架子与水平面平行，各激光测量头在架子上竖直安装，确保发射的激光束垂直于水平面，图2是本发明三维曲面生成示意图，如图2所示，架子上第1个激光点与第12个激光点之间的距离与车箱宽度相等，支架与车箱宽边平行，每个激光点投射的光线与车箱底面垂直。

激光控制装置，用于当所述粉料车的车箱驶入所述架子正下方时，控制所述激光测量装置中各所述激光测量头对所述粉体物料堆表面同时按照预设的相同频率进行扫描，用于获得各扫描点的高度数据，当所述粉料车的车箱驶出所述架子正下方时控制所述激光测量装置结束扫描，将所述各扫描点的高度数据发送给坐标记算软件模块；

例如，当有待检测的煤炭运输车进入该煤炭运输监控系统的工作区域时，12个激光头同时对煤堆自上而下进行照射扫描，并将激光在煤堆表面扫描点反馈的光信号通过模数转换器将光信号转换为数字信号，对所述数字信号进行编码和处理以获得粉体物料堆表面各扫描点的高度数据。当煤炭运输车的车尾经过激光头时，结束扫描。

坐标记算软件模块，用于接收所述激光控制装置发送的各扫描点的高度数据，利用计算机程序分别计算各扫描点在以粉料车纵向为X轴，以粉料车横向为Y轴，以垂直水平面的方向为Z轴的三维迪卡尔坐标系下的三维坐标，将所有扫描点的三维坐标发送给体积计算软件模块；

通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

例如，如图2所示，各激光测量头第1次扫描的扫描点的X坐标均为0；各激光测量头第2次扫描的扫描点的X坐标均为第2次扫描距第1次扫描时车箱移动的距离；各激光测量头第3次扫描的扫描点的X坐标均为第3次扫描距第1次扫描时车箱移动的距离；依次类推，各激光测量头第最后一次扫描的扫描点的X坐标均为该次扫描距第1次扫描时车箱移动的距离，也为车箱的长度值。

第1个激光测量头所测得的所有扫描点的Y坐标均为0，第2个激光测量头所测得的所有扫描点的Y坐标均为第2个激光测量头到第1个激光测量头的距离；第3个激光测量头所测得的所有扫描点的Y坐标均为第3个激光测量头到第1个激光测量头的距离；依次类推，第12个激光测量头所测得的所有扫描点的Y坐标均第12个激光测量头到第1个激光测量头的距离，也为车箱的宽度值。

各扫描点的Z值均由步骤S102所获得各扫描点的高度数据确定，需要说明的是，测量所得的直接高度数据是指激光测量头照射到的点到该激光测量头之间的高度，该扫描点在上述三维迪卡尔坐标系下的Z坐标需要将激光测量头到车箱底部的距离减去上述直接高度数据获得。

体积计算软件模块，用于接收所述坐标记算软件模块发送的所有扫描点的三维坐标，利用计算机程序分别对每个激光测量头扫描到的所有扫描点进行样条插值形成光滑曲线，分别对每个扫描时刻所有激光测量头扫描到的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

其中，样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差，这样就避免了使用高阶多项式所出现的龙格现象，样条插值具体方法有很多，例如，包括线性样条插值、二次样条插值、三次样条插值。

本实施例以二次样条插值为例，使用Matlab函数库来实现,主要用到的Matlab中的griddata函数、mesh函数、linspace函数。例如，对10个扫描点以二次样条插值为例，所述10个扫描点坐标是

A=[x1,y1,z1;x2,y2,z2;........;x10,y10,z10];

x=A(:,1);

y=A(:,2);

z=A(:,3);

N=100;//插值点数，实际使用时需确定，暂写100

x0=linspace(min(x),max(x),N);

y0=linspace(min(y),max(y),N);

[X,Y,Z]=griddata(x,y,z,x0’,y0,'v4');//插值拟合曲面。

figure;

mesh(X,Y,Z);//体积的近似值用下代码

L=sum(sum(abs(Z)));

S=abs((y0(2)-y0(1))*(x0(2)-x0(1)));

V=L*S;

所述V的值就是体积值。

依据该算法计算出粉料车内粉体物料堆的体积，进一步地，通过地磅测量出的所述粉体物料堆的质量数据，使用密度计算公式：密度=质量/体积，可求出密度。

其中，对各种不同的多重积分进行介绍和分析参见实施例一所述的对各种不同的多重积分进行介绍和分析部分，在此不作赘述。

本发明能准确真实地获取煤炭、粉料的现场数据，能节省大量的检验用人力物力和资金，可以在应用层解决用户的部分运输数据对比与监控问题，且可以让用户在第一时间得到准确真实的现场数据，为用户在进行详细标准正确的数据对比分析工作提供重要数据的辅助决策。本发明能极大的提升煤炭、粉料相关工作的检验效率，能客观地监督运输方使用合法、安全、正规的手段进行煤炭、粉料的运输操作。

作为变形，本发明除了用于运输数据对比与监控方面以外，也可应用于其他很多方面，例如，用于以粉体物料为源料的生产或对粉体物料产品进行包装等应用中，用于对传送带上粉体物料的体积测量，密度测量等。

以上实施例提供的技术方案中的全部或部分内容可以通过软件编程实现，其软件程序存储在可读取的存储介质中，存储介质例如：计算机中的硬盘、光盘或软盘。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种粉体物料堆的体积测量方法，其特征在于，包括：

S1、撘建激光测量装置，将N个激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，将两个最远的激光测量头之间的距离调整到等于粉料车的宽度，N为大于等于2的正整数；

S2、当粉料车的车箱驶入所述架子正下方时，控制各所述激光测量头对所述粉体物料堆表面同时按照预设的相同频率进行扫描，获得各扫描点的高度数据，当所述粉料车的车箱驶出所述架子正下方时结束扫描；

S3、利用计算机程序以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴构建三维迪卡尔坐标系，计算各扫描点在所述三维迪卡尔坐标系下的三维坐标；

S4、利用计算机程序分别对每行和每列的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

2.如权利要求1所述的粉体物料堆的体积测量方法，其特征在于，步骤S3中所述计算各扫描点在所述三维迪卡尔坐标系下的三维坐标具体包括：通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

3.如权利要求1或2所述的粉体物料堆的体积测量方法，其特征在于，步骤S4中所述样条插值方法为二次样条插值。

4.如权利要求3所述的粉体物料堆的体积测量方法，其特征在于，步骤S1中所述激光测量头个数为12。

5.一种粉体物料堆的密度测量方法，其特征在于，包括，通过地磅测量出的所述粉体物料堆的质量数据，通过如权利要求1-4之一所述的粉体物料堆的体积测量方法测量所述粉体物料堆的体积数据，通过所述质量数据除以所述体积数据获取所述粉体物料堆的密度数据。

6.一种粉体物料堆的体积测量系统，其特征在于，包括：

激光测量装置，所述装置包括N个激光测量头，所述激光测量头并排固定于粉料车通过处的上方架子上，测量之前确保各所述激光测量头发射的激光束垂直于水平面，确保两个最远的激光测量头之间的距离调整到等于粉料车的宽度，其中所述N为大于等于2的正整数；

激光控制装置，用于当所述粉料车的车箱驶入所述架子正下方时，控制所述激光测量装置中各所述激光测量头对所述粉体物料堆表面同时按照预设的相同频率进行扫描，用于获得各扫描点的高度数据，当所述粉料车的车箱驶出所述架子正下方时控制所述激光测量装置结束扫描，将所述各扫描点的高度数据发送给坐标记算软件模块；

坐标记算软件模块，用于接收所述激光控制装置发送的各扫描点的高度数据，利用计算机程序分别计算各扫描点在以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴的三维迪卡尔坐标系下的三维坐标，将所有扫描点的三维坐标发送给体积计算软件模块；

体积计算软件模块，用于接收所述坐标记算软件模块发送的所有扫描点的三维坐标，利用计算机程序分别对每行和每列的扫描点进行样条插值形成光滑曲线，对所述光滑曲线组成的曲面进行二重积分获取所述粉料车中粉体物料堆的体积。

7.如权利要求6所述的粉体物料堆的体积测量系统，其特征在于，所述坐标记算软件模块中所述利用计算机程序分别计算各扫描点在以平行于粉料车行驶方向的方向为X轴，以垂直于粉料车行驶方向的方向为Y轴，以垂直于水平面的方向为Z轴的三维迪卡尔坐标系下的三维坐标具体包括：通过激光测量头在第一次扫描和对某扫描点扫描期间所述粉料车行驶的距离确定该扫描点的X坐标；通过各激光测量头所在的位置的Y坐标确定各激光测量头扫描到的所有扫描点Y坐标；通过各扫描点的高度数据确定各扫描点的Z坐标。

8.如权利要求6或7所述的粉体物料堆的体积测量系统，其特征在于，体积计算软件模块中所述样条插值方法为二次样条插值。

9.如权利要求8所述的粉体物料堆的体积测量系统，其特征在于，所述激光测量装置中所述激光测量头个数为12。

10.一种粉体物料堆的密度测量系统，其特征在于，包括

质量测定模块，用于通过地磅测量出的所述粉体物料堆的质量数据；

体积测定模块，用于通过如权利要求6-9之一所述的粉体物料堆的体积测量系统测量所述粉体物料堆的体积数据；

密度计算模块，用于通过所述质量数据除以所述体积数据获取所述粉体物料堆的密度数据。

Images