# Research on the Method of Detecting the Spreading Rate of the Simultaneous Crushed Stone Sealing Layer Based on Machine Vision

Authors: Xin Rong; Honghai Liu; Xinmin Gao; Qinghua Bian
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Synchronous chip seal is an advanced road constructing technology, and the gravel coverage rate is an important indicator of the construction quality. The traditional method to measure the gravel coverage rate usually depends on observation by human eyes, which is rough and inefficient. In this paper, a detection method of gravel coverage based on improved wavelet algorithm is proposed. By decomposing the image with two-dimensional discrete wavelet, the high-frequency and low-frequency coefficients are extracted. The noise of the high-frequency coefficients in the image is removed by improving the threshold function, and the contrast of the gravel target in the low-frequency coefficients is improved by the multiscale Retinex algorithm, and then two-dimensional wavelet reconstruction is carried out. Finally, the gravel target is segmented by the block threshold method, and the pixel ratio of the gravel is calculated to complete the detection of the gravel coverage. The experimental results show that the proposed method can effectively segment the gravel target and reduce the influence of environmental factors on the detection accuracy. The detection accuracy error is within ±2%, which can meet the detection requirements. The improved wavelet algorithm improves the signal-to-noise ratio of the denoised image, reduces the mean square error, and achieves a relatively good denoising effect.


---

## Body

## 1. Introduction

Synchronous gravel seal construction is an advanced asphalt pavement construction method in recent years. The basic method is as follows: during construction, the synchronous gravel seal vehicles drive over the construction road surface, synchronously spread asphalt, and gravel on the ground. The asphalt and gravel are bonded on the ground and rolled by tires or steel wheel rollers to form a single layer of asphalt and gravel wear layer. Compared with other methods, this construction method has extremely high construction efficiency, and it can reduce the traffic control of the construction section and complete the construction on the road where vehicles pass. The gravel coverage is the percentage of the projected area covered by the gravel to the entire pavement and is an important construction index in the construction of the synchronous gravel seal. In the synchronous gravel seal construction, improper gravel spreading rate or uneven gravel spreading will cause disease (too many gravels will cause threshing, and insufficient gravel will cause oil panning). Generally speaking, the upper seal layer of crushed stones should be in the size of 9.5-13.2 mm, the distribution rate of crushed stones should be 80%∼90%, the lower seal layer should be in the size of 9.5 mm-16 mm (highway) or 4.75-9.5 mm (other Grade), the distribution rate of crushed stone should be 60%∼80%, the stress absorbing layer should use 9.5 mm∼16 mm gauge, and the distribution rate of crushed stone should be 60%∼80% [1]. At present, in the actual construction process, only the construction personnel can use the method of visual inspection to estimate the rubble distribution rate, which is subjective and arbitrary. Therefore, it has a great significance for the machine vision-based method for detecting the distribution rate of the crushed stone in the synchronous gravel seal.Many scholars have conducted research on the use of machine vision to solve problems in production and life. Chowdhury et al. of Texas A&M University in the United States used Hough transform and other image analysis methods to quantify the angularity of fine aggregate particles and distinguish the quality of the aggregate in the hot-mix asphalt mixture [2]. Song et al. used digital image processing to evaluate the uniformity of the asphalt pavement surface texture structure distribution [3]. Song and Wang studied the application of nuclear-free densitometer and digital image segregation evaluation method in asphalt pavement construction quality monitoring and conducted correlation analysis of the two methods [4]. Wang et al. studied the detection method of the coverage of the rubble spreading of the synchronous gravel seal layer. First, the Retinex algorithm was used to restore the image, remove the shadows, and enhance the image, and then the image binarization method was used to obtain the rough spreading coverage. The method fails to treat particulate matter as an independent object, so it is difficult to obtain accurate results [5]. Pan and Tutumluer and others in the United States proposed a three-view-based method for measuring the surface area of irregular aggregate particles and compared the results of the method with the results of the three-dimensional laser scanner to verify the accuracy of the method. The method can also be used for the calculation of asphalt film thickness [6]. Browne et al. designed and realized the ability to dynamically collect images of aggregate particles scattered in a light-assisted background environment and analyze the size characteristics of the particles [7]. Abdullah et al. of Jordan University of Science and Technology used digital image analysis method to quantify the porosity of mineral material in asphalt mixture and compared it with traditional detection methods [8]. In addition, AMIS developed by Texas A&M University, UIAIA developed by the University of Illinois at Urbana-Champaign, and WipShape developed by the University of Missouri all use image analysis technology to achieve automated grading analysis [9, 10]. Tajeripour and Fekri-Ershad proposed detecting abnormalities in stone textures based on one-dimensional local binary patterns, and the proposed approach is fully automatic and all of the necessary parameters can be tuned [11]. Al-Utaibi et al. and Basheera M. Mahmmod used Krawtchouk and Hahn polynomials to reconstruct images and analyzed the influence of different parameters on the quality of reconstructed images to complete the detection of the targets [12, 13]. Currently, widely used image segmentation algorithms include threshold segmentation, watershed algorithm, cluster segmentation, genetic algorithm, and so on. These algorithms can measure the target and background segmentation well. However, after studying the impact on the accuracy of rubble spreading detection, it is found that environmental factors such as light and noise in actual construction often cause the final treatment result to fail to meet the requirements.Aiming at the abovementioned problems, this paper proposes a detection method of gravel coverage based on improved wavelet algorithm. By decomposing the image with two-dimensional discrete wavelet, the high-frequency and low-frequency coefficients are extracted, the high-frequency coefficients are removed by improving the threshold function, and the low-frequency coefficients are improved by the multiscale Retinex algorithm to improve the contrast of the gravel target, and then the two-dimensional wavelet reconstruction is carried out. Finally, the crushed stone target is segmented by the block threshold method, the pixel proportion of the crushed stone is calculated, the detection of the crushed stone coverage is completed, and the relationship between the spreading rate and the spreading amount is established. In the laboratory, the spreading conditions of 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm crushed stones under different spreading rates were tested and verified in actual road engineering.


## 2. Experimental Programme

### 2.1. Image Segmentation Based on Machine Vision

In this study, an industrial camera was used to capture images of scattered gravel. The picture is processed by MATLAB software. According to the result of segmentation, the distribution rate of crushed stone of the synchronous gravel sealer can be obtained, and then the relationship between the distribution rate and the amount of distribution can be established, so as to adjust the subsequent distribution section to ensure the quality of the distribution. First, the impact of crushed stone specifications on the accuracy of spreading rate detection was tested in the laboratory, and secondly, the relationship between the spreading rate of different specifications of crushed stone and the amount of spreading was established. The laboratory test equipment used for testing is shown in Figure1. In this experiment, the crushed stone type used is AC-16, and the asphalt type is SBS emulsified asphalt. The flow chart of the laboratory test is shown in Figure 2.Figure 1
Laboratory test device for the segmentation.Figure 2
Flow chart of the laboratory test.In the process of detecting the distribution rate of the synchronous gravel seal layer, after preprocessing operations such as distortion correction and image enhancement, the coverage of the gravel on the scattered road section can be detected. This paper adopts the method of image binary segmentation to detect the spreading rate of the spreading road. However, in the process of image acquisition, the target and the background cannot be completely separated due to environmental factors such as light and noise, resulting in wrong segmentation, and the results cannot meet the requirements. Therefore, in response to this problem, this paper uses discrete wavelet to extract the high-frequency and low-frequency coefficients of the image for processing to eliminate the influence of environmental factors such as light and noise, and then perform the binary segmentation processing on the wavelet-reconstructed image. Carry out the detection of the rubble spreading rate and establish the relationship between the spreading rate and the amount of spreading.Take the method of image binary segmentation to calculate the spreading rate. Since the stones and the uncovered asphalt have significant color differences, the original image can be segmented using a certain segmentation algorithm to segment the gravel from the background (asphalt) [14], which is as follows:(1)gx,y=1fx,y∈gravel0fx,y∈asphalt,,Then, the spread rate can be calculated:(2)r=∑x∑ygx,y∑x∑y1.It can be seen from the above two formulas that the accuracy of the detection method is determined by the accuracy of image segmentation.The specific steps are as follows:Step 1.
Select the right industrial camera and lens so that the picture accuracy can reach 0.1 mm, that is, the 1.18 mm size gravel can occupy at least 10 pixels. After comparison, the CCD camera model MV-EM200 M and the lens of BT-23C1214MP5 were finally selected.Step 2.
Adjust the camera distance according to the camera focus position. LED lights should be installed in front of the test equipment to reduce the impact of uneven natural light on shooting. In this experiment, the focal length of the camera is 1.4 mm to 16 mm, the frame rate is 20 fps, the area of the imaging part is 1/1.8'', and the number of effective pixels is 1920000.Step 3.
Simulate the spreading process of the synchronous gravel sealer in the laboratory. Simulate the spreading scenes of gravel in different specifications and spreading rates. The camera can sample the road surface after spreading. The picture can be processed by MATLAB software.(1)
Use the imread function to read the selected picture.(2)
Use the imcrop function to crop the image and remove the unrelated parts.(3)
Convert the color image into a grayscale image by the rgb2gray function.(4)
Use discrete wavelet algorithm to decompose the gray image in two dimensions.(5)
Extract discrete wavelet high-frequency and low-frequency coefficients.(6)
Use improved threshold function to process high-frequency coefficients.(7)
Use improved Retinex algorithm to process low-frequency coefficients.(8)
Perform wavelet reconstruction on the image.(9)
Use block Otsu threshold segmentation to convert grayscale images into binary images.(10)
Use hole filling to optimize the morphology of the binary image.(11)
Detect the spread rate of gravel.(12)
Establish the relationship between spreading rate and spreading amount.Step 4.
Calculate the actual distribution rate of gravel using existing methods [15].(1)
Sprinkle a layer of 1∼2 mm thick cementing material evenly on the enamel plate or oil felt with a known area, and weigh the combined weight of the oil felt and the cementm1.(2)
Heat the clean aggregate (single particle size) to be spread to 120-140°C.(3)
Spread the heated aggregate on an enamel pan or linoleum felt.(4)
After the aggregate and cement are cooled to room temperature, weigh the total weight of the enamel pan or felt, cement, and aggregatem2.(5)
Calculate the spreading rate of full-paved aggregate and design aggregate spreading rate.(3)Pm=m2−m1S1,Pd=Pm×r,where Pm is the spreading amount of aggregate in full spread, kg/m2; Pd is the designed amount of aggregate, kg/m2; m1 is the combined weight of enamel plate or oil felt and cementing material, kg; m2 is the combined weight of enamel plate or linoleum felt, cement, and aggregate, kg; r is the aggregate design coverage rate, %; and S1 is the enamel plate or oil felt area, m2.
Compare the actual spread rate with the detected spread rate, analyze the error rate of the algorithm and the cause of the error, and make improvements.Step 5.
Calculate the actual amount of spreading using existing methods.(1)
Place the enamel tray or felt with area S1 on the place where the distributor passes.(2)
Take out the enamel tray or linoleum felt immediately after the car with distributor passes through.(3)
Shovel off the aggregate and part of the cement with a scraper, and it is better to shovel off all the aggregates on the porcelain plate or linoleum.(4)
The weight of aggregatem1 is obtained by extraction method or trichloroethylene soaking and washing method.(5)
Calculate the aggregate spreading amount:(4)Ps=m1S1,where Ps is the amount of aggregate spreading, kg/m2; and m1 is the weight of aggregate, kg.
Compare the detected value with the actual value, analyze the algorithm error rate and the cause of the error, and make improvements.Step 6.
Repeat steps 3-8 to obtain the detected and actual values of the spreading rate and spreading amount of crushed stones of different specifications.Step 7.
After the detection algorithm is corrected, the actual gravel seal pavement is tested.
For laboratory tests, according to actual working conditions, aggregates of 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm are divided. The 0-1.18 mm aggregate has a small particle size, and a camera with a higher resolution should be used to shoot clearly. Therefore, the sieving of the 0-1.18 mm aggregate is not considered in this study.


### 2.2. Noise Reduction Algorithm Based on Discrete Wavelet

#### 2.2.1. Wavelet Decomposition and Reconstruction

The wavelet decomposition method is used to obtain low-frequency coefficients and high-frequency coefficients. The low-frequency coefficients mainly include the global information of the image, including the contour information of the image, and the high-frequency coefficients mainly include the local information of the image, including the edges, details, and noise of the image. The expression of two-dimensional discrete wavelet decomposition is as follows [16]:(5)Wφj0,m,n=1/MN∑x=0M−1∑y=0N−1fx,yφj0,m,nx,yWψij,m,n=1/MN∑x=0M−1∑y=0N−1fx,yψj,m,nix,yi=H,V,D,where Wφj0,m,n is the low-frequency coefficient after decomposition; Wψij,m,n is the high-frequency coefficient after decomposition, means H, V, and D; j0 is any starting scale, usually let it be 0; f(x,y) is the discrete function, x and y represent discrete variables; m and n are expressed as relative offsets; M×N represents the image is composed of pixels; φj0,m,nx,y is a two-dimensional scaling function, and H,V,D is the wavelet function corresponding to i in the horizontal, vertical, and diagonal directions.Process the low-frequency coefficientWφj0,m,n and high-frequency coefficient Wψij,m,n obtained after decomposition. The reconstructed image is obtained through the inverse discrete wavelet transform, and the wavelet reconstruction expression is as follows. The schematic diagram of the two-layer wavelet decomposition is shown in Figure 3 [17].(6)fx,y=1MN∑m∑nWφj0,m,nφj0,m,nx,y+1MN∑i=H,V,D∑j=j0∞∑m∑nWψij,m,nψj,m,nix,y,where S represents the decomposed image, cA1 and cA2, respectively, represent the low-frequency components of the two-layer wavelet decomposition, that is, low-frequency images, and cD1 and cD2 represent the high-frequency components of the two-level wavelet decomposition, that is, the horizontal (H), the vertical (V), and the diagonal (D) components of each layer. The core processing flow is as follows:(1)
Image wavelet decomposition: including selecting a certain wavelet base function and determining the optimal wavelet decomposition level N, perform N-layer wavelet decomposition on the noisy imagef(x, y) to obtain high-frequency wavelet coefficients ωj,k and low-frequency wavelet image S(x, y).(2)
Threshold quantization of high-frequency coefficients of wavelet decomposition: according to a threshold quantization criterion, an appropriate threshold is determined for each high-frequency coefficient of the first to Nth layers. By thresholdingωj,k, the estimated wavelet coefficients ω∧j,k are obtained, and ω∧j,k−ωj,k is as small as possible.(3)
Retinex enhancement of the low-frequency coefficients of wavelet decomposition: according to the principle of the Retinex algorithm, the low-frequency image is converted to the logarithmic domain to obtain the reflection image, and an illumination adjustment parameterk is added in the process of eliminating the illumination image, so that the obtained reflection image is more natural. After Gaussian filtering, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(4)
Wavelet reconstruction: inverse wavelet transform is performed on the high-frequency coefficientsωj,k of the 1st to Nth layers after threshold quantization and the low-frequency wavelet coefficients of the Nth layer to obtain the estimated image.Figure 3
Schematic diagram of wavelet decomposition.


#### 2.2.2. Improved Threshold Function Wavelet Denoising Method

(1) Selection of Wavelet Threshold. An improved threshold function method is used to denoise the high-frequency coefficients after wavelet decomposition. In order to effectively separate the noise components in the high-frequency coefficients, the function threshold must first be determined. When the threshold function is too small, the denoising effect of the image is not obvious, and the denoising result image is not much different from the original image, which loses the meaning of threshold denoising. When the threshold function is too large, the degree of image noise removal is too large. Although the noise is removed more thoroughly, the image is too smooth and the image is blurred, which cannot reflect the essential information of the image. Therefore, choosing a suitable threshold is of great significance to the results of threshold denoising [18].The fixed threshold estimation method has a more thorough denoising effect when the noise is more distributed at high-frequency coefficients. Therefore, the fixed threshold estimation method is used to estimate the function threshold, and the expression is as follows [19]:(7)λ=σm2logM⋅Nσm=medianabsk/0.6745,where λ is the determined fixed threshold; and σm is the standard deviation of noise.(2) Selection of Threshold Function. The traditional soft threshold (ST) and hard threshold (HT) function methods are less computationally intensive and easy to implement, but the method will produce constant errors when the high-frequency coefficient is greater than the fixed threshold, resulting in blurred edge details. Through the analysis of ST function in Figure 4, we can see that the continuity of ST function at ±λ is better, but when ωj,k>λ, the estimated wavelet coefficient ωj,k∧ obtained by the ST function has a constant deviation λ from the original threshold coefficient ωj,k. The constant deviation will cause the loss of high-frequency information and cause distortion. From the analysis of the HT function in Figure 5, it can be seen that the HT method solves the constant error problem of the ST method, but there is a discontinuity point ±λ. The HT function is not continuous at the discontinuity point ±λ, and it is easy to cause oscillation when reconstructing the image.(1)
Soft threshold function is as follows:(8)ωj,k∧=sgnωj,k⋅ωj,k−λ0ωj,k>λωj,k≤λ.(2)
Hard threshold function is as follows:(9)ωj,k∧=ωj,k0ωj,k>λωj,k≤λ.(3)
This paper improves the threshold function:(10)ωj,k∧=sgnωj,k⋅ωj,k−λeωj,k−λ0ωj,k≥λωj,k<λ.Figure 4
Soft threshold function.Figure 5
Hard threshold function.From the analysis of equation (13), we can see that in terms of function continuity, when the wavelet coefficients tend to a fixed threshold, the estimated wavelet coefficient limit is 0, so the function value is continuous at the fixed threshold, avoiding the shock phenomenon when the hard threshold function has a discontinuity point. In terms of function asymptoticity, when the wavelet coefficients tend to positive or negative infinity, the limit of the ratio of the estimated wavelet coefficients to the original wavelet coefficients is 1, which proves that the improved threshold function takes ωj,k∧=ωj,k as the asymptote; thus, eliminating the constant error of the ST function problem. In terms of function deviation, the difference between the estimated wavelet coefficient and the original wavelet coefficient is calculated with the wavelet coefficient tending to plus or minus infinity as the limit, and the difference is 0, so when the wavelet coefficient tends to infinity, the estimated new wavelet coefficient is infinite closed to the original wavelet coefficients to avoid the problem of deviation.


#### 2.2.3. Image Enhancement Algorithm Based on Improved Retinex

An improved Retinex image enhancement algorithm is used to estimate the illumination image of the low-frequency coefficients after wavelet decomposition, and the reflection image is obtained through logarithmic domain conversion. The improved Retinex algorithm in this paper adds a light adjustment parameter based on the single-scale filter Retinex (SSR) algorithm, which makes the reflected image more natural, avoids the halo phenomenon in effect, and has a better edge retention effect.(1) Illumination Image Estimation. According to the principle of Retinex algorithm, the illumination image is the convolution of the original image and the Gaussian kernel function. In order to obtain the reflection image and eliminate the influence of the illumination on the image, the function is converted to the logarithmic domain, and then the reflection image after the logarithmic domain conversion is obtained. The number field conversion expression is as follows [20]:(11)Rx,y=logaSx,y−logaSx,y∗Fx,y.However, through experiments on the SSR algorithm, it is found that if the illumination estimation part in the original image is completely eliminated, the image will often be unnatural. Therefore, an illumination adjustment parameter is added in the process of eliminating the illumination image to make the reflected image more natural. The expression of the reflection image after introducing the adjustment factor is as follows:(12)Rx,y=logaSx,y−klogaSx,y∗Fx,y,where k is the light adjustment parameter; F(x,y) is the Gaussian kernel function; Rx,i is the illumination image; and a is the number of image channels, usually 3.After Gaussian filtering is performed, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(13)Rx,y=∑k=1NωklogaSx,y−klogaSx,y∗Fkx,y,where ωk is the weight coefficient.


#### 2.2.4. Reflection Image Stretch

The reflection image obtained after the image is transformed in the logarithmic domain often cannot reach the range of 0∼255 in the image value domain, which may cause the image to be gray and white, which affects the visual effect and subsequent processing. Therefore, the reflection image needs to be stretched. Firstly, by traversing the pixel values, find the maximum and minimum pixel values in the reflection image, and then stretch the image by the following formula to stretch the reflection image to the range of [0,255] [21].(14)dx,y=Rx,y−RminRmax−Rmin×255−0,where d(x,y) represents the reflected image after stretching; R(x,y) represents the reflected image before stretching; and Rmax and Rmin represent the maximum and minimum gray values of the reflected image before stretching, respectively.


#### 2.2.5. Piecewise Linear Transformation

Aiming at the pan-gray phenomenon in the reconstructed image, this paper uses a three-segment piecewise linear transformation method to improve the image contrast. The three-stage piecewise linear transformation expression is as follows:(15)fi,j=k1×di,j0≤di,j<ab+k2×di,ja≤di,j≤cd+k3×di,jc<di,j≤255,where f(i,j) is the output image after contrast enhancement; d(i,j) is the input image; and k1, k2, and k3 represent the slope of the three-stage transformation. and the expression is as follows [22]:(16)k1=ab,k2=d−bc−a,k3=255−d255−c,where (a,b) and (c,d) indicate the point where the slope of the piecewise function changes.The segmentation function can set different segmentation points according to the image characteristics, enhance the details of the region of interest, and can also suppress the region of interest.


#### 2.2.6. Detection of Spreading Rate and Spreading Amount Based on Segmentation Threshold Segmentation Method

By dividing the image into several blocks and performing threshold segmentation, the uneven effect caused by illumination or reflection can be solved to a certain extent. The blocks are chosen to be small enough so that the lighting of each block is approximately uniform, so that when auto-thresholding, high-threshold segmentation is used in high-gray areas and low-threshold segmentation is used in low-gray areas. Since the size of a single aggregate is too small, if a single aggregate is divided into blocks, the amount of calculation will be greatly increased. Therefore, in this paper, the size of 5 aggregates is selected to be divided inton × n small images of equal length and width. After the blocks are divided, the global threshold method is performed according to the block, and the threshold T of each image block is calculated. By calculating the average gray difference Δd between classes of each image block, the interference of pure background or pure object is excluded, and combined with Otsu algorithm, the image block is binarized based on this condition. Then, the binarized images of each image block are stitched together to form a complete crushed stone binarized image.(17)Δd=m1k−m2k,m1k=∑i=1kiPiC1,m2k=∑i=k+1LiPiC2,where C1 is the pixels whose gray level is [1,2, ⋯, k]; C2 is the pixels whose gray level is [k+1, k+2, ⋯, L]; Pi is the frequency distribution of gray levels; m(k) is the average gray value of gray levels from 1 to k; and Δd is the average grayscale difference between classes.The average grayscale difference Δd between classes calculated according to the formula is combined with the threshold value T of each image block to realize the binarization operation on each image block. When the interclass grayscale difference Δd of the image block is less than the set value d0, it is determined that the pixels in the image block are of the same class, and the image block is not binarized; otherwise, the image block is binarized. The detection method in this paper divides the image into 4 × 4 image blocks (i.e., n = 4), the accuracy of the detection result is 98.1%, and the average processing speed of a single image is less than 85 ms. In this paper, comparisons of single aggregate, 5 aggregates, and 10 aggregates as a block size are listed in Table 1.Table 1
Comparison of single aggregate, 5 aggregates, and 10 aggregates.
Aggregate quantityAccuracy (%)OH (ms)1 aggregate98.21095 aggregates98.18510 aggregates96.283It can be seen from Table1 that when a single aggregate and 5 aggregates are selected as the size of the block, the recognition accuracy is basically the same, but the operation speed of a single aggregate is much slower than that of 5 aggregates. Too many blocks will slow down the recognition speed; the operation speed of 5 aggregates and 10 aggregates is roughly the same, but the accuracy of 5 aggregates is better than 10 aggregates, which is because the block is too large, the thresholds of the targets in the block are sometimes quite different, and the entire target cannot be accurately and completely segmented, resulting in a decrease in the accuracy rate. Therefore, 5 aggregates are selected as one block in this paper.Under a fixed shooting area, after obtaining a certain specification of the crushed stone distribution rate and the weight of the used crushed stone, the crushed stone weight and distribution amount of a certain specification of crushed stone in this area can be obtained under the condition of spreading. Furthermore, it is possible to detect the weight and the amount of crushed stones of different specifications of crushed stones under different spreading rates.(18)mi+2=mi−1/ri−1+mi/ri+mi+1/ri+1300×ri+2,(19)Ps=mi+2s,where ri-1, ri, and ri+1 are the known distributing rates of crushed stones at different moments; mi-1, mi, and mi+1 are the known crushed rock weights at different distributing rates of crushed stones; ri+2 is the distributing rate of crushed stones at any time (≤100%); mi+2 is the weight of the gravel at the spreading rate ri+2; and s is the shooting area.


## 2.1. Image Segmentation Based on Machine Vision

In this study, an industrial camera was used to capture images of scattered gravel. The picture is processed by MATLAB software. According to the result of segmentation, the distribution rate of crushed stone of the synchronous gravel sealer can be obtained, and then the relationship between the distribution rate and the amount of distribution can be established, so as to adjust the subsequent distribution section to ensure the quality of the distribution. First, the impact of crushed stone specifications on the accuracy of spreading rate detection was tested in the laboratory, and secondly, the relationship between the spreading rate of different specifications of crushed stone and the amount of spreading was established. The laboratory test equipment used for testing is shown in Figure1. In this experiment, the crushed stone type used is AC-16, and the asphalt type is SBS emulsified asphalt. The flow chart of the laboratory test is shown in Figure 2.Figure 1
Laboratory test device for the segmentation.Figure 2
Flow chart of the laboratory test.In the process of detecting the distribution rate of the synchronous gravel seal layer, after preprocessing operations such as distortion correction and image enhancement, the coverage of the gravel on the scattered road section can be detected. This paper adopts the method of image binary segmentation to detect the spreading rate of the spreading road. However, in the process of image acquisition, the target and the background cannot be completely separated due to environmental factors such as light and noise, resulting in wrong segmentation, and the results cannot meet the requirements. Therefore, in response to this problem, this paper uses discrete wavelet to extract the high-frequency and low-frequency coefficients of the image for processing to eliminate the influence of environmental factors such as light and noise, and then perform the binary segmentation processing on the wavelet-reconstructed image. Carry out the detection of the rubble spreading rate and establish the relationship between the spreading rate and the amount of spreading.Take the method of image binary segmentation to calculate the spreading rate. Since the stones and the uncovered asphalt have significant color differences, the original image can be segmented using a certain segmentation algorithm to segment the gravel from the background (asphalt) [14], which is as follows:(1)gx,y=1fx,y∈gravel0fx,y∈asphalt,,Then, the spread rate can be calculated:(2)r=∑x∑ygx,y∑x∑y1.It can be seen from the above two formulas that the accuracy of the detection method is determined by the accuracy of image segmentation.The specific steps are as follows:Step 1.
Select the right industrial camera and lens so that the picture accuracy can reach 0.1 mm, that is, the 1.18 mm size gravel can occupy at least 10 pixels. After comparison, the CCD camera model MV-EM200 M and the lens of BT-23C1214MP5 were finally selected.Step 2.
Adjust the camera distance according to the camera focus position. LED lights should be installed in front of the test equipment to reduce the impact of uneven natural light on shooting. In this experiment, the focal length of the camera is 1.4 mm to 16 mm, the frame rate is 20 fps, the area of the imaging part is 1/1.8'', and the number of effective pixels is 1920000.Step 3.
Simulate the spreading process of the synchronous gravel sealer in the laboratory. Simulate the spreading scenes of gravel in different specifications and spreading rates. The camera can sample the road surface after spreading. The picture can be processed by MATLAB software.(1)
Use the imread function to read the selected picture.(2)
Use the imcrop function to crop the image and remove the unrelated parts.(3)
Convert the color image into a grayscale image by the rgb2gray function.(4)
Use discrete wavelet algorithm to decompose the gray image in two dimensions.(5)
Extract discrete wavelet high-frequency and low-frequency coefficients.(6)
Use improved threshold function to process high-frequency coefficients.(7)
Use improved Retinex algorithm to process low-frequency coefficients.(8)
Perform wavelet reconstruction on the image.(9)
Use block Otsu threshold segmentation to convert grayscale images into binary images.(10)
Use hole filling to optimize the morphology of the binary image.(11)
Detect the spread rate of gravel.(12)
Establish the relationship between spreading rate and spreading amount.Step 4.
Calculate the actual distribution rate of gravel using existing methods [15].(1)
Sprinkle a layer of 1∼2 mm thick cementing material evenly on the enamel plate or oil felt with a known area, and weigh the combined weight of the oil felt and the cementm1.(2)
Heat the clean aggregate (single particle size) to be spread to 120-140°C.(3)
Spread the heated aggregate on an enamel pan or linoleum felt.(4)
After the aggregate and cement are cooled to room temperature, weigh the total weight of the enamel pan or felt, cement, and aggregatem2.(5)
Calculate the spreading rate of full-paved aggregate and design aggregate spreading rate.(3)Pm=m2−m1S1,Pd=Pm×r,where Pm is the spreading amount of aggregate in full spread, kg/m2; Pd is the designed amount of aggregate, kg/m2; m1 is the combined weight of enamel plate or oil felt and cementing material, kg; m2 is the combined weight of enamel plate or linoleum felt, cement, and aggregate, kg; r is the aggregate design coverage rate, %; and S1 is the enamel plate or oil felt area, m2.
Compare the actual spread rate with the detected spread rate, analyze the error rate of the algorithm and the cause of the error, and make improvements.Step 5.
Calculate the actual amount of spreading using existing methods.(1)
Place the enamel tray or felt with area S1 on the place where the distributor passes.(2)
Take out the enamel tray or linoleum felt immediately after the car with distributor passes through.(3)
Shovel off the aggregate and part of the cement with a scraper, and it is better to shovel off all the aggregates on the porcelain plate or linoleum.(4)
The weight of aggregatem1 is obtained by extraction method or trichloroethylene soaking and washing method.(5)
Calculate the aggregate spreading amount:(4)Ps=m1S1,where Ps is the amount of aggregate spreading, kg/m2; and m1 is the weight of aggregate, kg.
Compare the detected value with the actual value, analyze the algorithm error rate and the cause of the error, and make improvements.Step 6.
Repeat steps 3-8 to obtain the detected and actual values of the spreading rate and spreading amount of crushed stones of different specifications.Step 7.
After the detection algorithm is corrected, the actual gravel seal pavement is tested.
For laboratory tests, according to actual working conditions, aggregates of 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm are divided. The 0-1.18 mm aggregate has a small particle size, and a camera with a higher resolution should be used to shoot clearly. Therefore, the sieving of the 0-1.18 mm aggregate is not considered in this study.


## 2.2. Noise Reduction Algorithm Based on Discrete Wavelet

### 2.2.1. Wavelet Decomposition and Reconstruction

The wavelet decomposition method is used to obtain low-frequency coefficients and high-frequency coefficients. The low-frequency coefficients mainly include the global information of the image, including the contour information of the image, and the high-frequency coefficients mainly include the local information of the image, including the edges, details, and noise of the image. The expression of two-dimensional discrete wavelet decomposition is as follows [16]:(5)Wφj0,m,n=1/MN∑x=0M−1∑y=0N−1fx,yφj0,m,nx,yWψij,m,n=1/MN∑x=0M−1∑y=0N−1fx,yψj,m,nix,yi=H,V,D,where Wφj0,m,n is the low-frequency coefficient after decomposition; Wψij,m,n is the high-frequency coefficient after decomposition, means H, V, and D; j0 is any starting scale, usually let it be 0; f(x,y) is the discrete function, x and y represent discrete variables; m and n are expressed as relative offsets; M×N represents the image is composed of pixels; φj0,m,nx,y is a two-dimensional scaling function, and H,V,D is the wavelet function corresponding to i in the horizontal, vertical, and diagonal directions.Process the low-frequency coefficientWφj0,m,n and high-frequency coefficient Wψij,m,n obtained after decomposition. The reconstructed image is obtained through the inverse discrete wavelet transform, and the wavelet reconstruction expression is as follows. The schematic diagram of the two-layer wavelet decomposition is shown in Figure 3 [17].(6)fx,y=1MN∑m∑nWφj0,m,nφj0,m,nx,y+1MN∑i=H,V,D∑j=j0∞∑m∑nWψij,m,nψj,m,nix,y,where S represents the decomposed image, cA1 and cA2, respectively, represent the low-frequency components of the two-layer wavelet decomposition, that is, low-frequency images, and cD1 and cD2 represent the high-frequency components of the two-level wavelet decomposition, that is, the horizontal (H), the vertical (V), and the diagonal (D) components of each layer. The core processing flow is as follows:(1)
Image wavelet decomposition: including selecting a certain wavelet base function and determining the optimal wavelet decomposition level N, perform N-layer wavelet decomposition on the noisy imagef(x, y) to obtain high-frequency wavelet coefficients ωj,k and low-frequency wavelet image S(x, y).(2)
Threshold quantization of high-frequency coefficients of wavelet decomposition: according to a threshold quantization criterion, an appropriate threshold is determined for each high-frequency coefficient of the first to Nth layers. By thresholdingωj,k, the estimated wavelet coefficients ω∧j,k are obtained, and ω∧j,k−ωj,k is as small as possible.(3)
Retinex enhancement of the low-frequency coefficients of wavelet decomposition: according to the principle of the Retinex algorithm, the low-frequency image is converted to the logarithmic domain to obtain the reflection image, and an illumination adjustment parameterk is added in the process of eliminating the illumination image, so that the obtained reflection image is more natural. After Gaussian filtering, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(4)
Wavelet reconstruction: inverse wavelet transform is performed on the high-frequency coefficientsωj,k of the 1st to Nth layers after threshold quantization and the low-frequency wavelet coefficients of the Nth layer to obtain the estimated image.Figure 3
Schematic diagram of wavelet decomposition.


### 2.2.2. Improved Threshold Function Wavelet Denoising Method

(1) Selection of Wavelet Threshold. An improved threshold function method is used to denoise the high-frequency coefficients after wavelet decomposition. In order to effectively separate the noise components in the high-frequency coefficients, the function threshold must first be determined. When the threshold function is too small, the denoising effect of the image is not obvious, and the denoising result image is not much different from the original image, which loses the meaning of threshold denoising. When the threshold function is too large, the degree of image noise removal is too large. Although the noise is removed more thoroughly, the image is too smooth and the image is blurred, which cannot reflect the essential information of the image. Therefore, choosing a suitable threshold is of great significance to the results of threshold denoising [18].The fixed threshold estimation method has a more thorough denoising effect when the noise is more distributed at high-frequency coefficients. Therefore, the fixed threshold estimation method is used to estimate the function threshold, and the expression is as follows [19]:(7)λ=σm2logM⋅Nσm=medianabsk/0.6745,where λ is the determined fixed threshold; and σm is the standard deviation of noise.(2) Selection of Threshold Function. The traditional soft threshold (ST) and hard threshold (HT) function methods are less computationally intensive and easy to implement, but the method will produce constant errors when the high-frequency coefficient is greater than the fixed threshold, resulting in blurred edge details. Through the analysis of ST function in Figure 4, we can see that the continuity of ST function at ±λ is better, but when ωj,k>λ, the estimated wavelet coefficient ωj,k∧ obtained by the ST function has a constant deviation λ from the original threshold coefficient ωj,k. The constant deviation will cause the loss of high-frequency information and cause distortion. From the analysis of the HT function in Figure 5, it can be seen that the HT method solves the constant error problem of the ST method, but there is a discontinuity point ±λ. The HT function is not continuous at the discontinuity point ±λ, and it is easy to cause oscillation when reconstructing the image.(1)
Soft threshold function is as follows:(8)ωj,k∧=sgnωj,k⋅ωj,k−λ0ωj,k>λωj,k≤λ.(2)
Hard threshold function is as follows:(9)ωj,k∧=ωj,k0ωj,k>λωj,k≤λ.(3)
This paper improves the threshold function:(10)ωj,k∧=sgnωj,k⋅ωj,k−λeωj,k−λ0ωj,k≥λωj,k<λ.Figure 4
Soft threshold function.Figure 5
Hard threshold function.From the analysis of equation (13), we can see that in terms of function continuity, when the wavelet coefficients tend to a fixed threshold, the estimated wavelet coefficient limit is 0, so the function value is continuous at the fixed threshold, avoiding the shock phenomenon when the hard threshold function has a discontinuity point. In terms of function asymptoticity, when the wavelet coefficients tend to positive or negative infinity, the limit of the ratio of the estimated wavelet coefficients to the original wavelet coefficients is 1, which proves that the improved threshold function takes ωj,k∧=ωj,k as the asymptote; thus, eliminating the constant error of the ST function problem. In terms of function deviation, the difference between the estimated wavelet coefficient and the original wavelet coefficient is calculated with the wavelet coefficient tending to plus or minus infinity as the limit, and the difference is 0, so when the wavelet coefficient tends to infinity, the estimated new wavelet coefficient is infinite closed to the original wavelet coefficients to avoid the problem of deviation.


### 2.2.3. Image Enhancement Algorithm Based on Improved Retinex

An improved Retinex image enhancement algorithm is used to estimate the illumination image of the low-frequency coefficients after wavelet decomposition, and the reflection image is obtained through logarithmic domain conversion. The improved Retinex algorithm in this paper adds a light adjustment parameter based on the single-scale filter Retinex (SSR) algorithm, which makes the reflected image more natural, avoids the halo phenomenon in effect, and has a better edge retention effect.(1) Illumination Image Estimation. According to the principle of Retinex algorithm, the illumination image is the convolution of the original image and the Gaussian kernel function. In order to obtain the reflection image and eliminate the influence of the illumination on the image, the function is converted to the logarithmic domain, and then the reflection image after the logarithmic domain conversion is obtained. The number field conversion expression is as follows [20]:(11)Rx,y=logaSx,y−logaSx,y∗Fx,y.However, through experiments on the SSR algorithm, it is found that if the illumination estimation part in the original image is completely eliminated, the image will often be unnatural. Therefore, an illumination adjustment parameter is added in the process of eliminating the illumination image to make the reflected image more natural. The expression of the reflection image after introducing the adjustment factor is as follows:(12)Rx,y=logaSx,y−klogaSx,y∗Fx,y,where k is the light adjustment parameter; F(x,y) is the Gaussian kernel function; Rx,i is the illumination image; and a is the number of image channels, usually 3.After Gaussian filtering is performed, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(13)Rx,y=∑k=1NωklogaSx,y−klogaSx,y∗Fkx,y,where ωk is the weight coefficient.


### 2.2.4. Reflection Image Stretch

The reflection image obtained after the image is transformed in the logarithmic domain often cannot reach the range of 0∼255 in the image value domain, which may cause the image to be gray and white, which affects the visual effect and subsequent processing. Therefore, the reflection image needs to be stretched. Firstly, by traversing the pixel values, find the maximum and minimum pixel values in the reflection image, and then stretch the image by the following formula to stretch the reflection image to the range of [0,255] [21].(14)dx,y=Rx,y−RminRmax−Rmin×255−0,where d(x,y) represents the reflected image after stretching; R(x,y) represents the reflected image before stretching; and Rmax and Rmin represent the maximum and minimum gray values of the reflected image before stretching, respectively.


### 2.2.5. Piecewise Linear Transformation

Aiming at the pan-gray phenomenon in the reconstructed image, this paper uses a three-segment piecewise linear transformation method to improve the image contrast. The three-stage piecewise linear transformation expression is as follows:(15)fi,j=k1×di,j0≤di,j<ab+k2×di,ja≤di,j≤cd+k3×di,jc<di,j≤255,where f(i,j) is the output image after contrast enhancement; d(i,j) is the input image; and k1, k2, and k3 represent the slope of the three-stage transformation. and the expression is as follows [22]:(16)k1=ab,k2=d−bc−a,k3=255−d255−c,where (a,b) and (c,d) indicate the point where the slope of the piecewise function changes.The segmentation function can set different segmentation points according to the image characteristics, enhance the details of the region of interest, and can also suppress the region of interest.


### 2.2.6. Detection of Spreading Rate and Spreading Amount Based on Segmentation Threshold Segmentation Method

By dividing the image into several blocks and performing threshold segmentation, the uneven effect caused by illumination or reflection can be solved to a certain extent. The blocks are chosen to be small enough so that the lighting of each block is approximately uniform, so that when auto-thresholding, high-threshold segmentation is used in high-gray areas and low-threshold segmentation is used in low-gray areas. Since the size of a single aggregate is too small, if a single aggregate is divided into blocks, the amount of calculation will be greatly increased. Therefore, in this paper, the size of 5 aggregates is selected to be divided inton × n small images of equal length and width. After the blocks are divided, the global threshold method is performed according to the block, and the threshold T of each image block is calculated. By calculating the average gray difference Δd between classes of each image block, the interference of pure background or pure object is excluded, and combined with Otsu algorithm, the image block is binarized based on this condition. Then, the binarized images of each image block are stitched together to form a complete crushed stone binarized image.(17)Δd=m1k−m2k,m1k=∑i=1kiPiC1,m2k=∑i=k+1LiPiC2,where C1 is the pixels whose gray level is [1,2, ⋯, k]; C2 is the pixels whose gray level is [k+1, k+2, ⋯, L]; Pi is the frequency distribution of gray levels; m(k) is the average gray value of gray levels from 1 to k; and Δd is the average grayscale difference between classes.The average grayscale difference Δd between classes calculated according to the formula is combined with the threshold value T of each image block to realize the binarization operation on each image block. When the interclass grayscale difference Δd of the image block is less than the set value d0, it is determined that the pixels in the image block are of the same class, and the image block is not binarized; otherwise, the image block is binarized. The detection method in this paper divides the image into 4 × 4 image blocks (i.e., n = 4), the accuracy of the detection result is 98.1%, and the average processing speed of a single image is less than 85 ms. In this paper, comparisons of single aggregate, 5 aggregates, and 10 aggregates as a block size are listed in Table 1.Table 1
Comparison of single aggregate, 5 aggregates, and 10 aggregates.
Aggregate quantityAccuracy (%)OH (ms)1 aggregate98.21095 aggregates98.18510 aggregates96.283It can be seen from Table1 that when a single aggregate and 5 aggregates are selected as the size of the block, the recognition accuracy is basically the same, but the operation speed of a single aggregate is much slower than that of 5 aggregates. Too many blocks will slow down the recognition speed; the operation speed of 5 aggregates and 10 aggregates is roughly the same, but the accuracy of 5 aggregates is better than 10 aggregates, which is because the block is too large, the thresholds of the targets in the block are sometimes quite different, and the entire target cannot be accurately and completely segmented, resulting in a decrease in the accuracy rate. Therefore, 5 aggregates are selected as one block in this paper.Under a fixed shooting area, after obtaining a certain specification of the crushed stone distribution rate and the weight of the used crushed stone, the crushed stone weight and distribution amount of a certain specification of crushed stone in this area can be obtained under the condition of spreading. Furthermore, it is possible to detect the weight and the amount of crushed stones of different specifications of crushed stones under different spreading rates.(18)mi+2=mi−1/ri−1+mi/ri+mi+1/ri+1300×ri+2,(19)Ps=mi+2s,where ri-1, ri, and ri+1 are the known distributing rates of crushed stones at different moments; mi-1, mi, and mi+1 are the known crushed rock weights at different distributing rates of crushed stones; ri+2 is the distributing rate of crushed stones at any time (≤100%); mi+2 is the weight of the gravel at the spreading rate ri+2; and s is the shooting area.


## 2.2.1. Wavelet Decomposition and Reconstruction

The wavelet decomposition method is used to obtain low-frequency coefficients and high-frequency coefficients. The low-frequency coefficients mainly include the global information of the image, including the contour information of the image, and the high-frequency coefficients mainly include the local information of the image, including the edges, details, and noise of the image. The expression of two-dimensional discrete wavelet decomposition is as follows [16]:(5)Wφj0,m,n=1/MN∑x=0M−1∑y=0N−1fx,yφj0,m,nx,yWψij,m,n=1/MN∑x=0M−1∑y=0N−1fx,yψj,m,nix,yi=H,V,D,where Wφj0,m,n is the low-frequency coefficient after decomposition; Wψij,m,n is the high-frequency coefficient after decomposition, means H, V, and D; j0 is any starting scale, usually let it be 0; f(x,y) is the discrete function, x and y represent discrete variables; m and n are expressed as relative offsets; M×N represents the image is composed of pixels; φj0,m,nx,y is a two-dimensional scaling function, and H,V,D is the wavelet function corresponding to i in the horizontal, vertical, and diagonal directions.Process the low-frequency coefficientWφj0,m,n and high-frequency coefficient Wψij,m,n obtained after decomposition. The reconstructed image is obtained through the inverse discrete wavelet transform, and the wavelet reconstruction expression is as follows. The schematic diagram of the two-layer wavelet decomposition is shown in Figure 3 [17].(6)fx,y=1MN∑m∑nWφj0,m,nφj0,m,nx,y+1MN∑i=H,V,D∑j=j0∞∑m∑nWψij,m,nψj,m,nix,y,where S represents the decomposed image, cA1 and cA2, respectively, represent the low-frequency components of the two-layer wavelet decomposition, that is, low-frequency images, and cD1 and cD2 represent the high-frequency components of the two-level wavelet decomposition, that is, the horizontal (H), the vertical (V), and the diagonal (D) components of each layer. The core processing flow is as follows:(1)
Image wavelet decomposition: including selecting a certain wavelet base function and determining the optimal wavelet decomposition level N, perform N-layer wavelet decomposition on the noisy imagef(x, y) to obtain high-frequency wavelet coefficients ωj,k and low-frequency wavelet image S(x, y).(2)
Threshold quantization of high-frequency coefficients of wavelet decomposition: according to a threshold quantization criterion, an appropriate threshold is determined for each high-frequency coefficient of the first to Nth layers. By thresholdingωj,k, the estimated wavelet coefficients ω∧j,k are obtained, and ω∧j,k−ωj,k is as small as possible.(3)
Retinex enhancement of the low-frequency coefficients of wavelet decomposition: according to the principle of the Retinex algorithm, the low-frequency image is converted to the logarithmic domain to obtain the reflection image, and an illumination adjustment parameterk is added in the process of eliminating the illumination image, so that the obtained reflection image is more natural. After Gaussian filtering, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(4)
Wavelet reconstruction: inverse wavelet transform is performed on the high-frequency coefficientsωj,k of the 1st to Nth layers after threshold quantization and the low-frequency wavelet coefficients of the Nth layer to obtain the estimated image.Figure 3
Schematic diagram of wavelet decomposition.


## 2.2.2. Improved Threshold Function Wavelet Denoising Method

(1) Selection of Wavelet Threshold. An improved threshold function method is used to denoise the high-frequency coefficients after wavelet decomposition. In order to effectively separate the noise components in the high-frequency coefficients, the function threshold must first be determined. When the threshold function is too small, the denoising effect of the image is not obvious, and the denoising result image is not much different from the original image, which loses the meaning of threshold denoising. When the threshold function is too large, the degree of image noise removal is too large. Although the noise is removed more thoroughly, the image is too smooth and the image is blurred, which cannot reflect the essential information of the image. Therefore, choosing a suitable threshold is of great significance to the results of threshold denoising [18].The fixed threshold estimation method has a more thorough denoising effect when the noise is more distributed at high-frequency coefficients. Therefore, the fixed threshold estimation method is used to estimate the function threshold, and the expression is as follows [19]:(7)λ=σm2logM⋅Nσm=medianabsk/0.6745,where λ is the determined fixed threshold; and σm is the standard deviation of noise.(2) Selection of Threshold Function. The traditional soft threshold (ST) and hard threshold (HT) function methods are less computationally intensive and easy to implement, but the method will produce constant errors when the high-frequency coefficient is greater than the fixed threshold, resulting in blurred edge details. Through the analysis of ST function in Figure 4, we can see that the continuity of ST function at ±λ is better, but when ωj,k>λ, the estimated wavelet coefficient ωj,k∧ obtained by the ST function has a constant deviation λ from the original threshold coefficient ωj,k. The constant deviation will cause the loss of high-frequency information and cause distortion. From the analysis of the HT function in Figure 5, it can be seen that the HT method solves the constant error problem of the ST method, but there is a discontinuity point ±λ. The HT function is not continuous at the discontinuity point ±λ, and it is easy to cause oscillation when reconstructing the image.(1)
Soft threshold function is as follows:(8)ωj,k∧=sgnωj,k⋅ωj,k−λ0ωj,k>λωj,k≤λ.(2)
Hard threshold function is as follows:(9)ωj,k∧=ωj,k0ωj,k>λωj,k≤λ.(3)
This paper improves the threshold function:(10)ωj,k∧=sgnωj,k⋅ωj,k−λeωj,k−λ0ωj,k≥λωj,k<λ.Figure 4
Soft threshold function.Figure 5
Hard threshold function.From the analysis of equation (13), we can see that in terms of function continuity, when the wavelet coefficients tend to a fixed threshold, the estimated wavelet coefficient limit is 0, so the function value is continuous at the fixed threshold, avoiding the shock phenomenon when the hard threshold function has a discontinuity point. In terms of function asymptoticity, when the wavelet coefficients tend to positive or negative infinity, the limit of the ratio of the estimated wavelet coefficients to the original wavelet coefficients is 1, which proves that the improved threshold function takes ωj,k∧=ωj,k as the asymptote; thus, eliminating the constant error of the ST function problem. In terms of function deviation, the difference between the estimated wavelet coefficient and the original wavelet coefficient is calculated with the wavelet coefficient tending to plus or minus infinity as the limit, and the difference is 0, so when the wavelet coefficient tends to infinity, the estimated new wavelet coefficient is infinite closed to the original wavelet coefficients to avoid the problem of deviation.


## 2.2.3. Image Enhancement Algorithm Based on Improved Retinex

An improved Retinex image enhancement algorithm is used to estimate the illumination image of the low-frequency coefficients after wavelet decomposition, and the reflection image is obtained through logarithmic domain conversion. The improved Retinex algorithm in this paper adds a light adjustment parameter based on the single-scale filter Retinex (SSR) algorithm, which makes the reflected image more natural, avoids the halo phenomenon in effect, and has a better edge retention effect.(1) Illumination Image Estimation. According to the principle of Retinex algorithm, the illumination image is the convolution of the original image and the Gaussian kernel function. In order to obtain the reflection image and eliminate the influence of the illumination on the image, the function is converted to the logarithmic domain, and then the reflection image after the logarithmic domain conversion is obtained. The number field conversion expression is as follows [20]:(11)Rx,y=logaSx,y−logaSx,y∗Fx,y.However, through experiments on the SSR algorithm, it is found that if the illumination estimation part in the original image is completely eliminated, the image will often be unnatural. Therefore, an illumination adjustment parameter is added in the process of eliminating the illumination image to make the reflected image more natural. The expression of the reflection image after introducing the adjustment factor is as follows:(12)Rx,y=logaSx,y−klogaSx,y∗Fx,y,where k is the light adjustment parameter; F(x,y) is the Gaussian kernel function; Rx,i is the illumination image; and a is the number of image channels, usually 3.After Gaussian filtering is performed, the filtering results on different scales are averagely weighted to obtain the estimated illuminance image.(13)Rx,y=∑k=1NωklogaSx,y−klogaSx,y∗Fkx,y,where ωk is the weight coefficient.


## 2.2.4. Reflection Image Stretch

The reflection image obtained after the image is transformed in the logarithmic domain often cannot reach the range of 0∼255 in the image value domain, which may cause the image to be gray and white, which affects the visual effect and subsequent processing. Therefore, the reflection image needs to be stretched. Firstly, by traversing the pixel values, find the maximum and minimum pixel values in the reflection image, and then stretch the image by the following formula to stretch the reflection image to the range of [0,255] [21].(14)dx,y=Rx,y−RminRmax−Rmin×255−0,where d(x,y) represents the reflected image after stretching; R(x,y) represents the reflected image before stretching; and Rmax and Rmin represent the maximum and minimum gray values of the reflected image before stretching, respectively.


## 2.2.5. Piecewise Linear Transformation

Aiming at the pan-gray phenomenon in the reconstructed image, this paper uses a three-segment piecewise linear transformation method to improve the image contrast. The three-stage piecewise linear transformation expression is as follows:(15)fi,j=k1×di,j0≤di,j<ab+k2×di,ja≤di,j≤cd+k3×di,jc<di,j≤255,where f(i,j) is the output image after contrast enhancement; d(i,j) is the input image; and k1, k2, and k3 represent the slope of the three-stage transformation. and the expression is as follows [22]:(16)k1=ab,k2=d−bc−a,k3=255−d255−c,where (a,b) and (c,d) indicate the point where the slope of the piecewise function changes.The segmentation function can set different segmentation points according to the image characteristics, enhance the details of the region of interest, and can also suppress the region of interest.


## 2.2.6. Detection of Spreading Rate and Spreading Amount Based on Segmentation Threshold Segmentation Method

By dividing the image into several blocks and performing threshold segmentation, the uneven effect caused by illumination or reflection can be solved to a certain extent. The blocks are chosen to be small enough so that the lighting of each block is approximately uniform, so that when auto-thresholding, high-threshold segmentation is used in high-gray areas and low-threshold segmentation is used in low-gray areas. Since the size of a single aggregate is too small, if a single aggregate is divided into blocks, the amount of calculation will be greatly increased. Therefore, in this paper, the size of 5 aggregates is selected to be divided inton × n small images of equal length and width. After the blocks are divided, the global threshold method is performed according to the block, and the threshold T of each image block is calculated. By calculating the average gray difference Δd between classes of each image block, the interference of pure background or pure object is excluded, and combined with Otsu algorithm, the image block is binarized based on this condition. Then, the binarized images of each image block are stitched together to form a complete crushed stone binarized image.(17)Δd=m1k−m2k,m1k=∑i=1kiPiC1,m2k=∑i=k+1LiPiC2,where C1 is the pixels whose gray level is [1,2, ⋯, k]; C2 is the pixels whose gray level is [k+1, k+2, ⋯, L]; Pi is the frequency distribution of gray levels; m(k) is the average gray value of gray levels from 1 to k; and Δd is the average grayscale difference between classes.The average grayscale difference Δd between classes calculated according to the formula is combined with the threshold value T of each image block to realize the binarization operation on each image block. When the interclass grayscale difference Δd of the image block is less than the set value d0, it is determined that the pixels in the image block are of the same class, and the image block is not binarized; otherwise, the image block is binarized. The detection method in this paper divides the image into 4 × 4 image blocks (i.e., n = 4), the accuracy of the detection result is 98.1%, and the average processing speed of a single image is less than 85 ms. In this paper, comparisons of single aggregate, 5 aggregates, and 10 aggregates as a block size are listed in Table 1.Table 1
Comparison of single aggregate, 5 aggregates, and 10 aggregates.
Aggregate quantityAccuracy (%)OH (ms)1 aggregate98.21095 aggregates98.18510 aggregates96.283It can be seen from Table1 that when a single aggregate and 5 aggregates are selected as the size of the block, the recognition accuracy is basically the same, but the operation speed of a single aggregate is much slower than that of 5 aggregates. Too many blocks will slow down the recognition speed; the operation speed of 5 aggregates and 10 aggregates is roughly the same, but the accuracy of 5 aggregates is better than 10 aggregates, which is because the block is too large, the thresholds of the targets in the block are sometimes quite different, and the entire target cannot be accurately and completely segmented, resulting in a decrease in the accuracy rate. Therefore, 5 aggregates are selected as one block in this paper.Under a fixed shooting area, after obtaining a certain specification of the crushed stone distribution rate and the weight of the used crushed stone, the crushed stone weight and distribution amount of a certain specification of crushed stone in this area can be obtained under the condition of spreading. Furthermore, it is possible to detect the weight and the amount of crushed stones of different specifications of crushed stones under different spreading rates.(18)mi+2=mi−1/ri−1+mi/ri+mi+1/ri+1300×ri+2,(19)Ps=mi+2s,where ri-1, ri, and ri+1 are the known distributing rates of crushed stones at different moments; mi-1, mi, and mi+1 are the known crushed rock weights at different distributing rates of crushed stones; ri+2 is the distributing rate of crushed stones at any time (≤100%); mi+2 is the weight of the gravel at the spreading rate ri+2; and s is the shooting area.


## 3. Results and Discussions

### 3.1. Processing Result of a Single Image

In this study, using 9.5 mm-13.2 mm gravel particles as an example, a picture is randomly selected for segmentation, as shown in Figure6. The processing of each phase is showed in Figures 7-10.Figure 6
9.5-13.2 mm indoor spreading image.Figure 7
Improve the threshold function to process high-frequency coefficients.Figure 8
Improved Retinex algorithm to process low-frequency coefficients.Figure 9
Wavelet reconstruction.Figure 10
Binary segmentation image.The test value of the spreading rate of the scattered rubble is 54.11%, the actual value is 54.18%, and the difference is 0.07%, within the error range±2%, indicating that this algorithm can meet the detection requirements.As shown in Table2, comparing the proposed method with hard threshold, soft threshold, and multiscale Retinex, it can be found that the peak signal-to-noise ratio and mean square error of the proposed method are better than other single methods, and the detection speed is basically the same.Table 2
Comparison of proposed method with a single detection method.
Evaluation indicatorOriginal imageHard thresholdSoft thresholdMsrProposedPSNR18.984225.617328.361720.452628.3943MSE0.02670.00470.00290.00270.0020OH—0.64710.67160.75430.7556The proposed method is compared with the image reconstructed by Krawtchouk and Hahn polynomials. As shown in Table3, when the mean square error is smaller than the original image, the orders of these two polynomials are 10 and 14, respectively, but the peak signal-to-noise ratio and the mean square error are still inferior to the proposed method. Also, the calculation speed is much lower than that of the proposed method. If the orders of Krawtchouk and Hahn polynomials continue to increase, the detection accuracy can be better than the proposed method, but the calculation speed is too slow which cannot meet the real-time requirements.Table 3
Comparison of different reconstruction algorithms.
Evaluation indicatorOriginal imageKrawtchouk (p1 = 0.5, p2 = 0.5, 10th order)Hahn (a = 10, b = 10, 14th order)ProposedPSNR18.984219.013526.536228.3943MSE0.02670.01290.00780.0020OH—1.35661.43280.7556


### 3.2. Different Specifications of Crushed Stone and Testing Test of Spreading Rate

Using the proposed method, a series of experiments have been carried out on the spreading conditions of different specifications of crushed stones under different spreading rates and have been tested. Figure11 is an indoor weighing picture of indoor gravel.Figure 11
9.5-13.2 mm indoor gravel weighing.Table 4
Existing method indoor gravel scattering test gravel scattering rate (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)115.620.929.4233.931.731.4350.939.340.6Table 5
Detected value of indoor gravel distribution rate in this method (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)115.221.828.8234.431.931.4350.138.339.9Take the 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm crushed stones, respectively, for test verification, as shown in Tables4 to 7.Table 6
Indoor gravel spread gravel weight (g).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)150.6183.1303.62109.7277.6324.93164.9344.2419.7Table 7
Difference between the existing method and the proposed method (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)10.40.90.620.50.20.030.81.00.7Average0.60.70.4The results of the gravel treatment of 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm are shown in Figures12-14.Figure 12
16-19 mm original image and inspection image.
(a)(b)(c)(d)Figure 13
9.5-13.2 mm original image and inspection image.
(a)(b)(c)(d)Figure 14
2.36-4.75 mm original image and inspection image.
(a)(b)(c)(d)This paper carries out testing tests for different specifications of crushed stones and spreading rates. Compared with the existing method, the average value of the difference in the detection accuracy of the three types of gravel is 0.6%, 0.7%, and 0.4% respectively. From the detection results, it can be seen that the proposed method can meet the requirements of detection accuracy.


### 3.3. Detection of the Amount of Scattered Stones of Different Specifications

#### 3.3.1. Calculation of Spreading Amount and Full Spreading Amount of Different Specifications of Gravel

According to (18) and (19), the weight of each kind of crushed stone can be calculated under the full spread, and then the scattering quality and amount of each kind of stone can be calculated under different spreading rates, as shown in Table 8.Table 8
Spreading amount and full coverage of different specifications of gravel.
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)Full paved gravel weight/g323.98875.871033.7Amount of rubble/(kg/m2)1.591.601.50


#### 3.3.2. Calculation Result Verification

The spreading pictures of two kinds of stone materials randomly collected are shown in Figures15-16. Thespreading rates detected by the proposed method are 19.5% and 33.1%, respectively, and the weights of crushed stones are calculated to be 63.18 g and 342.15 g, respectively. The actual weighed weight is 60.3 g and 340.6 g, and the difference is 2.88 g and 1.55 g, respectively. It is within the allowable range of error, indicating that the proposed method can basically meet the requirements of the rubble spreading detection.Figure 15
16-19 mm original image and inspection image.
(a)(b)Figure 16
2.36-4.75 mm original image and inspection image.
(a)(b)


## 3.1. Processing Result of a Single Image

In this study, using 9.5 mm-13.2 mm gravel particles as an example, a picture is randomly selected for segmentation, as shown in Figure6. The processing of each phase is showed in Figures 7-10.Figure 6
9.5-13.2 mm indoor spreading image.Figure 7
Improve the threshold function to process high-frequency coefficients.Figure 8
Improved Retinex algorithm to process low-frequency coefficients.Figure 9
Wavelet reconstruction.Figure 10
Binary segmentation image.The test value of the spreading rate of the scattered rubble is 54.11%, the actual value is 54.18%, and the difference is 0.07%, within the error range±2%, indicating that this algorithm can meet the detection requirements.As shown in Table2, comparing the proposed method with hard threshold, soft threshold, and multiscale Retinex, it can be found that the peak signal-to-noise ratio and mean square error of the proposed method are better than other single methods, and the detection speed is basically the same.Table 2
Comparison of proposed method with a single detection method.
Evaluation indicatorOriginal imageHard thresholdSoft thresholdMsrProposedPSNR18.984225.617328.361720.452628.3943MSE0.02670.00470.00290.00270.0020OH—0.64710.67160.75430.7556The proposed method is compared with the image reconstructed by Krawtchouk and Hahn polynomials. As shown in Table3, when the mean square error is smaller than the original image, the orders of these two polynomials are 10 and 14, respectively, but the peak signal-to-noise ratio and the mean square error are still inferior to the proposed method. Also, the calculation speed is much lower than that of the proposed method. If the orders of Krawtchouk and Hahn polynomials continue to increase, the detection accuracy can be better than the proposed method, but the calculation speed is too slow which cannot meet the real-time requirements.Table 3
Comparison of different reconstruction algorithms.
Evaluation indicatorOriginal imageKrawtchouk (p1 = 0.5, p2 = 0.5, 10th order)Hahn (a = 10, b = 10, 14th order)ProposedPSNR18.984219.013526.536228.3943MSE0.02670.01290.00780.0020OH—1.35661.43280.7556


## 3.2. Different Specifications of Crushed Stone and Testing Test of Spreading Rate

Using the proposed method, a series of experiments have been carried out on the spreading conditions of different specifications of crushed stones under different spreading rates and have been tested. Figure11 is an indoor weighing picture of indoor gravel.Figure 11
9.5-13.2 mm indoor gravel weighing.Table 4
Existing method indoor gravel scattering test gravel scattering rate (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)115.620.929.4233.931.731.4350.939.340.6Table 5
Detected value of indoor gravel distribution rate in this method (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)115.221.828.8234.431.931.4350.138.339.9Take the 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm crushed stones, respectively, for test verification, as shown in Tables4 to 7.Table 6
Indoor gravel spread gravel weight (g).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)150.6183.1303.62109.7277.6324.93164.9344.2419.7Table 7
Difference between the existing method and the proposed method (%).
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)10.40.90.620.50.20.030.81.00.7Average0.60.70.4The results of the gravel treatment of 2.36-4.75 mm, 9.5-13.2 mm, and 16-19 mm are shown in Figures12-14.Figure 12
16-19 mm original image and inspection image.
(a)(b)(c)(d)Figure 13
9.5-13.2 mm original image and inspection image.
(a)(b)(c)(d)Figure 14
2.36-4.75 mm original image and inspection image.
(a)(b)(c)(d)This paper carries out testing tests for different specifications of crushed stones and spreading rates. Compared with the existing method, the average value of the difference in the detection accuracy of the three types of gravel is 0.6%, 0.7%, and 0.4% respectively. From the detection results, it can be seen that the proposed method can meet the requirements of detection accuracy.


## 3.3. Detection of the Amount of Scattered Stones of Different Specifications

### 3.3.1. Calculation of Spreading Amount and Full Spreading Amount of Different Specifications of Gravel

According to (18) and (19), the weight of each kind of crushed stone can be calculated under the full spread, and then the scattering quality and amount of each kind of stone can be calculated under different spreading rates, as shown in Table 8.Table 8
Spreading amount and full coverage of different specifications of gravel.
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)Full paved gravel weight/g323.98875.871033.7Amount of rubble/(kg/m2)1.591.601.50


### 3.3.2. Calculation Result Verification

The spreading pictures of two kinds of stone materials randomly collected are shown in Figures15-16. Thespreading rates detected by the proposed method are 19.5% and 33.1%, respectively, and the weights of crushed stones are calculated to be 63.18 g and 342.15 g, respectively. The actual weighed weight is 60.3 g and 340.6 g, and the difference is 2.88 g and 1.55 g, respectively. It is within the allowable range of error, indicating that the proposed method can basically meet the requirements of the rubble spreading detection.Figure 15
16-19 mm original image and inspection image.
(a)(b)Figure 16
2.36-4.75 mm original image and inspection image.
(a)(b)


## 3.3.1. Calculation of Spreading Amount and Full Spreading Amount of Different Specifications of Gravel

According to (18) and (19), the weight of each kind of crushed stone can be calculated under the full spread, and then the scattering quality and amount of each kind of stone can be calculated under different spreading rates, as shown in Table 8.Table 8
Spreading amount and full coverage of different specifications of gravel.
Size2.36-4.75 (mm)9.5-13.2 (mm)16-19 (mm)Full paved gravel weight/g323.98875.871033.7Amount of rubble/(kg/m2)1.591.601.50


## 3.3.2. Calculation Result Verification

The spreading pictures of two kinds of stone materials randomly collected are shown in Figures15-16. Thespreading rates detected by the proposed method are 19.5% and 33.1%, respectively, and the weights of crushed stones are calculated to be 63.18 g and 342.15 g, respectively. The actual weighed weight is 60.3 g and 340.6 g, and the difference is 2.88 g and 1.55 g, respectively. It is within the allowable range of error, indicating that the proposed method can basically meet the requirements of the rubble spreading detection.Figure 15
16-19 mm original image and inspection image.
(a)(b)Figure 16
2.36-4.75 mm original image and inspection image.
(a)(b)


## 4. Actual Engineering Verification

This paper collects images of the road surface after the synchronous gravel sealer has been spread in actual engineering to verify the effectiveness of the algorithm. The collection location was a high-speed construction bid section in Inner Mongolia, and the road surface after the synchronous gravel seal truck was spread. The working process of the synchronous gravel sealing truck is shown in Figure17. The collected pictures are shown in Figure 18.Figure 17
Spreading process of synchronous gravel seal car.Figure 18
Original and inspection images of the actual spreading road.
(a)(b)(c)(d)(e)(f)The pictures collected this time are the spreading construction pictures of the lower seal layer, and the spreading gravel specification is 9.5-13.2 mm. According to the processing results of the algorithm, the processing results of the three collected images were 71.08%, 71.69%, and 68.95%, respectively, with an average value of 70.57%. The average of the actual spreading rate test result is 69.81%, and the difference is 0.76%. Furthermore, the amount of spreading can be calculated to be 1.52 kg/m2. The test results of spreading rate and spreading amount are both within the allowable range of error.


## 5. Conclusions

For the pavement spread by the synchronous gravel sealer, environmental factors such as light and noise, as well as gravel of different specifications and spreading rates, have an impact on the detection accuracy. This paper processes the high-frequency and low-frequency coefficients of the image after two-dimensional wavelet decomposition and then reconstructs the image to eliminate the influence of environmental factors such as lighting. Next, the segmented Otsu threshold method is used to separate the asphalt and gravel to detect the spreading rate of gravel and to further establish the relationship between the spreading amount and the spreading rate.(i)
After extracting the high-frequency and low-frequency coefficients of the image by using two-dimensional wavelet decomposition, the high-frequency coefficients are processed by the improved threshold method, the low-frequency coefficients are processed by the improved Retinex method, and then the image is reconstructed, which can quickly, accurately, and effectively detect the spreading rate of gravel.(ii)
Through the segmented Otsu threshold processing of the image, the regional threshold can be set in each block, which can further eliminate the influence of environmental factors such as light, so as to accurately segment the asphalt and gravel, detect the rubble spread rate, and establish the relationship between the rubble spread rate and the amount of spread.(iii)
By collecting and testing crushed stones of different specifications and spreading rates, and comparing them with actual testing results, it can be found that the proposed method can meet the requirements of actual construction testing and has good reproducibility.(iv)
Through the proposed method, the poor quality of the spreading part can be found in time in the spreading process, so as to make real-time adjustments to ensure the quality of road construction.



---

*Source: 4017071-2022-10-19.xml*