A method for three-dimensional soil reconstruction

Abstract

A soil three-dimensional reconstruction method comprises the following steps: step 1: determining a preset sample point and a reconstruction sample point of a site, collecting soil samples and taking photos, and respectively obtaining soil parameters of the preset sample points and the reconstruction sample points; step 2: performing CT tomography scanning of the reconstructed soil sample, extracting cross-sectional image features, and preprocessing sample point parameters, taking the preset sample point parameters as a training set and the reconstruction sample point parameters as a test set, and using a composite machine learning model to obtain a target parameter correction value; step 3: based on the cross-sectional image features of the reconstructed soil sample, selecting a feature cross-sectional subset, combining the porosity correction value, and using an algorithm to generate and optimize a soil three-dimensional reconstruction result; step 4: obtaining a soil organic matter occurrence reconstruction model, and then obtaining a soil three-dimensional reconstruction model; step 5: improving and correcting the soil three-dimensional reconstruction model to obtain a final reconstruction result.

Description

Technical Field

[0001] The present invention relates to the field of environmental analysis technology, and in particular to a technology for accurately reconstructing soil samples in three dimensions and achieving accurate restoration of soil pores. Background Art

[0002] At present, under the guarantee of hydraulic direct push collection technology, the soil samples obtained under normal conditions have high fidelity and can be used for subsequent soil analysis. However, in complex site environments, due to the combined influence of various changing conditions, the high fidelity of soil samples is still insufficient. How to accurately obtain the in-situ physical properties of soil still requires further technical assistance.

[0003] In response to the above industry needs, technicians in related fields have made some improvements:

[0004] For example, patent application number CN 104933760 A proposes "a method for reconstructing soil CT images in three dimensions and searching for soil pores". Although the patent of this invention can reconstruct a three-dimensional stereogram of soil, visually observe the spatial morphological distribution and connectivity of pores, and search for the number and volume of soil pore clusters, the technical solution in the application document still has deficiencies. For example, the reconstruction of soil information depends entirely on the CT scanning of soil samples, which requires continuous layered scanning of soil samples, which is too expensive. In addition, the reconstruction of soil pores between two adjacent scanning sections uses a simple interpolation method, and the image is binarized without considering the organic matter in the soil, and the soil restoration degree has not been significantly improved.

[0005] In addition, patent application number CN 109959602 A proposes "a method for quantitatively determining the spatial structure of particulate organic matter in soil". Although the present invention patent extracts particulate organic matter in soil through physical and chemical means, and realizes the quantification of soil organic morphological characteristics, organic matter quantity, volume ratio and volume distribution, and pore size distribution through CT image scanning-artifact removal-threshold calculation-image segmentation-three-dimensional reconstruction and other process analysis, the technical solution in the application document is still insufficient. For example, the reconstruction of organic matter in the soil after wet screening and classification will cause disturbance to the original soil sample, resulting in changes in the soil pore structure. Therefore, the three-dimensional image obtained by subsequent reconstruction cannot accurately restore the original soil morphology. In addition, in order to accurately reconstruct the three-dimensional distribution of soil organic matter, it is necessary to perform a full-range CT scan under the condition of uniform rotation of the soil, and reconstruct with thousands of scanned images. In addition to considering the monetary cost, the operation is cumbersome and time-consuming, which is not conducive to large-scale promotion in actual engineering.

[0006] In view of the above shortcomings, the applicant proposed a soil three-dimensional reconstruction method. Summary of the invention

[0007] The purpose of the present invention is to solve the problems in existing soil reconstruction methods that soil mineral particles and organic matter components cannot be comprehensively inverted and the experimental costs are too high, and to propose a low-cost, high-efficiency, and scalable soil three-dimensional reconstruction technology that can simultaneously invert the soil three-dimensional pore structure and organic matter content.

[0008] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0009] A soil three-dimensional reconstruction method comprises the following steps:

[0010] Step 1: Determine the preset sampling points and reconstructed sampling points of the site, collect soil samples and take photos, and obtain the soil parameters of the preset sampling points and reconstructed sampling points respectively;

[0011] Step 2: Perform CT tomography scanning of the reconstructed soil sample, extract cross-sectional image features, and preprocess the sample point parameters. The preset sample point parameters are used as the training set and the reconstructed sample point parameters are used as the test set. The composite machine learning model is used to obtain the target parameter correction value.

[0012] Step 3: Based on the reconstructed soil sample cross-section image features, a feature cross-section subset is selected, and the porosity correction value is combined to generate and optimize the soil pore reconstruction result;

[0013] Step 4: Analyze the pore deformation during the sampling process to obtain the in-situ soil pore reconstruction model; combine the organic matter volume fraction correction value to obtain the soil organic matter occurrence reconstruction model, and then obtain the soil three-dimensional reconstruction model;

[0014] Step 5: Combined with the soil pore connectivity assessment and considering the dynamic disturbance process of sampling, the soil 3D reconstruction model is improved and corrected to obtain the final reconstruction result.

[0015] In step 1, the following sub-steps are included:

[0016] Step 1-1: define the sample points for the required three-dimensional soil reconstruction as reconstruction sample points, the sample points for collecting soil parameters as preset sample points, determine the preset sample points, and select any point as the reconstruction sample point;

[0017] Step 1-2: Use sampling equipment to collect soil samples from all sampling points and obtain soil resistance parameters during the uniform speed drilling sampling process;

[0018] Step 1-3: After the collected soil samples are air-dried, all soil samples are photographed in a dark box with the same shooting parameters, exposure, and distance to obtain soil sample photo image information;

[0019] It also includes steps 1-4: obtaining the measured values ​​of the target parameters of the preset soil sample (ie, the soil sample collected at the preset sampling point), that is, measuring the porosity and the organic matter volume fraction.

[0020] In step 2, it includes extracting features of the cross-sectional image of the reconstructed soil sample (i.e., the soil sample collected at the reconstructed sampling point) and obtaining the deviation correction value of the target parameter;

[0021] The steps of reconstructing the feature extraction of the soil sample cross-section image include:

[0022] Step 2-1-1) Use a spiral CT scanner (maximum resolution of 50 μm), put the reconstructed soil sample (i.e., the soil sample collected at the reconstructed sampling point) into the CT scanner, scan a cross section from the top to the bottom at a fixed distance (0.5-1.0 mm), and scan at equal intervals to obtain M cross-sectional images of the soil sample;

[0023] Step 2-1-2) After removing artifacts from the obtained soil sample cross-sectional image, a grayscale histogram is established according to the different grayscale values ​​of each component; the original X-CT image has many noise points, which greatly interferes with the extraction and analysis of the pore structure. Before threshold segmentation of the two-dimensional image, the non-local means method (None Local Means) is used to process the noise of all soil sample cross-sectional images to remove a large number of poor noise points that are too small to avoid interference with the pore structure calculation;

[0024] The CT scan cross-sectional image obtained after Gaussian high-pass filtering is selected as the boundary condition required for soil three-dimensional reconstruction calculation. Its mathematical expression is:

[0025]

[0026] Where Q is the filtered image at the image coordinate (i, j); f(i, j) is the original pixel value at the image coordinate (i, j); k and l are the parameters of the matrix function in the x and y directions; K is the filter kernel matrix function, in Gaussian high-pass filtering:

[0027]

[0028] In the filtered images, the pores and soil skeleton are well segmented, and the microstructural information is preserved to a large extent;

[0029] Step 2-1-3) In the grayscale value histogram of each cross-sectional image of the reconstructed soil sample, the grayscale value of the middle trough between two peaks is selected as the threshold, and the image is segmented into three parts: pores, organic matter, and soil skeleton, so as to extract the features of the reconstructed cross-sectional image of the soil sample. The pores and soil skeleton parts in the image are separated to obtain the M pore cross-sections required for the subsequent three-dimensional reconstruction;

[0030] The acquisition of the correction values ​​of the target parameters (porosity and organic matter volume fraction) includes the following steps:

[0031] Step 2-2-1) Preprocessing the soil sample parameter set in the same way for all sample points, including feature extraction of multi-color space of soil sample image and normalization of multi-type parameters;

[0032] Step 2-2-2) Perform composite machine learning model training, hyperparameter adjustment and optimization output;

[0033] Step 2-2-3) Use the composite machine learning model to predict the target parameters of the reconstructed sample points, perform secondary correction based on the correction model and measured data, and output the correction value.

[0034] In step 2-2-1, the steps of soil sample parameter set preprocessing are as follows:

[0035] (1) Use the OpenCV software library to convert soil sample photo images into RGB, XYZ, HSV, HSL, Lab, and Luv color spaces, and extract feature information in multiple color spaces of the image;

[0036] (2) principal component analysis is used to reduce the dimension of the obtained multi-color space feature information, that is, the correlation between the color space parameters is calculated, and the color parameters with correlation greater than 0.7 are removed to obtain the sub-color space containing the main color space feature information;

[0037] (3) Using the Min-max normalization method, the sub-color space and soil resistance parameters are linearly transformed to map the range of all types of parameter sets to the interval [0,1].

[0038] In step 2-2-2, the steps of model training, hyperparameter adjustment and optimization output are as follows:

[0039] (1) The preset sample point sub-color space and soil resistance parameters obtained after normalization in step 2-2-1 are used as input sets, wherein 80% of the randomly divided input sets are used as training sets to complete the basic training of the composite machine learning model; the remaining 20% ​​are used as validation sets to achieve hyperparameter tuning of the composite machine learning model;

[0040] (2) The mean square error (MSE) function is used as the loss function for model training. The formula is:

[0041]

[0042] Where Y i , Y i are the i-th value in the true value and predicted value vector respectively, and n is the vector dimension.

[0043] (3) The training set is input into the composite machine learning model, and an output item is obtained through the first-layer support vector machine regression (SVR) operation; after the validation set is input into the SVR for hyperparameter tuning, the first-layer SVR operation is performed again using the same training set to obtain an optimized output item;

[0044] (4) The output of the first optimization is used as the input of the second-layer minimum absolute shrinkage and selection operator (Lasso) regression algorithm to obtain the secondary output. Similar to the first layer, the validation set is input into the Lasso regression for hyperparameter tuning, and the output of the first optimization is used again to perform the second-layer Lasso regression operation to obtain the secondary optimization output. At this point, the training of the composite machine learning model is completed, and the secondary optimization output is the model output value.

[0045] In step 2-2-3, the steps of reconstructing the prediction and correction of the sample point target parameters are as follows:

[0046] (1) The reconstructed sample point sub-color space and soil resistance parameters obtained after normalization in step 2-2-1 are used as the test set and input into the composite machine learning model; under the condition that the hyperparameters in step 2-2-2 have been tuned, the test set is subjected to the first layer SVR operation to obtain a first optimization output item; this item is used as the input of the second layer Lasso regression, and the second optimization output item obtained by the operation is the predicted value of the target parameter of the reconstructed sample point;

[0047] (2) Use a convolutional neural network (CNN) model to correct the target parameters. Input the test set (i.e., the reconstructed sample point parameters) into the pre-trained CNN model (the model training set comes from a multi-source heterogeneous soil database, which is a data set of soil image feature parameters and measured values ​​of target parameters obtained by integrating data from channels such as geological survey departments, literature journals, and construction units). Perform a weighted linear combination of the obtained CNN model output item and the secondary optimization output item (prediction value) to obtain a primary correction value. It can be expressed as:

[0048] Where Z1 is the output item of secondary optimization; Z CNN is the output item of CNN model; ω1 is the weight of CNN model output item; ω1 is set to range from 0.1 to 0.3, and the training sample size in soil multi-source heterogeneous database and the similarity between the properties of training sample soil samples and reconstructed soil samples are used as the basis for weight calibration. The larger the training sample size and the closer the properties between soil samples are, the greater the weight given to the correction value;

[0049] (3) Correction using the measured values ​​of preset sample points; after completing steps 1-4 to obtain the measured values ​​of the target parameters (i.e., the actual porosity and organic matter volume fraction) of the preset soil sample (i.e., the soil sample collected at the preset sample point), select n preset sample points closest to the reconstructed sample point as the correction sample point group, and use the correction value obtained in step 2-2-1 as the correction point group. Perform secondary deviation correction;

[0050] Based on the inverse distance weighting method, the weight calibration is performed on the correction effect of the preset sample points, that is, it is assumed that the correction effect of the preset sample points on the reconstructed sample points weakens with the increase of distance; the definition formula of the inverse distance weighting method is: Where Z2 is the measured deviation correction value; n is the number of points in the deviation correction sample group; Z i is the measured value of the target parameter of the i-th sample point in the correction sample point group; i is the weight coefficient of the sample point measurement value, which is given by the following formula:

[0051]

[0052] In the formula, Point to i The bth power. i is the distance between the i-th preset sample point and the reconstructed sample point, and b is the distance attenuation coefficient;

[0053] The first correction value Perform a weighted linear combination with the measured correction value Z2 to output the secondary correction value of the reconstructed sample point target parameter (porosity or organic matter volume fraction) It can be expressed by the formula:

[0054] In the formula, is the one-time correction value; Z2 is the measured correction value; ω2 is the weight of the measured correction value; the range of ω2 is set to 0.1-0.3, and the number of points n of the selected correction sample group and the similarity between the properties of the correction sample group and the reconstructed soil sample are used as the basis for weight calibration. The larger the n is, the greater the similarity of the physical and chemical properties between the soil samples is, and the greater the weight given to the correction value;

[0055] (4) Output the primary or secondary correction value of the target parameter of the reconstructed sample point, where p * Indicates the porosity correction value, Indicates the correction value of organic matter volume fraction.

[0056] In step 3, the following sub-steps are included:

[0057] Step 3-1: For the M reconstructed soil sample pore sections obtained in step 2-1-3, select m characteristic sections by a permutation and combination method based on certain rules, that is, Form a characteristic section subset;

[0058] Step 3-2: Evaluate the pore connectivity of the above-mentioned characteristic section subsets in order to carry out connectivity correction of the subsequent soil three-dimensional reconstruction model;

[0059] Step 3-3: Use the marching cube algorithm to generate the pore structure between characteristic sections (two adjacent characteristic sections divide a soil unit), and replace the generated pores with equivalent diameters to obtain the soil pore reconstruction result generated based on the characteristic sections;

[0060] Step 3-4: Iteratively generate different characteristic section subsets according to the rules described in step 3-1 to obtain a soil pore reconstruction result set generated by different characteristic sections;

[0061] Step 3-5: Combined with the porosity correction value obtained in step 2-2-3, the soil pore reconstruction result set is optimized to obtain the soil pore reconstruction optimization result.

[0062] In step 3-1, the specific rules for feature section subset selection and mutual iteration are as follows:

[0063] (1) The characteristic section subset must include the top and bottom sections of the soil sample;

[0064] (2) In the characteristic section subset, the distance between any group of adjacent sections does not exceed 5 times the CT scan spacing;

[0065] (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different;

[0066] (4) Define the section position P as the ratio of the distance between the section and the top of the soil mass to the height of the soil mass. It is stipulated that in the characteristic section subset that meets the requirements, in addition to the sections at the top and bottom of the soil sample, the four intervals of P = 0-0.25, 0.25-0.50, 0.50-0.75, and 0.75-1 must contain at least one section.

[0067] In step 3-2, the pore connectivity of the characteristic cross section is evaluated as follows:

[0068] (1) For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices and serve as the basis for the correction of pore connectivity in the subsequent step 5; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula:

[0069]

[0070] Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices.

[0071] When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected.

[0072] In step 3-3, the steps for generating soil pore structure between characteristic cross sections are as follows;

[0073] (1) For the first characteristic section subset obtained in step 3-2 The calculation is performed with the top section of the reconstructed soil sample (i.e., the first pore section) as the initial condition, and the calculation ends at the characteristic section adjacent to the top section, that is, the calculation of the first soil unit is completed.

[0074] The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces. The relationship can be expressed as:

[0075] F i,j,k =F(x i ,y j ,z k )

[0076] (i=1,2,…,N x ; j = 1, 2, ..., N y ; k = 1, 2, ..., N z )

[0077] In the formula, x i ,y j 、z k N is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are the x, y, z coordinates of pixels.

[0078] (2) Use software such as Avizo to count the pores generated by soil units, sort them by pore volume, and classify all pores into three categories based on the quantile classification method. The volume mean of each type of pore can be directly obtained from the software statistical results. The equivalent diameter of each type of pore is calculated (where is the volume mean of all pores under the i-th type of pores), and all the generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage.

[0079] (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit (if there are m characteristic sections, m-1 units are divided) from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil. The calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset. The three-dimensional visualization results can be drawn using image processing software such as Avizo.

[0080] The total volume of all pores generated in this reconstruction result (i.e., the sum of the pore volumes of each soil unit) is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now:

[0081]

[0082] Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section.

[0083] In step 3-4, the specific process is as follows:

[0084] Based on the specific rules for selecting feature section subsets in step 3-1, a combination of permutation and combination iterative algorithms and rule constraints are used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step 3-2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step 3-3, and the iteration ends;

[0085] In steps 3-5, the specific optimization process for the soil pore reconstruction result set is as follows:

[0086] (1) Using the porosity correction value p obtained in step 2-2-3 * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entropy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as:

[0087] Loss = -p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,...N

[0088] With the goal of minimizing the loss function, the optimal result of soil pore reconstruction is obtained. Let the porosity of the optimal result be p b .

[0089] In step 4, the following sub-steps are specifically included:

[0090] 4-1) Carry out compression test on reconstructed soil samples, measure the porosity of soil samples under different loads, and draw compression relationship curve between load and porosity;

[0091] 4-2) Estimate the load level of the reconstructed soil sample at the depth of the stratum, calculate the change in porosity of the soil sample during the process of transferring from the underground in-situ to the ground state according to the compression relationship curve, and use the porosity p of the optimal result to calculate the porosity of the soil sample. b Inversely infer the porosity p of the soil sample in the original state of the soil b ′, defines the compression factor

[0092] 4-3) Multiply all pore volumes of the soil pore reconstruction optimization results in step 3-5 by S, that is, compress the pores equally by S times, so as to restore the pore structure of the soil sample before the pore deformation occurs during the sampling process, and obtain the in-situ soil pore reconstruction model. The porosity of this model is p b ′;

[0093] (4) Considering that there may be overlap between pores and organic matter during CT scanning, the porosity p of the in situ soil pore reconstruction model is combined with the known b ′, the organic matter volume fraction correction value obtained in step 2.2-3 Corrected to The formula is:

[0094] Where p b ′ is the porosity of the in-situ soil pore reconstruction model, is the correction value of organic matter volume fraction.

[0095] 4-5) Considering the actual storage form of soil organic matter, a spherical shell model of soil organic matter storage is proposed: the storage form of organic matter is simplified to a spherical shell of a certain thickness outside a sphere equivalent to the pore structure.

[0096] Under the condition that the equivalent size and number of all porous spheres are known, an equivalent thickness d is preset so that a spherical shell with a thickness of d is attached to the outside of each porous sphere to represent the organic matter in the reconstructed model. Then the volume of organic matter on the sphere under a certain equivalent diameter is: The volume fraction of organic matter in the total volume of the soil sample on all equivalent spheres is: Among them, N iis the number of spheres contained in a certain equivalent size range, and V0 is the total volume of the reconstructed soil.

[0097] 4-6) Increase or decrease d within a specific range with a specific step size to obtain the soil organic matter storage reconstruction result set under different equivalent thicknesses. Iterative calculation is performed to obtain each equivalent thickness d i The corresponding The corrected organic matter volume fraction is the expected value, calculated using the cross entropy loss function and expected value The difference between them, the loss function can be expressed as:

[0098]

[0099] With the goal of minimizing the loss function, the optimal result of soil organic matter storage reconstruction is obtained. Assume that the equivalent thickness of the optimal result is d b The spherical shell of equivalent thickness is the reconstruction model of soil organic matter storage.

[0100] The in-situ soil pore reconstruction model (pore sphere) and the soil organic matter storage reconstruction model (organic matter spherical shell) were combined to obtain a three-dimensional soil reconstruction model that included pore structure and organic matter storage.

[0101] In step 5, the following sub-steps are included:

[0102] 5-1) Combined with the pore connectivity evaluation results in step 3-2, the spherical morphology of the pores and organic matter in the soil three-dimensional reconstruction model obtained in step 4 is artificially and volume-fixed longitudinally appropriately stretched, that is, the pore connectivity in the reconstruction model is corrected, so that the pore spaces that are far apart and separated from each other are brought closer to each other, and the soil pore connectivity structure and organic matter occurrence status in the actual scenario are as close as possible;

[0103] 5-2) Considering the influence and change of sampling disturbance on soil pore structure, the rationality of the soil three-dimensional reconstruction model after pore stretching is analyzed. The analysis of the dynamic process of sampling includes the disturbance of soil when the sampling equipment drills the sample, and the creep that occurs when the soil is transferred from the underground burial position to the ground. Based on this comprehensive expert opinion and prior knowledge base, the soil three-dimensional reconstruction model obtained after pore longitudinal stretching is artificially disturbed and corrected, and the corrected model is used as the final reconstruction result.

[0104] The present invention also includes a method for generating and optimizing soil pore reconstruction results, comprising the following steps:

[0105] Step S1: For M reconstructed soil sample cross-section images, m characteristic cross-sections are selected by a permutation and combination method based on certain rules, namely Form a characteristic section subset;

[0106] Step S2: evaluating the pore connectivity of a subset of characteristic sections;

[0107] Step S3: using the marching cube algorithm to generate the pore structure between the characteristic sections, performing equivalent replacement of the generated pores with equivalent diameters, and then obtaining the soil pore reconstruction result generated based on the characteristic sections;

[0108] Step S4: Iteratively generate different characteristic section subsets according to the rules described in step S1 to obtain a soil pore reconstruction result set generated by different characteristic sections;

[0109] Step S5: Optimize the soil pore reconstruction result set in combination with the obtained porosity correction value to obtain the soil pore reconstruction optimization result.

[0110] In step S1, the specific rules for feature section subset selection and mutual iteration are as follows:

[0111] (1) The characteristic section subset must include the top and bottom sections of the soil sample;

[0112] (2) In the characteristic section subset, the distance between any group of adjacent sections does not exceed 5 times the CT scan spacing;

[0113] (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different;

[0114] (4) Define the cross-sectional position P as the ratio of the distance between the cross-sectional area and the top of the soil mass to the height of the soil mass;

[0115] In step S2, the pore connectivity of the characteristic cross section is evaluated as follows:

[0116] (1) For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula:

[0117]

[0118] Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices;

[0119] When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected.

[0120] In step S3, the steps of generating the soil pore structure between characteristic sections are as follows;

[0121] (1) For the first feature section subset obtained in step S2 Based on the pore cross section obtained after gray value segmentation, the marching cubes algorithm is used to perform calculations with the top cross section of the reconstructed soil sample as the initial condition, and the calculation ends at the characteristic cross section adjacent to the top cross section, that is, the calculation of the first soil unit is completed;

[0122] The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces; its relationship can be expressed as:

[0123] F i,j,k =F(x i ,y j ,z k )

[0124] (i=1,2,…,N x ; j = 1, 2, ..., N y ; k = 1, 2, ..., N z )

[0125] In the formula, x i ,y j 、z k N is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are x, y, z coordinate pixels;

[0126] (2) The pores generated by the soil unit are counted using software such as Avizo, sorted by pore volume, and classified into three categories according to the quantile classification method. The volume mean of each type of pore can be directly obtained from the statistical results of the software; The equivalent diameter of each type of pore is calculated, where is the volume mean of all pores under the i-th type of pores, and all generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage;

[0127] (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil; the calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset, and the three-dimensional visualization results can be drawn by image processing software such as Avizo;

[0128] The total volume of all pores generated in this reconstruction result is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now:

[0129]

[0130] Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section;

[0131] In step S4, the specific process is as follows:

[0132] Based on the specific rules for selecting the feature section subset in step S1, a combination of permutation and combination iterative algorithm and rule constraints is used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step S2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step S3, and the iteration ends;

[0133] In step S5, the specific optimization process for the soil pore reconstruction result set is as follows:

[0134] The porosity correction value p * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entroy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as:

[0135] Loss = -p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,…N

[0136] Taking the minimization of the loss function as the goal, the optimal result of soil pore reconstruction is obtained; let the porosity of the optimal result be p b .

[0137] Compared with the prior art, the present invention has the following technical effects:

[0138] 1) The present invention collects a large amount of relevant information and data of soil samples from contaminated sites, and uses a variety of technical methods such as big data analysis and multimodal machine learning inversion to accurately reconstruct soil samples in three dimensions, and can accurately restore soil pores (including organic matter volume distribution information), thereby fully improving the authenticity and reliability of soil survey results;

[0139] 2) In view of the randomness problem of the three-dimensional reconstruction algorithm, the present invention introduces a machine learning method based on soil parameters, and proposes a set of optimization correction schemes for target parameter prediction. The loss function is constructed with the correction values ​​of soil porosity and organic matter volume fraction to achieve the optimal output of soil three-dimensional reconstruction results;

[0140] 3) The present invention improves and enhances the existing soil three-dimensional reconstruction method that only considers a single factor, realizes the multivariate reconstruction of the pore characteristics and organic matter occurrence conditions in the soil, and dynamically corrects the reconstruction results according to the disturbance and compression of soil sampling, making the final reconstruction model more realistic;

[0141] 4) The present invention combines multiple types of soil information parameters and adopts a composite machine learning model that interacts with heterogeneous databases and measured data to achieve the prediction and correction of soil porosity and organic matter volume fraction at any point in the research site, which can effectively eliminate the process of measuring soil porosity and organic matter at each sample point during the sampling process. BRIEF DESCRIPTION OF THE DRAWINGS

[0142] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0143] Figure 1 It is the main flow chart of the present invention;

[0144] Figure 2 A schematic diagram of the process of obtaining the correction value of the target parameter;

[0145] Figure 3 Generate a schematic flow chart for soil pore and organic matter reconstruction model;

[0146] Figure 4 It is a schematic diagram of the logical relationship of the soil three-dimensional reconstruction model and its sub-concepts;

[0147] Figure 5 To reconstruct two pore cross sections of soil sample #3 after CT tomography and image segmentation, the black part is the soil pores;

[0148] Figure 6 This is the result of soil 3D reconstruction based on reconstructed soil sample #3. The white part contains soil pores and organic matter, and the black part is the soil skeleton. DETAILED DESCRIPTION

[0149] like Figure 1 As shown, a soil three-dimensional reconstruction method comprises the following steps:

[0150] Step 1: Determine the preset sampling points and reconstructed sampling points of the site, collect soil samples and take photos, and obtain the soil parameters of the preset sampling points and reconstructed sampling points respectively;

[0151] Step 2: Perform CT tomography scanning of the reconstructed soil sample, extract cross-sectional image features, and preprocess the sample point parameters. Use the preset sample point parameters as the training set and the reconstructed sample point parameters as the test set, and use the composite machine learning model to obtain the target parameter correction value, such as Figure 2 As shown;

[0152] Step 3: Based on the reconstructed soil sample cross-section image features, select the characteristic cross-section subset, combine the porosity correction value, and generate and optimize the soil pore reconstruction results, such as Figure 3 , 4 As shown;

[0153] Step 4: Considering the pore deformation during the sampling process, the in-situ soil pore reconstruction model is obtained; combined with the organic matter volume fraction correction value, the soil organic matter occurrence reconstruction model is obtained, and then the soil three-dimensional reconstruction model is obtained, such as Figure 3 , 4 As shown;

[0154] Step 5: Combined with the soil pore connectivity assessment and considering the dynamic disturbance process of sampling, the soil 3D reconstruction model is improved and corrected to obtain the final reconstruction result.

[0155] In step 1, the following sub-steps are included:

[0156] Step 1-1: Define the sample points for the required soil 3D reconstruction as reconstruction sample points, and the sample points for collecting soil parameters as preset sample points. Determine the preset sample points and select any point as the reconstruction sample point;

[0157] Step 1-2: Use the EP2000s direct-push hydraulic sampling drill to drill the soil downward at a constant speed of 10 m / h to obtain soil samples at all sampling points at the target depth, and use a pressure sensor to obtain the horizontal tangential soil resistance parameters during the sampling process;

[0158] Step 1-3: Collect soil samples under a light intensity of 1.5×10 4 Lux, air-dried for 5 hours under the condition of wind speed of 16m / s, placed in a dark box, and photographed all soil samples with the same shooting parameters, exposure, and distance to obtain soil sample photo image information. The dark box is a rectangular imaging box with a circular window on the top for placing a camera to take pictures of soil samples; the soil samples are placed and fixed at the bottom of the box; the inner wall of the box is painted black, and LED strips are attached to the inner wall of the box to provide fill light when taking pictures of soil samples.

[0159] Step 1-4: Obtain the measured values ​​of the target parameters of the preset soil sample (i.e., the soil sample collected at the preset sampling point), that is, measure the porosity and organic matter volume fraction.

[0160] In step 2, the main contents include feature extraction of cross-sectional images of reconstructed soil samples (i.e., soil samples collected at reconstructed sampling points) and acquisition of deviation correction values ​​of target parameters;

[0161] The steps of reconstructing the feature extraction of the soil sample cross-section image include:

[0162] Step 2-1-1) Use a SOMATOM Definition Flash spiral CT scanner (maximum effective resolution is about 50 μm), place the reconstructed soil sample into the CT scanner, scan a cross section from the top to the bottom at a fixed distance (0.5-1.0 mm), and scan at equal intervals to obtain M DICOM format soil sample cross-sectional images. The feature extraction process of the soil sample cross-sectional image, the preprocessing of the soil sample parameter set, and the visualization of the soil pore and organic matter structure are all based on the same computer. The computer processor used is 4.00 GHz Intel Core i7-4790, the memory is 16 GB, the operating system is Windows 10, and the configuration environment is Visual Studio 2021 platform. To facilitate subsequent computer processing, the original soil CT image in DICOM format is converted into a Jpeg or bmp format file and saved;

[0163] Step 2-1-2) After removing artifacts from the obtained soil sample cross-sectional image, a grayscale histogram is established using the 3D visualization tool VTK image processing library according to the different grayscale values ​​of each component; the original X-CT image has many noise points, which greatly interferes with the extraction and analysis of the pore structure. Before threshold segmentation of the two-dimensional image, the non-local mean method (None Local Means) is used to process the noise of all soil sample cross-sectional images to remove a large number of poor noise points that are too small to avoid interference with the pore structure calculation;

[0164] The CT scan cross-sectional image obtained after Gaussian high-pass filtering is selected as the boundary condition required for soil three-dimensional reconstruction calculation. Its mathematical expression is:

[0165]

[0166] Where Q is the filtered image at the image coordinate (i, j); f(i, j) is the original pixel value at the image coordinate (i, j); k and l are the parameters of the matrix function in the x and y directions; K is the filter kernel matrix function in Gaussian high-pass filtering;

[0167]

[0168] In the filtered images, the pores and soil skeleton are well segmented, and the microstructural information is preserved to a large extent;

[0169] Step 2-1-3) In the grayscale value histogram of each cross-sectional image of the reconstructed soil sample, the grayscale value of the middle trough between two peaks is selected as the threshold, and the image is segmented into three parts: pores, organic matter, and soil skeleton, so as to extract the features of the reconstructed cross-sectional image of the soil sample. The pores and soil skeleton parts in the image are separated to obtain the M pore cross-sections required for the subsequent three-dimensional reconstruction;

[0170] The acquisition of the correction values ​​of the target parameters (porosity and organic matter volume fraction) includes the following steps:

[0171] Step 2-2-1) Preprocessing the soil sample parameter set in the same way for all sample points, including feature extraction of multi-color space of soil sample image and normalization of multi-type parameters;

[0172] Step 2-2-2) Perform composite machine learning model training, hyperparameter adjustment and optimization output;

[0173] Step 2-2-3) Use the composite machine learning model to predict the target parameters of the reconstructed sample points, perform secondary correction based on the correction model and measured data, and output the correction value.

[0174] In step 2-2-1, the steps of soil sample parameter set preprocessing are as follows:

[0175] (1) Use the OpenCV software library to convert soil sample photo images into RGB, XYZ, HSV, HSL, Lab, and Luv color spaces, and extract feature information in multiple color spaces of the image;

[0176] (2) principal component analysis is used to reduce the dimension of the obtained multi-color space feature information, that is, the correlation between the color space parameters is calculated, and the color parameters with correlation greater than 0.7 are removed to obtain the sub-color space containing the main color space feature information;

[0177] (3) Using the Min-max normalization method, the sub-color space and soil resistance parameters are linearly transformed to map the range of all types of parameter sets to the interval [0,1].

[0178] In step 2-2-2, the steps of model training, hyperparameter adjustment and optimization output are as follows:

[0179] (1) The preset sample point sub-color space and soil resistance parameters obtained after normalization in step 2.2-1 are used as input sets, wherein 80% of the randomly divided input sets are used as training sets to complete the basic training of the composite machine learning model; the remaining 20% ​​are used as validation sets to achieve hyperparameter tuning of the composite machine learning model;

[0180] (2) The mean square error (MSE) function is used as the loss function for model training. The formula is:

[0181]

[0182] Where Y i , are the i-th value in the true value and predicted value vector respectively, and n is the vector dimension.

[0183] (3) The training set is input into the composite machine learning model, and an output item is obtained through the first-layer support vector machine regression (SVR) operation; after the validation set is input into the SVR for hyperparameter tuning, the first-layer SVR operation is performed again using the same training set to obtain an optimized output item;

[0184] (4) The output of the first optimization is used as the input of the second-layer minimum absolute shrinkage and selection operator (Lasso) regression to obtain the secondary output. Similar to the first layer, the validation set is input into the Lasso regression for hyperparameter tuning, and the output of the first optimization is used again to perform the second-layer Lasso regression operation to obtain the secondary optimization output. At this point, the training of the composite machine learning model is completed, and the secondary optimization output is the model output value.

[0185] In step 2-2-3, the steps of reconstructing the prediction and correction of the sample point target parameters are as follows:

[0186] (1) The reconstructed sample point sub-color space and soil resistance parameters obtained after normalization in step 2-2-1 are used as the test set and input into the composite machine learning model; under the condition that the hyperparameters in step 2-2-2 have been tuned, the test set is subjected to the first layer SVR operation to obtain a first optimization output item; this item is used as the input of the second layer Lasso regression, and the second optimization output item obtained by the operation is the predicted value of the target parameter of the reconstructed sample point;

[0187] (2) Use the convolutional neural network (CNN) model to correct the target parameters; input the test set (i.e., the reconstructed sample point parameters) into the pre-trained CNN model (the model training set comes from the soil multi-source heterogeneous database, which is a data set of soil image feature parameters and target parameter measured values ​​obtained by integrating data from channels such as geological survey departments, literature journals, and construction units), and perform a weighted linear combination of the obtained CNN model output items and the secondary optimization output items (prediction values) to obtain a primary correction value. It can be expressed as:

[0188] Where Z1 is the output item of secondary optimization; Z CNN is the output item of CNN model; ω1 is the weight of CNN model output item; ω1 is set to range from 0.1 to 0.3, and the training sample size in soil multi-source heterogeneous database and the similarity between the properties of training sample soil samples and reconstructed soil samples are used as the basis for weight calibration. The larger the training sample size and the closer the properties between soil samples are, the greater the weight given to the correction value;

[0189] (3) Correction using the measured values ​​of preset sample points; after completing steps 1-4 to obtain the measured values ​​of the target parameters (i.e., the actual porosity and organic matter volume fraction) of the preset soil sample (i.e., the soil sample collected at the preset sample point), select n preset sample points closest to the reconstructed sample point as the correction sample point group, and use the correction value obtained in step 2-2-1 as the correction point group. Perform secondary deviation correction;

[0190] Based on the inverse distance weighting method, the weight calibration is performed on the correction effect of the preset sample points, that is, it is assumed that the correction effect of the preset sample points on the reconstructed sample points weakens with the increase of distance; the definition formula of the inverse distance weighting method is: Where Z2 is the measured deviation correction value; n is the number of points in the deviation correction sample group; Z i is the measured value of the target parameter of the i-th sample point in the correction sample point group; i is the weight coefficient of the sample point measurement value, which is given by the following formula:

[0191]

[0192] In the formula, Point to i The bth power. i is the distance between the i-th preset sample point and the reconstructed sample point, and b is the distance attenuation coefficient;

[0193] The first correction value Perform a weighted linear combination with the measured correction value Z2 to output the secondary correction value of the reconstructed sample point target parameter (porosity or organic matter volume fraction) It can be expressed by the formula:

[0194] In the formula, is the one-time correction value; Z2 is the measured correction value; ω2 is the weight of the measured correction value; the range of ω2 is set to 0.1-0.3, and the number of points n of the selected correction sample group and the similarity between the properties of the correction sample group and the reconstructed soil sample are used as the basis for weight calibration. The larger the n is, the greater the similarity of the physical and chemical properties between the soil samples is, and the greater the weight given to the correction value;

[0195] (4) Output the primary or secondary correction value of the target parameter of the reconstructed sample point, where p * Indicates the porosity correction value, Indicates the correction value of organic matter volume fraction.

[0196] In step 3, the following sub-steps are included:

[0197] Step 3-1: For the M reconstructed soil sample pore sections obtained in step 2-1-3, select m characteristic sections by a permutation and combination method based on certain rules, that is, Form a characteristic section subset;

[0198] Step 3-2: Evaluate the pore connectivity of the above-mentioned characteristic section subsets in order to carry out connectivity correction of the subsequent soil three-dimensional reconstruction model;

[0199] Step 3-3: Use the marching cube algorithm to generate the pore structure between characteristic sections (two adjacent characteristic sections divide a soil unit), and replace the generated pores with equivalent diameters to obtain the soil pore reconstruction result generated based on the characteristic sections;

[0200] Step 3-4: Iteratively generate different characteristic section subsets according to the rules described in step 3-1 to obtain a soil pore reconstruction result set generated by different characteristic sections;

[0201] Step 3-5: Combined with the porosity correction value obtained in step 2-2-3, the soil pore reconstruction result set is optimized to obtain the soil pore reconstruction optimization result.

[0202] In step 3-1, the specific rules for feature section subset selection and mutual iteration are as follows:

[0203] (1) The characteristic section subset must include the top and bottom sections of the soil sample;

[0204] (2) In the characteristic section subset, the distance between any group of adjacent sections does not exceed 5 times the CT scan spacing;

[0205] (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different;

[0206] (4) Define the section position P as the ratio of the distance between the section and the top of the soil mass to the height of the soil mass. It is stipulated that in the characteristic section subset that meets the requirements, in addition to the sections at the top and bottom of the soil sample, the four intervals of P = 0-0.25, 0.25-0.50, 0.50-0.75, and 0.75-1 must contain at least one section.

[0207] In step 3-2, the pore connectivity of the characteristic cross section is evaluated as follows:

[0208] (1) For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices and serve as the basis for the correction of pore connectivity in the subsequent step 5; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula:

[0209]

[0210] Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices.

[0211] When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected.

[0212] In step 3-3, the steps for generating soil pore structure between characteristic cross sections are as follows;

[0213] (1) For the first characteristic section subset obtained in step 3-2 The marching cubes algorithm is used to perform calculations with the top section of the reconstructed soil sample (i.e., the first pore section) as the initial condition, and the calculations are completed to the characteristic section adjacent to the top section, that is, the calculations for the first soil unit are completed.

[0214] The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces. The relationship can be expressed as:

[0215] F i,j,k =F(x i ,y j ,z k )

[0216] (i=1,2,…,N x ; j = 1, 2, ..., N y ; k = 1, 2, ..., Nz )

[0217] In the formula, x i ,y j 、z k N is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are the x, y, z coordinates of pixels.

[0218] (2) Use software such as Avizo to count the pores generated by soil units, sort them by pore volume, and classify all pores into three categories based on the quantile classification method. The volume mean of each type of pore can be directly obtained from the software statistical results. The equivalent diameter of each type of pore is calculated (where V i is the volume mean of all pores under the i-th type of pores), and all the generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage.

[0219] (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit (if there are m characteristic sections, m-1 units are divided) from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil. The calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset. The three-dimensional visualization results can be drawn using image processing software such as Avizo.

[0220] The total volume of all pores generated in this reconstruction result (i.e., the sum of the pore volumes of each soil unit) is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now:

[0221]

[0222] Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section.

[0223] In step 3-4, the specific process is as follows:

[0224] Based on the specific rules for selecting feature section subsets in step 3-1, a combination of permutation and combination iterative algorithms and rule constraints are used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step 3-2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step 3-3, and the iteration ends;

[0225] In steps 3-5, the specific optimization process for the soil pore reconstruction result set is as follows:

[0226] (1) Using the porosity correction value p obtained in step 2-2-3 * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entropy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as:

[0227] Loss = -p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,…N

[0228] With the goal of minimizing the loss function, the optimal result of soil pore reconstruction is obtained. Let the porosity of the optimal result be p b .

[0229] In step 4, the following sub-steps are specifically included:

[0230] 4-1) Carry out compression test on reconstructed soil samples, measure the porosity of soil samples under different loads, and draw compression relationship curve between load and porosity;

[0231] 4-2) Estimate the load level of the reconstructed soil sample at the depth of the stratum, calculate the change in porosity of the soil sample during the process of transferring from the underground in-situ to the ground state according to the compression relationship curve, and use the porosity p of the optimal result to calculate the porosity of the soil sample. b Inversely infer the porosity p of the soil sample in the original state of the soil b ′, defines the compression factor

[0232] 4-3) Multiply all pore volumes of the soil pore reconstruction optimization results in step 3-5 by S, that is, compress the pores equally by S times, so as to restore the pore structure of the soil sample before the pore deformation occurs during the sampling process, and obtain the in-situ soil pore reconstruction model. The porosity of this model is p b ′;

[0233] (4) Considering that there may be overlap between pores and organic matter during CT scanning, the porosity p of the in situ soil pore reconstruction model is combined with the known b ′, the organic matter volume fraction correction value obtained in step 2.2-3 Corrected to The formula is:

[0234] Where p b ′ is the porosity of the in-situ soil pore reconstruction model, is the correction value of organic matter volume fraction.

[0235] 4-5) Considering the actual storage form of soil organic matter, a spherical shell model of soil organic matter storage is proposed: the storage form of organic matter is simplified to a spherical shell of a certain thickness outside a sphere equivalent to the pore structure.

[0236] Under the condition that the equivalent size and number of all porous spheres are known, an equivalent thickness d is preset so that a spherical shell with a thickness of d is attached to the outside of each porous sphere to represent the organic matter in the reconstructed model. Then the volume of organic matter on the sphere under a certain equivalent diameter is: The volume fraction of organic matter in the total volume of the soil sample on all equivalent spheres is: Among them, N i is the number of spheres contained in a certain equivalent size range, and V0 is the total volume of the reconstructed soil.

[0237] 4-6) Increase or decrease d within a specific range with a specific step size to obtain the soil organic matter storage reconstruction result set under different equivalent thicknesses. Iterative calculation is performed to obtain each equivalent thickness d i The corresponding The corrected organic matter volume fraction is the expected value, calculated using the cross entropy loss function and expected value The difference between them, the loss function can be expressed as:

[0238]

[0239] With the goal of minimizing the loss function, the optimal result of soil organic matter storage reconstruction is obtained. Assume that the equivalent thickness of the optimal result is d b The spherical shell of equivalent thickness is the reconstruction model of soil organic matter storage.

[0240] The in-situ soil pore reconstruction model (pore sphere) and the soil organic matter storage reconstruction model (organic matter spherical shell) were combined to obtain a three-dimensional soil reconstruction model that included pore structure and organic matter storage.

[0241] In step 5, the following sub-steps are included:

[0242] 5-1) Combined with the pore connectivity evaluation results in step 3-2, the spherical morphology of the pores and organic matter in the soil three-dimensional reconstruction model obtained in step 4 is artificially and volume-fixed longitudinally appropriately stretched, that is, the pore connectivity in the reconstruction model is corrected, so that the pore spaces that are far apart and separated from each other are brought closer to each other, and the soil pore connectivity structure and organic matter occurrence status in the actual scenario are as close as possible;

[0243] 5-2) Considering the influence and change of sampling disturbance on soil pore structure, the rationality of the soil three-dimensional reconstruction model after pore stretching is analyzed. The analysis of the dynamic process of sampling includes the disturbance of soil when the sampling equipment drills the sample, and the creep that occurs when the soil is transferred from the underground burial position to the ground. Based on this comprehensive expert opinion and prior knowledge base, the soil three-dimensional reconstruction model obtained after pore longitudinal stretching is artificially disturbed and corrected, and the corrected model is used as the final reconstruction result.

[0244] The present invention also includes a method for generating and optimizing soil pore reconstruction results, comprising the following steps:

[0245] Step S1: For the M reconstructed soil sample pore sections, m characteristic sections are selected by a permutation and combination method based on certain rules, namely Form a characteristic section subset;

[0246] Step S2: evaluating the pore connectivity of a subset of characteristic sections;

[0247] Step S3: using the marching cube algorithm to generate the pore structure between the characteristic sections, performing equivalent replacement of the generated pores with equivalent diameters, and then obtaining the soil pore reconstruction result generated based on the characteristic sections;

[0248] Step S4: Iteratively generate different characteristic section subsets according to the rules described in step S1 to obtain a soil pore reconstruction result set generated by different characteristic sections;

[0249] Step S5: Optimize the soil pore reconstruction result set in combination with the obtained porosity correction value to obtain the soil pore reconstruction optimization result.

[0250] In step S1, the specific rules for feature section subset selection and mutual iteration are as follows:

[0251] (1) The characteristic section subset must include the top and bottom sections of the soil sample;

[0252] (2) In the characteristic section subset, the distance between any group of adjacent sections does not exceed 5 times the CT scan spacing;

[0253] (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different;

[0254] (4) Define the cross-sectional position P as the ratio of the distance between the cross-sectional area and the top of the soil mass to the height of the soil mass;

[0255] In step S2, the pore connectivity of the characteristic cross section is evaluated as follows:

[0256] (1) For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula:

[0257]

[0258] Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices;

[0259] When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected.

[0260] In step S3, the steps of generating the soil pore structure between characteristic sections are as follows;

[0261] (1) For the first feature section subset obtained in step S2 Based on the pore cross section obtained after gray value segmentation, the marching cubes algorithm is used to perform calculations with the top cross section of the reconstructed soil sample as the initial condition, and the calculation ends at the characteristic cross section adjacent to the top cross section, that is, the calculation of the first soil unit is completed;

[0262] The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces; its relationship can be expressed as:

[0263] F i,j,k =F(x i ,y j ,z k )

[0264] (i=1,2,…,N x ; j = 1, 2, ..., N y ; k = 1, 2, ..., N z )

[0265] In the formula, x i ,y j 、z kN is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are x, y, z coordinate pixels;

[0266] (2) The pores generated by the soil unit are counted using software such as Avizo, sorted by pore volume, and classified into three categories according to the quantile classification method. The volume mean of each type of pore can be directly obtained from the statistical results of the software; The equivalent diameter of each type of pore is calculated, where is the volume mean of all pores under the i-th type of pores, and all generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage;

[0267] (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil; the calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset, and the three-dimensional visualization results can be drawn by image processing software such as Avizo;

[0268] The total volume of all pores generated in this reconstruction result is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now:

[0269]

[0270] Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section;

[0271] In step S4, the specific process is as follows:

[0272] Based on the specific rules for selecting the feature section subset in step S1, a combination of permutation and combination iterative algorithm and rule constraints is used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step S2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step S3, and the iteration ends;

[0273] In step S5, the specific optimization process for the soil pore reconstruction result set is as follows:

[0274] The porosity correction value p * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entroy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as:

[0275] Loss = -p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,...N

[0276] Taking the minimization of the loss function as the goal, the optimal result of soil pore reconstruction is obtained; let the porosity of the optimal result be p b .

[0277] Example:

[0278] Soil samples were collected from a pesticide-contaminated site in Binhai New Area. The contaminated soil layer was mainly clay-silt soil interaction. The pollutants in the site that exceeded the screening value for residential land included phenols and benzene series; the main contaminated layers were the first soil layer (0-1m) and the second soil layer (1-3.0m). Due to the long pollution time and the serious aging of the polluted organic matter, the method of the present invention was used to perform three-dimensional reconstruction of the collected sample soil.

[0279] Among them, the experimental equipment and specifications used are shown in Table 1:

[0280]

[0281]

[0282] The specific implementation steps are the same as the above specific implementation methods, namely:

[0283] Step 1: Determine the preset sampling points and reconstructed sampling points of the site, collect soil samples and take photos, and obtain the soil parameters of the preset sampling points and reconstructed sampling points respectively;

[0284] Step 2: Perform CT tomography scanning of the reconstructed soil sample, extract cross-sectional image features, and preprocess the sample point parameters. The preset sample point parameters are used as the training set and the reconstructed sample point parameters are used as the test set. The composite machine learning model is used to obtain the target parameter correction value.

[0285] Step 3: Based on the reconstructed soil sample cross-section image features, a feature cross-section subset is selected, and the porosity correction value is combined to generate and optimize the soil pore reconstruction result;

[0286] Step 4: Considering the pore deformation during the sampling process, the in-situ soil pore reconstruction model is obtained; combined with the organic matter volume fraction correction value, the soil organic matter occurrence reconstruction model is obtained, and then the soil three-dimensional reconstruction model is obtained;

[0287] Step 5: Combined with the soil pore connectivity assessment and considering the dynamic disturbance process of sampling, the soil 3D reconstruction model is improved and corrected to obtain the final reconstruction result.

[0288] Some experimental analysis results of the reconstructed soil sample #3 (soil layer depth 80-90 cm) collected from the contaminated site are shown in Table 2:

[0289]

[0290] As shown in Table 2, the measured porosity of the reconstructed soil sample #3 of the contaminated site is 0.47, and the measured organic matter volume fraction is 5.5%. After the soil is reconstructed in three dimensions by the method of the present invention, the porosity of the reconstructed result is 0.45, and the organic matter volume fraction is 5.2%, which is consistent with the measured values, proving that the model has a high reduction accuracy for soil pores and organic matter parameters.

Claims

1. A soil three-dimensional reconstruction method, characterized in that: The following steps are involved: Step 1: Determine the preset sampling points and reconstructed sampling points of the site, collect soil samples and take photos, and obtain the soil parameters of the preset sampling points and reconstructed sampling points respectively; Step 2: Perform CT tomography scanning of the reconstructed soil sample, extract cross-sectional image features, and preprocess the sample point parameters. The preset sample point parameters are used as the training set and the reconstructed sample point parameters are used as the test set. The composite machine learning model is used to obtain the target parameter correction value. Step 3: Based on the reconstructed soil sample cross-section image features, a feature cross-section subset is selected, and the porosity correction value is combined to generate and optimize the soil pore reconstruction result; Step 4: Analyze the pore deformation during the sampling process to obtain the in-situ soil pore reconstruction model; combine the organic matter volume fraction correction value to obtain the soil organic matter occurrence reconstruction model, and then obtain the soil three-dimensional reconstruction model; Step 5: Combined with the soil pore connectivity assessment and considering the dynamic disturbance process of sampling, the soil 3D reconstruction model is improved and corrected to obtain the final reconstruction result.

2. The method according to claim 1, characterized in that In step 1, the following sub-steps are included: Step 1-1: define the sample points for the required three-dimensional soil reconstruction as reconstruction sample points, the sample points for collecting soil parameters as preset sample points, determine the preset sample points, and select any point as the reconstruction sample point; Step 1-2: Use sampling equipment to collect soil samples from all sampling points and obtain soil resistance parameters during the uniform speed drilling sampling process; Step 1-3: After the collected soil samples are air-dried, all soil samples are photographed in a dark box with the same shooting parameters, exposure, and distance to obtain soil sample photo image information; It also includes steps 1-4: obtaining the measured values ​​of the preset soil sample target parameters, that is, measuring the porosity and organic matter volume fraction.

3. The method according to claim 1, characterized in that In step 2, it includes extracting features of the reconstructed soil sample cross-section image and obtaining the target parameter correction value; The steps of reconstructing the feature extraction of the soil sample cross-section image include: Step 2-1-1) Using a spiral CT scanner, place the reconstructed soil sample into the CT scanner, scan a cross section from the top to the bottom at a fixed distance, and scan at equal intervals to obtain M cross-sectional images of the soil sample; Step 2-1-2) After removing artifacts from the obtained soil sample cross-sectional image, a grayscale histogram is established according to the different grayscale values ​​of each component; the original X-CT image has many noise points, which greatly interferes with the extraction and analysis of the pore structure. Before threshold segmentation of the two-dimensional image, all soil sample cross-sectional images are processed with the non-local mean method to remove a large number of poor noise points that are too small to avoid interference with the pore structure calculation; The CT scan cross-sectional image obtained after Gaussian high-pass filtering is selected as the boundary condition required for soil three-dimensional reconstruction calculation. Its mathematical expression is: Where Q is the filtered image at the image coordinate (i, j); f(i, j) is the original pixel value at the image coordinate (i, j); k and l are the parameters of the matrix function in the x and y directions; K(k, l) is the filter kernel matrix function, in Gaussian high-pass filtering: In the filtered images, the pores and soil skeleton are well segmented, and the microstructural information is preserved to a large extent; Step 2-1-3) In the gray value histogram of each cross-sectional image of the reconstructed soil sample, the gray value of the middle trough between two peaks is selected as the threshold, and the image is segmented into three parts: pores, organic matter, and soil skeleton, so as to extract the features of the reconstructed soil sample cross-sectional image, separate the pores and soil skeleton parts in the image, and obtain M pore cross sections required for subsequent three-dimensional reconstruction; The acquisition of the target parameter correction value includes the following steps: Step 2-2-1) Preprocessing the soil sample parameter set in the same way for all sample points, including feature extraction of multi-color space of soil sample image and normalization of multi-type parameters; Step 2-2-2) Perform composite machine learning model training, hyperparameter adjustment and optimization output; Step 2-2-3) Use the composite machine learning model to predict the target parameters of the reconstructed sample points, perform secondary correction based on the correction model and measured data, and output the correction value.

4. The method according to claim 3, characterized in that In step 2-2-1, the steps of soil sample parameter set preprocessing are as follows: (1) Use the OpenCV software library to convert soil sample photo images into RGB, XYZ, HSV, HSL, Lab, and Luv color spaces, and extract feature information in multiple color spaces of the image; (2) principal component analysis is used to reduce the dimension of the obtained multi-color space feature information, that is, the correlation between the color space parameters is calculated, and the color parameters with correlation greater than 0.7 are removed to obtain the sub-color space containing the main color space feature information; (3) Using the Min-max normalization method, the sub-color space and soil resistance parameters are linearly transformed to map the range of all types of parameter sets to the interval [0,1].

5. The method according to claim 4, characterized in that In step 2-2-2, the steps of model training, hyperparameter adjustment and optimization output are as follows: (1) The preset sample point sub-color space and soil resistance parameters obtained after normalization in step 2-2-1 are used as input sets, wherein a certain proportion of the input set is randomly divided as a training set to complete the basic training of the composite machine learning model; the remaining proportion is used as a validation set to achieve hyperparameter tuning of the composite machine learning model; (2) The mean square error (MSE) function is used as the loss function for model training. The formula is: Where Y i , are the i-th value in the true value and predicted value vector respectively, and n is the vector dimension; (3) The training set is input into the composite machine learning model, and an output item is obtained through the first-layer support vector machine regression (SVR) operation; after the validation set is input into the SVR for hyperparameter tuning, the first-layer SVR operation is performed again using the same training set to obtain an optimized output item; (4) The output of the first optimization is used as the input of the second-layer minimum absolute shrinkage and selection operator Lasso regression algorithm to obtain the secondary output. Similar to the first layer, the validation set is input into the Lasso regression for hyperparameter tuning, and the output of the first optimization is used again to perform the second-layer Lasso regression operation to obtain the secondary optimization output. At this point, the training of the composite machine learning model is completed, and the secondary optimization output is the model output value.

6. The method according to claim 5, characterized in that In step 2-2-3, the steps of reconstructing the prediction and correction of the sample point target parameters are as follows: (1) The reconstructed sample point sub-color space and soil resistance parameters obtained after normalization in step 2-2-1 are used as the test set and input into the composite machine learning model; under the condition that the hyperparameters in step 2-2-2 have been tuned, the test set is subjected to the first layer SVR operation to obtain a first optimization output item; this item is used as the input of the second layer Lasso regression, and the second optimization output item obtained by the operation is the predicted value of the target parameter of the reconstructed sample point; (2) Use a convolutional neural network (CNN) model to correct the target parameters; input the test set into the pre-trained CNN model, perform a weighted linear combination of the CNN model output and the secondary optimization output to obtain a correction value. It can be expressed as: Where Z1 is the output item of secondary optimization; Z CNN is the output item of CNN model; ω1 is the output item weight of CNN model; The range of ω1 is set to 0.1-0.3, and the size of the training sample in the soil multi-source heterogeneous database and the similarity between the properties of the training sample soil sample and the reconstructed soil sample are used as the basis for weight calibration. The larger the size of the training sample and the greater the degree of similarity between the properties of the soil samples, the greater the weight given to the correction value. (3) Correction using the measured values ​​of preset sample points; After completing the measured values ​​of the preset soil sample target parameters in step 1-4, select n preset sample points closest to the reconstructed sample points as the correction sample point group, and correct the correction value obtained in step 2-2-1. Perform secondary deviation correction; Based on the inverse distance weighting method, the correction effect of the preset sample points is calibrated, that is, it is assumed that the correction effect of the preset sample points on the reconstructed sample points weakens with the increase of distance; The definition formula of the inverse distance weighted method is: Where Z2 is the measured correction value; n is the number of points in the correction sample group; Z i is the measured value of the target parameter of the i-th sample point in the correction sample point group; i is the weight coefficient of the sample point measurement value, which is given by the following formula: In the formula, Point to i b to the power of d i is the distance between the i-th preset sample point and the reconstructed sample point, and b is the distance attenuation coefficient; The first correction value Perform a weighted linear combination with the measured deviation correction value Z2 to output the secondary deviation correction value of the reconstructed sample point target parameter It can be expressed by the formula: In the formula, is the one-time correction value; Z2 is the measured correction value; ω2 is the weight of the measured correction value; the range of ω2 is set to 0.1-0.3, and the number of points n of the selected correction sample group and the similarity between the properties of the correction sample group and the reconstructed soil sample are used as the basis for weight calibration. The larger the n is, the greater the similarity of the physical and chemical properties between the soil samples is, and the greater the weight given to the correction value; (4) Output the primary or secondary correction value of the target parameter of the reconstructed sample point, where p * Indicates the porosity correction value, Indicates the correction value of organic matter volume fraction.

7. The method according to claim 1, characterized in that In step 3, the following sub-steps are included: Step 3-1: For the M reconstructed soil sample pore sections obtained in step 2-1-3, select m characteristic sections by a permutation and combination method based on certain rules, that is, Form a characteristic section subset; Step 3-2: Evaluate the pore connectivity of the characteristic cross-section subset to facilitate connectivity correction of the subsequent soil 3D reconstruction model; Step 3-3: Use the moving cube algorithm to generate the pore structure between the characteristic sections, and replace the generated pores with equivalent diameters to obtain the soil pore reconstruction result based on the characteristic sections; Step 3-4: Iteratively generate different characteristic section subsets according to the rules described in step 3-1 to obtain a soil pore reconstruction result set generated by different characteristic sections; Step 3-5: Combined with the porosity correction value obtained in step 2-2-3, the soil pore reconstruction result set is optimized to obtain the soil pore reconstruction optimization result.

8. The method according to claim 7, characterized in that In step 3-1, the specific rules for feature section subset selection and mutual iteration are as follows: (1) The characteristic section subset must include the top and bottom sections of the soil sample; (2) In the characteristic section subset, the distance between any set of adjacent sections does not exceed a certain multiple of the CT scan spacing; (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different; (4) Define the cross-sectional position P as the ratio of the distance between the cross-sectional area and the top of the soil mass to the height of the soil mass; In step 3-2, the pore connectivity of the characteristic cross section is evaluated as follows: For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices and serve as the basis for the correction of pore connectivity in the subsequent step 5; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula: Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices; When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected; In step 3-3, the steps for generating soil pore structure between characteristic cross sections are as follows; (1) For the first characteristic section subset obtained in step 3-2 The marching cubes algorithm is used to perform calculations with the top section of the reconstructed soil sample as the initial condition, and the calculation ends at the characteristic section adjacent to the top section, that is, the calculation of the first soil unit is completed; The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces; its relationship can be expressed as: F i,j,k =F(x i ,y j ,z k ) (i=1,2,…,N x ;j=1,2,…,N y ;k=1,2,…,N z ) In the formula, x i ,y j 、z k N is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are x, y, z coordinate pixels; (2) The pores generated by the soil unit were counted using Avizo software, sorted by pore volume, and classified into three categories according to the quantile classification method. The volume mean of each type of pore can be directly obtained from the software statistical results; By formula The equivalent diameter of each type of pore is calculated, where is the volume mean of all pores under the i-th type of pores, and all generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage; (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil; the calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset, and the three-dimensional visualization results can be drawn by image processing software such as Avizo; The total volume of all pores generated in this reconstruction result is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now: Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section; In step 3-4, the specific process is as follows: Based on the specific rules for selecting feature section subsets in step 3-1, a combination of permutation and combination iterative algorithms and rule constraints are used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step 3-2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step 3-3, and the iteration ends; In steps 3-5, the specific optimization process for the soil pore reconstruction result set is as follows: Take the porosity correction value p obtained in step 2-2-3 * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entroy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as: Loss=-p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,…N Taking the minimization of the loss function as the goal, the optimal result of soil pore reconstruction is obtained; let the porosity of the optimal result be p b .

9. The method according to claim 8, characterized in that In step 4, the following sub-steps are specifically included: 4-1) Carry out compression test on reconstructed soil samples, measure the porosity of soil samples under different loads, and draw compression relationship curve between load and porosity; 4-2) Estimate the load level of the reconstructed soil sample at the depth of the stratum, calculate the change in porosity of the soil sample during the process of transferring from the underground in-situ to the ground state according to the compression relationship curve, and use the porosity p of the optimal result to calculate the porosity of the soil sample. b Inversely infer the porosity p of the soil sample in the original state of the soil b ′, defines the compression factor 4-3) Multiply all pore volumes of the soil pore reconstruction optimization results in step 3-5 by S, that is, compress the pores equally by S times, so as to restore the pore structure of the soil sample before the pore deformation occurs during the sampling process, and obtain the in-situ soil pore reconstruction model. The porosity of this model is p b ′; (4) Considering that there may be overlap between pores and organic matter during CT scanning, the porosity p of the in situ soil pore reconstruction model is combined with the known b ′, the organic matter volume fraction correction value obtained in step 2.2-3 Corrected to The formula is: Where p b ′ is the porosity of the in-situ soil pore reconstruction model, is the correction value of organic matter volume fraction; 4-5) Considering the actual occurrence form of soil organic matter, a spherical shell model of soil organic matter occurrence is proposed: the occurrence form of organic matter is simplified to a spherical shell of a certain thickness outside a sphere equivalent to the pore structure; Under the condition that the equivalent size and number of all porous spheres are known, an equivalent thickness d is preset so that a spherical shell with a thickness of d is attached to the outside of each porous sphere to represent the organic matter in the reconstructed model; then the volume of organic matter on the sphere under a certain equivalent diameter is: The volume fraction of organic matter in the total volume of the soil sample on all equivalent spheres is: Among them, N i is the number of spheres contained in a certain equivalent size range, V0 is the total volume of the reconstructed soil; 4-6) Increase or decrease d within a specific range with a specific step size to obtain a set of soil organic matter storage reconstruction results under different equivalent thicknesses; iterative calculation is performed to obtain the soil organic matter storage reconstruction result set under each equivalent thickness d i The corresponding The corrected organic matter volume fraction is the expected value, calculated using the cross entropy loss function and expected value The difference between them, the loss function can be expressed as: Taking the minimization of the loss function as the goal, the optimal result of soil organic matter storage reconstruction is obtained; let the equivalent thickness of the optimal result be d b , the spherical shell of equivalent thickness is the soil organic matter storage reconstruction model; The in-situ soil pore reconstruction model and the soil organic matter storage reconstruction model were combined to obtain a three-dimensional soil reconstruction model that included pore structure and organic matter storage.

10. The method according to any one of claims 1 to 9, characterized in that In step 5, the following sub-steps are included: 5-1) Combined with the pore connectivity evaluation results in step 3-2, the spherical morphology of the pores and organic matter in the soil three-dimensional reconstruction model obtained in step 4 is artificially and volume-fixed longitudinally appropriately stretched, that is, the pore connectivity in the reconstruction model is corrected, so that the pore spaces that are far apart and separated from each other are brought closer to each other, and the soil pore connectivity structure and organic matter occurrence status in the actual scenario are as close as possible; 5-2) Considering the impact and changes of sampling disturbance on soil pore structure, the rationality of the soil three-dimensional reconstruction model after pore stretching is analyzed; considering the analysis of the dynamic process of sampling, including the disturbance of the soil when the sampling equipment drills the sample, and the creep that occurs when the soil is transferred from the underground burial position to the ground; in this way, based on the expert opinions and prior knowledge base, the soil three-dimensional reconstruction model obtained after the longitudinal stretching of the pores is artificially disturbed and corrected, and the corrected model is used as the final reconstruction result.

11. A method for generating and optimizing soil pore reconstruction results, characterized in that: The following steps are involved: Step S1: For the M reconstructed soil sample pore sections, m characteristic sections are selected by a permutation and combination method based on certain rules, namely Form a characteristic section subset; Step S2: evaluating the pore connectivity of a subset of characteristic sections; Step S3: using the marching cube algorithm to generate the pore structure between the characteristic sections, performing equivalent replacement of the generated pores with equivalent diameters, and then obtaining the soil pore reconstruction result generated based on the characteristic sections; Step S4: Iteratively generate different characteristic section subsets according to the rules described in step S1 to obtain a soil pore reconstruction result set generated by different characteristic sections; Step S5: Optimize the soil pore reconstruction result set in combination with the obtained porosity correction value to obtain the soil pore reconstruction optimization result.

12. The method according to claim 11, characterized in that In step S1, the specific rules for feature section subset selection and mutual iteration are as follows: (1) The characteristic section subset must include the top and bottom sections of the soil sample; (2) In the characteristic section subset, the distance between any group of adjacent sections does not exceed 5 times the CT scan spacing; (3) If any of the cross sections in the two sets of characteristic cross section subsets are different, the two cross section subsets are considered to be different; (4) Define the cross-sectional position P as the ratio of the distance between the cross-sectional area and the top of the soil mass to the height of the soil mass; In step S2, the pore connectivity of the characteristic cross section is evaluated as follows: (1) For each group of adjacent characteristic sections in the characteristic section subset, the double-section method is used to analyze and obtain the Euler characteristic number corresponding to the pore space, so as to evaluate the connectivity of the three-dimensional space composed of the two slices; the Euler characteristic number of the upper and lower sections X1 and X2 and their intersection image can be obtained by the following formula: Where E2 is the two-dimensional Euler characteristic number, which is obtained by subtracting the number of interconnected pores from the number of completely isolated pores in the two-dimensional image; V is the volume of the three-dimensional space formed by the two slices; When E v When E is less than 0, the soil pore space is considered to be a continuum. v When it is equal to 0, the pore space is considered to be completely disconnected; In step S3, the steps of generating the soil pore structure between characteristic sections are as follows; (1) For the first feature section subset obtained in step S2 Based on the pore cross section obtained after gray value segmentation, the marching cubes algorithm is used to perform calculations with the top cross section of the reconstructed soil sample as the initial condition, and the calculation ends at the characteristic cross section adjacent to the top cross section, that is, the calculation of the first soil unit is completed; The marching cube algorithm uses each subunit formed by the existing upper and lower sections as a small cubic grid, determines the internal isosurface of the voxel according to the value of each vertex of the grid voxel, and constructs a pore model by connecting the isosurfaces; its relationship can be expressed as: F i,j,k =F(x i ,y j ,z k ) (i=1,2,…,N x ;j=1,2,…,N y ;k=1,2,…,N z ) In the formula, x i ,y j 、z k N is the coordinate value of the x, y, and z axes of the three-dimensional data field; x 、N y 、N z are x, y, z coordinate pixels; (2) The pores generated by the soil unit were counted using Avizo software, sorted by pore volume, and classified into three categories according to the quantile classification method. The volume mean of each type of pore can be directly obtained from the software statistical results; By formula The equivalent diameter of each type of pore is calculated, where is the volume mean of all pores under the i-th type of pores, and all generated pores are replaced by pores of corresponding equivalent diameters to facilitate the subsequent three-dimensional reconstruction of soil organic matter storage; (3) For the entire reconstructed soil sample, the above-mentioned moving cube algorithm is used to iteratively perform the reconstruction operation of the pore structure of each soil unit from top to bottom, and the generated pores of each soil unit are expressed as pore equivalent diameters. The iterative process ends at the bottom section of the soil; the calculation results of continuous soil units are stacked to generate the soil pore reconstruction results under the characteristic section subset, and the three-dimensional visualization results can be drawn by image processing software such as Avizo; The total volume of all pores generated in this reconstruction result is divided by the volume of the soil sample to obtain the characteristic section subset The porosity p of the reconstructed result is i ,Right now: Where V por,i is the pore volume of the i-th soil unit, V0 is the total volume of the reconstructed soil sample, and m-1 is the number of soil units divided by the characteristic section; In step S4, the specific process is as follows: Based on the specific rules for selecting the feature section subset in step S1, a combination of permutation and combination iterative algorithm and rule constraints is used to find all feature section subsets that meet the rules and generate new feature section subsets. Iterate and solve each subset according to step S2 to obtain the corresponding soil pore reconstruction result, and calculate the porosity according to the formula in step S3, and the iteration ends; In step S5, the specific optimization process for the soil pore reconstruction result set is as follows: The porosity correction value p * is the expected value, for each characteristic section subset The porosity p of the reconstructed result is i ,i∈[1,2,3,…,N], and use the cross entroy loss function to calculate p i With the expected value p * The difference between them, the loss function can be expressed as: Loss=-p * lg(p i )-(1-p * )lg(1-p i ),i∈1,2,…N Taking the minimization of the loss function as the goal, the optimal result of soil pore reconstruction is obtained; let the porosity of the optimal result be p b .

Images