基于图像处理的铁路路基设计土石方计算方法

Earthwork Calculation Method for Railway Subgrade Design Based on Image Processing

摘要: 本文讨论了铁路路基设计中土石方计算的挑战以及现有方法的局限性，并提出了一种基于图像处理的新型计算方法，该方法将地质剖面图中地质界线进行一定的修正，再将其按一定比例系数放大并绘制到特定位图中，进而利用图像处理技术快速识别出各地层区域轮廓，并提取出每个地层的区域边界范围，最后通过预定义的填挖设计线构造填挖区域，利用区域之间的布尔运算，实现了铁路路基设计中断面填挖土石方的分类计算，极大提高了计算效率和准确性。

Abstract: This paper discusses the challenges of earth and rock calculation in railway subgrade design and the limitations of existing methods, and proposes a new calculation method based on image processing, which corrects the geological boundary in geological profile to a certain extent, amplifies it according to a certain scale factor and draws it into a specific bitmap, and then uses image processing technology to quickly identify the contours of various layers. The boundary range of each stratum is extracted, and the filling area is constructed by the pre-defined cut-fill design line, and the classification calculation of the section filling and excavating earth and rock in the railway subgrade design is realized by using the Boolean operation between the areas, which greatly improves the calculation efficiency and accuracy.

1. 引言

在铁路路基设计中，需要计算出工程设计中形成的填挖土石方量[1] [2]，用于土石方调配[3]，以便后期寻找取弃土场，或采购填料。在计算土石方量时，由于人工或填料价格会随土石类别变化，因此不仅需要计算土石方量，对于挖方还需要根据地层土的类型详细区分每一个类型土的开挖方量，对于填方也要根据填料的不同位置区分不同的填料类型。

目前路基设计中土石方计算主要是基于地质横断面[4] [5] (也称地质剖面)，但是由于地质横断面可能错综复杂(如包含岩脉、透镜体、断层等)，难以数字化解析，且目前商业软件基本只能计算横断面总体土石方量[6] [7]，无法按地层分类计算，因此传统的土石方计算方法主要为手动量取，即人工量取每个填挖区域面积，再分类统计，该方法需要大量手动操作，效率极低。

苏伟、白青波等[8] [9]提出并公开了一种地质横断面数字化解析方法，可以根据地质剖面图的地层线和地层名及其位置，解析出各个地层线的上下部地层名，再通过条分法计算各个地层的开挖面积，实现了土石方计算程序化，提高了设计计算效率。但该方法复杂度较高，且对地层信息解析程度较低，对于岩脉、断层等特殊地质构造仍存在一定的局限性。

土石方计算的最终目的是需要计算出每一个地层的开挖量，只要能计算出每一个地层的区域范围边界，其与开挖区域的交集即为该地层的开挖区域，进而得到各个地层的开挖量。此时，关键就在于计算各个地层的区域范围边界。因此，本文提出一种基于图像处理的方法，通过将地质剖面绘制到特定位图中，利用图像处理技术识别区域轮廓，将每个地层对应的地层区域范围边界提取出来，进而利用地层所在区域范围边界，根据多边形的布尔运算方便快捷地计算出每个地层的开挖面积，实现土石方的分类计算。

2. 修正地质界线

在地质剖面中，任意建立一个局部平面坐标系，则可得知地质界线(地面线、地层分界线)的各个关键点坐标、地层名及其位置坐标。根据地质界线的关键点坐标和地层名的位置坐标，计算出这些坐标的矩形包围盒，即该矩形左侧范围为所有点坐标的X值最小值，右侧范围为所有点坐标的X值最大值，顶部范围为所有点坐标的Y值最大值，底部范围为所有点坐标的Y值最小值，如图1。

遍历所有地质界线，获取其起点和终点判断该地质界线的起点和终点与其他地质界线有无交点：

如果起点无交点，则说明起点处并不是该分界线的结束点，将地质界线沿起点处方向延伸至与矩形包围盒相交(或超出)；如果终点无交点，则说明终点处也不是该分界线的结束点，将地质界线沿终点处切线方向延伸至与矩形包围盒相交(或超出)。修正完的地质剖面如图2所示。

Figure 1. Schematic diagram of rectangular bounding box

图1. 矩形包围盒示意图

Figure 2. Schematic diagram of geological boundary correction

图2. 地质界线修正示意图

3. 创建地质剖面位图

以矩形包围盒的高$\text{H}$ 、宽$\text{W}$ 为尺寸，放大系数为$\text{k}$ ，创建尺寸高$\text{k}\cdot \text{H}$ 、宽$\text{k}\cdot \text{W}$ 的白色背景的位图(Bitmap，使用像素点来描述图像，也称为点阵图像)，并将矩形包围盒和地层分界线的关键点坐标乘以放大系数$\text{k}$ 后，以黑色线(像素点拟合的线)将其绘制到位图中，如图3。

其中，k为单位距离所包含的像素点个数。由于地质剖面中，点坐标单位可能是米m (dm、cm等)，其数值放到位图中单位则是像素，即1个像素点表示1 m，可能造成精度降低，因此为保证计算精度，需对地质界线坐标进行适当缩放。例如，地质界线单位为m，则可以取$\text{k}=100$ ，即位图中1000个像素点表示1 m。相应地，放大系数$\text{k}$ 值越大，计算精度越高，计算量也会越大。特别地，当不需要缩放时，可以取$\text{k}=1$ 。

可以看出，所有白色像素区域即为地层区域(或无效区域)，白色像素区域的边界即为地层区域的边界。利用OpenCV等图像处理库对地质剖面位图进行二值化处理[10] [11]，使地层线(黑色)与其他区域(白色)形成明显对比，随后利用其轮廓查找技术，准确定位出图中所有白色像素区域的边界，获取多个边界像素坐标点集。

Figure 3. Schematic diagram of geological profile bitmap

图3. 地质剖面位图示意图

遍历所得边界，在每个边界内部填充不同的颜色，并在上述位图中绘制出来，如图4，即将边界区域内的像素点赋予颜色，同时将每个边界与颜色值建立对应关系，便于通过颜色值找到对应的边界点集。

Figure 4. Schematic diagram of regional boundary filling

图4. 区域边界填充示意图

4. 获取地层区域名称

在获取所有区域边界(点集)后，此时区域边界对应的地层还未知，因此需提取该区域所对应的地层名。遍历所有区域边界，对于某一个区域边界：

遍历所有地层名，提取该地层名的位置点坐标${\text{p}}_{\text{i}}$ ，将其通过放大系数$\text{k}$ 进行缩放得到新坐标${\text{p}}_{\text{ki}}$ ，在上述位图中提取坐标${\text{p}}_{\text{ki}}$ 处的像素点的颜色值，根据颜色值即可找到与之对应的区域边界点集。将区域边界点集除以放大系数$\text{k}$ ，可得到与实际地层坐标相对应的地层区域边界。此时，将地层区域边界(点集)与该地层名建立对应关系，并存储到地层列表中。

遍历完成后，地层列表中所有建立对应关系的地层名和地层区域边界则是最终有效的地层区域，没有建立对应关系的区域则放弃(如地面线以上无效区域，则找不到对应的地层名)。最终获取地层名称后的区域，如图5。

5. 定义结构填挖设计线

根据李伯根等[12]提出的链式路基设计方法，在定义结构的同时，增加各个结构填挖设计线的定义。此处所谓的填挖设计线，是指结构在设计时，对应位置的开挖线(或回填线)，通常为结构某两个关键点(如前后连接点)之间与地层相接触的轮廓线(点集)，如图6。将路基工程中涉及的结构，如侧沟、路堑边坡、挡土墙、锚固桩、桩板墙等，全部定义出填挖设计线。

Figure 5. Schematic diagram of stratum-named regions

图5. 获取地层名称后的区域示意图

Figure 6. Schematic diagram of cut-fill design line (black) in roadbed structure

图6. 路基工程结构填挖设计线(黑色线)示意图

6. 更新填挖设计线

在放置、修改或删除路基工程结构时，各个路基工程结构都会根据其前连接点更新自身的图形关键点、连接点和属性信息[12]。此时根据更新后的结构关键点和属性信息，计算更新结构的填挖设计线。

7. 连接填挖设计线

在设计完成后，各个结构的填挖设计线都已计算更新完成。在进行填挖土石方计算时，根据左侧结构列表中工程结构顺序，依次提取出各个结构的填挖设计线点集，并存储到设计线点集列表中。由于路基结构工程列表是按照从路基面往两侧为顺序进行的存储，因此左侧填挖设计线点集也是同样的顺序，此时将列表中的点集按照相反的顺序重新排列。路基面的填挖设计线为基床表层的底面线，将其点集添加到设计线点集列表中。再根据右侧结构列表中工程结构顺序，依次提取出各个结构的填挖设计线点集，存储到设计线点集列表中。

此时，填挖设计线点集列表中的点均为从左往右的顺序(也可从右往左)，连接填挖设计线点集列表中的点所形成的线，即为整个设计横断面的填挖设计线，如图7，据此可进行整个横断面填挖土石方的计算。

8. 计算挖方区域

根据提取的填挖设计线点集，计算出设计线的Y坐标最大值${\text{topY}}_1$ ，并获得该设计线的$\text{X}$ 坐标范围$\left[\text{minX},\text{maxX}\right]$ ，同时提取地质横断面中的地面线，根据其点集计算出地面线的Y坐标最大值${\text{topY}}_2$ 。此时设计区域Y坐标最大值$\text{topY}$ 即为两者较大值，即$\text{topY}=\text{Max}\left\{{\text{topY}}_1,{\text{topY}}_2\right\}$ 。为计算方便，可将Y坐标最大值往上延伸一定距离$\text{exDist}$ ，即$\text{topY}=\text{topY}+\text{exDist}$ ，其中延伸距离满足$\text{exDist}>0$ 。

在填挖设计线点集列表中添加两个点$\text{p}1\left(\text{maxX},\text{topY}\right)$ 、$\text{p}2\left(\text{minX},\text{topY}\right)$ ，此时设计线所形成封闭多边形区域如图8。

Figure 7. Schematic diagram of cross-sectional cut-fill design line (red) in roadbed

图7. 路基横断面填挖设计线(红色线)示意图

Figure 8. Formation area of cut-fill design line

图8. 填挖设计线形成区域

在路基设计中，不仅路基工程结构会形成土石方的填挖，对地基的挖除换填，同样会造成填挖，因此还需要考虑地基的挖除换填形成的土石方。提取挖除换填线的坐标点集，并存储到挖除换填点集列表中，获得该点集列表的$\text{X}$ 坐标范围$\left[\text{min<msup> X <mo>′</mo> </msup>},\text{max<msup> X <mo>′</mo> </msup>}\right]$ ，再同样在列表中添加两个点$\text{p}3\left(\text{max<msup> X <mo>′</mo> </msup>},\text{topY}\right)$ 、$\text{p}4\left(\text{min<msup> X <mo>′</mo> </msup>},\text{topY}\right)$ 。此时，挖除换填线所形成封闭多边形区域如图9。

Figure 9. Formation area of excavation-replacement lines

图9. 挖除换填线形成区域

将填挖设计线形成区域和挖除换填线形成区域合并起来(并集)，即为挖方计算区域，如图10。

Figure 10. Computational domain for excavation quantities

图10. 挖方计算区域

将待计算地质横断面中的所有地层区域边界提取出来，如图11。此时，遍历所有地层区域，依次将地层区域与挖方计算区域求交集，如果没有交集则表示该地层未开挖，否则所得交集区域即为该地层的开挖区域，计算其区域面积即可得到该地层的开挖方量，如图12。

Figure 11. Schematic diagram of stratum region extraction

图11. 提取地层区域

Figure 12. Schematic diagram of excavation calculation results

图12. 挖方计算结果示意图

9. 计算填方区域

与挖方区域计算类似，提取填挖设计线点集，计算出设计线的Y坐标最小值$\text{botY}$ 。为计算方便，可将Y坐标最小值往下延伸一定距离$\text{exDist}$ ，即$\text{botY}=\text{botY}-\text{exDist}$ 。在填挖设计线点集列表中添加两个点$\text{p}5\left(\text{maxX},\text{botY}\right)$ 、$\text{p}6\left(\text{minX},\text{botY}\right)$ ，此时设计线所形成封闭多边形区域如图13。

Figure 13. Formation area of fill design lines

图13. 填方设计线形成区域

提取地面线点集，并存储到地面线点集列表中，获得该点集列表的$\text{X}$ 坐标范围$\left[\text{min<msup> X <mo>″</mo> </msup>},\text{max<msup> X <mo>″</mo> </msup>}\right]$ ，在列表中添加两个点$\text{p}7\left(\text{max<msup> X <mo>″</mo> </msup>},\text{topY}\right)$ 、$\text{p}8\left(\text{min<msup> X <mo>″</mo> </msup>},\text{topY}\right)$ 。此时，地面线所形成封闭多边形区域如图14。

Figure 14. Formation area of natural ground profile

图14. 地面线形成区域

同样，仍需要考虑地基的挖除换填形成的土石方。提取挖除换填线的坐标点集，并存储到挖除换填点集列表中，再同样在列表中添加两个点$\text{p}3$ 、$\text{p}4$ 。此时，挖除换填线所形成封闭多边形区域如图9。

将地面线形成区域和挖除换填线形成区域合并起来(并集)，即为填方计算区域，如图15。此时，将填挖设计线形成区域，与填方计算区域求交集，所得交集区域即为所求填方区域，计算其区域面积即可得到该地质横断面的回填方量，如图16。

Figure 15. Computational domain for fill quantities

图15. 填方计算区域

Figure 16. Schematic diagram of fill calculation results

图16. 填方计算结果示意图

10. 实例分析

以某铁路工程路基横断面为例，如图17。其土石方数量计算按照传统计算方法，可利用辅助软件求取设计断面与自然地形断面的闭合区域面积，即为该断面总挖方量，如图18(a)。若在地质条件单一的情况下，尚可满足设计需求，然而，当遇到地质条件复杂的情况，例如地层层数较多或含有透镜体时，现有辅助软件无法区分不同土石类别，因此，仍需设计人员手动在cad中绘制出透镜体与开挖线形成的封闭区域，通过计算其面积来获取各类土石方量，如图18(b)。这一过程不仅耗时费力，还容易因人为操作失误导致结果偏差。

Figure 17. Roadbed cross-section schematic diagram

图17. 路基横断面示意图

(a) 断面总挖方量自动计算

(b) 断面不同土石类别手动框算

Figure 18. Schematic of traditional earthwork calculation

图18. 传统土石方计算示意

本文提出的基于图像处理的土石方计算方法，通过精准识别不同地层区域，实现了不同土石类别的分类计算，如图19。由图可知，该方法能够直接识别并计算各类土石方量，无需人工进行二次处理，极大提高了复杂地质条件下土石方计算效率。从对比结果来看，该方法计算得到的土石方量与手动框算值的差值约为0.2‰，这表明其在精度上与传统手动方法相当，但优势在于显著减少了人工干预，避免了人为失误，提升了计算效率。

Figure 19. Schematic of image-based earthwork calculation

图19. 基于图像处理的土石方计算示意

11. 结论

本文提出的基于图像处理的土石方计算方法，以地质剖面图为基础，通过将地质界线绘制到特定位图中，利用图像处理技术识别区域轮廓，将每个地层对应的地层区域范围边界提取出来，再在所得地层区域范围边界的基础上，构造填挖土石方计算区域，利用地层区域与填挖区域边界所形成多边形区域间的布尔运算，实现了对填挖土石方的分类计算，克服了传统方法的局限性，可极大提高铁路路基工程设计中土石方计算效率与准确性。

该方法不仅适用于简单地层，对于包含断层、岩脉、透镜体等特殊地质构造的地层同样具有很好的效果。这一方法为铁路路基工程设计提供了可靠的土石方计算方案，可以大大提升设计效率和结果准确性。