Artificial intelligence (ai) assisted three-dimensional design method and system for foam soil filled subgrade

Abstract

This invention relates to the field of roadbed engineering technology, and in particular to an AI-assisted three-dimensional design method and system for foamed soil filler roadbeds. By acquiring roadbed geometric parameters, traffic load spectra, and foundation soil layer distribution, the method generates a foundation bearing capacity map and calculates the stress distribution on the roadbed surface to divide treatment sections. For the foamed soil treatment sections, a three-dimensional voxel mesh is constructed and assigned density values ​​to generate a three-dimensional density field. The optimal density field configuration is obtained through numerical solution and Bayesian optimization iteration, and finally transformed into construction drawings. This invention achieves intelligent optimization design of the three-dimensional distribution of foamed soil filler, improving the accuracy and efficiency of roadbed design.

Description

Technical Field

[0001] This invention relates to the field of roadbed engineering technology, and in particular to an AI-assisted three-dimensional design method and system for foamed soil filler roadbeds. Background Technology

[0002] In the design of foamed soil subgrades, current conventional practices typically rely on geological survey reports and engineering experience, using simplified mechanical models to determine the thickness and density of the foamed soil fill. Designers, based on the distribution of soft soil layers, use recommended bearing capacity formulas or semi-empirical algorithms to calculate the settlement and stability of the foundation under embankment loads, thereby identifying sections requiring replacement or treatment. For these sections, the foamed soil filler is often configured with a single density value or a fixed density for different blocks, and is filled in layers during construction, with the thickness and compaction degree of each layer set based on empirical data. The entire design process is based on static parameters, lacking dynamic analysis of the relationship between load distribution and foundation bearing capacity in the three-dimensional space of the subgrade.

[0003] Conventional methods have significant drawbacks. On the one hand, the stress distribution on the subgrade surface is influenced by traffic loads, subgrade geometry, and soil properties. However, conventional designs simplify calculations using only average or maximum loads, neglecting the non-uniformity of stress along the transverse and longitudinal directions. This leads to inaccurate division of foamed soil treatment zones, potentially causing over-excavation or under-excavation in certain areas, resulting in material waste or insufficient bearing capacity. On the other hand, the density distribution of foamed soil filler plays a decisive role in subgrade settlement and stability. However, existing methods often employ homogeneous assumptions or empirical zoning assignments, failing to achieve adaptive matching between density values ​​and local stresses and foundation bearing capacity. This makes it difficult to achieve an optimal balance between controlling total settlement and ensuring slope stability, often requiring subsequent field tests or repeated adjustments to meet engineering requirements. Summary of the Invention

[0004] This invention provides an AI-assisted three-dimensional design method and system for foamed soil filler roadbeds, which can solve the problems in the prior art.

[0005] A first aspect of the present invention provides an AI-assisted three-dimensional design method for foamed soil filler roadbeds, comprising:

[0006] Obtain the geometric parameters of the roadbed, traffic load spectrum, and distribution of the foundation soil layers for the roadbed engineering;

[0007] The features of the foundation soil layer distribution are extracted to generate a foundation bearing capacity map. The roadbed surface stress distribution is calculated based on the traffic load spectrum. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, the foamed soil filling treatment section and the conventional filling section are divided.

[0008] For the foamed soil filler treatment section, a three-dimensional voxel grid is constructed in combination with the roadbed geometric parameters. Based on the foundation bearing capacity map, density values ​​are assigned to each voxel node, and a three-dimensional density field distribution of the foamed soil filler is generated through a spatial interpolation algorithm.

[0009] Boundary constraints and loads are applied to the three-dimensional density field distribution, and the settlement value and stability coefficient are obtained by numerical solution. The deviation is calculated by comparing with the target design value. The three-dimensional density field distribution is iteratively adjusted by Bayesian optimization algorithm until the convergence condition is met, and the optimal density field configuration is obtained.

[0010] A three-dimensional spatial distribution model of foamed soil filler is generated based on the optimal density field configuration, and the optimal density field configuration is converted into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and the construction drawings are output.

[0011] Feature extraction is performed on the distribution of the foundation soil layers to generate a foundation bearing capacity map. The roadbed surface stress distribution is calculated based on the traffic load spectrum. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, foamed soil filling treatment sections and conventional filling sections are divided, including:

[0012] The distribution of the foundation soil layers is processed by layering and slicing, and the rock and soil type identification and thickness data of each soil layer are extracted to construct a vertical profile model containing soil layer sequence and burial depth information;

[0013] Bearing characteristic parameters are assigned to each soil layer in the vertical profile model, the ultimate bearing capacity and allowable settlement of each soil layer are calculated, and the distribution of total bearing capacity and total settlement potential of the foundation are generated through vertical cumulative calculation, and integrated into a foundation bearing performance map.

[0014] The dominant frequency load component is extracted by performing spectral analysis on the traffic load spectrum. The dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at the interfaces of each soil layer is calculated using the elastic layered system theory to obtain the stress distribution on the roadbed surface.

[0015] The difference between the total bearing capacity distribution in the foundation bearing capacity performance map and the stress distribution on the roadbed surface is calculated to generate a bearing surplus distribution.

[0016] Identify the areas with negative values ​​in the load-bearing surplus distribution, classify the areas with insufficient load-bearing capacity according to the amount of load loss, and determine the treatment intensity corresponding to each level. Areas with treatment intensity greater than the treatment threshold are divided into foamed soil filling treatment sections, and the remaining areas are divided into conventional filling sections.

[0017] The dominant frequency load component is extracted by spectral analysis of the traffic load spectrum. This dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at each soil layer interface is calculated using elastic layered system theory to obtain the stress distribution on the roadbed surface, including:

[0018] The traffic load spectrum is segmented according to a preset time window and subjected to spectrum transformation to obtain a local spectrum. The frequency components that recur in all local spectra are counted, and the frequency component with the most recurrences is determined as the stable main frequency. The load amplitude corresponding to the stable main frequency is extracted as the main frequency load component.

[0019] The dominant frequency load component is divided into distributed load nodes according to the transverse width and longitudinal length of the roadbed top surface, and the load intensity is determined according to the distance of each load node from the load action center.

[0020] Extract the compression modulus and lateral compression coefficient of each soil layer in the vertical profile model, calculate the strain response coefficient of each soil layer under vertical load, and determine the damping characteristics of each soil layer.

[0021] Select the load nodes on the top surface of the subgrade and use the corresponding load intensity as the initial stress input to the first soil layer. Calculate the stress attenuation within the layer and its transmission to the interface based on the damping characteristics of the first soil layer. Determine the degree of interface abrupt change by the difference in soil properties on both sides of the interface.

[0022] When the degree of interface mutation exceeds the interface sensitivity threshold, an interface reflection loss coefficient is introduced. Combined with the stress at the interface, the stress continues to propagate to the next soil layer. The stress attenuation and interface treatment are repeated until the top surface of the foundation is reached. The stresses propagated from all load nodes to the top surface of the foundation are spatially superimposed to obtain the stress distribution of the roadbed surface.

[0023] For the foamed soil filler section, a three-dimensional voxel grid is constructed based on the subgrade geometric parameters. Density values ​​are assigned to each voxel node based on the foundation bearing capacity map. A three-dimensional density field distribution of the foamed soil filler is generated using a spatial interpolation algorithm, including:

[0024] Determine the start and end positions of the foamed soil filling treatment section in the longitudinal direction of the route. Combine the top width of the subgrade and the slope ratio in the subgrade geometric parameters to calculate the outline boundary of each cross section of the treatment section. Set multiple cross section slices at mileage intervals along the longitudinal direction. Generate planar grid points in each cross section slice according to the grid step size. Connect them in space to form a three-dimensional voxel grid.

[0025] Based on the total bearing capacity distribution in the foundation bearing capacity spectrum, weak areas are identified with the bearing critical threshold as the boundary. The reinforcement density corresponding to the weak areas and the reference density corresponding to the non-weak areas are set respectively. Each voxel node in the three-dimensional voxel grid is traversed and it is determined whether it belongs to a weak area, and a differentiated density value is assigned.

[0026] The density adjustment amount is propagated by the voxel nodes on the boundary of the weak area as the diffusion source. The density decay rate is determined according to the bearing capacity gradient between voxel nodes. A density value that gradually transitions from the reinforced density to the reference density is formed in the boundary area. The adjusted density values ​​of each voxel node are organized into a three-dimensional data array according to spatial coordinates to generate the three-dimensional density field distribution of the foamed soil filler.

[0027] Density adjustment is propagated using voxel nodes at the boundary of the weak region as diffusion sources. The density decay rate is determined based on the load-bearing capacity gradient between voxel nodes, forming a density value that gradually transitions from reinforced density to reference density in the boundary region, including:

[0028] Identify the spatial boundaries between weak and non-weak regions in a 3D voxel mesh, mark the voxel nodes located in the weak regions on the spatial boundaries as diffusion source nodes, and determine the connection paths extending from the diffusion source nodes to the non-weak regions on the spatial boundaries.

[0029] The total bearing capacity value of each voxel node on the connection path is obtained from the foundation bearing capacity performance map. The cumulative change of the total bearing capacity is calculated along the connection path. The cumulative change of the total bearing capacity and the spatial length of the connection path are mapped together to the density decay rate of each voxel node.

[0030] The density adjustment amount is propagated from the diffusion source node along the connection path. Each voxel node adjusts the amplitude of the incoming density adjustment amount according to the corresponding density decay rate and transmits it to the downstream node, forming a density propagation sequence that decreases along the path.

[0031] For voxel nodes located on multiple connection paths, the propagation priority of each connection path is determined based on the degree of matching between the cumulative change in the total carrying capacity of each connection path and the current total carrying capacity value of the node. The density adjustment amount transmitted by each connection path is then fused according to the propagation priority to form the gradually transitioning density value of each voxel node on the spatial boundary.

[0032] Boundary constraints and loads are applied to the three-dimensional density field distribution, and settlement values ​​and stability coefficients are obtained through numerical solutions. These values ​​are then compared with the target design values ​​to calculate the deviations. The three-dimensional density field distribution is iteratively adjusted using a Bayesian optimization algorithm until the convergence condition is met, resulting in the optimal density field configuration, including:

[0033] The displacement degrees of freedom of the boundary nodes of the three-dimensional density field distribution are fixed according to the foundation boundary conditions, and the load is converted into nodal forces acting on the nodes in the top region according to the engineering load distribution.

[0034] The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution, and the settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results.

[0035] The settlement value is compared with the target settlement value to obtain the settlement deviation, and the stability coefficient is compared with the target stability coefficient to obtain the stability deviation. A deviation function is then constructed.

[0036] The deviation function is used as the objective function of the Bayesian optimization algorithm. The Bayesian optimization algorithm predicts the adjustment direction of the three-dimensional density field distribution based on the current deviation function value, and modifies the density value of each voxel node in the three-dimensional density field distribution along the adjustment direction to generate a new three-dimensional density field distribution.

[0037] The new three-dimensional density field distribution is re-input into the mechanical response model for numerical solution and calculation of the new deviation function value. It is then determined whether the new deviation function value meets the convergence condition. If the convergence condition is met, the current three-dimensional density field distribution is determined as the optimal density field configuration. If the convergence condition is not met, the three-dimensional density field distribution is further adjusted based on the new deviation function value.

[0038] The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution. The settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results, including:

[0039] The density values ​​of each voxel node are extracted from the three-dimensional density field distribution, and the corresponding elastic modulus and Poisson's ratio are calculated to construct the distribution of foundation material properties. The stress balance control equation is established by combining the spatial topological relationship of the three-dimensional voxel mesh.

[0040] The displacement constraints of the boundary nodes and the nodal forces of the nodes in the top region are substituted into the stress balance control equation as boundary conditions. The stress balance control equation is discretized to form a linear equation system and numerically solved to obtain the displacement field distribution of each voxel node. The vertical displacement component of the monitoring node on the top surface of the foundation is extracted from the displacement field distribution as the settlement value.

[0041] The strain field distribution of each voxel node is calculated based on the displacement field distribution, and the stress field distribution of each voxel node is calculated in combination with the distribution of the foundation material properties. Regions where the stress gradient exceeds the stress threshold are identified as potential failure regions. Within the potential failure regions, the stress transmission path is traced along the principal stress direction to form a slip surface. The anti-slip moment and sliding moment are calculated along the slip surface. All slip surfaces are traversed and the minimum stability coefficient is extracted as the overall stability coefficient of the foundation.

[0042] A second aspect of the present invention provides an AI-assisted three-dimensional design system for foamed soil filler roadbeds, comprising:

[0043] The parameter acquisition unit is used to acquire the roadbed geometric parameters, traffic load spectrum, and foundation soil layer distribution of the roadbed project.

[0044] The section division unit is used to extract features of the foundation soil layer distribution, generate a foundation bearing capacity map, calculate the roadbed surface stress distribution based on the traffic load spectrum, and divide the foamed soil filling treatment section and the conventional filling section according to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution.

[0045] The density field construction unit is used to construct a three-dimensional voxel grid for the foamed soil filler treatment section, in combination with the roadbed geometric parameters, assign density values ​​to each voxel node based on the foundation bearing capacity map, and generate the three-dimensional density field distribution of the foamed soil filler through a spatial interpolation algorithm.

[0046] The optimization solution unit is used to apply boundary constraints and loads to the three-dimensional density field distribution, perform numerical solutions to obtain settlement values ​​and stability coefficients, compare them with the target design values ​​to calculate the deviations, and iteratively adjust the three-dimensional density field distribution through a Bayesian optimization algorithm until the convergence condition is met to obtain the optimal density field configuration.

[0047] The construction output unit is used to generate a three-dimensional spatial distribution model of foamed soil filler according to the optimal density field configuration, and convert the optimal density field configuration into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and output the construction drawings.

[0048] A third aspect of the present invention provides an electronic device, comprising:

[0049] processor;

[0050] Memory used to store processor-executable instructions;

[0051] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0052] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0053] This method enables intelligent optimization of foamed soil fill subgrade design, significantly improving design accuracy and efficiency. By automatically identifying the matching relationship between foundation bearing capacity and pavement stress distribution, it accurately delineates foamed soil and conventional fill sections, avoiding material waste or uneven settlement. Based on a density field generated by three-dimensional voxel mesh and spatial interpolation, it can adaptively adjust the density distribution of foamed soil, ensuring that the filler performance highly matches the foundation bearing capacity requirements, fundamentally reducing the risk of differential subgrade deformation.

[0054] The iterative mechanism of numerical solution and Bayesian optimization automatically converges the settlement value and stability coefficient to the target design value, eliminating the need for repeated manual calculations. This process significantly reduces reliance on professional experience, shortening the traditional design cycle from weeks to days, while ensuring the long-term service performance of the roadbed. The final output of the three-dimensional spatial distribution model and layered construction parameters directly guides on-site filling, avoiding subjective errors and rework costs during construction.

[0055] This method also enables precise prediction and control of foamed soil usage. Through global optimization of the density field, it minimizes the amount of high-density foamed soil used while meeting load-bearing and stability requirements, directly reducing material costs. The dynamic adjustment mechanism allows the design to flexibly respond to complex foundation changes, making it particularly suitable for special conditions such as soft soil and high embankments, significantly improving the settlement resistance and overall durability of subgrade engineering. Attached Figure Description

[0056] Figure 1 This is a flowchart illustrating the AI-assisted three-dimensional design method for foamed soil filler subgrade according to an embodiment of the present invention;

[0057] Figure 2 This is a flowchart illustrating the method for dividing the foamed soil filler treatment section and the conventional filling section in an embodiment of the present invention. Detailed Implementation

[0058] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0059] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0060] Figure 1 This is a flowchart illustrating the AI-assisted 3D design method for foamed soil filler roadbed according to an embodiment of the present invention. The present invention provides an AI-assisted 3D design method for foamed soil filler roadbed, comprising:

[0061] Obtain the geometric parameters of the roadbed, traffic load spectrum, and distribution of the foundation soil layers for the roadbed engineering;

[0062] The features of the foundation soil layer distribution are extracted to generate a foundation bearing capacity map. The roadbed surface stress distribution is calculated based on the traffic load spectrum. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, the foamed soil filling treatment section and the conventional filling section are divided.

[0063] For the foamed soil filler treatment section, a three-dimensional voxel grid is constructed in combination with the roadbed geometric parameters. Based on the foundation bearing capacity map, density values ​​are assigned to each voxel node, and a three-dimensional density field distribution of the foamed soil filler is generated through a spatial interpolation algorithm.

[0064] Boundary constraints and loads are applied to the three-dimensional density field distribution, and the settlement value and stability coefficient are obtained by numerical solution. The deviation is calculated by comparing with the target design value. The three-dimensional density field distribution is iteratively adjusted by Bayesian optimization algorithm until the convergence condition is met, and the optimal density field configuration is obtained.

[0065] A three-dimensional spatial distribution model of foamed soil filler is generated based on the optimal density field configuration, and the optimal density field configuration is converted into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and the construction drawings are output.

[0066] Figure 2 This is a flowchart illustrating the method for dividing the foamed soil fill treatment section and the conventional filling section according to an embodiment of the present invention. The method involves extracting features from the foundation soil layer distribution to generate a foundation bearing capacity map. Based on the traffic load spectrum, the roadbed surface stress distribution is calculated. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, the foamed soil fill treatment section and the conventional filling section are divided, including:

[0067] The distribution of the foundation soil layers is processed by layering and slicing, and the rock and soil type identification and thickness data of each soil layer are extracted to construct a vertical profile model containing soil layer sequence and burial depth information;

[0068] Bearing characteristic parameters are assigned to each soil layer in the vertical profile model, the ultimate bearing capacity and allowable settlement of each soil layer are calculated, and the distribution of total bearing capacity and total settlement potential of the foundation are generated through vertical cumulative calculation, and integrated into a foundation bearing performance map.

[0069] The dominant frequency load component is extracted by performing spectral analysis on the traffic load spectrum. The dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at the interfaces of each soil layer is calculated using the elastic layered system theory to obtain the stress distribution on the roadbed surface.

[0070] The difference between the total bearing capacity distribution in the foundation bearing capacity performance map and the stress distribution on the roadbed surface is calculated to generate a bearing surplus distribution.

[0071] Identify the areas with negative values ​​in the load-bearing surplus distribution, classify the areas with insufficient load-bearing capacity according to the amount of load loss, and determine the treatment intensity corresponding to each level. Areas with treatment intensity greater than the treatment threshold are divided into foamed soil filling treatment sections, and the remaining areas are divided into conventional filling sections.

[0072] When performing layered slicing of the foundation soil distribution, the foundation profile is divided into several continuous horizontal thin layers along the vertical direction according to the soil layer interfaces revealed by the actual boreholes. Each layer corresponds to a soil type. For each thin layer, its soil type identifier (such as soft clay, silty clay, gravel layer, bedrock, etc.) and corresponding layer thickness data are extracted. At the same time, the burial depth coordinates of the top and bottom surfaces of each layer are recorded, thereby constructing a vertical profile model containing soil layer sequence and burial depth information. This model expresses the heterogeneous characteristics of the foundation in the vertical direction in a discretized form, providing a structured data foundation for the subsequent assignment of bearing capacity parameters. In the horizontal direction, based on the planar distribution of the exploration points, the vertical profile model of each borehole is spatially expanded using Kriging interpolation or inverse distance weighting methods to form a three-dimensional stratigraphic model covering the entire roadbed area, ensuring that the spatial variability of the stratigraphy is fully expressed.

[0073] When assigning bearing capacity parameters to each soil layer in the vertical profile model, key parameters such as internal friction angle, cohesion, compression modulus, and void ratio are configured for each soil type based on laboratory and in-situ test results. The ultimate bearing capacity of each soil layer is calculated layer by layer based on Terzaghi's ultimate bearing capacity theory or modified standard formulas. subscript Indicates the first One soil layer. Allowable settlement. The bearing capacity is obtained by integrating the compressive modulus and additional stress within the depth of the layer, reflecting the deformation potential of the layer under the design load level. Through vertical cumulative calculation, the ultimate bearing capacity of each layer is converted to the top surface of the foundation according to the stress diffusion angle, and the total bearing capacity distribution of the foundation is obtained after comprehensive superposition. The allowable settlement of each layer By summing the results layer by layer, the total settlement potential distribution of the foundation is obtained. .Will and The data are integrated into a foundation bearing capacity map, which is stored in the form of a planar grid. Each grid node carries the total bearing capacity value and total settlement potential value of the corresponding planar location, intuitively reflecting the spatial bearing capacity differences of the foundation.

[0074] When performing spectral analysis on traffic load spectra, measured or statistically obtained axle load time history data are input into a Fast Fourier Transform (FFT) to decompose the amplitude and phase information corresponding to each frequency component. Several dominant frequency load components with the most concentrated energy are identified from the spectrum; these components represent the dynamic load modes that have the most significant impact on the subgrade structure. The extracted dominant frequency load components are then equivalent to static uniformly distributed loads or moving loads, acting on the top surface of the subgrade. Based on the theory of elastic layered systems, the subgrade and foundation are considered as a multi-layered elastic half-space system, utilizing the elastic modulus of each layer. Compared to Poisson (subscript) Indicates the first A continuity equation is established between each elastic layer, and the vertical, horizontal, and shear stress components at each layer interface are solved by numerical integration. The stress decreases with increasing depth during downward transmission, and the vertical stress value is finally extracted at the subgrade surface, forming the subgrade surface stress distribution. This distribution is stored in the same planar grid coordinate system as the foundation bearing capacity map to ensure the spatial correspondence of subsequent calculations.

[0075] Distribution of total bearing capacity in the foundation bearing capacity map Stress distribution on roadbed surface Perform grid-by-grid interpolation to obtain the load surplus distribution. The calculation relationship is . When the value is positive, it indicates that the bearing capacity of the foundation at that location is sufficient and no special treatment is required; A negative value indicates that the bearing capacity of the foundation at that location is insufficient, and the stress transmitted to the roadbed surface exceeds the limit that the foundation can withstand, posing a risk of instability or excessive settlement. The bearing capacity surplus distribution map visually presents the spatial distribution pattern of the foundation's strength and weakness areas on a plane, providing a quantitative basis for subsequent section division.

[0076] Identify the distribution of carrying capacity surplus The absolute value of these undercapacitated areas As a measure of missing data, it is classified and processed according to the magnitude of the missing data. For example, when... When the temperature is in the range [0,1), it is classified as a light treatment level, and low-density foamed soil filler can meet the weight reduction and reinforcement requirements; when... When the range is [1,2), it is classified as a moderate treatment level, requiring the use of medium-density foamed soil filler and an appropriate expansion of the treatment range; when When classified as a severe treatment level, it requires the use of low-density foamed soil filler with high weight reduction effect, combined with a widened treatment range. The threshold values ​​for each grade are determined based on engineering experience and specification requirements, and can be flexibly adjusted according to foundation conditions and design safety factors for different projects. For grid cells with a bearing capacity deficiency exceeding the threshold, they are classified into foamed soil filler treatment sections, and the corresponding treatment intensity level is marked; for... In areas where the foundation bearing capacity meets the design requirements, these areas are designated as conventional filling sections, and can be constructed using ordinary fill material and standard processes.

[0077] The boundary between the foamed soil filling treatment section and the conventional filling section is represented by a vector polygon on the plan view. The boundary line is fitted by connecting grid nodes with positive and negative conversions of adjacent bearing surplus values, and smoothing is performed when necessary to facilitate construction layout. The graded boundaries within the treatment section are also represented by superimposed vectors to form a multi-level section division result map. After being overlaid with the subgrade plan design drawing, this result map can clearly indicate the spatial range of each treatment level on the cross-section and longitudinal section of the subgrade, providing accurate spatial boundary basis for the subsequent construction of three-dimensional voxel mesh and density field assignment. Through the complete process described above, from extracting stratum features, calculating bearing capacity, analyzing load transfer, to calculating bearing surplus difference, a quantitative assessment of the subgrade foundation conditions is achieved, ensuring that the delineation of the foamed soil filling treatment section has sufficient mechanical basis and avoiding the problem of conservative or aggressive treatment ranges caused by the reliance on engineers' experience judgment in traditional methods.

[0078] The dominant frequency load component is extracted by spectral analysis of the traffic load spectrum. This dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at each soil layer interface is calculated using elastic layered system theory to obtain the stress distribution on the roadbed surface, including:

[0079] The traffic load spectrum is segmented according to a preset time window and subjected to spectrum transformation to obtain a local spectrum. The frequency components that recur in all local spectra are counted, and the frequency component with the most recurrences is determined as the stable main frequency. The load amplitude corresponding to the stable main frequency is extracted as the main frequency load component.

[0080] The dominant frequency load component is divided into distributed load nodes according to the transverse width and longitudinal length of the roadbed top surface, and the load intensity is determined according to the distance of each load node from the load action center.

[0081] Extract the compression modulus and lateral compression coefficient of each soil layer in the vertical profile model, calculate the strain response coefficient of each soil layer under vertical load, and determine the damping characteristics of each soil layer.

[0082] Select the load nodes on the top surface of the subgrade and use the corresponding load intensity as the initial stress input to the first soil layer. Calculate the stress attenuation within the layer and its transmission to the interface based on the damping characteristics of the first soil layer. Determine the degree of interface abrupt change by the difference in soil properties on both sides of the interface.

[0083] When the degree of interface mutation exceeds the interface sensitivity threshold, an interface reflection loss coefficient is introduced. Combined with the stress at the interface, the stress continues to propagate to the next soil layer. The stress attenuation and interface treatment are repeated until the top surface of the foundation is reached. The stresses propagated from all load nodes to the top surface of the foundation are spatially superimposed to obtain the stress distribution of the roadbed surface.

[0084] Traffic load spectra are typically composed of measured axle load data or simulation data, and their time series contains both random fluctuations and periodic dominant components. To extract representative dominant frequency load components, the load spectrum needs to be segmented according to a preset time window. The selection of the time window length should consider both the periodicity and statistical stability of the load, generally taking an integer multiple of the duration of a single vehicle passage to avoid the truncation effect from biasing the spectrum analysis results. A Fast Fourier Transform is applied to the load time history data within each time window to obtain the corresponding local spectrum, i.e., the amplitude distribution of each frequency component. The entire local spectrum is traversed, and the frequency component with the most repetitions is counted across all windows. The frequency component with the most repetitions is determined as the stable dominant frequency. Extracting a stable core frequency The average of the amplitudes corresponding to each local frequency spectrum is taken as the amplitude of the dominant frequency load component. This amplitude represents the core excitation intensity that generates a sustained dynamic response from traffic load in the roadbed structure.

[0085] The main frequency load component According to the transverse width of the top surface of the roadbed and longitudinal length Spatial discretization is performed to establish a load node mesh on the top surface of the roadbed. The spacing of the load nodes is determined based on the roadbed width and calculation accuracy requirements, typically with a lateral spacing of 0.5m to 1.0m and a longitudinal spacing of 1.0m to 2.0m. The load intensity at each load node is then determined. The distance of the node from the center of the load. Correspondingly, the load intensity decreases with distance from the center. The load intensity is distributed according to the distance attenuation relationship, satisfying that the integral of the load intensity at all nodes equals the total value of the dominant frequency load components. The constraints are used to ensure the global balance of the load. The load center is usually taken as the intersection of the transverse centerline of the top surface of the roadbed and the longitudinal calculation section. For a two-lane roadbed, two load centers can be set respectively, corresponding to the wheel track positions of the left and right lanes.

[0086] In the vertical profile model, the compression modulus of each soil layer and confined compressibility Obtained through indoor or in-situ testing, among which Reflects the vertical compressive stiffness of the soil layer under lateral constraints. Characterizes the volumetric strain caused by a unit stress increment. Based on the above parameters, calculate the... Strain response coefficient of soil layer under vertical load This coefficient comprehensively reflects the soil layer's efficiency in transmitting vertical stress. Simultaneously, the damping characteristics of each soil layer are determined based on its dynamic shear modulus and damping ratio. The energy dissipation capacity of a soil layer under dynamic load is determined by resonant column tests or dynamic triaxial tests. Damping characteristics determine the amplitude attenuation rate of stress waves propagating within the soil layer and are the core parameter for calculating stress attenuation within the layer.

[0087] Stress propagation calculations are performed layer by layer downwards, starting from the first soil layer on the top surface of the roadbed. A specific load node on the top surface of the roadbed is selected. The corresponding load strength The initial stress is input to the top surface of the first soil layer. The stress attenuation within the first soil layer is calculated using an exponential attenuation model, with the attenuation amount varying with the layer thickness. Damping ratio The stress at the bottom surface of the first soil layer (i.e., the first soil layer interface) after in-layer attenuation is denoted as... At the interface, the degree of abrupt change in interface abruptness needs to be determined based on the ratio of the compression moduli of the soil layers on both sides of the interface. An interface abrupt change index is defined. It is the absolute value of the deviation between the ratio of the compression modulus of two adjacent soil layers and 1, i.e. ,when Exceeding the preset interface sensitivity threshold At that time, it was considered that there was a significant abrupt change in the mechanical properties of the interface, and it was necessary to introduce an interface reflection loss coefficient. Interface reflection loss coefficient Determined by the transmission coefficient in wave theory, its value ranges from 0 to 1, reflecting the energy loss caused by impedance mismatch when a stress wave crosses an interface. The actual downward propagating stress at the interface is... When the degree of interface mutation does not exceed the threshold, Taking 1 means that the stress passes through the interface without loss.

[0088] Will As the top surface input stress of the second soil layer, the in-layer attenuation calculation is repeated to obtain the bottom surface stress of the second soil layer. Then determine the degree of abrupt change at the second interface and determine the corresponding interface reflection loss coefficient. This process continues layer by layer downwards. With each soil layer, the stress undergoes two stages: in-layer damping attenuation and interfacial transmission loss, until it reaches the top surface of the foundation, i.e., the bottom surface of the last soil layer, where the load node is obtained. Vertical stress component contributed at the top surface of the foundation Repeat the above propagation calculation process for all load nodes on the top surface of the subgrade to obtain the stress contribution value of each node on the top surface of the subgrade.

[0089] The total stress at a point on the top surface of the foundation is obtained by superimposing the stress contributions from all load nodes at that point. When superimposing, the spatial diffusion effect of stress must be considered, i.e., the stress contributions from all load nodes. At its horizontal distance The stress component generated at the top surface of the foundation varies with The stress decreases as the load increases. Spatial superposition is corrected using Boussinesq stress diffusion theory, which weights and sums the stress contributions of each load node on the top surface of the foundation according to spatial location, ultimately obtaining the stress distribution of the roadbed surface across the entire top surface of the foundation. ,in and These represent the horizontal and vertical coordinates in the plane coordinate system of the foundation top surface, respectively. The stress distribution is stored in a two-dimensional grid, with the grid resolution consistent with the spacing between the load nodes on the roadbed top surface, facilitating subsequent spatial matching calculations with the foundation bearing capacity map.

[0090] Stress distribution on roadbed surface The calculation results intuitively reflect the spatial non-uniformity of traffic loads after being transferred from the roadbed structure to the top surface of the foundation. In the cross-sectional direction, the stress distribution exhibits a pattern of high stress in the center and low stress on both sides, corresponding to the location of the driving lanes; in the longitudinal direction, the stress distribution fluctuates with changes in the thickness and properties of the foundation soil layers. By comparing the distribution of ultimate bearing capacity and allowable settlement with the distribution of foundation bearing capacity in the foundation bearing capacity spectrum, areas with insufficient bearing capacity can be identified, providing a quantitative basis for the subsequent division of foamed soil filling treatment sections and ensuring that the spatial accuracy of the section division matches the actual bearing state of the foundation.

[0091] For the foamed soil filler section, a three-dimensional voxel grid is constructed based on the subgrade geometric parameters. Density values ​​are assigned to each voxel node based on the foundation bearing capacity map. A three-dimensional density field distribution of the foamed soil filler is generated using a spatial interpolation algorithm, including:

[0092] Determine the start and end positions of the foamed soil filling treatment section in the longitudinal direction of the route. Combine the top width of the subgrade and the slope ratio in the subgrade geometric parameters to calculate the outline boundary of each cross section of the treatment section. Set multiple cross section slices at mileage intervals along the longitudinal direction. Generate planar grid points in each cross section slice according to the grid step size. Connect them in space to form a three-dimensional voxel grid.

[0093] Based on the total bearing capacity distribution in the foundation bearing capacity spectrum, weak areas are identified with the bearing critical threshold as the boundary. The reinforcement density corresponding to the weak areas and the reference density corresponding to the non-weak areas are set respectively. Each voxel node in the three-dimensional voxel grid is traversed and it is determined whether it belongs to a weak area, and a differentiated density value is assigned.

[0094] The density adjustment amount is propagated by the voxel nodes on the boundary of the weak area as the diffusion source. The density decay rate is determined according to the bearing capacity gradient between voxel nodes. A density value that gradually transitions from the reinforced density to the reference density is formed in the boundary area. The adjusted density values ​​of each voxel node are organized into a three-dimensional data array according to spatial coordinates to generate the three-dimensional density field distribution of the foamed soil filler.

[0095] After determining the spatial extent of the foamed soil filling treatment section, this section needs to be discretized into a three-dimensional voxel mesh suitable for numerical calculations. Along the longitudinal direction of the route, the longitudinal coverage of the mesh is determined based on the start and end station numbers of the treatment section. Combining the roadbed top width and slope ratio from the roadbed geometry parameters, the outline boundary of each cross-section is calculated. Specifically, the roadbed top width determines the horizontal extent of the top of the cross-section, and the slope ratio and filling height together determine the toe position, thus determining the outer contour polygon of each cross-section. Multiple cross-section slices are set at equal intervals along the longitudinal direction. The selection of the mileage interval should comprehensively consider the longitudinal variation of the foundation soil layers and the required calculation accuracy. In sections with drastic changes in geological conditions, the slice spacing can be appropriately increased. Within each cross-section slice, planar grid points are evenly distributed within the outline boundary according to a preset grid step size. The selection of the grid step size balances computational efficiency and spatial resolution, and independent step size parameters are usually set in both the horizontal and vertical directions. The planar grid points at corresponding positions on adjacent cross-sectional slices are connected sequentially in the longitudinal direction to form a three-dimensional voxel grid covering the entire processing section. Each voxel unit is enclosed by eight nodes, and the spatial coordinates of the nodes are determined by the lateral position, vertical elevation, and longitudinal mileage.

[0096] After constructing the 3D voxel mesh, initial density values ​​need to be assigned to each voxel node based on the foundation bearing capacity map. The foundation bearing capacity map records the distribution of total bearing capacity, reflecting the differences in bearing capacity of the foundation on the plane. Using the critical bearing capacity threshold as a criterion, the distribution of total bearing capacity is binarized for identification: areas with a total bearing capacity below this threshold are identified as weak areas, and the remaining areas are considered non-weak areas. The determination of the critical bearing capacity threshold comprehensively considers the stress distribution of the subgrade surface and the safety reserve requirements of the foundation bearing capacity. Typically, the product of the maximum stress on the subgrade surface under design load and a certain safety factor is taken as the reference benchmark. For weak areas, a higher reinforcement density value is set to give the foamed soil filler greater self-weight compensation and deformation coordination capabilities in these areas; for non-weak areas, a relatively lower benchmark density value is set to reduce the cost of the filler and the overall self-weight of the subgrade while meeting the bearing capacity requirements. All voxel nodes in the 3D voxel mesh are traversed, and the planar coordinates of each node are... The system queries the foundation bearing capacity map to determine whether the location is a weak area, and assigns the corresponding reinforcement density or reference density accordingly, thus completing the global assignment of the initial density value.

[0097] After the initial density assignment, a sudden density transition occurs between weak and non-weak regions. This transition is difficult to accurately achieve in actual construction and introduces stress concentration effects into the numerical solution, affecting the accuracy of the calculation results. Therefore, a density adjustment is propagated outward from the voxel nodes at the boundary of the weak region, forming a smoothly transitioning density gradient zone at the interface. During the diffusion process, the decay rate of the density adjustment is determined by the bearing capacity gradient between the voxel nodes: locations with a larger bearing capacity gradient indicate drastic changes in foundation performance over a short distance, corresponding to a faster density decay rate and a narrower transition zone; locations with a smaller bearing capacity gradient indicate gradual changes in foundation performance, corresponding to a slower density decay rate and a wider transition zone. Let the density adjustment at the diffusion source node be... The spatial distance of a certain node from the diffusion source is The magnitude of the carrying capacity gradient at this node is The density adjustment amount at that node Determine as follows:

[0098]

[0099] in This is the attenuation control coefficient, and its value is calibrated according to the smoothness requirements of the transition band. When When the density is less than the preset attenuation cutoff threshold, the node is considered to have exceeded the transition zone range, and no further density adjustment is applied. The density adjustment amount is superimposed on the initial density value of each node, so that the node density near the boundary of the weak region gradually transitions from the reinforcement density value to the reference density value, forming a continuous and smooth density distribution.

[0100] After the density adjustment of the transition zone is completed, the physical rationality of the density values ​​of all voxel nodes is verified. The density value of the foamed soil filler should be within its engineering-producible range, and the lower limit is generally not lower than [the specified value]. (Corresponding to the minimum foaming rate formula ratio), the upper limit shall not exceed (Corresponding to the density limit of plain concrete), for nodes that exceed the range, their density values ​​are truncated to the corresponding boundary values, and the truncation position is recorded for key attention during subsequent optimization iterations.

[0101] After assigning and validating the density values ​​of all nodes, the density values ​​of each voxel node are organized in an ordered manner according to their spatial coordinates to construct a three-dimensional data array. The index dimensions of the three-dimensional data array correspond to the horizontal grid numbers. Vertical grid numbering and longitudinal cross-sectional slice numbering Each index position stores the spatial coordinates of the corresponding node. With density value This three-dimensional data array represents the three-dimensional density field distribution of the foamed soil filler, serving as the fundamental input for subsequent numerical solutions and Bayesian optimization iterations. The three-dimensional density field distribution comprehensively describes the density configuration of the foamed soil filler at various spatial locations within the treated area, encompassing the reinforced density core zone in weak areas, the baseline density zone in non-weak areas, and the gradual transition zone between the two. It effectively reflects the impact of spatial differences in foundation bearing capacity on the filler density configuration, providing a spatially resolved material parameter field for subsequent settlement analysis and stability verification.

[0102] In practical engineering applications, the appropriate selection of the mesh step size and cross-sectional slice spacing of the 3D voxel mesh has a significant impact on both computational accuracy and efficiency. For engineering scenarios where the foundation soil layers are uniformly distributed laterally and change slowly longitudinally, a larger mesh step size and slice spacing can be used to improve computational efficiency. For scenarios where the foundation has weak interlayers, locally ultra-soft soil areas, or complex strata, the mesh step size and slice spacing should be appropriately reduced to ensure that the identification accuracy of weak areas and the spatial resolution of the density transition zone meet design requirements. The value of the bearing capacity threshold should also be adjusted according to the design safety level of the specific project. For high-grade roadbed projects such as highways and urban expressways, the threshold should be appropriately increased to expand the identification range of weak areas, thereby obtaining a more conservative reinforcement density configuration area and ensuring the long-term operational stability of the roadbed.

[0103] Density adjustment is propagated using voxel nodes at the boundary of the weak region as diffusion sources. The density decay rate is determined based on the load-bearing capacity gradient between voxel nodes, forming a density value that gradually transitions from reinforced density to reference density in the boundary region, including:

[0104] Identify the spatial boundaries between weak and non-weak regions in a 3D voxel mesh, mark the voxel nodes located in the weak regions on the spatial boundaries as diffusion source nodes, and determine the connection paths extending from the diffusion source nodes to the non-weak regions on the spatial boundaries.

[0105] The total bearing capacity value of each voxel node on the connection path is obtained from the foundation bearing capacity performance map. The cumulative change of the total bearing capacity is calculated along the connection path. The cumulative change of the total bearing capacity and the spatial length of the connection path are mapped together to the density decay rate of each voxel node.

[0106] The density adjustment amount is propagated from the diffusion source node along the connection path. Each voxel node adjusts the amplitude of the incoming density adjustment amount according to the corresponding density decay rate and transmits it to the downstream node, forming a density propagation sequence that decreases along the path.

[0107] For voxel nodes located on multiple connection paths, the propagation priority of each connection path is determined based on the degree of matching between the cumulative change in the total carrying capacity of each connection path and the current total carrying capacity value of the node. The density adjustment amount transmitted by each connection path is then fused according to the propagation priority to form the gradually transitioning density value of each voxel node on the spatial boundary.

[0108] After assigning initial density field values ​​to the 3D voxel mesh, abrupt density changes often occur at the boundary between weak and non-weak regions. Directly dividing the density values ​​of the two types of regions with hard boundaries can lead to stress concentration and uneven settlement within the roadbed. Therefore, it is necessary to construct a smoothly transitioning density gradient zone at the boundary.

[0109] When identifying the spatial boundaries between weak and non-weak regions in a 3D voxel mesh, the neighborhood of each voxel node is scanned sequentially. If a voxel node belongs to a weak region and at least one of its adjacent voxel nodes belongs to a non-weak region, the voxel node is determined to be located on the spatial boundary and marked as a diffusion source node. The marking process for diffusion source nodes covers all six adjacency directions of the 3D voxel mesh, ensuring that diagonally adjacent boundary nodes can also be completely identified. After marking the diffusion source nodes, starting from each diffusion source node, the search extends along one side of the non-weak region to search for a sequence of non-weak region voxel nodes that are directly connected to the diffusion source node in space, forming a connection path. The search direction for the connection path preferentially extends along the direction of the fastest increase in bearing capacity gradient. The path termination condition is that the path length exceeds a preset penetration depth threshold, or the total bearing capacity value of the node at the end of the path has reached the baseline level of the non-weak region.

[0110] Extract the total bearing capacity value at each voxel node on the connection path from the foundation bearing capacity map, and denote the i-th voxel node as the i-th voxel node. The first connection path The total carrying capacity at each node is [value]. The total carrying capacity at the path starting point (i.e., the diffusion source node) is... The total carrying capacity at the end node of the path is The cumulative change in total carrying capacity is calculated along the connecting path, defined as the change from the path origin to the [missing information]. The cumulative change of each node is This value reflects the degree of recovery of bearing capacity from the weak boundary to the interior of the non-weak area. Simultaneously, the spatial length of the connecting path is recorded, calculated by summing the Euclidean distances between adjacent nodes on the path, from the path start point to the [missing value]. Path arc length of each node .Will and Co-mapping is the first Density decay rate at each node The mapping relationship is ,in To normalize the reference path length, To prevent extremely small positive numbers with a denominator of zero. When the cumulative change in bearing capacity is large and the path arc length is short, the attenuation rate is large, indicating that the density transition needs to be completed within a shorter space; when the cumulative change in bearing capacity is small and the path arc length is long, the attenuation rate is small, indicating that the density can transition slowly over a longer space.

[0111] The initial density adjustment at the diffusion source node as the density adjustment propagates downstream along the connection path from the diffusion source node. The density is determined by the difference between the reinforcement density and the reference density in the weak region where the node is located. (The last part, "on the path," appears to be a typo and can be left as is.) Density adjustment received by each node It is obtained by multiplying the density adjustment amount from the upstream node by the attenuation factor at that node. The attenuation factor is defined as follows: ,in For the first The node and the first The spatial step size between nodes is determined. Consequently, the density adjustment amount decreases node by node along the path, forming a density propagation sequence. The density adjustment amount at the end of the path approaches zero, smoothly connecting with the baseline density of the non-weak region. After receiving the density adjustment amount, each voxel node superimposes it onto the baseline density value of that node to obtain the actual density value of the transition region. Simultaneously, the superposition result is truncated to ensure that the final density value remains consistent with the baseline density of the transition region. and Within the defined physical reasonable range.

[0112] In practical 3D voxel meshes, voxel nodes located within boundary regions often lie on multiple connection paths simultaneously, and the density adjustments transmitted through different paths need to be fused. For voxel nodes simultaneously located on... For a voxel node on a connection path, read the actual total bearing capacity value at that node from the foundation bearing capacity map. And extract the range of cumulative changes in the total carrying capacity of each connection path. .Will Matching the carrying capacity variation range of each connection path, and calculating... Falling in The degree of matching within the range of carrying capacity variation of the path The higher the matching degree, the higher the propagation priority of that path to the current node. Specifically, the matching degree... Defined as ,in For the first The reference carrying capacity value corresponding to the current node for each path is obtained by linear interpolation of the total carrying capacity values ​​of the adjacent nodes on the path. For all The matching degree of each path is normalized to obtain the fusion weight of each path. The density adjustment amounts transmitted to the current node from each path are weighted and summed according to the fusion weight to obtain the final density adjustment amount for that node. .

[0113] The fused density adjustment is superimposed on the baseline density value of the current node to form the gradually transitioning density values ​​of each voxel node on the spatial boundary. Because multi-path fusion fully considers the bearing capacity characteristics of each connecting path and the actual bearing capacity state of the current node, the density distribution within the transition area is spatially continuous and smooth, avoiding local density anomalies caused by single-path propagation. After the entire density gradual transition process is completed, the voxel node density values ​​in the boundary area start from the reinforced density at the weak area boundary, smoothly decay along each connecting path, and finally recover to the baseline density level in the non-weak area, forming a spatially continuous density field that conforms to the foundation bearing capacity characteristics, providing a physically reasonable initial density field configuration for subsequent numerical solutions.

[0114] Boundary constraints and loads are applied to the three-dimensional density field distribution, and settlement values ​​and stability coefficients are obtained through numerical solutions. These values ​​are then compared with the target design values ​​to calculate the deviations. The three-dimensional density field distribution is iteratively adjusted using a Bayesian optimization algorithm until the convergence condition is met, resulting in the optimal density field configuration, including:

[0115] The displacement degrees of freedom of the boundary nodes of the three-dimensional density field distribution are fixed according to the foundation boundary conditions, and the load is converted into nodal forces acting on the nodes in the top region according to the engineering load distribution.

[0116] The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution, and the settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results.

[0117] The settlement value is compared with the target settlement value to obtain the settlement deviation, and the stability coefficient is compared with the target stability coefficient to obtain the stability deviation. A deviation function is then constructed.

[0118] The deviation function is used as the objective function of the Bayesian optimization algorithm. The Bayesian optimization algorithm predicts the adjustment direction of the three-dimensional density field distribution based on the current deviation function value, and modifies the density value of each voxel node in the three-dimensional density field distribution along the adjustment direction to generate a new three-dimensional density field distribution.

[0119] The new three-dimensional density field distribution is re-input into the mechanical response model for numerical solution and calculation of the new deviation function value. It is then determined whether the new deviation function value meets the convergence condition. If the convergence condition is met, the current three-dimensional density field distribution is determined as the optimal density field configuration. If the convergence condition is not met, the three-dimensional density field distribution is further adjusted based on the new deviation function value.

[0120] After constructing the three-dimensional density field distribution, boundary constraints and loads need to be applied. Mechanical response indices are obtained through numerical solutions, and the density field is iteratively adjusted using a Bayesian optimization algorithm until the design objectives are met. The method of applying boundary constraints directly affects the accuracy of the numerical solution; therefore, the boundary nodes of the three-dimensional density field distribution need to be processed according to the actual foundation boundary conditions. For nodes located on the foundation surface, their vertical displacement degrees of freedom are fixed; for nodes located on the foundation side, their horizontal displacement degrees of freedom are fixed; for internal nodes, all displacement degrees of freedom are retained, allowing them to deform freely under load. This method of applying displacement constraints in different regions enables the numerical model to accurately reflect the deformation boundary conditions of the roadbed foundation under actual working conditions, avoiding distortion of the solution results due to improper boundary treatment.

[0121] The transformation of the engineering load distribution is another crucial preliminary step in numerical solution. The distributed loads described by the traffic load spectrum are discretized into concentrated nodal forces acting on each top node, based on the spatial correspondence between the load application area and the nodes in the top region. The magnitude of the nodal forces is determined by the product of the area of ​​influence controlled by each node and the load intensity of that area, ensuring that the discretized nodal force system remains consistent with the original distributed load in terms of resultant force and resultant moment. For dynamic load components, they are converted into static nodal forces according to the principle of equivalent static force, and the conversion factor is determined based on the standard of dynamic amplification factor for roadbed engineering, thereby incorporating the dynamic load effect into the static numerical solution framework.

[0122] Apply the above boundary constraints and nodal forces A mechanical response model based on a three-dimensional density field distribution is input and numerically solved. The mechanical response model is established using the finite element method, and the elastic modulus of each voxel element is calculated. Compared to Poisson The density value at the voxel node is determined through material constitutive relation mapping. There is a power function relationship between the density and elastic modulus of foamed soil filler; higher density results in a higher elastic modulus, stronger stiffness, and higher bearing capacity. After numerical solution, the vertical displacement of each node on the top surface of the foundation is extracted from the solution results, and the maximum value is taken as the representative settlement value. Simultaneously, the strength reduction method is used to conduct overall stability analysis of the foundation, gradually reducing the shear strength parameters of the foamed soil filler and the foundation soil until a through sliding surface appears in the foundation. The reduction factor at this point is the stability coefficient. .

[0123] Will With the target settlement value By comparison, the settlement deviation was obtained. ;Will With the target stability coefficient By comparison, the stability deviation was obtained. To unify the two types of deviations into the same optimization objective, a comprehensive deviation function is constructed. Its expression is ,in This is the weighting coefficient for settlement deviation. This is the stability deviation weighting coefficient, and the sum of the two is 1. It is assigned a value based on the importance level of the project and the design focus. When the roadbed has stricter requirements for settlement control, it should be appropriately increased. The value of ; when stability control is the dominant design constraint, appropriately increase The value of . Deviation function The smaller the value, the closer the current three-dimensional density field distribution is to the design target.

[0124] Bayesian optimization algorithm with bias function As the objective function, a surrogate model is constructed to guide the adjustment direction of the 3D density field distribution. The surrogate model is established using Gaussian process regression, using the parameterized representation of the 3D density field distribution and the corresponding deviation function values ​​from each iteration as training samples to fit the response surface of the objective function. At the beginning of each iteration, the Bayesian optimization algorithm calculates the acquisition function based on the current surrogate model. The acquisition function comprehensively considers the expected improvement in the currently known region and the uncertainty in the unexplored region, achieving a balance between the two, thereby determining the direction and magnitude of the next density field adjustment.

[0125] Following the adjustment direction predicted by the Bayesian optimization algorithm, the density values ​​of each voxel node in the three-dimensional density field distribution are modified. Specifically, for areas with large settlement deviations, the density values ​​of voxel nodes in these areas are appropriately increased to enhance local stiffness and reduce settlement. For areas near potential sliding surfaces with insufficient stability coefficients, the density values ​​of voxel nodes in these areas are appropriately increased to enhance shear strength and improve overall stability. The adjustment amount of the density values ​​is determined by the acquisition function output of the Bayesian optimization algorithm and is subject to upper and lower density limits. and The constraints ensure that the adjusted density value remains within the physically feasible range. After adjustment, a new three-dimensional density field distribution is generated, and it is re-input into the mechanical response model for numerical solution to calculate the new settlement value. Stability coefficient and the new deviation function value .

[0126] The new deviation function value With convergence threshold The comparison is performed to determine whether the convergence condition is met. The determination of the convergence condition involves two aspects: first, whether the deviation function value itself is sufficiently small, i.e. Secondly, is the change in the deviation function value between two adjacent iterations small enough? ,in The deviation function value from the previous iteration. This is the threshold for iterative convergence accuracy. Both conditions must be met simultaneously for the algorithm to be considered convergent. If the convergence condition is met, the current 3D density field distribution is determined as the optimal density field configuration and output to the subsequent construction drawing generation stage. If the convergence condition is not met, the new deviation function value and the corresponding density field parameters are used as new training samples to supplement the surrogate model, the Gaussian process regression model is updated, and the acquisition function is recalculated based on the updated surrogate model to determine the next adjustment direction, continuing iterative optimization.

[0127] To prevent the algorithm from getting trapped in local optima, a maximum number of iterations is set during the Bayesian optimization process. When the number of iterations exceeds If the convergence condition is still not met, the three-dimensional density field distribution with the smallest deviation function value is selected from the results of each iteration as the optimal density field configuration output. The settlement value and stability coefficient of this configuration are marked on the construction drawings for engineers to review. The entire optimization process fully utilizes the advantages of Bayesian optimization algorithms in efficient search with small samples. Compared with traditional exhaustive methods or gradient descent methods, it can find the optimal density field configuration that meets the engineering design requirements with fewer numerical solutions, significantly reducing computational costs and improving the efficiency of three-dimensional design.

[0128] The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution. The settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results, including:

[0129] The density values ​​of each voxel node are extracted from the three-dimensional density field distribution, and the corresponding elastic modulus and Poisson's ratio are calculated to construct the distribution of foundation material properties. The stress balance control equation is established by combining the spatial topological relationship of the three-dimensional voxel mesh.

[0130] The displacement constraints of the boundary nodes and the nodal forces of the nodes in the top region are substituted into the stress balance control equation as boundary conditions. The stress balance control equation is discretized to form a linear equation system and numerically solved to obtain the displacement field distribution of each voxel node. The vertical displacement component of the monitoring node on the top surface of the foundation is extracted from the displacement field distribution as the settlement value.

[0131] The strain field distribution of each voxel node is calculated based on the displacement field distribution, and the stress field distribution of each voxel node is calculated in combination with the distribution of the foundation material properties. Regions where the stress gradient exceeds the stress threshold are identified as potential failure regions. Within the potential failure regions, the stress transmission path is traced along the principal stress direction to form a slip surface. The anti-slip moment and sliding moment are calculated along the slip surface. All slip surfaces are traversed and the minimum stability coefficient is extracted as the overall stability coefficient of the foundation.

[0132] After extracting the density values ​​of each voxel node from the three-dimensional density field distribution, the density values ​​are converted into the corresponding elastic modulus and Poisson's ratio based on the density-mechanical property relationship curve of the foamed soil material. Specifically, the elastic modulus of the foamed soil filler increases nonlinearly with increasing density, and the conversion can be accomplished through a pre-calibrated material constitutive relation mapping function; the Poisson's ratio remains approximately stable within a certain density range. For areas of lightweight foamed soil with lower density, the Poisson's ratio is set relatively small, while for areas with higher density, the Poisson's ratio is appropriately increased. The elastic modulus and Poisson's ratio of each voxel node are combined into a foundation material property distribution matrix, so that each node in the entire three-dimensional voxel mesh carries unique material parameter information, thereby forming a material description system for heterogeneous foundations.

[0133] When establishing the stress balance governing equations, the spatial topology of the three-dimensional voxel mesh is considered, and each voxel element is discretized using an eight-node hexahedral isoparametric element. For any voxel element, its stiffness matrix is ​​determined by the element's elastic modulus. Compared to Poisson It was jointly decided that the element stiffness matrix would be formed by numerical integration within each element using a Gaussian integration scheme, and then assembled into a global stiffness matrix according to node numbers. The global stiffness matrix reflects the mechanical coupling relationship between all elements in the 3D voxel mesh, and its size is directly related to the total number of nodes in the voxel mesh. For large-scale 3D voxel meshes, the global stiffness matrix has sparse, symmetric, and positive definite characteristics, and can be efficiently stored and computed using a compressed storage format.

[0134] Boundary conditions are applied in two ways: displacement constraint boundary conditions and load boundary conditions. Displacement constraint boundary conditions are applied to the bottom and lateral boundary nodes. The bottom boundary nodes are subject to triaxial displacement constraints, meaning they are fixed in the vertical and two horizontal directions. The lateral boundary nodes are subject to horizontal displacement constraints, allowing free vertical deformation to simulate the lateral constraint effect of a real foundation. The nodal forces at the top region nodes... Substituting the load boundary conditions into the right-hand side of the equation, we get the equivalent concentrated force transmitted from the roadbed surface to the top surface of the subgrade. After substituting the above displacement constraints and nodal force boundary conditions into the stress equilibrium governing equations, we discretize the equations to form a linear equation system with nodal displacement vectors as unknowns.

[0135] The numerical solution of the linear equation system employs the preconditioned conjugate gradient method, utilizing incomplete Cholesky decomposition as a preconditioner to improve the iterative convergence speed. After solving, the three-dimensional displacement field distribution of all voxel nodes is obtained, including lateral, longitudinal, and vertical displacement components. The vertical displacement components of the foundation top surface monitoring nodes are extracted from the displacement field distribution, representing the settlement values ​​at the corresponding locations. In practical engineering, foundation top surface monitoring nodes are typically located below the roadbed centerline and at characteristic points on both sides. By extracting the vertical displacements of these key nodes, settlement values ​​representing the overall settlement state of the foundation are obtained. Used for subsequent comparison with target settlement values Compare them.

[0136] Based on the solved displacement field distribution, the strain field distribution at each voxel node is further calculated. The strain field is obtained by spatial differentiation of the displacement field. For each voxel element, the strain components at each Gaussian integration point within the element are obtained by multiplying the spatial derivative matrix of the shape function (i.e., the strain-displacement matrix) with the element node displacement vector, including three normal strain components and three shear strain components. After mapping the strain components to the element nodes, the stress field distribution at each voxel node is calculated according to the generalized Hooke's law, based on the distribution of foundation material properties, resulting in six independent stress components. The stress field distribution is then post-processed to calculate the stress gradient amplitude at each node, i.e., the combined value of the rate of change of each component of the stress tensor in the spatial direction.

[0137] The stress gradient magnitude is compared node-by-node with a preset stress threshold, and nodes with stress gradients exceeding the threshold are marked as potential failure zones. Potential failure zones are typically concentrated at stress concentration points within the foundation, near interfaces between materials of different densities, and below areas of concentrated load. After identifying potential failure zones, the stress transfer path is traced along the principal stress direction within these zones to form candidate slip surfaces. The principal stress direction is determined by the eigenvalue decomposition of the stress tensor, with the direction of maximum principal stress representing the direction in which the material is most prone to shear failure. The tracing process starts from the boundary node of the potential failure zone and extends gradually along the direction of maximum principal stress until the path traverses the entire potential failure zone and extends to the foundation boundary, forming a complete slip surface trajectory. The above tracing process is repeated for all starting nodes within the potential failure zone that meet the path connectivity condition, generating multiple candidate slip surfaces.

[0138] For each candidate slip surface, calculate the resisting moment and the sliding moment. The resisting moment is determined by the normal stress at each node on the slip surface, the corresponding internal friction angle, and cohesion of the material. The total resisting moment for that slip surface is obtained by integrating and summing the shear strength contributions at each node along the slip surface, and then multiplying by the corresponding lever arm length. The sliding moment is obtained by integrating and summing the self-weight of the soil surrounding the slip surface and the component of the external load in the slip direction, multiplied by the lever arm length. The stability coefficient of a single slip surface is defined as the ratio of the resisting moment to the sliding moment. By iterating through all candidate slip surfaces and calculating the stability coefficient for each surface, the minimum value is extracted as the overall stability coefficient of the foundation. This minimum value corresponds to the most dangerous slip surface, representing the most unfavorable stability state of the foundation under the current density field configuration, and is used in conjunction with the target stability coefficient. The comparison drives the subsequent Bayesian optimization iteration process.

[0139] In practical engineering applications, for sections treated with foamed soil fill, the overall stress transfer path and failure mode of the foundation differ significantly from those of conventional roadbeds because the density of foamed soil is much lower than that of ordinary fill. The elastic modulus of the low-density foamed soil region... While exhibiting lower density and greater deformation under load, its lightweight properties effectively reduce additional stress on the foundation surface, thereby minimizing the potential damage zone. The numerical solution process described above accurately captures the comprehensive impact of the heterogeneous density field of the foamed soil filler on foundation settlement and stability, providing a reliable mechanical response evaluation basis for the Bayesian optimization algorithm and ensuring that the final optimal density field configuration simultaneously meets settlement control and stability requirements.

[0140] A second aspect of the present invention provides an AI-assisted three-dimensional design system for foamed soil filler roadbeds, comprising:

[0141] The parameter acquisition unit is used to acquire the roadbed geometric parameters, traffic load spectrum, and foundation soil layer distribution of the roadbed project.

[0142] The section division unit is used to extract features of the foundation soil layer distribution, generate a foundation bearing capacity map, calculate the roadbed surface stress distribution based on the traffic load spectrum, and divide the foamed soil filling treatment section and the conventional filling section according to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution.

[0143] The density field construction unit is used to construct a three-dimensional voxel grid for the foamed soil filler treatment section, in combination with the roadbed geometric parameters, assign density values ​​to each voxel node based on the foundation bearing capacity map, and generate the three-dimensional density field distribution of the foamed soil filler through a spatial interpolation algorithm.

[0144] The optimization solution unit is used to apply boundary constraints and loads to the three-dimensional density field distribution, perform numerical solutions to obtain settlement values ​​and stability coefficients, compare them with the target design values ​​to calculate the deviations, and iteratively adjust the three-dimensional density field distribution through a Bayesian optimization algorithm until the convergence condition is met to obtain the optimal density field configuration.

[0145] The construction output unit is used to generate a three-dimensional spatial distribution model of foamed soil filler according to the optimal density field configuration, and convert the optimal density field configuration into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and output the construction drawings.

[0146] A third aspect of the present invention provides an electronic device, comprising:

[0147] processor;

[0148] Memory used to store processor-executable instructions;

[0149] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0150] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0151] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0152] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. An AI-assisted three-dimensional design method for foamed soil filler subgrade, characterized in that, include: Obtain the roadbed geometric parameters, traffic load spectrum, and foundation soil layer distribution of the roadbed engineering; The features of the foundation soil layer distribution are extracted to generate a foundation bearing capacity map. The roadbed surface stress distribution is calculated based on the traffic load spectrum. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, the foamed soil filling treatment section and the conventional filling section are divided. For the foamed soil filler treatment section, a three-dimensional voxel grid is constructed in combination with the roadbed geometric parameters. Based on the foundation bearing capacity map, density values ​​are assigned to each voxel node, and a three-dimensional density field distribution of the foamed soil filler is generated through a spatial interpolation algorithm. Boundary constraints and loads are applied to the three-dimensional density field distribution, and the settlement value and stability coefficient are obtained by numerical solution. The deviation is calculated by comparing with the target design value. The three-dimensional density field distribution is iteratively adjusted by Bayesian optimization algorithm until the convergence condition is met, and the optimal density field configuration is obtained. A three-dimensional spatial distribution model of foamed soil filler is generated based on the optimal density field configuration, and the optimal density field configuration is converted into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and the construction drawings are output.

2. The method according to claim 1, characterized in that, Feature extraction is performed on the distribution of the foundation soil layers to generate a foundation bearing capacity map. The roadbed surface stress distribution is calculated based on the traffic load spectrum. According to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution, foamed soil filling treatment sections and conventional filling sections are divided, including: The distribution of the foundation soil layers is processed by layering and slicing, and the rock and soil type identification and thickness data of each soil layer are extracted to construct a vertical profile model containing soil layer sequence and burial depth information; Bearing characteristic parameters are assigned to each soil layer in the vertical profile model, the ultimate bearing capacity and allowable settlement of each soil layer are calculated, and the distribution of total bearing capacity and total settlement potential of the foundation are generated through vertical cumulative calculation, and integrated into a foundation bearing performance map. The dominant frequency load component is extracted by performing spectral analysis on the traffic load spectrum. The dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at the interfaces of each soil layer is calculated using the elastic layered system theory to obtain the stress distribution on the roadbed surface. The difference between the total bearing capacity distribution in the foundation bearing capacity performance map and the stress distribution on the roadbed surface is calculated to generate a bearing surplus distribution. Identify the areas with negative values ​​in the load-bearing surplus distribution, classify the areas with insufficient load-bearing capacity according to the amount of load loss, and determine the treatment intensity corresponding to each level. Areas with treatment intensity greater than the treatment threshold are divided into foamed soil filling treatment sections, and the remaining areas are divided into conventional filling sections.

3. The method according to claim 2, characterized in that, The dominant frequency load component is extracted by spectral analysis of the traffic load spectrum. This dominant frequency load component is then applied to the top surface of the roadbed, and the stress transmission and attenuation process at each soil layer interface is calculated using elastic layered system theory to obtain the stress distribution on the roadbed surface, including: The traffic load spectrum is segmented according to a preset time window and subjected to spectrum transformation to obtain a local spectrum. The frequency components that recur in all local spectra are counted, and the frequency component with the most recurrences is determined as the stable main frequency. The load amplitude corresponding to the stable main frequency is extracted as the main frequency load component. The dominant frequency load component is divided into distributed load nodes according to the transverse width and longitudinal length of the roadbed top surface, and the load intensity is determined according to the distance of each load node from the load action center. Extract the compression modulus and lateral compression coefficient of each soil layer in the vertical profile model, calculate the strain response coefficient of each soil layer under vertical load, and determine the damping characteristics of each soil layer. Select the load nodes on the top surface of the subgrade and use the corresponding load intensity as the initial stress input to the first soil layer. Calculate the stress attenuation within the layer and its transmission to the interface based on the damping characteristics of the first soil layer. Determine the degree of interface abrupt change by the difference in soil properties on both sides of the interface. When the degree of interface mutation exceeds the interface sensitivity threshold, an interface reflection loss coefficient is introduced. Combined with the stress at the interface, the stress continues to propagate to the next soil layer. The stress attenuation and interface treatment are repeated until the top surface of the foundation is reached. The stresses propagated from all load nodes to the top surface of the foundation are spatially superimposed to obtain the stress distribution of the roadbed surface.

4. The method according to claim 1, characterized in that, For the foamed soil filler section, a three-dimensional voxel grid is constructed based on the subgrade geometric parameters. Density values ​​are assigned to each voxel node based on the foundation bearing capacity map. A three-dimensional density field distribution of the foamed soil filler is generated using a spatial interpolation algorithm, including: Determine the start and end positions of the foamed soil filling treatment section in the longitudinal direction of the route. Combine the top width of the subgrade and the slope ratio in the subgrade geometric parameters to calculate the outline boundary of each cross section of the treatment section. Set multiple cross section slices at mileage intervals along the longitudinal direction. Generate planar grid points in each cross section slice according to the grid step size. Connect them in space to form a three-dimensional voxel grid. Based on the total bearing capacity distribution in the foundation bearing capacity spectrum, weak areas are identified with the bearing critical threshold as the boundary. The reinforcement density corresponding to the weak areas and the reference density corresponding to the non-weak areas are set respectively. Each voxel node in the three-dimensional voxel grid is traversed and it is determined whether it belongs to a weak area, and a differentiated density value is assigned. The density adjustment amount is propagated by the voxel nodes on the boundary of the weak area as the diffusion source. The density decay rate is determined according to the bearing capacity gradient between voxel nodes. A density value that gradually transitions from the reinforced density to the reference density is formed in the boundary area. The adjusted density values ​​of each voxel node are organized into a three-dimensional data array according to spatial coordinates to generate the three-dimensional density field distribution of the foamed soil filler.

5. The method according to claim 4, characterized in that, Density adjustment is propagated using voxel nodes at the boundary of the weak region as diffusion sources. The density decay rate is determined based on the load-bearing capacity gradient between voxel nodes, forming a density value that gradually transitions from reinforced density to reference density in the boundary region, including: Identify the spatial boundaries between weak and non-weak regions in a 3D voxel mesh, mark the voxel nodes in the weak regions located on the spatial boundaries as diffusion source nodes, and determine the connection paths extending from the diffusion source nodes to the non-weak regions on the spatial boundaries. The total bearing capacity value of each voxel node on the connection path is obtained from the foundation bearing capacity performance map. The cumulative change of total bearing capacity is calculated along the connection path. The cumulative change of total bearing capacity and the spatial length of the connection path are mapped together to the density decay rate of each voxel node. The density adjustment amount is propagated from the diffusion source node along the connection path. Each voxel node adjusts the amplitude of the incoming density adjustment amount according to the corresponding density decay rate and passes it to the downstream node, forming a density propagation sequence that decreases along the path. For voxel nodes located on multiple connection paths, the propagation priority of each connection path is determined based on the degree of matching between the cumulative change in the total carrying capacity of each connection path and the current total carrying capacity value of the node. The density adjustment amount transmitted by each connection path is then fused according to the propagation priority to form the gradually transitioning density value of each voxel node on the spatial boundary.

6. The method according to claim 1, characterized in that, Boundary constraints and loads are applied to the three-dimensional density field distribution, and settlement values ​​and stability coefficients are obtained through numerical solutions. These values ​​are then compared with the target design values ​​to calculate the deviations. The three-dimensional density field distribution is iteratively adjusted using a Bayesian optimization algorithm until the convergence condition is met, resulting in the optimal density field configuration, including: The displacement degrees of freedom of the boundary nodes of the three-dimensional density field distribution are fixed according to the foundation boundary conditions, and the load is converted into nodal forces acting on the nodes in the top region according to the engineering load distribution. The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution, and the settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results. The settlement value is compared with the target settlement value to obtain the settlement deviation, and the stability coefficient is compared with the target stability coefficient to obtain the stability deviation. A deviation function is then constructed. The deviation function is used as the objective function of the Bayesian optimization algorithm. The Bayesian optimization algorithm predicts the adjustment direction of the three-dimensional density field distribution based on the current deviation function value, and modifies the density value of each voxel node in the three-dimensional density field distribution along the adjustment direction to generate a new three-dimensional density field distribution. The new three-dimensional density field distribution is re-input into the mechanical response model for numerical solution and calculation of the new deviation function value. It is then determined whether the new deviation function value meets the convergence condition. If the convergence condition is met, the current three-dimensional density field distribution is determined as the optimal density field configuration. If the convergence condition is not met, the three-dimensional density field distribution is further adjusted based on the new deviation function value.

7. The method according to claim 6, characterized in that, The nodal forces are input into the mechanical response model constructed based on the three-dimensional density field distribution for numerical solution. The settlement value of the top surface of the foundation and the overall stability coefficient of the foundation are extracted from the solution results, including: The density values ​​of each voxel node are extracted from the three-dimensional density field distribution, and the corresponding elastic modulus and Poisson's ratio are calculated to construct the distribution of foundation material properties. The stress balance control equation is established by combining the spatial topological relationship of the three-dimensional voxel mesh. The displacement constraints of the boundary nodes and the nodal forces of the nodes in the top region are substituted into the stress balance control equation as boundary conditions. The stress balance control equation is discretized to form a linear equation system and numerically solved to obtain the displacement field distribution of each voxel node. The vertical displacement component of the monitoring node on the top surface of the foundation is extracted from the displacement field distribution as the settlement value. The strain field distribution of each voxel node is calculated based on the displacement field distribution, and the stress field distribution of each voxel node is calculated in combination with the distribution of the foundation material properties. Regions where the stress gradient exceeds the stress threshold are identified as potential failure regions. Within the potential failure regions, the stress transmission path is traced along the principal stress direction to form a slip surface. The anti-slip moment and sliding moment are calculated along the slip surface. All slip surfaces are traversed and the minimum stability coefficient is extracted as the overall stability coefficient of the foundation.

8. An AI-assisted three-dimensional design system for foamed soil filler subgrade, used to implement the method as described in any one of claims 1-7, characterized in that, include: The parameter acquisition unit is used to acquire the roadbed geometric parameters, traffic load spectrum, and foundation soil layer distribution of the roadbed project. The section division unit is used to extract features of the foundation soil layer distribution, generate a foundation bearing capacity map, calculate the roadbed surface stress distribution based on the traffic load spectrum, and divide the foamed soil filling treatment section and the conventional filling section according to the matching relationship between the foundation bearing capacity map and the roadbed surface stress distribution. The density field construction unit is used to construct a three-dimensional voxel grid for the foamed soil filler treatment section, in combination with the roadbed geometric parameters, assign density values ​​to each voxel node based on the foundation bearing capacity map, and generate the three-dimensional density field distribution of the foamed soil filler through a spatial interpolation algorithm. The optimization solution unit is used to apply boundary constraints and loads to the three-dimensional density field distribution, perform numerical solutions to obtain settlement values ​​and stability coefficients, compare them with the target design values ​​to calculate the deviations, and iteratively adjust the three-dimensional density field distribution through a Bayesian optimization algorithm until the convergence condition is met to obtain the optimal density field configuration. The construction output unit is used to generate a three-dimensional spatial distribution model of foamed soil filler according to the optimal density field configuration, and convert the optimal density field configuration into the filling range, density requirements and thickness parameters of each layer according to the construction layers, and output the construction drawings.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images