Intelligent design of loess water-blocking pile filling ratio and method for predicting impermeability

By mapping the proportion of loess water-blocking pile filler to an unconstrained Euclidean space through closed-loop computation and equidistant logarithmic ratio transformation, and combining it with genetic programming-driven symbolic regression, the problems of prediction bias and high experimental cost in the proportion design of loess water-blocking pile filler were solved, and accurate prediction of flowability and permeability coefficient was achieved.

CN122287404APending Publication Date: 2026-06-26LANZHOU INST OF TECH +1

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
LANZHOU INST OF TECH
Filing Date
2026-05-28
Publication Date
2026-06-26

Smart Images

  • Figure CN122287404A_ABST
    Figure CN122287404A_ABST
Patent Text Reader

Abstract

This invention discloses an intelligent design method for the proportioning and anti-seepage performance prediction of loess water-blocking pile filler, belonging to the field of data-driven modeling technology. The method includes the following steps: Step 1, collecting multiple sets of proportioning test records for loess water-blocking pile filler, and generating a perturbation sample template based on equidistant logarithmic ratio coordinate components; Step 2, using the input feature vector as input, constructing a dual-response analytical model for filler performance; Step 3, jointly generating a set of candidate proportioning points within the normalized solid composition space and the engineering allowable range of the water-solid ratio, and outputting the optimal proportioning scheme set for loess water-blocking pile filler under a preset total solid mass benchmark after inverse equidistant logarithmic ratio transformation. This invention achieves the establishment of an interpretable explicit mapping relationship between proportioning parameters and performance indicators, ensures predictive stability under raw material fluctuation conditions, and achieves multi-objective collaborative optimization of fluidity and anti-seepage performance.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to computer-aided engineering, and more particularly to the field of data-driven modeling and processing, specifically to a method for intelligent design of the mix proportion and prediction of the impermeability of loess water-blocking pile filler. Background Technology

[0002] Fluctuations in groundwater levels and surface water infiltration in loess regions often lead to increased soil moisture content and decreased bearing capacity. Loess water-blocking piles, as an effective underground water-blocking structure, form continuous or intermittent water-blocking curtains by setting low-permeability piles in the foundation, blocking the spread of water along horizontal or vertical migration paths. The filler material for loess water-blocking piles typically uses loess as the main matrix, mixed with industrial solid waste such as coal gangue and carbide slag as auxiliary cementitious or filling components, and chemically cured with cement-based or lime-based curing agents. A suitable amount of water is added and mixed to prepare a fluid or semi-fluid solidified filler material, which is then poured into pre-formed pile holes and cured to form a pile structure with certain strength and low permeability. The proportioning design of the filler material directly determines two core construction and service performance characteristics of the pile: first, fluidity, i.e., the flowability of the filler material in its unhardened state, affecting the compaction degree of grouting and the operability of construction; and second, permeability coefficient, i.e., the water permeability of the hardened pile under the action of a water head difference, determining the long-term water-blocking effect of the pile.

[0003] Currently, the mix design of loess water-blocking pile fill mainly relies on indoor orthogonal experiments or uniform design experiments. These methods involve conducting systematic mix design experiments by setting several factor level combinations, measuring the flowability and permeability coefficient under each combination, and then using range analysis or variance analysis to screen for the optimal level combination of each factor. While these methods are logically clear and widely used in engineering practice, they have significant limitations. First, orthogonal experiments or uniform design experiments can only cover a limited number of factor level combinations, searching only on a pre-set discrete level grid, and cannot conduct a refined exploration of the continuous mix design space, easily overlooking potential optimal mix proportions located between discrete levels. Second, range analysis and variance analysis essentially assume that the influence of each factor on performance indicators is independent and additive, making it difficult to capture complex nonlinear interactions between multiple components, such as the coupled effect of the synergistic pozzolanic reaction of coal gangue and carbide slag in an alkaline solidification environment on pore structure and permeability characteristics. Summary of the Invention

[0004] The purpose of this invention is to provide a method for intelligent design of the mix proportion and prediction of the impermeability of loess water-blocking pile filler, comprising the following steps: Step 1: Collect multiple sets of mix design test records for loess water-blocking pile filler. Each set of mix design test records includes four solid component variables: loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction, as well as the water-solid ratio and the corresponding fluidity and permeability coefficient values. Take the common logarithm of the permeability coefficient value to obtain the logarithmic permeability coefficient value. Perform zero-value substitution preprocessing, closure operation, and equidistant logarithmic ratio transformation on the four solid component variables in sequence to obtain the equidistant logarithmic ratio coordinate vector. Concatenate the equidistant logarithmic ratio coordinate vector with the water-solid ratio to form the input feature vector and stack them row by row to form the input feature matrix. Generate a perturbation sample template based on the equidistant logarithmic ratio coordinate components. Step 2: Using the input feature vector as input and the flowability value and log permeability coefficient value as output, the candidate expression tree is evaluated simultaneously by genetic programming-driven symbolic regression, which evaluates the fitting error, complexity index and component perturbation robustness index. After non-dominated sorting-assisted iterative evolution and screening from the non-dominated set, the explicit prediction formulas for flowability and log permeability coefficient are determined, thus forming a dual-response analytical model for packing performance. Step 3: Within the normalized solid composition space and the engineering allowable range of water-solid ratio, a set of candidate mix proportions is jointly generated. After equidistant logarithmic ratio transformation, the mix proportions are input into the dual-response analytical model of the filler performance to obtain the predicted values ​​of fluidity and logarithmic permeability coefficient. Candidate mix proportions with predicted fluidity values ​​exceeding the preset allowable fluidity range are eliminated. For the remaining candidate mix proportions, a non-dominated sorting is performed with the dual objectives of maximizing the predicted fluidity value and minimizing the predicted logarithmic permeability value to extract the Pareto optimal mix proportion set. After equidistant logarithmic ratio inverse transformation, the optimal mix proportion scheme set of loess water-blocking pile filler is output under the preset total solid mass benchmark.

[0005] Furthermore, the test records for each group of proportions were obtained under the same curing period, the same specimen preparation method, and the same penetration test conditions.

[0006] Furthermore, in step 1, the zero-value substitution preprocessing process is as follows: when the value of any solid component variable is less than the preset positive lower limit, it is substituted with the preset positive lower limit; the closed operation process is as follows: the four solid component variables after the zero-value substitution preprocessing are arranged into a 4-dimensional component vector in a fixed order of loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction; the four solid component variables are divided by the sum of the four solid component variables respectively, so that the sum of the four solid component variables after the closed operation is equal to 1, and the normalized solid component vector is obtained.

[0007] Further, in step 1, the process of the equidistant logarithmic ratio transformation is as follows: Construct an orthogonal basis matrix with 3 rows and 4 columns; the first row of the orthogonal basis matrix contains 4 elements, where the first element is positive, the second to fourth elements are equal and all negative, the sum of the squares of all elements in the first row is equal to 1 and the sum of all elements is equal to 0; the first element of the second row is 0, the second element is positive, the third to fourth elements are equal and all negative, the sum of the squares of all elements in the second row is equal to 1 and the sum of all elements is equal to 0; the first to second elements of the third row are both 0, the third element is positive, the fourth element is negative, the sum of the squares of all elements in the third row is equal to 1 and the sum of all elements is equal to 0; take the natural logarithm of each element in the normalized solid component vector to obtain a logarithmic vector, and perform matrix multiplication between the orthogonal basis matrix and the logarithmic vector to obtain a 3D equidistant logarithmic ratio coordinate vector.

[0008] Furthermore, in step 1, the process of generating the perturbation sample template is as follows: for the logarithmic ratio coordinate components of each sample in the input feature matrix, a preset number of perturbation vectors are generated using a preset fixed random seed. Each component of each perturbation vector independently follows a normal distribution with a mean of 0 and a standard deviation of a preset perturbation amplitude. Each perturbation vector is superimposed on the logarithmic ratio coordinate components of the corresponding sample while keeping the water-solid ratio component unchanged, thus forming a fixed perturbation sample set for the corresponding sample. The fixed perturbation sample sets of all samples are then summarized into a perturbation sample template.

[0009] Furthermore, in step 2, the candidate operator set includes five operators: addition, subtraction, multiplication, protected division, and protected logarithm operation. Protected division sets the output value to 1 when the absolute value of the divisor is less than a preset minimum threshold. Protected logarithm operation replaces the input value with the preset positive threshold before performing the natural logarithm operation when the input value is less than or equal to a preset positive threshold. The candidate terminal set includes equidistant logarithmic ratio coordinate components, water-solid ratio variables, and random constants.

[0010] Furthermore, in step 2, the fitting error is the sum of the absolute values ​​of the differences between the expression output value of each sample and the corresponding output variable after substituting all samples in the input feature matrix into the candidate expression tree; the complexity index is the total number of all nodes in the candidate expression tree; the calculation process of the component perturbation robustness index is as follows: take out the fixed perturbation sample set of each sample from the perturbation sample template, substituting each fixed perturbation sample in the fixed perturbation sample set into the candidate expression tree to calculate the output value, and statistically analyze the variance of all perturbation output values ​​of the same sample as the single sample perturbation variance, and take the mean of the single sample perturbation variance of all samples.

[0011] Furthermore, in step 2, the process of non-dominated sorting-assisted iterative evolution is as follows: non-dominated sorting is performed on the current candidate expression population with three objectives: minimizing fitting error, minimizing complexity index, and minimizing component perturbation robustness index, and a non-dominated level is assigned to each candidate expression tree; each time, a fixed number of candidate expression trees are randomly selected from the candidate expression population to form a tournament group, and the candidate expression tree with the best non-dominated level is selected from the tournament group as the parent expression tree. If there are multiple candidate expression trees with the same non-dominated level in the tournament group, the candidate expression tree with the smallest fitting error is selected as the parent expression tree.

[0012] Furthermore, after obtaining two parent expression trees, a subtree crossover operation is performed. In the first parent expression tree, a non-root node and a subtree rooted at the selected non-root node are randomly selected. In the second parent expression tree, a non-root node and a subtree rooted at the selected non-root node are also randomly selected. The two subtrees are then interchanged and grafted onto each other's parent expression trees to generate two descendant expression trees. Point mutation operations are performed on the descendant expression trees with a preset mutation probability. One node is randomly selected from the descendant expression trees. If the selected node is an operator node, it is replaced with another operator from the candidate operator set. If the selected node is a terminal node, it is replaced with another terminal from the candidate terminal set. After the iteration, the candidate expression tree with the smallest fitting error and a complexity index not exceeding the preset node limit is selected from the non-dominated set. The analytical expression corresponding to the selected candidate expression tree is determined as the explicit prediction formula.

[0013] Furthermore, in step 3, the process of generating the candidate mix proportion set is as follows: For the first three of the four solid component variables, set upper and lower bounds for normalized values ​​respectively. Insert multiple sampling nodes at equal intervals between the upper and lower bounds of each variable and perform a full combination. In each combination, the normalized value of the fourth solid component variable is determined by subtracting the sum of the normalized values ​​of the first three variables from 1. Eliminate combinations where the normalized value of the fourth solid component variable is negative or exceeds the allowable range for engineering purposes to obtain a normalized solid component candidate set. Combine the normalized solid component candidate set with the water-to-solid ratio sampling nodes to generate the candidate mix proportion set. In the non-dominated sorting, compare any two candidate mix proportion points one by one. If the predicted flowability value of the first candidate mix proportion point is larger... If the predicted flowability of the first candidate mix point is equal to the predicted flowability of the second candidate mix point, and the predicted logarithmic permeability coefficient of the first candidate mix point is less than or equal to the predicted logarithmic permeability coefficient of the second candidate mix point, and at least one target is strictly better than the second candidate mix point, then the first candidate mix point dominates the second candidate mix point. The candidate mix points that are not dominated by any other candidate mix points are extracted as the Pareto optimal mix set. The process of the inverse transformation of the equidistant logarithmic ratio is as follows: the transpose of the orthogonal basis matrix is ​​multiplied by the equidistant logarithmic ratio coordinate vector to obtain a 4-dimensional logarithmic vector. The natural exponent is taken for each element in the 4-dimensional logarithmic vector, and all elements are divided by the sum of all elements to complete the closed normalization. Each normalized value is multiplied by the preset total solid mass benchmark to obtain the mass fraction.

[0014] The intelligent design method for the proportioning and anti-seepage performance prediction of loess water-blocking pile fill material of the present invention has the following beneficial effects: First, this invention treats the four solid component variables—loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction—as closed-loop component data. Through closed-loop operations and equidistant logarithmic ratio transformation, these variables are mapped to an unconstrained Euclidean space. This eliminates the spurious correlations between components caused by the constraint of constant total sum, enabling subsequent modeling to accurately identify the true independent effects of each solid component on fluidity and impermeability. This avoids the prediction bias caused by closed-loop interference when using the original mass fractions directly for modeling in traditional methods.

[0015] Secondly, this invention uses genetic programming-driven symbolic regression to establish explicit prediction formulas for mobility and log permeability coefficients. The output results are analytical expressions that can be directly written and interpreted. Engineering technicians can complete performance prediction and mix design verification by manual calculation on the construction site without relying on professional calculation software. This overcomes the limitations of black box models such as neural networks and support vector machines, which are opaque, uninterpretable, and difficult to apply on site.

[0016] Third, this invention introduces a component perturbation robustness index into the fitness evaluation of genetic programming. Together with the fitting error and complexity index, it forms a three-criteria non-dominated sorting screening mechanism. In the model screening stage, engineering uncertainties such as raw material weighing error and batch fluctuation are taken into consideration, ensuring that the final selected explicit prediction formula can still give stable and reliable prediction results when the input components are slightly shifted, thus improving the transfer reliability of the model from the laboratory to the construction site.

[0017] Fourth, this invention performs non-dominated sorting with the dual objectives of maximizing the predicted flowability value and minimizing the predicted logarithmic permeability coefficient value. At the same time, it sets a preset flowability allowable interval to constrain the flowability. The output Pareto optimal mix design set systematically reveals the trade-off boundary between flowability and impermeability. Engineering technicians can flexibly select from it according to the water blocking requirements and pumping conditions of the specific construction site. A set of optional mix design schemes that are not mutually dominant in the two performance dimensions can be obtained without repeated experiments, which significantly reduces the experimental cost and time consumption of mix design optimization. Attached Figure Description

[0018] Figure 1 This is a schematic diagram illustrating the principle of equidistant logarithmic ratio transformation of closed-type component data provided in an embodiment of the present invention. Figure 2 This is a schematic diagram of the genetic programming expression tree structure and subtree crossover operation provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of three-criteria non-dominated sorting and expression filtering provided in an embodiment of the present invention. Detailed Implementation

[0019] A method for intelligent design of loess water-blocking pile fill material ratio and prediction of impermeability includes the following steps: Step 1: Collect multiple sets of mix design test records for loess water-blocking pile filler. Each set of mix design test records includes four solid component variables: loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction, as well as the water-solid ratio and the corresponding fluidity and permeability coefficient values. Take the common logarithm of the permeability coefficient value to obtain the logarithmic permeability coefficient value. Perform zero-value substitution preprocessing, closure operation, and equidistant logarithmic ratio transformation on the four solid component variables in sequence to obtain the equidistant logarithmic ratio coordinate vector. Concatenate the equidistant logarithmic ratio coordinate vector with the water-solid ratio to form the input feature vector and stack them row by row to form the input feature matrix. Generate a perturbation sample template based on the equidistant logarithmic ratio coordinate components. Step 2: Using the input feature vector as input and the flowability value and log permeability coefficient value as output, the candidate expression tree is evaluated simultaneously by genetic programming-driven symbolic regression, which evaluates the fitting error, complexity index and component perturbation robustness index. After non-dominated sorting-assisted iterative evolution and screening from the non-dominated set, the explicit prediction formulas for flowability and log permeability coefficient are determined, thus forming a dual-response analytical model for packing performance. Step 3: Within the normalized solid composition space and the engineering allowable range of water-solid ratio, a set of candidate mix proportions is jointly generated. After equidistant logarithmic ratio transformation, the mix proportions are input into the dual-response analytical model of the filler performance to obtain the predicted values ​​of fluidity and logarithmic permeability coefficient. Candidate mix proportions with predicted fluidity values ​​exceeding the preset allowable fluidity range are eliminated. For the remaining candidate mix proportions, a non-dominated sorting is performed with the dual objectives of maximizing the predicted fluidity value and minimizing the predicted logarithmic permeability value to extract the Pareto optimal mix proportion set. After equidistant logarithmic ratio inverse transformation, the optimal mix proportion scheme set of loess water-blocking pile filler is output under the preset total solid mass benchmark.

[0020] The technical details of step 1 are explained below.

[0021] In the mix design of loess water-blocking pile filler, the filler is typically prepared by mixing four solid raw materials—loess, coal gangue, carbide slag, and curing agent—in a specific mass ratio, with the addition of an appropriate amount of water. The combination of different raw material proportions directly determines the flowability and impermeability of the filler after hardening; therefore, it is necessary to obtain basic data through systematic indoor mix design tests. When collecting multiple sets of mix design test records for loess water-blocking pile filler, each set of test records was obtained under the same curing age, the same specimen preparation method, and the same permeability test conditions to ensure that the performance differences between different mix proportions are caused only by the mix proportions themselves, rather than by fluctuations in test conditions. Each set of test records includes four solid component variables: the mass fractions of loess, coal gangue, carbide slag, and curing agent, the water-to-solid ratio, and the corresponding flowability and permeability coefficient values. The water-to-solid ratio, which is the ratio of the mass of added water to the total mass of the four solid component variables, is an independent process variable and does not have a total closed constraint relationship with the four solid component variables. The fluidity value was determined according to the current cement mortar fluidity test method, and the permeability coefficient value was determined according to the current geotechnical test method using variable head or constant head permeability tests. In one specific embodiment, a total of 60 sets of mix proportion test records were collected. The mass fractions of loess ranged from 30 to 70, coal gangue from 10 to 40, carbide slag from 5 to 25, curing agent from 3 to 15, and water-to-solid ratio from 0.4 to 0.8. All mass fractions are expressed in mass units. In another optional embodiment, the number of mix proportion test records can be any integer between 30 and 200, with the specific number determined according to project requirements and test conditions.

[0022] The logarithmic permeability coefficient is obtained by taking the common logarithm of the permeability coefficient values ​​in each set of test records. The permeability coefficient often varies across multiple orders of magnitude under different mixing conditions; for example, from... centimeters per second to If the raw permeability coefficient value is directly used as the target variable for subsequent modeling (centimeters per second), the numerical distribution will exhibit a severe right skewness. Large numerical samples will dominate the model fitting process, leading to a significant decrease in prediction accuracy in smaller numerical intervals. Taking the common logarithm maps the raw permeability coefficient value to an approximately uniformly distributed numerical interval; for example, the logarithmic permeability coefficient value corresponding to the above range is... arrive The numerical spacing between samples tends to be balanced, which is beneficial for subsequent symbolic regression to obtain an explicit prediction formula with stable fitting accuracy across the entire ratio domain. Specifically, the log penetration coefficient value... The calculation method is as follows ,in For the first The permeability coefficient value recorded in the mix design test. This represents common logarithmic operations with base 10. This is the serial number of the mixing ratio test record.

[0023] Subsequently, zero-value substitution preprocessing was performed sequentially on the four solid component variables. In actual mix proportioning experiments, some combinations may set the dosage of a certain solid raw material to 0 or an extremely small value; for example, the mass fraction of the curing agent may be 0 when exploring the mix proportion of pure loess. However, the subsequent closed-loop operation and equidistant logarithmic ratio transformation both involve logarithmic operations. The logarithmic function requires the input value to be strictly positive. When any solid component variable is 0 or close to 0, the logarithmic operation will fail or produce extreme values. Therefore, when the value of any solid component variable is less than a preset positive lower limit, it is replaced by the preset positive lower limit. In one specific embodiment, the preset positive lower limit is 0.001, that is, when the value of a solid component variable is less than 0.001, it is replaced with 0.001. In another optional embodiment, the preset positive lower limit can be 0.01 or 0.0001. The specific value is determined according to the engineering precision of the raw material dosage, as long as the value after substitution is much less than the lower limit of the normal dosage range, it will not cause substantial deviations in subsequent modeling.

[0024] After zero-value substitution pretreatment, a closed-loop operation was performed on the four solid component variables. The four solid component variables—loess, coal gangue, carbide slag, and solidifying agent—together constitute the composition system of the filler's solid portion. There is a natural summation constraint relationship among them; that is, the core information of interest is the relative proportion of each component in the total solids, rather than its absolute mass value. Different test batches may use different total solids baselines; for example, one batch might use 100 grams of solids as the baseline, while another uses 200 grams. However, as long as the mass proportions of the four components are the same, the performance of the prepared filler will be consistent. The purpose of the closed-loop operation is to eliminate the influence of differences in the total solids baseline and unify all proportion test records to the same relative proportion expression. The four solid component variables after zero-value replacement preprocessing are arranged into a 4-dimensional component vector in a fixed order: loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction. Each of the four solid component variables is then divided by the sum of the four solid component variables, ensuring that the sum of the four solid component variables after the closed-loop operation equals 1, resulting in the normalized solid component vector. (The text then repeats the first two sentences.) Taking the group ratio test record as an example, let the zero-value substitution be the four solid component variables after pretreatment as follows: , , , The components of the normalized solid composition vector obtained after closed-loop calculation are respectively the mass fractions of loess, coal gangue, carbide slag, and solidifying agent. ,in It is the sum of four solid component variables. Take an integer from 1 to 4. After the closing operation, Each component Indicates the first The mass proportion of a certain solid raw material in the total solids.

[0025] refer to Figure 1 , Figure 1 The diagram shows the principle of logarithmic ratio transformation for closed-type component data. Figure 1 The left side illustrates the distribution of the normalized solid component vector after the four solid component variables are closed and constrained within the three-dimensional simplex. Due to the closure constraint that the sum equals 1, there are non-independent coupling relationships between the components. Figure 1 The right side illustrates that after performing an equidistant logarithmic ratio transformation on the normalized solid component vector, the original 4-dimensional vector, which was subject to closure constraints, is mapped to a 3-dimensional real vector without closure constraints. Each coordinate component takes independent values ​​in the real domain, making it suitable as an input variable for subsequent symbolic regression modeling.

[0026] While normalized solid component vectors unify the proportional representation, they remain closed-form component data, meaning each component is constrained within a simplex, and their values ​​are not independent. This closed constraint can lead to spurious correlations in traditional statistical analysis; for example, when the proportion of one component increases, the proportions of other components are inevitably compressed, even if they have no real causal relationship in engineering. To eliminate spurious correlations caused by closed constraints, subsequent modeling is performed in unconstrained Euclidean space, applying an isometric logarithmic ratio transformation to each normalized solid component vector. The isometric logarithmic ratio transformation is a standard mapping method in component data analysis. Its core principle is to project the component data on the simplex into a lower-dimensional real space using a set of orthogonal bases. The projected coordinate components are independent of each other while preserving the geometric distance relationships of the original component data.

[0027] The specific process of the equidistant logarithmic ratio transformation is as follows. Construct an orthogonal basis matrix with 3 rows and 4 columns, denoted as . The first row of an orthogonal basis matrix contains four elements: the first element is positive, and the second through fourth elements are all negative. The sum of the squares of all elements in the first row is 1, and the sum of all elements is 0. The first element of the second row is 0, the second element is positive, and the third and fourth elements are all negative. The sum of the squares of all elements in the second row is 1, and the sum of all elements is 0. The first and second elements of the third row are both 0, the third element is positive, and the fourth element is negative. The sum of the squares of all elements in the third row is 1, and the sum of all elements is 0. These construction rules ensure that the orthogonal basis matrix... The inner product between any two rows is 0, ensuring the orthogonality between the coordinate components after transformation. In a specific implementation, the orthogonal basis matrix... The elements in the first row are as follows , , , The above constraints are not satisfied, therefore the actual Helmert-type orthogonal basis construction method used is: The... The first line The elements are , No. The first to the fourth element are all The first to the second All elements are 0, among which Take an integer from 1 to 3. This Helmert-type construction method automatically satisfies all constraints: the sum of squares in each row equals 1, the sum of elements in each row equals 0, and the rows are orthogonal. In another optional implementation, the orthogonal basis matrix can also be constructed using other orthogonal matrix construction methods that satisfy the above constraints, as long as the number of rows is 3, the number of columns is 4, the sum of squares in each row is 1, the sum of elements in each row is 0, and the rows are orthogonal.

[0028] Taking the natural logarithm of each element in the normalized solid composition vector yields the logarithmic vector. Taking a group as an example, the logarithmic vector is ,in Represents the natural logarithm operation. This represents the transpose. Since each component in the normalized solid component vector is strictly positive after zero-value substitution preprocessing, the natural logarithmic operation is always meaningful. This relates to the orthogonal basis matrix. With logarithmic vector Performing matrix multiplication yields a 3-dimensional logarithmic coordinate vector. ,in It contains 3 components, denoted as follows: , , After the equidistant logarithmic ratio transformation, the originally closed-loop 4-dimensional normalized solid component vector is mapped to a 3-dimensional unconstrained real number vector. The three components can take values ​​freely in the entire real number domain and are independent of each other, making them suitable as input variables for subsequent symbolic regression.

[0029] The 3D logarithmic ratio coordinate vector recorded in each batch of experiments is concatenated with the corresponding water-to-solid ratio to form a 4D input feature vector. The water-to-solid ratio, as an independent process variable, does not participate in closed-loop operations or logarithmic ratio transformations; instead, it is directly appended to the logarithmic ratio coordinate vector in its original numerical form. Taking the first... Taking the group as an example, let the water-to-solid ratio be... Then the 4-dimensional input feature vector is The first three components are from the output of the logarithmic ratio transformation, and the fourth component is the original value of the water-to-solid ratio. The four-dimensional input feature vectors from all the proportioning test records are stacked row-wise to form the input feature matrix, denoted as... Input feature matrix The number of rows equals the total number of sets of matching test records, and the number of columns equals 4. In the specific implementation described above, the input feature matrix... It is a matrix with 60 rows and 4 columns.

[0030] After constructing the input feature matrix, a perturbation sample template is generated based on the equidistant logarithmic ratio coordinate components. The purpose of the perturbation sample template is to evaluate the sensitivity of each candidate expression tree to small fluctuations in the input components during subsequent genetic programming. In engineering practice, there is always a certain deviation between the actual dosage of raw materials and the designed ratio. For example, weighing errors, uneven mixing, and batch fluctuations of raw materials can all cause the actual solid component ratio to deviate from the design value. If a prediction formula is too sensitive to small changes in components, even if it has a high fitting accuracy on the training data, it may produce unreliable prediction results in actual engineering applications due to fluctuations in raw materials. Therefore, local stability needs to be considered during the model selection stage, and the perturbation sample template serves this purpose.

[0031] The generation process of the perturbation sample template is as follows. For the logarithmic ratio coordinate components of each sample in the input feature matrix, a preset number of perturbation vectors are generated using a preset fixed random seed. The purpose of using a preset fixed random seed is to ensure that the perturbation vector corresponding to each sample remains unchanged throughout the entire genetic programming iteration process, avoiding different robustness evaluations of the same candidate expression tree in different iteration generations due to the randomness of the perturbation, thereby ensuring the stability and repeatability of the fitness index. Each component of each perturbation vector independently follows a normal distribution with a mean of 0 and a standard deviation of a preset perturbation amplitude. In one specific implementation, the preset perturbation amplitude is set to 0.05, and the preset number is set to 20, that is, 20 perturbation vectors are generated for each sample. In another optional implementation, the preset perturbation amplitude can be a value between 0.01 and 0.1, and the preset number can be an integer between 10 and 50. The larger the preset perturbation amplitude, the wider the range of fluctuation of the components under consideration; the more preset numbers, the better the statistical stability of the robustness evaluation, but the computational cost also increases accordingly. Each perturbation vector is superimposed onto the equidistant logarithmic ratio coordinate components of the corresponding sample, while keeping the water-to-solid ratio component unchanged, thus forming a fixed perturbation sample set for the corresponding sample. Taking the first... The first group of samples perturbation vector For example, the perturbed 4D input feature vector is ,in , , The first The perturbation values ​​of the perturbation vector on the three equally spaced logarithmic ratio coordinate components. Take an integer from 1 to a preset number. The perturbation is applied only to the logarithmic ratio component while keeping the water-to-solid ratio component unchanged. This is because the logarithmic ratio component reflects the relative proportions between solid components, and raw material weighing errors and uneven mixing primarily affect this part. The water-to-solid ratio is usually controlled by metering water addition, which has relatively high precision and belongs to a different dimension of process parameters. The fixed perturbation sample set of all samples is compiled into a perturbation sample template. This template remains unchanged throughout the entire genetic programming iteration process and can be repeatedly used. In the specific implementation described above, each of the 60 sample groups generates 20 perturbation samples, and the perturbation sample template contains a total of 1200 perturbed 4D input feature vectors.

[0032] The technical details of step 2 are explained below.

[0033] After completing step 1, the input feature matrix has been obtained. The goal of step 2 is to find an explicit mathematical expression between the input feature vector and the performance indicators, rather than training a black box model with fixed parameters but an opaque structure. The core consideration for using an explicit prediction formula instead of a black box model such as a neural network or support vector machine is that the mix design of loess water-blocking piles ultimately needs to be used by on-site engineering technicians. These technicians need to be able to directly see the functional relationship between the mix parameters and performance, thereby determining how changes in the amount of a certain raw material will affect the flowability or impermeability, and to quickly verify this through manual calculations on the construction site where computing equipment is lacking. A black box model cannot meet this requirement, while an explicit prediction formula can be directly written into the construction technical specifications, facilitating project implementation and third-party review.

[0034] Genetic programming-driven symbolic regression is a method that represents mathematical expressions using a tree structure and searches for the optimal expression by simulating the natural evolutionary process. Unlike traditional regression analysis, which predefines the function form, symbolic regression makes no assumptions about the function structure. Instead, it allows the algorithm to freely combine, compete, and eliminate elements within an expression space composed of basic operators and variables, ultimately resulting in an analytical formula with high fitting accuracy and a concise structure. This "let the data speak for itself" approach is particularly suitable for multi-component mixed material systems such as loess water-blocking pile fillers, because the interactions between the components are complex and there is a lack of mature theoretical models to predefine the function form.

[0035] Using the input feature vector as input, we first establish an explicit prediction formula for mobility with the mobility value as output. The specific process is as follows.

[0036] The candidate operator set includes five operators: addition, subtraction, multiplication, protected division, and protected logarithm. The meanings of addition, subtraction, and multiplication are consistent with regular mathematical operations. Protected division sets the output value to 1 when the absolute value of the divisor is less than a preset minimum threshold. The reason for introducing this protection mechanism is that genetic programming randomly generates various expression structures during the search process, some of which may produce divisors approaching 0. Without protection, the division operation would produce infinity or overflow, causing fitness evaluation to collapse and interrupting the entire evolutionary process. Setting the output to 1 instead of 0 when the absolute value of the divisor is less than the preset minimum threshold is because 1, as the multiplication unit, does not amplify or diminish the numerical propagation to other branches of the expression tree, thus limiting the impact of abnormal situations to a minimum. In one specific implementation, the preset minimum threshold is set to... In another alternative implementation, the preset minimum threshold can be set to a value of [value missing]. to The protected logarithm operation replaces the input value with a preset positive threshold before performing the natural logarithm operation if the input value is less than or equal to a preset positive threshold. The natural logarithm function requires the input to be strictly positive, while the coordinate components of the logarithm ratio can take negative values ​​or values ​​close to 0. Intermediate results during the genetic programming search process may also be non-positive. Directly performing the natural logarithm operation on non-positive numbers will produce undefined results. The protected logarithm operation avoids this problem by replacing non-positive inputs with a preset positive threshold. In one specific implementation, the preset positive threshold is set to... In another alternative implementation, the preset positive threshold can be set to a value of [value missing]. to The values ​​between these two values. These five operators were chosen from the candidate operator set because addition, subtraction, multiplication, and division constitute the complete set of basic arithmetic operations, while the natural logarithm operation can express power-law relationships and logarithmic decay, which are common nonlinear laws in materials science. The combination of these five operators is sufficient to cover the function structure of most engineering empirical formulas.

[0037] The candidate terminal set is defined as including log-ratio coordinate components, a water-to-solid ratio variable, and a random constant. There are three log-ratio coordinate components, corresponding to the first, second, and third components of the 3D log-ratio coordinate vector obtained after the log-ratio transformation in step 1. The water-to-solid ratio variable is the fourth component of the input feature vector. The random constant is a real value randomly generated within a preset range, randomly assigned during expression tree initialization and updated through mutation operations during subsequent evolution. In one specific implementation, the generation range of the random constant is... to A uniform distribution between them. In another alternative implementation, the range of the generated random constants can be adjusted according to the magnitude of the target variable.

[0038] Candidate expressions are represented using a tree structure. In this tree structure, internal nodes store operators from the candidate operator set, and leaf nodes store terminals from the candidate terminal set. For example, if the root node of an expression tree is the addition operator, the left subtree is the multiplication operator connecting the first component of the equidistant logarithmic ratio coordinate vector and the random constant 2.3, and the right subtree is the water-solid ratio variable, then the mathematical expression represented by this expression tree is: ,in The first component of the logarithmic ratio coordinate vector is the equidistant logarithmic ratio. Let be the water-to-solid ratio variable. More complex expression trees can nest multiple levels of operators and terminals, thus representing analytical formulas of arbitrary complexity.

[0039] refer to Figure 2 , Figure 2 The diagram shows the structure of the genetic programming expression tree and the subtree crossover operation. Figure 2 This illustrates the process of selecting a subtree rooted at a non-root node from the first parent expression tree and another subtree rooted at a non-root node from the second parent expression tree, then swapping and grafting these two subtrees into the other parent expression tree to generate a descendant expression tree. The subtree crossover operation recombines local structural fragments from two different expression trees, resulting in a descendant expression tree with new function forms not found in either parent expression tree. This is beneficial for expanding the search range of the expression space in genetic programming.

[0040] A first-generation candidate expression population is generated through random initialization. During initialization, each candidate expression tree is randomly constructed using either the growth method or the full-tree method. The growth method selects an operator or a terminal at each node with a certain probability, resulting in trees of varying depths and diverse structures. The full-tree method selects an operator at all nodes before reaching a specified depth, selecting a terminal only at the maximum depth level, thus generating a more regular tree structure. Combining the two methods contributes to the diversity of the initial population structure. In one specific implementation, the candidate expression population consists of 500 expression trees with an initial maximum depth of 6 levels, with the growth method and the full-tree method each accounting for 50%. In another optional implementation, the population size can be an integer between 200 and 2000, and the initial maximum depth can be an integer between 4 and 8.

[0041] After the first generation of candidate expression population is generated, three fitness indicators are calculated simultaneously for each candidate expression tree.

[0042] The first term is the fitting error. Each sample in the input feature matrix is ​​substituted into the candidate expression tree. For each sample, the components of its 4-dimensional input feature vector replace the corresponding terminal variables in the expression tree. Operator operations are performed layer by layer along the tree structure from the leaf nodes to the root node, finally obtaining the expression output value for that sample at the root node. The fitting error is defined as the sum of the absolute values ​​of the differences between the expression output values ​​of all samples and their corresponding mobility values, denoted as . ,in This represents the total number of groups recorded in the mixing experiment. For the first The expression output value obtained after substituting the group of samples into the candidate expression tree. For the first The mobility values ​​corresponding to the sample groups. A smaller fitting error indicates a better fit of the candidate expression tree to the training data. The sum of absolute values, rather than the sum of squares, is used as the error metric because absolute error is less sensitive to individual extreme bias samples than square error, which is beneficial for obtaining robust expressions across the entire proportion domain. In the specific implementation described above, It equals 60.

[0043] The second item is the complexity metric. The total number of nodes in the candidate expression tree is used as the complexity metric, denoted as . The total number of nodes includes the sum of internal operator nodes and leaf terminal nodes. The purpose of introducing a complexity metric is to control the complexity of the expression and avoid genetic programming searching for expressions with extremely low fitting errors but excessively large structures. Overly complex expressions are not only difficult for engineers to understand and verify, but also prone to overfitting, i.e., performing well on training data but failing to predict accurately on new ratios. In one specific implementation, the preset upper limit for the number of nodes is set to 30, meaning that only candidate expression trees with a total number of nodes not exceeding 30 are considered in the final screening stage. In another optional implementation, the preset upper limit for the number of nodes can be an integer between 20 and 50.

[0044] The third item is the component perturbation robustness index. This index measures the sensitivity of the candidate expression tree's output to small fluctuations in the isochronous logarithmic ratio coordinate components. A fixed perturbation sample set is taken from the perturbation sample template for each sample, and each fixed perturbation sample in the fixed perturbation sample set is substituted into the candidate expression tree one by one to calculate the output value. For the ... A set of samples, whose fixed perturbation sample set contains The perturbated 4D input feature vector will be the first... Substituting a fixed perturbation sample into the candidate expression tree yields the perturbation output value. The variance of all perturbation output values ​​for the same sample is used as the single-sample perturbation variance. ,in For the first The mean of all perturbation output values ​​in the sample group. The preset number is the fixed number of perturbation samples for each sample. The mean of the variances of the single-sample perturbations of all samples is taken as the component perturbation robustness index of the candidate expression tree, denoted as . The smaller the robustness index for component perturbations, the less sensitive the candidate expression tree is to small perturbations in the logarithmic ratio coordinate components, and the stronger the prediction stability under raw material fluctuation conditions. In the specific implementation described above, Equals 20, It equals 60.

[0045] It is worth noting that the introduction of the component perturbation robustness index means that the expression selection process not only focuses on fitting accuracy and formula simplicity, but also incorporates the unavoidable raw material fluctuation factors in engineering practice into the model selection criteria. This design differs substantially from pure symbolic regression methods: traditional symbolic regression only compromises between fitting error and complexity, and the resulting expression may fit well mathematically, but in practical applications, it will produce drastic prediction fluctuations due to input perturbations; while this method, by applying perturbations in the equidistant logarithmic ratio coordinate space and evaluating the output variance, directly embeds the local stability of the expression into the evolutionary selection pressure, so that the final selected expression can still provide reliable performance predictions when the component ratios shift slightly.

[0046] After calculating the three fitness metrics, a non-dominated ranking is performed on all candidate expression trees in the current candidate expression population, with the three objectives of minimizing fitting error, minimizing complexity, and minimizing component perturbation robustness. The specific process of non-dominated ranking is as follows: Any two candidate expression trees in the population are compared one by one. If the fitness metrics of the first candidate expression tree are all no worse than the corresponding metrics of the second candidate expression tree, and at least one metric is strictly better than the second candidate expression tree, then the first candidate expression tree is determined to dominate the second candidate expression tree. Candidate expression trees that are not dominated by any other candidate expression tree are assigned to the first non-dominated level and removed from the population. Then, non-dominated individuals are extracted from the remaining candidate expression trees and assigned to the second non-dominated level. This process is repeated until all candidate expression trees are assigned a non-dominated level. The smaller the non-dominated level value, the better the overall performance of the candidate expression tree across the three objectives.

[0047] After assigning a non-dominated level to each candidate expression tree, a tournament selection process is performed within the current candidate expression population to determine the parent expression trees to participate in subsequent crossover and mutation operations. Each time, a fixed number of candidate expression trees are randomly selected from the candidate expression population to form a tournament group. In one specific implementation, the fixed number of the tournament group is 7. In another optional implementation, the fixed number of the tournament group can be an integer between 3 and 10; a larger value indicates stronger selection pressure and a higher probability of selecting superior individuals, but also a faster decline in population diversity. The candidate expression tree with the optimal non-dominated level from the tournament group is selected as the parent expression tree. If multiple candidate expression trees with the same non-dominated level exist in the tournament group, the candidate expression tree with the smallest fitting error is selected as the parent expression tree. Using the non-dominated level as the primary selection criterion ensures that the direction of selection pressure during the evolutionary process remains consistent with the final selection criterion, achieving a balanced drive between fitting accuracy, formula simplicity, and component perturbation stability, rather than a single bias towards fitting accuracy.

[0048] After obtaining two parent expression trees, a subtree crossover operation is performed. In the first parent expression tree, a non-root node and its root subtree are randomly selected. In the second parent expression tree, a non-root node and its root subtree are also randomly selected. These two subtrees are then swapped and grafted into the other parent expression tree to generate two descendant expression trees. The purpose of the subtree crossover operation is to recombine local structural fragments in the two different expression trees, exploring new function forms. For example, if the first parent expression tree contains a substructure representing multiplicative combinations, and the second parent expression tree contains a substructure representing logarithmic transformations, the crossover may produce a new structure that nests the logarithmic transformation within a multiplicative combination—a structure that does not exist in either parent expression tree. Through the accumulation of numerous crossover operations, the expression structure in the population is continuously enriched and reorganized, increasing the chance of finding the globally optimal expression.

[0049] Point mutation operations are performed on the descendant expression tree with a preset mutation probability. A node is randomly selected from the descendant expression tree; if the selected node is an operator node, it is replaced with another operator from the candidate operator set; if the selected node is a terminal node, it is replaced with another terminal from the candidate terminal set. In one specific implementation, the preset mutation probability is 0.1, meaning each descendant expression tree has a 10% probability of point mutation. In another optional implementation, the preset mutation probability can be a value between 0.05 and 0.2. The purpose of point mutation operations is to introduce small, random changes locally into the expression tree, preventing the population from prematurely converging to a specific function structure and losing its exploration ability. Crossover operations focus on the reorganization of the global structure, while mutation operations focus on the perturbation of local details; the two work together to ensure a balance between exploration and utilization in the search process.

[0050] After merging all descendant expression trees to form the next generation of candidate expression trees, the processes of fitness evaluation, non-dominated sorting, tournament selection, subtree crossover, and point mutation are repeated until the preset maximum number of iterations is reached. In one specific implementation, the preset maximum number of iterations is 200. In another optional implementation, the preset maximum number of iterations can be an integer between 100 and 500. The choice of the number of iterations needs to balance search sufficiency and computation time. Too few generations may cause the search to terminate before the population has fully evolved, while too many generations will increase unnecessary computational overhead. In the specific implementation above, the population size is 500, the preset maximum number of iterations is 200, and approximately 100,000 candidate expression trees are evaluated in the entire genetic programming process.

[0051] After iteration, non-dominated ranking is performed on all candidate expression trees generated throughout the iteration process, with the three selection objectives of minimizing fitting error, minimizing complexity index, and minimizing component perturbation robustness index. The scope of non-dominated ranking here is not limited to the last generation but covers all candidate expression trees that have appeared in all iterations, ensuring that excellent expressions emerging in early iterations are not lost due to random drift in subsequent evolution. A set of non-dominated expression trees is retained where none of the three indices are simultaneously surpassed by any other candidate expression tree. From this set, candidate expression trees with the smallest fitting error and a complexity index not exceeding a preset node limit are selected, and the analytical expression corresponding to the selected candidate expression tree is determined as the explicit prediction formula for mobility. The logic of this selection strategy is: first, expressions that are too complex or have poor robustness are eliminated through non-dominated ranking and complexity limit constraints; then, the one with the highest fitting accuracy is selected from the remaining candidates, thus achieving the best trade-off between interpretability, stability, and accuracy.

[0052] refer to Figure 3 , Figure 3 The diagram shows a three-criterion non-dominated sorting and expression filtering. Figure 3 The diagram illustrates the non-dominated ranking of the candidate expression population in the three-dimensional target space composed of fitting error, complexity index, and component perturbation robustness index. The candidate expression tree of the first non-dominated level constitutes the Pareto front, reflecting the trade-off between fitting accuracy, expression complexity, and perturbation robustness. The position in the Pareto front where the complexity index does not exceed the preset upper limit of the number of nodes and the fitting error is the candidate expression tree corresponding to the finally selected explicit prediction formula.

[0053] In the specific implementation described above, the explicit prediction formula for mobility determined after 200 generations of genetic programming iterations from the non-dominated set in the three-dimensional target space is: yf = 186.42 + 14.27 × w − 8.45 × z1 + 6.83 × ln(z2 + 3.5) − 4.27 × z3, where yf is the predicted mobility value in millimeters, z1, z2, and z3 are the 1st, 2nd, and 3rd components of the equidistant logarithmic ratio coordinate vector, w is the water-to-solid ratio, and the five constants 186.42, 14.27, −8.45, 6.83, −4.27, and the constant offset 3.5 within the logarithmic operation are all real constants automatically determined during the genetic programming process. The complexity metric of this formula is the total number of nodes, which is 13, within the preset upper limit of 30. The fitting error, the sum of the absolute values ​​of the differences between the predicted and measured mobility values ​​of all training samples, is 52.7 mm, corresponding to an average prediction bias of approximately 0.88 mm per sample. The robustness metric for component perturbation is the mean variance of the output after applying a perturbation of 0.05 standard deviation to all training samples in the log-ratio coordinate space, which is 0.43 square millimeters, corresponding to an output standard deviation of approximately 0.66 mm. It should be noted that the specific constant values ​​and function structure of this formula are related to the training dataset and the random seed used in genetic programming. Re-running the genetic programming process under different training datasets or different random seeds may yield explicit mobility prediction formulas with different structures and constant values, but all will satisfy the upper limit constraint on complexity and the non-dominated sorting selection condition.

[0054] After determining the explicit prediction formula for mobility, the output variable is replaced with the log-permeability coefficient value, and the explicit prediction formula for the log-permeability coefficient is obtained using the same genetic programming process. Specifically, using the 4-dimensional input feature vector of each row in the input feature matrix as input, and the log-permeability coefficient value corresponding to each set of matching experiment records as output, a complete genetic programming process is run independently once using the same set of candidate operators, candidate terminal set, population size, initial maximum depth, tournament group size, preset mutation probability, and preset maximum number of iterations. In fitness evaluation, the method for calculating the fitting error becomes... ,in For the first The output value of the expression obtained by substituting the group of samples into the current candidate expression tree. For the first The logarithmic penetration coefficient values ​​corresponding to the sample groups. The calculation methods for complexity indices and component perturbation robustness indices are completely consistent with those in the mobility modeling stage, still using the same perturbation sample template. After non-dominated sorting and screening, the explicit prediction formula for the logarithmic penetration coefficient is determined.

[0055] In the specific implementation described above, the explicit prediction formula for the logarithmic permeability coefficient, determined after 200 generations of genetic programming iterations from the non-dominated set in the three-dimensional target space, is: ys = -6.82 - 1.05 × w + 0.47 × z1 - 0.28 × z2 × z3, where ys is the predicted value of the logarithmic permeability coefficient (dimensionless, corresponding to the permeability coefficient in centimeters per second, taking the commonly used logarithm), z1, z2, and z3 are the 1st, 2nd, and 3rd components of the equidistant logarithmic ratio coordinate vector, w is the water-to-solid ratio, and the four constants in the formula -6.82, -1.05, 0.47, and -0.28 are all real constants automatically determined during the genetic programming process. The complexity metric of this formula is the total number of nodes, which is 11, within the preset upper limit of 30. The fitting error, the sum of the absolute values ​​of the differences between the predicted values ​​and the measured log-permeability coefficient values ​​of all training samples, is 3.84, corresponding to an average prediction bias of approximately 0.064 per sample. The component perturbation robustness metric is the mean variance of the output after applying a perturbation of 0.05 to all training samples in the equidistant log-ratio coordinate space, which is 0.0018, corresponding to an output standard deviation of approximately 0.042. It should also be noted that the specific constant values ​​and function structure of this formula are related to the training dataset and the random seed used in genetic programming.

[0056] The explicit prediction formulas for flowability and logarithmic permeability are combined to form a dual-response analytical model for packing performance. This model receives a 4-dimensional input feature vector and outputs a predicted flowability value and a predicted logarithmic permeability value. The two explicit prediction formulas are established independently, each with different expression structures and constant parameters, but they share the same input variable space and the same equidistant logarithmic ratio transformation framework. In step 3, they will be used together to achieve a synergistic optimization design of flowability and impermeability.

[0057] The technical details of step 3 are explained below.

[0058] After step 2 is completed, a dual-response analytical model for packing performance, consisting of explicit prediction formulas for fluidity and logarithmic permeability coefficient, has been obtained. The goal of step 3 is to use this dual-response analytical model to systematically search for the optimal mix design within the feasible range of mix proportions, and output a set of recommended mix proportions that achieve the best trade-off between fluidity and impermeability, for engineering technicians to make the final selection based on specific construction conditions.

[0059] First, a set of candidate mix proportions is jointly generated within the engineering allowable range of the normalized solid composition space and the water-to-solid ratio. The reason for sampling within the normalized solid composition space, rather than within the original absolute mass space of the four solid composition variables, is that a closure operation has already been performed on the four solid composition variables in step 1. After the closure operation, what truly matters is the proportional relationship of each component in the total solids, not the absolute mass value. If samples were directly and independently combined within the original mass space of the four solid composition variables, a large number of candidate mix proportions that are equivalent or approximately equivalent after the closure operation would be generated. For example, the two mix proportions of 40 parts by mass of loess, 20 parts by mass of coal gangue, 10 parts by mass of carbide slag, and 10 parts by mass of solidifying agent are different from those of 80 parts by mass of loess, 40 parts by mass of coal gangue, 20 parts by mass of carbide slag, and 20 parts by mass of solidifying agent. Although the absolute quantities are different, the normalized proportions after the closure operation are exactly the same, and the corresponding logarithmic ratio coordinate vectors are also exactly the same. Therefore, substituting them into the dual-response analytical model of filler performance will yield completely consistent prediction results. Direct sampling within the normalized solid composition space can fundamentally avoid this redundancy, ensuring that each candidate composition point corresponds to a unique composition ratio, significantly improving sampling efficiency and computational resource utilization.

[0060] The process of generating the candidate mix proportion set is as follows. A preset total solid mass benchmark is set. In one specific implementation, the preset total solid mass benchmark is set to 100, meaning that the mass fractions of each solid component variable in the final output are based on a total solid mass of 100 parts. In another optional implementation, the preset total solid mass benchmark can be set to 1 or 1000, with the specific value determined according to engineering practice, without affecting the relative proportions and performance prediction results.

[0061] For the first three of the four solid component variables, normalized upper and lower bounds are set respectively. The four solid component variables, in the fixed order determined in step 1, are loess mass fractions, coal gangue mass fractions, carbide slag mass fractions, and solidifying agent mass fractions, with the first three being loess mass fractions, coal gangue mass fractions, and carbide slag mass fractions. The normalized upper and lower bounds reflect the maximum and minimum allowable mass proportions of each component in the total solids, determined by engineering experience and construction conditions. In a specific implementation, the normalized lower bound for loess mass fractions is 0.30, and the normalized upper bound is 0.70; the normalized lower bound for coal gangue mass fractions is 0.10, and the normalized upper bound is 0.40; and the normalized lower bound for carbide slag mass fractions is 0.05, and the normalized upper bound is 0.25. In another optional implementation, the above-mentioned upper and lower bounds of the normalized values ​​can be adjusted according to the loess characteristics and raw material supply of different engineering sites.

[0062] Multiple sampling nodes are inserted at equal intervals between the upper and lower bounds of the normalized value for each item, and then fully combined. In one specific implementation, 20 equally spaced sampling nodes are inserted for each of the first three solid component variables. Taking loess mass fraction as an example, the lower bound of the normalized value is 0.30, the upper bound of the normalized value is 0.70, and the equal interval step size is... The 20 sampling nodes are successively 0.30, 0.321, 0.342, and so on up to 0.70. Sampling nodes for the mass fractions of coal gangue and carbide slag are generated in the same way according to their respective normalized upper and lower bounds. The sampling nodes for all three items are then combined to generate a total of... In another alternative implementation, the number of sampling nodes can be an integer between 10 and 50. The more nodes there are, the higher the coverage density of candidate matching points and the better the search accuracy, but the computational cost also increases accordingly.

[0063] The normalized value of the fourth solid component variable in each combination is determined by subtracting the sum of the normalized values ​​of the first three variables from 1. The fourth solid component variable is the mass fraction of the curing agent, and its normalized value is determined in a way that reflects the constraint of closed operations, i.e., the sum of the four components in the normalized solid component vector must equal 1. Let the normalized value of the loess mass fraction in a certain combination be... The normalized value of the mass fraction of coal gangue is The normalized value of the mass fraction of carbide slag is The normalized value of the mass fraction of the curing agent is... The reason for using this method, instead of sampling all four component variables independently, is that there is a constraint that the sum of the four normalized values ​​equals 1. In reality, only three of the four variables are independent; the fourth is entirely determined by the first three. If each of the four variables were sampled independently and then normalized, it would not only generate a large number of equivalent points but also lead to an uneven distribution of samples on the simplex, with some regions being oversampled while others are missed. By fixing the first three variables to be sampled independently and the fourth variable to be determined by the constraint, we can ensure that the uniform sampling in the 3D independent space maintains a relatively uniform coverage on the 4D simplex.

[0064] The normalized candidate set of solid components is obtained by eliminating combinations whose normalized values ​​for the fourth solid component variable are negative or exceed the allowable range in engineering. When the sum of the first three normalized values ​​exceeds 1, the normalized value of the fourth variable is... The result will be negative, which is physically meaningless because the mass ratio cannot be negative, and such combinations must be eliminated. Furthermore, even... For positive values, it is also necessary to verify whether they fall within the engineering allowable range for the mass fraction of the curing agent. In the specific embodiment described above, the engineering allowable range for the mass fraction of the curing agent is 0.03 to 0.15, meaning the normalized value should be between 0.03 and 0.15. Combinations outside this range are also eliminated. After elimination, in the specific embodiment described above, approximately 3200 valid combinations remain from the initial 8000 combinations to form the normalized solid component candidate set. The specific number depends on the size of the intersection between the sampling node distribution of the first three items and the engineering allowable range of the fourth item.

[0065] Multiple water-to-solid ratio sampling nodes are inserted at equal intervals between the upper and lower limits of the engineering allowable range for the water-to-solid ratio. In the specific implementation described above, the lower limit of the engineering allowable range for the water-to-solid ratio is 0.4, the upper limit is 0.8, and 10 equally spaced water-to-solid ratio sampling nodes are inserted with a step size of [missing value]. The 10 water-to-solid ratio sampling nodes are 0.4, 0.444, 0.489, and so on up to 0.8. In another optional embodiment, the number of water-to-solid ratio sampling nodes can be an integer between 5 and 20.

[0066] Each normalized combination in the normalized solids composition candidate set is fully combined with each water-to-solids ratio sampling node. In the specific implementation described above, 3200 normalized combinations are fully combined with 10 water-to-solids ratio sampling nodes to generate... There are candidate matching points. For each normalized combination, the same isometric logarithmic transformation as in step 1 is performed sequentially to obtain a 3D isometric logarithmic coordinate vector. Specifically, the natural logarithm of the four normalized values ​​in the normalized combination is taken to obtain a logarithmic vector, and the orthogonal basis matrix constructed in step 1 is... A matrix multiplication is performed with the logarithmic vector to obtain a 3D isochronous logarithmic ratio coordinate vector. This 3D isochronous logarithmic ratio coordinate vector is then concatenated with the corresponding water-to-solid ratio sample value to form a 4D input feature vector, generating a candidate ratio point set. In the specific implementation described above, the candidate ratio point set contains 32,000 candidate ratio points, each represented by a 4D input feature vector and its corresponding normalized solid component combination and water-to-solid ratio sample value.

[0067] The 4-dimensional input feature vector of each candidate mix design point in the candidate mix design set is substituted into the explicit prediction formulas for mobility and log-permeability coefficients, respectively, to calculate the predicted values ​​for mobility and log-permeability coefficients. Since both the explicit prediction formulas for mobility and log-permeability coefficients are analytical expressions, the substitution calculations only involve basic arithmetic and logarithmic operations, resulting in extremely fast computation speed. In the specific implementation described above, the dual-response prediction of 32,000 candidate mix design points can be completed within seconds on a typical desktop computer. This is another practical advantage of using explicit prediction formulas instead of a black-box model.

[0068] Subsequently, candidate mix proportions with predicted flowability values ​​below the preset lower limit or above the preset upper limit are eliminated. This elimination process reflects the engineering essence of flowability as a range constraint rather than a simple maximization objective. During construction, loess water-blocking pile filler needs to be injected into the pile holes via pumping or grouting. If the flowability is too low, the filler cannot flow smoothly into place, resulting in voids or interlayers in the pile hole; if the flowability is too high, solid particles in the filler are prone to sedimentation and separation, with coarse particles settling at the bottom while fine particles and water rise to the top, leading to severely uneven impermeability performance between the upper and lower parts of the hardened pile. Therefore, only mix proportions with flowability falling within the allowable construction window have practical engineering value. The preset lower limit and preset upper limit of flowability constitute the preset allowable flowability range. In one specific implementation, the preset lower limit of flowability is set at 160 mm, and the preset upper limit of flowability is set at 220 mm. In another optional implementation, the preset lower limit and preset upper limit of fluidity can be adjusted according to the requirements of specific construction processes. For example, when using high-pressure grouting, the preset lower limit of fluidity can be appropriately lowered, while when using gravity-flow grouting, the preset lower limit of fluidity needs to be raised. In the above specific implementation, after eliminating fluidity intervals, approximately 18,000 candidate mix proportions that meet the fluidity constraints are retained out of 32,000 candidate mix proportions.

[0069] For the retained candidate mix proportions, a non-dominated ranking was performed with two optimization objectives: maximizing the predicted flowability and minimizing the predicted logarithmic permeability. The reason for maximizing the predicted flowability, rather than using a fixed value, is that, under the constraint of the construction flowability window, higher flowability means better pumpability of the filler, more thorough grouting, and higher pile density. Conversely, a lower predicted logarithmic permeability means a lower permeability coefficient after hardening (because the logarithm is a monotonically increasing function, a smaller logarithmic permeability coefficient corresponds to a smaller original permeability coefficient), and stronger water-blocking performance. These two objectives are often contradictory: increasing flowability often requires increasing the water-to-solid ratio or reducing the amount of hardener, both of which can lead to increased porosity and permeability after hardening. Therefore, no single mix proportion can simultaneously achieve both objectives at their optimal values; rather, there exists a set of compromise mix proportions where improving the performance of one objective inevitably comes at the expense of the other. This set of compromise mix proportions is the Pareto optimal mix proportion set.

[0070] The specific process of non-dominated ordination is as follows: Compare any two candidate mixes one by one. If the predicted mobility of the first candidate mix is ​​greater than or equal to the predicted mobility of the second candidate mix, and the predicted log-permeability coefficient of the first candidate mix is ​​less than or equal to the predicted log-permeability coefficient of the second candidate mix, and at least one objective is strictly better than the second candidate mix, then the first candidate mix dominates the second candidate mix. "At least one objective is strictly better" means that the predicted mobility of the first candidate mix is ​​strictly greater than the predicted mobility of the second candidate mix, or the predicted log-permeability coefficient of the first candidate mix is ​​strictly less than the predicted log-permeability coefficient of the second candidate mix, or both. This definition ensures the strictness of the dominance relationship and avoids classifying two candidate mixes with identical performance as having a dominance relationship. Candidate mixes that are not dominated by any other candidate mix are extracted as the Pareto optimal mix set. For any candidate mix designation point in the Pareto optimal mix designation set, there is no other candidate mix designation point that is simultaneously non-inferior to it in both the predicted mobility and the predicted log permeability coefficient, and at least one candidate mix designation point is strictly superior to it in at least one objective. In the specific implementation described above, the Pareto optimal mix designation set contains approximately 50 to 200 candidate mix designation points, the specific number depending on the sampling density of the candidate mix designation set and the nonlinear characteristics of the dual-response analytical model for packing performance.

[0071] The candidate proportions in the Pareto optimal proportion set are currently stored as normalized solid component combinations and water-to-solid ratios, and need to be converted into a proportion expression directly usable by engineers through an inverse logarithmic ratio transformation. An inverse logarithmic ratio transformation is performed on each candidate proportion in the Pareto optimal proportion set. This inverse logarithmic ratio transformation is the reverse of the logarithmic ratio transformation in step 1, and its purpose is to convert the 3D logarithmic ratio coordinate vector back to a 4D normalized solid component vector. Specifically, the transpose of the orthogonal basis matrix is ​​used... The coordinate vector of the 3D isochronous logarithmic ratio of the candidate matching point Performing matrix multiplication yields a 4-dimensional logarithmic vector. ,in Orthogonal basis matrix The transpose of the matrix has 4 rows and 3 columns. It is a 3D logarithmic ratio coordinate vector. It is a 4-dimensional logarithmic vector. Orthogonal basis matrix. satisfy The orthogonality condition for a 3x3 identity matrix means that transpose matrix multiplication can accurately inversely map the isochronous logarithmic ratio coordinates back to logarithmic space. Taking the natural exponent for each element in the 4-dimensional logarithmic vector, i.e., for each element... The Each component calculate ,in This indicates the operation of the natural exponent. Taking integers from 1 to 4, we get four positive values. The sum of these four positive values ​​is generally not equal to 1, so we need to perform closed normalization: divide all elements by the sum of all elements so that the sum of the four normalized values ​​equals 1. Let the four values ​​after taking the natural exponent be... , , , The normalized value after closed normalization is ,in For the first The normalized value of the solid component variable. Take an integer from 1 to 4. Closed normalization ensures that the inverse transform result still satisfies the component constraint that the sum equals 1, which is completely consistent with the output format of the closed operation in step 1.

[0072] It is worth noting that in the process of generating the candidate ratio point set in step 3, the candidate ratio points themselves are directly sampled and generated in the normalized solid composition space, and their normalized combination is known. Therefore, the role of the inverse logarithmic ratio transformation here is mainly to verify consistency and provide a unified processing flow for subsequent engineering outputs. For candidate ratio points sampled directly from the normalized space, the result of the inverse transformation is consistent with the original normalized combination within the range of numerical precision. However, in more general application scenarios, such as the optimal point obtained by continuous optimization algorithms in the logarithmic ratio coordinate space, the inverse transformation is a necessary step to restore the normalized ratio. Therefore, the inverse logarithmic ratio transformation operation is uniformly retained in the method flow to ensure the integrity and universality of the method.

[0073] The mass fraction is obtained by multiplying each normalized value by a preset total solid mass benchmark. Taking a solid component variable as an example, its mass fraction under the preset total solid mass standard is: ,in After the inverse transformation, the first The normalized value of the solid component variable. To establish a baseline for the total mass of solids, For the first The mass fraction of each solid component variable under a preset total solid mass standard. Integers from 1 to 4 correspond to the mass fractions of loess, coal gangue, carbide slag, and solidifying agent, respectively. In the specific implementation method described above, a preset total solid mass benchmark is used. If the normalized value of a Pareto optimal candidate mix is ​​0.50 for loess, 0.25 for coal gangue, 0.15 for carbide slag, and 0.10 for solidifying agent, then the converted mass fractions of loess are 50, coal gangue is 25, carbide slag is 15, and solidifying agent is 10, all in mass units. Together with the corresponding water-to-solid ratio, the output is a set of optimal mix proportions for loess water-blocking pile fillers, consisting of all mass fractions and the water-to-solid ratio.

[0074] The optimal mix design scheme for loess water-blocking pile filler comprises multiple schemes. Each scheme consists of five parameters: the mass fractions of loess, coal gangue, carbide slag, and solidifying agent, and the water-to-solid ratio, along with corresponding predicted values ​​for fluidity and logarithmic permeability. These schemes are distributed along the Pareto front, exhibiting a systematic trade-off between the predicted fluidity and logarithmic permeability: moving along the Pareto front from one end to the other, the predicted fluidity gradually increases while the predicted logarithmic permeability also gradually increases, or vice versa. Engineers can select the appropriate scheme based on the specific site requirements. For example, for sites with high groundwater levels and strict water-blocking requirements, the scheme with the lowest predicted logarithmic permeability along the Pareto front is preferred; for sites with deep pile holes and long pumping distances, the scheme with the highest predicted fluidity along the Pareto front is preferred. The Pareto optimal set does not provide a single optimal solution, but a series of alternatives that do not dominate each other in two performance dimensions. This multi-option output provides ample flexibility for engineering decision-making.

[0075] In another alternative implementation, decision preferences can be further introduced based on the Pareto optimal mix design set. For example, engineers specify the expected target value for the predicted flowability and select the mix design with the predicted flowability value closest to the expected target value from the Pareto optimal mix design set as the final recommended solution. Alternatively, a comprehensive score can be calculated for each mix design in the Pareto optimal mix design set. The comprehensive score is a weighted combination of the normalized value of the predicted flowability and the normalized value of the predicted logarithmic permeability coefficient. The weights of flowability and impermeability in the weighted combination are set by engineers based on the construction site conditions, and the mix design with the highest comprehensive score is selected as the final recommended solution. These methods of introducing decision preferences do not change the generation process of the Pareto optimal mix design set itself; they only add a selection mechanism tailored to specific engineering scenarios on top of the Pareto frontier.

[0076] The present invention has been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of the invention. The descriptions of the embodiments above are merely for the purpose of helping to understand the method and core ideas of the present invention. It should be noted that those skilled in the art can make various improvements and modifications to the present invention without departing from its principles, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A method for intelligent design of loess water-blocking pile filler and prediction of impermeability, characterized in that, Includes the following steps: Step 1: Collect multiple sets of mix design test records for loess water-blocking pile filler. Each set of mix design test records includes four solid component variables: loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction, as well as the water-solid ratio and the corresponding fluidity and permeability coefficient values. Take the common logarithm of the permeability coefficient value to obtain the logarithmic permeability coefficient value. Perform zero-value substitution preprocessing, closure operation, and equidistant logarithmic ratio transformation on the four solid component variables in sequence to obtain the equidistant logarithmic ratio coordinate vector. Concatenate the equidistant logarithmic ratio coordinate vector with the water-solid ratio to form the input feature vector and stack them row by row to form the input feature matrix. Generate a perturbation sample template based on the equidistant logarithmic ratio coordinate components. Step 2: Using the input feature vector as input and the flowability value and log permeability coefficient value as output, the candidate expression tree is evaluated simultaneously by genetic programming-driven symbolic regression, which evaluates the fitting error, complexity index and component perturbation robustness index. After non-dominated sorting-assisted iterative evolution and screening from the non-dominated set, the explicit prediction formulas for flowability and log permeability coefficient are determined, thus forming a dual-response analytical model for packing performance. Step 3: Within the normalized solid composition space and the engineering allowable range of water-solid ratio, a set of candidate mix proportions is jointly generated. After equidistant logarithmic ratio transformation, the mix proportions are input into the dual-response analytical model of the filler performance to obtain the predicted values ​​of fluidity and logarithmic permeability coefficient. Candidate mix proportions with predicted fluidity values ​​exceeding the preset allowable fluidity range are eliminated. For the remaining candidate mix proportions, a non-dominated sorting is performed with the dual objectives of maximizing the predicted fluidity value and minimizing the predicted logarithmic permeability value to extract the Pareto optimal mix proportion set. After equidistant logarithmic ratio inverse transformation, the optimal mix proportion scheme set of loess water-blocking pile filler is output under the preset total solid mass benchmark.

2. The method of claim 1, wherein, Each set of mixing ratio test records was obtained under the same curing age, the same specimen preparation method, and the same penetration test conditions.

3. The method of claim 1, wherein, In step 1, the zero-value substitution preprocessing process is as follows: when the value of any solid component variable is less than the preset positive lower limit, it is substituted with the preset positive lower limit; the closed operation process is as follows: the four solid component variables after the zero-value substitution preprocessing are arranged into a 4-dimensional component vector in a fixed order of loess mass fraction, coal gangue mass fraction, carbide slag mass fraction, and solidifying agent mass fraction; the four solid component variables are divided by the sum of the four solid component variables respectively, so that the sum of the four solid component variables after the closed operation is equal to 1, and the normalized solid component vector is obtained.

4. The method of claim 3, wherein, In step 1, the process of the equidistant logarithmic ratio transformation is as follows: Construct an orthogonal basis matrix with 3 rows and 4 columns; the first row of the orthogonal basis matrix contains 4 elements, where the first element is positive, the second to fourth elements are equal and all negative, the sum of the squares of all elements in the first row is equal to 1 and the sum of all elements is equal to 0; the first element of the second row is 0, the second element is positive, the third and fourth elements are equal and all negative, the sum of the squares of all elements in the second row is equal to 1 and the sum of all elements is equal to 0; the first and second elements of the third row are both 0, the third element is positive, the fourth element is negative, the sum of the squares of all elements in the third row is equal to 1 and the sum of all elements is equal to 0; take the natural logarithm of each element in the normalized solid component vector to obtain a logarithmic vector, and perform matrix multiplication between the orthogonal basis matrix and the logarithmic vector to obtain a 3D equidistant logarithmic ratio coordinate vector.

5. The method of claim 1, wherein, In step 1, the process of generating the perturbation sample template is as follows: For the logarithmic ratio coordinate components of each sample in the input feature matrix, a preset number of perturbation vectors are generated using a preset fixed random seed. Each component of each perturbation vector independently follows a normal distribution with a mean of 0 and a standard deviation of a preset perturbation amplitude. Each perturbation vector is superimposed on the logarithmic ratio coordinate components of the corresponding sample while keeping the water-solid ratio component unchanged, thus forming a fixed perturbation sample set for the corresponding sample. The fixed perturbation sample sets of all samples are then summarized into a perturbation sample template.

6. The method of claim 1, wherein, In step 2, the candidate operator set includes five operators: addition, subtraction, multiplication, protected division, and protected logarithm operation. Protected division sets the output value to 1 when the absolute value of the divisor is less than a preset minimum threshold. Protected logarithm operation replaces the input value with the preset positive threshold before performing the natural logarithm operation when the input value is less than or equal to a preset positive threshold. The candidate terminal set includes equidistant logarithmic ratio coordinate components, water-solid ratio variables, and random constants.

7. The method of claim 1, wherein, In step 2, the fitting error is the sum of the absolute values ​​of the differences between the expression output value of each sample and the corresponding output variable after substituting all samples in the input feature matrix into the candidate expression tree; the complexity index is the total number of all nodes in the candidate expression tree; the calculation process of the component perturbation robustness index is as follows: take out the fixed perturbation sample set of each sample from the perturbation sample template, substituting each fixed perturbation sample in the fixed perturbation sample set into the candidate expression tree to calculate the output value, and statistically analyze the variance of all perturbation output values ​​of the same sample as the single sample perturbation variance, and take the mean of the single sample perturbation variance of all samples.

8. The method of claim 1, wherein, In step 2, the process of non-dominated sorting-assisted iterative evolution is as follows: non-dominated sorting is performed on the current candidate expression population with three objectives: minimizing fitting error, minimizing complexity index, and minimizing component perturbation robustness index, and a non-dominated level is assigned to each candidate expression tree; each time, a fixed number of candidate expression trees are randomly selected from the candidate expression population to form a tournament group, and the candidate expression tree with the best non-dominated level is selected from the tournament group as the parent expression tree. If there are multiple candidate expression trees with the same non-dominated level in the tournament group, the candidate expression tree with the smallest fitting error is selected as the parent expression tree.

9. The method according to claim 8, characterized in that, After obtaining two parent expression trees, perform a subtree crossover operation. In the first parent expression tree, randomly select one non-root node and a subtree rooted at the selected non-root node. In the second parent expression tree, randomly select one non-root node and a subtree rooted at the selected non-root node. Then, swap the two subtrees and graft them onto the other parent expression tree to generate two offspring expression trees. Perform point mutation operations on the descendant expression tree with a preset mutation probability. Randomly select one node in the descendant expression tree. If the selected node is an operator node, replace it with another operator in the candidate operator set. If the selected node is a terminal node, replace it with another terminal in the candidate terminal set. After the iteration, select the candidate expression tree with the smallest fitting error and a complexity index that does not exceed the preset node number limit from the non-dominated set. Determine the analytical expression corresponding to the selected candidate expression tree as the explicit prediction formula.

10. The method according to claim 1, characterized in that, In step 3, the process of generating the candidate mix proportion set is as follows: For the first three of the four solid component variables, set upper and lower bounds for normalized values ​​respectively. Insert multiple sampling nodes at equal intervals between the upper and lower bounds of normalized values ​​for each variable and perform a full combination. In each combination, the normalized value of the fourth solid component variable is determined by subtracting the sum of the normalized values ​​of the first three variables from 1. Eliminate combinations where the normalized value of the fourth solid component variable is negative or exceeds the allowable range for engineering purposes to obtain a normalized solid component candidate set. Combine the normalized solid component candidate set with the water-to-solid ratio sampling nodes to generate the candidate mix proportion set. In the non-dominated sorting, compare any two candidate mix proportion points one by one. If the predicted flowability value of the first candidate mix proportion point is greater than the expected value, the candidate mix proportion is considered a candidate mix proportion. If the predicted flowability value of the second candidate mix point and the predicted log permeability coefficient value of the first candidate mix point are less than or equal to the predicted log permeability coefficient value of the second candidate mix point, and at least one target is strictly better than the second candidate mix point, then the first candidate mix point dominates the second candidate mix point. The candidate mix points that are not dominated by any other candidate mix points are extracted as the Pareto optimal mix set. The process of the inverse transformation of the equidistant logarithmic ratio is as follows: perform matrix multiplication between the transpose of the orthogonal basis matrix and the equidistant logarithmic ratio coordinate vector to obtain a 4-dimensional logarithmic vector. Take the natural exponent for each element in the 4-dimensional logarithmic vector and divide all elements by the sum of all elements to complete the closed normalization. Multiply each normalized value by the preset total solid mass benchmark to obtain the mass fraction.