2.1. Overview

The proposed method considers the density of the target soil to estimate its properties and perform autonomous excavation based on that information. Density significantly influences the penetration resistance of soil^13,14, which is critical in planning an excavation. In this study, we proposes a method to estimate soil density from 3D point-cloud data and generate optimal excavation path parameters based on the estimated density. An overview of the proposed method is presented in Fig. 2. The flow of the main process is as follows.

1. 3D point-cloud data of the ground before and after excavation are obtained.

2. The volume of the soil pile before and after excavation is calculated using the 3D point cloud data.

3. The bulking factor of soil is calculated from the volume.

4. The soil density is estimated from the bulking factor of the soil using a relational expression.

5. Soil with the estimated density is replicated in a physical simulator.

6. Optimal excavation path parameters for the soil are obtained using the simulator.

7. An automated excavator excavates the ground along the excavation path defined by the obtained path parameters.

Note that this method must be applied in a repeated excavation cycle because it requires 3D point-cloud data about the ground before and after excavation. Moreover, the method assumes that the current excavation point is close to the previous one and that their soil properties are nearly identical. Although the proposed method includes searching for an excavation path using a simulator, that would be too time-consuming in practical applications. Realistically, we consider that searching for candidate optimal excavation paths in advance for multiple expected soil conditions using the simulator as a preliminary process and then selecting the most appropriate path based on the estimated soil properties during excavation would be more practical.

Items 1–4 correspond to soil property estimation, whereas items 5–7 correspond to optimal path generation and autonomous excavation. In the estimation process, the input is a 3D point-cloud data of the target ground, and the output is the estimated soil density. In the path generation process, the input is the estimated soil density and the output is the optimal parameters representing an excavation path. Details of each process are provided in the following sections.

2.2. Soil Density Estimation Based on the Bulking Factor Using 3D Point-Cloud Data

2.2.1. Basic Principle of Soil Density Estimation

As mentioned above, soil density affects penetration resistance, which is critical in excavation planning. In general, direct measurement methods such as the core method are used to determine soil density at construction sites. However, these methods are not applicable during the excavation process because they are time-consuming and destructive ¹⁵. Therefore, we propose a new method here to measure or estimate soil density during the excavation process.

To estimate soil density, we focused on the bulking factor of the soil. Bulking factor is an indicator of the increase in soil volume that occurs during the excavation process ¹⁶. Soil volume increases after excavation due to changes in the structure of the soil that create additional void space. Experimental works have established that the bulking factor and soil density exhibit a strong positive correlation ¹⁷. Qualitatively, a large bulking factor indicates that the soil is dense and hard. Moreover, the soil volumes before and after excavation, which are necessary to calculate the bulking factor, can be derived from 3D point-cloud data acquired using sensors that are relatively easy to install such as LiDAR or RGB-D sensors. In particular, many studies have explored the use of RGB-D sensors on construction sites ^18,19,20. Thus, in this study, we calculated the bulking factor from 3D point-cloud data and estimated soil density based on the bulking factor. The procedure used to estimate the density of the soil is summarized in Fig. 3. The details of each step are described below.

2.2.2. Calculation of Soil Volume and Bulking Factor Using 3D Point-Cloud Data

First, the soil volume and bulking factor are calculated using 3D point-cloud data.

The bulking factor quantifies the increase in soil volume. Let \(k\) denote the bulking factor. It is defined as follows.

\begin{align} k = \frac{V_{\rm{after}}}{V_{\rm{before}}}, \label{eq:eq15} \end{align}

where \(V_{\rm{before}}\) and \(V_{\rm{after}}\) denote the soil volumes before and after excavation, respectively.

The system shown in Fig. 4 is used to calculate the volume of the soil. As shown in the figure, two areas are measured, including the excavation and dumping areas. 3D point-cloud data are obtained from each area.

• \(P_{\rm{before}}\): 3D point-cloud data of the excavation area before excavation.

• \(P_{\rm{after}}\): 3D point-cloud data of the excavation area after excavation.

• \(p_{\rm{before}}\): 3D point-cloud data of the dumping area before excavation.

• \(p_{\rm{after}}\): 3D point-cloud data of the dumping area after excavation.

The proposed method calculates the following three soil volumes from the 3D point-cloud data, as illustrated in Fig. 4.

• \(V_{\rm{up}}\): The volume of spilled soil.

• \(V_{\rm{down}}\): The volume of excavated soil.

• \(v_{\rm{up}}\): The volume of dumped soil.

The point cloud is processed for this calculation. First, changes in the 3D point-cloud data before and after excavation are detected. Subsequently, outliers in the point cloud are removed. Then, a voxel grid filter is applied to the point cloud. After applying the voxel grid filter, let the \(i\)-th point of the processed point cloud be \(p_i = (x_i, y_i, z_i)\). The soil volumes are then calculated as

\begin{align} V(\Omega, s) = \sum_{p_i \in \Omega} \bigl( s \cdot z_i \bigr)^+ \, l^2, \end{align}

where \((x)^+ = \max(x,0)\) denotes the operator that extracts the positive part, \(\Omega\) denotes the target region (excavation area \(\Omega_{\text{exca}}\) or dumping area \(\Omega_{\text{dump}}\)), and \(l\) is the horizontal grid spacing, i.e., the side length of the voxel in the base plane. From this formula, the volumes are obtained as \(V_{\text{up}} = V(\Omega_{\text{exca}}, +1)\), \(V_{\text{down}} = V(\Omega_{\text{exca}}, -1)\), \(v_{\text{up}} = V(\Omega_{\text{dump}}, +1)\). Finally, the volumes before and after excavation, \(V_{\rm{before}}\) and \(V_{\rm{after}}\), are defined as follows.

\begin{align} V_{\rm{before}} & = V_{\rm{down}}, \label{eq:eq19} \\ \end{align}

\begin{align} V_{\rm{after}} & = \ V_{\rm{up}} + v_{\rm{up}}. \label{eq:eq20} \end{align}

The bulking factor can be calculated by substituting \(V_{\rm{before}}\) and \(V_{\rm{after}}\) into Eq. \(\eqref{eq:eq15}\).

2.3. Generating Efficient Excavation Paths Based on Estimated Soil Density

2.3.2. Representation of Excavation Path

The proposed method applies an excavation path represented by a set of parameters. We consider an excavation path represented by three parameters as shown in Fig. 5. This path is based on excavation operations performed by experienced operators. Several studies have investigated excavation operations performed by experienced operators ^21,22. In general, such operations can be divided into three phases, including penetrating, dragging, and scooping. In the penetration phase, operators insert the bucket into the target soil until it reaches the bottom line of the desired excavation area while keeping the bottom of the bucket parallel to the penetration line to minimize resistance force. Subsequently, in the dragging phase, the tip of the bucket is dragged horizontally along the bottom line while maintaining a fixed angle between the bucket and the arm. Finally, in the scooping phase, the operators rotate only the bucket until it is oriented horizontally. To model these operations, we consider the following three parameters.

• \(\theta_{\rm{penetrate}}\): The angle between the horizontal surface and the penetration line.

• \(l_{\rm{penetrate}}\): The length of the penetration line.

• \(l_{\rm{drag}}\): The length of the dragging line.

These parameters determine the bucket posture at a given moment, which in turn determines the joint angles of the boom, arm, and bucket. Excavation along the planned path can be carried out by controlling each joint according to these parameters.

In this study, we made two assumptions regarding the excavation model. First, we assumed that the excavator does not move laterally while the bucket is in contact with the ground. This implies that the excavation model considers only the \(XZ\) plane. Second, we assumed that the position of the main body of the excavator is fixed.

2.3.3. Generation of Optimal Excavation Path Using a GA

Finally, an efficient excavation path is generated for the target soil. Thus, our approach finds the values of three parameters that represent an excavation path with the highest efficiency for soil of a given density.

To that end, the proposed method searches for an excavation path by applying a GA ²³. The GA is an evolutionary optimization algorithm inspired by natural selection. At each generation \(t\), a population of candidate solutions \(\{ x_1^t, x_2^t, \dots, x_N^t \}\) is maintained. The quality of each individual \(x_i^t\) is evaluated using a fitness function \(f(x_i^t)\). Individuals with higher fitness values are more likely to be selected as parents and new individuals are generated through crossover and mutation to form the next generation \(\{ x_1^{t+1}, x_2^{t+1}, \dots, x_N^{t+1} \}\). Through this iterative process, the population converges toward an optimal or nearly optimal solution without requiring gradient information. The GA is suitable for this study because the gradient of the excavation efficiency cannot be calculated.

Figure 6 illustrates the concept of searching for an optimal excavation path using a GA. In this method, the GA searches for the best excavation path by testing various paths over a pre-defined number of generations using a physics simulator that replicates the hydraulic excavator and the target soil. The fitness function of each individual is defined by the excavation efficiency \(f\) as described above. Therefore, the GA is expected to return an efficient excavation path that minimizes energy consumption while maximizing the amount of excavated soil.

In this study, each individual in the GA population was represented as a vector of the following three parameters. \((\theta_{\text{penetrate}}, l_{\text{penetrate}}, \textrm{and}~l_{\text{drag}})\), respectively represent the angle and length of penetration and the dragging length of the excavation path. In the initialization step, these parameters are assigned random values within specified ranges such that each individual corresponds to a distinct excavation motion. The fitness of each individual is then calculated using the excavation efficiency \(f\) defined above. To generate the next generation, the elitism and tournament selection methods are combined. In elitism, the best individual of the current generation is directly preserved in the next generation. For tournament selection, a small group of individuals is randomly sampled from the remaining population and the one with the highest fitness among them is copied into the next generation. Selected parents undergo crossover using a two-point recombination scheme in which contiguous segments of their parameter vectors are exchanged. After crossover, mutation is applied by randomly altering one or more parameters within their predefined ranges to maintain population diversity and prevent premature convergence. After a predefined number of generations, the individual with the highest fitness is taken as the optimal excavation path.

As mentioned above, it is considered practical to precompute candidate excavation paths for multiple expected soil conditions using this method. In actual operations, excavation can be performed without delay by selecting an optimal excavation path from the precomputed candidates based on the estimated soil density.

3.1. Experimental Setting

Simulations were conducted to evaluate the effectiveness of the proposed method.

In the experiments, we used the dynamic simulator Vortex Studio developed by CM Labs Simulations. Vortex Studio enables the replication of soil–tool interactions during the excavation process, which allows the mass of excavated soil and energy consumption to be evaluated. This simulator has been evaluated in a previous study in which soil behavior such as the angle of repose was reproduced ²⁴. In the simulator, an excavator model performed autonomous excavation operations. The model was based on the ZAXIS120 developed by Hitachi Construction Machinery, and was identical to the excavator used in the field experiments. Three soil density conditions were prepared, including soft, medium, and hard soil, as shown in Table 1. The soil type was set to loam, which was the same as that used in the field experiments. In addition, the initial ground surface was set to be flat. Using these settings, the excavation efficiency of the proposed method, which considers the bulking factor of the soil, was compared with that of excavation using fixed parameters. For the fixed parameters, \([\theta_{\rm{penetrate}}, l_{\rm{penetrate}}, l_{\rm{drag}}]\) were set to \([\pi/4, 0.35, 1.5]\). These values were determined empirically to ensure that a sufficient excavation volume was obtained to calculate the bulking factor.

The relationship between the bulking factor and soil density was determined in advance in the dynamic simulator. Fig. 7 shows the results of the preliminary experiment along with the regression line. The parameters in Eq. \(\eqref{eq:density}\) were determined as \(a = 161.57\) and \(b = -152.87\). Given that the coefficient of determination was \(R^2 = 0.9811\), estimating the soil density from the bulking factor using a linear function did not pose a significant problem.

Candidates for efficient excavation paths were precomputed using the GA in the simulator. The conditions for the GA are listed in Table 2 and the lower and upper bounds of the genes for each individual are given in Table 3. The number of generations was set to be relatively small due to computational time constraints and to obtain reasonable results within the available resources. Under these settings, candidates for efficient excavation paths were obtained for all density patterns ranging from 0% to 100% in continuous increments of 5%.

The detailed procedure for the simulation experiments is as follows.

1. The target ground is excavated once using fixed parameters.

2. The bulking factor is calculated using the height map of the target ground obtained directly from the simulator, instead of point cloud data obtained from RGB-D sensors or similar devices.

3. An efficient excavation path is selected based on the soil density estimated from the bulking factor.

4. The target ground is excavated again using the selected efficient excavation path.

3.2. Experimental Results

Figure 8 shows the experimental results for excavation efficiency under the three conditions. For each condition, the left bar represents the excavation efficiency using fixed parameters, while the right bar represents the excavation efficiency using the proposed method. As shown in Fig. 8, the proposed method improved excavation efficiency under all conditions. These results confirm the validity of the proposed method in the simulation environment.

For further discussion, Fig. 9 presents a comparison of the excavation paths. In each figure, the solid line represents the path of the bucket tip with fixed parameters, while the dashed line represents the path of the bucket tip generated by the proposed method. Here, \(x = 0\) is defined as the \(x\)-coordinate of the joint between the base and boom links. Similarly, \(z = 0\) is defined as the \(z\)-coordinate of the initial ground surface. As shown in Fig. 9, the excavation paths generated by the proposed method are adapted to each condition. In particular, it may be observed that the excavation area became larger as the target soil became harder. This result may be related to the definition of excavation efficiency used in this work, in which both the excavated volume and energy consumption are considered. For hard soil, the optimization may have selected a deeper and longer path to compensate for the larger energy required in penetration by increasing the excavated volume. This phenomenon will be addressed in more detail in future work.

4.1. Experimental Setup

Field experiments were conducted to evaluate the proposed method for actual excavators. Fig. 10 shows the experimental setup used in the field experiments.

In this experiment, we used a ZAXIS120 hydraulic excavator developed by Hitachi Construction Machinery. It was equipped with three rotary encoders, including a GNSS, a GNSS compass, a fuel flow meter, and a wireless LAN. The control input to operate the excavator can be transmitted via a wireless connection. The electromagnetic proportional valves of the hydraulic cylinders were opened or closed based on the input. Therefore, the excavator was automatically controlled by transmitting the appropriate input via a wireless connection. This system was implemented using robot operating system (ROS) ^a. Using this system, an excavator can excavate a predetermined path. For external sensing, we used an Intel RealSense depth camera D435i, developed by Intel, as an RGB-D sensor.

The soil type was loam and the initial ground surface was leveled to be almost flat. The experiment was conducted under two soil conditions, including soft and hard soil. The hardness of the soil was verified in advance using the coefficient of subgrade reaction. Table 4 shows the coefficients of subgrade reaction for each condition measured using a falling weight deflectometer. A larger coefficient of subgrade reaction indicates harder soil. In addition, the same relationship between the bulking factor and soil density described in Section 3 was applied here.

Two types of excavation paths were used, including a fixed-parameter path and an efficient path generated using the proposed method. The fixed-parameter path was used to estimate the soil density and for comparison with the proposed method. For the fixed parameters, \([\theta_{\rm{penetrate}}, l_{\rm{penetrate}}, l_{\rm{drag}}]\) were set to \([1.023, 0.496, 2.487]\), respectively. The parameters were changed from those used in the simulation experiments to ensure a sufficient volume of excavated soil to calculate the bulking factor. When using the same fixed parameters as in the simulation experiments, the volume of the excavated soil was sometimes insufficient due to the roughness of the ground surface and other factors. The same candidate excavation paths used in Section 3 were applied in this experiment. When the estimation of the soil density was sufficiently accurate, an efficient excavation path was selected based on the estimated density. An excavation was performed in an adjacent area where the soil conditions were confirmed to match those of the area used for the fixed-parameter path.

Although the mass of the soil and the energy consumption are necessary to calculate efficiency as shown in Eq. \(\eqref{eq:eq21}\), measuring them directly can be difficult. Thus, the volume of soil accumulated in the bucket after excavation was used instead of the mass of the soil. The volume corresponds to the mass by multiplying it by the density, which was measured using the RGB-D sensor. Similarly, fuel consumption was used instead of energy consumption. Energy consumption can be calculated from fuel consumption and energy density. A fuel flow meter was mounted on the excavator to measure fuel consumption during the excavation process. The total fuel consumption was calculated by integrating the instantaneous flow rate over time because the meter was able to measure the instantaneous rate of fuel flow.

The detailed procedure for the field experiments is as follows.

1. The excavator executes an excavation path using fixed parameters.

2. The bulking factor is calculated from the 3D point-cloud data obtained by the RGB-D sensors.

3. An efficient excavation path is selected based on the soil density estimated from the bulking factor.

4. The target ground is excavated again using the selected excavation path.

5. The total energy consumption during each excavation process and the soil volume accumulated in the bucket are measured, and the excavation efficiency is calculated based on these data.

4.2. Experimental Results on Soil Density Estimation

First, the experimental results under the soft soil conditions are presented. Fig. 11 shows the 3D point-cloud data measured by the RGB-D sensors during excavation with fixed parameters. In this figure, (a) and (b) show the excavation areas before and after excavation, respectively, and (c) and (d) show the dumping areas before and after dumping. Table 5 lists the bulking factor and the estimated soil density calculated from the 3D point cloud data. As shown in Table 5, a reasonable value was obtained for soft soil.

Here, we consider the experimental results for the hard soil. Fig. 12 shows the 3D point-cloud data measured by the RGB-D sensors during excavation with fixed parameters. Table 5 lists the bulking factor and the estimated soil density calculated from the 3D point-cloud data mentioned above. As shown in Table 5, the estimated density exceeded 100%, although it was theoretically defined to range from 0% to 100%. Based on observations from this experiment, we hypothesized that the presence of compacted soil in the target ground caused this issue. Fig. 13 shows the excavated area of the hard soil. As shown in Fig. 13, a significant amount of compacted soil remained in block-like forms. Given that the compacted soil did not break apart even after excavation, large voids remained in the soil. Because the proposed method assumes that there are no large voids in the soil in the volume calculation, the experimental results overestimated the post-excavation soil volume and underestimated the pre-excavation volume. Consequently, the estimated density exceeded 100%.

4.3. Experimental Results on Efficient Excavation

Based on the results of our estimation of soil density, an efficient excavation was conducted on the soft soil where a valid estimation result was obtained.

Table 6 shows the soil volume accumulated in the bucket after excavation (\(V_{\rm{d}}\)), the total energy consumption during the excavation process (\(E\)), and the excavation efficiency (\(V_{\rm{d}}/E\)). As shown in Table 6, the excavation efficiency using the proposed method was 0.212 \(\rm{m^3}/\rm{mL}\), whereas that obtained using fixed parameters was 0.166 \(\rm{m^3}/\rm{mL}\). This indicates that the proposed method achieved 27.7% higher excavation efficiency than the fixed-parameter approach. Specifically, the proposed method significantly reduced the energy consumption while excavating a similar volume of soil.

Figure 14 compares the fuel flow rates during the excavation process using fixed parameters and the proposed method. The gray line represents the fuel flow rate with fixed parameters, and the black line represents the fuel flow rate with the proposed method. As shown in Fig. 14, the proposed method required less energy during the excavation process. This occurred primarily because the proposed method shortened the dragging phase, which corresponds to the significant gap observed between 5,000 ms and 10,000 ms in Fig. 14.