Dynamic parametric modeling method for tunnel operation area

Abstract

The invention belongs to the technical field of parametric modeling, and discloses a dynamic parametric modeling method for a tunnel operation area. Comprising the steps of collecting actual measurement data of a tunnel operation area; based on the collected actual measurement data, an initial model of the tunnel operation area is established; dividing the initial model of the tunnel operation area into a plurality of model areas, and constructing an operation area undirected graph by taking the model areas as nodes; colors are distributed to nodes in the operation area undirected graph; for each model area, determining a corresponding parameter group according to the allocated color, and updating based on the initial model of the tunnel operation area and the parameter group corresponding to each model area to obtain an overall parameterized model; the dynamic parametric modeling result of the whole tunnel operation area provides a scientific decision basis for design optimization, construction scheme formulation, risk assessment, operation management and the like of tunnel engineering.

Description

Technical Field

[0001] The present invention relates to the technical field of parametric modeling, and more specifically, to a method for dynamic parametric modeling of a tunnel operation area. Background Art

[0002] As an important part of infrastructure construction, tunnel engineering is of great significance to national economic development and social progress. However, due to the complexity and variability of the tunnel operation area environment, the existing tunnel operation area modeling technology has many shortcomings, which seriously restricts the design, construction and operation management level of tunnel engineering. First, there are significant spatial differences in the multi-source parameters of the tunnel operation area, and it is difficult to accurately reflect the actual situation by using a unified parameter modeling. Taking the tunnel engineering in Tibet as an example, due to the complex geological structure, large differences in rock and soil properties, and significant impact of the plateau environment in the region, if a unified parameter combination is used for modeling, it will inevitably lead to a large deviation between the model and the actual situation, which will bring major risks to engineering design and construction. Secondly, most of the existing tunnel operation area modeling methods are based on simplified assumptions and empirical formulas, lack of full use of measured data, low modeling accuracy and reliability, and designers can only rely on limited empirical data for modeling, resulting in a large deviation between the model and the actual situation. Furthermore, traditional modeling methods often regard the tunnel operation area as a whole, and cannot adopt different parameter combinations for different regions, making it difficult to reflect the characteristics of local areas. Finally, if the parameter transition is not handled properly, it may cause discontinuity or parameter mutation in some areas of the model, affecting the accuracy and reliability of the model.

[0003] In view of this, the present invention proposes a tunnel operation area dynamic parameterized modeling method to solve the above problems. Summary of the invention

[0004] In order to overcome the above-mentioned defects of the prior art and to achieve the above-mentioned purpose, the present invention provides the following technical solutions: A method for dynamic parameterized modeling of a tunnel operation area, comprising: step 1, collecting measured data of the tunnel operation area;

[0005] Step 2: Based on the collected measured data, an initial model of the tunnel operation area is established;

[0006] Step 3: Divide the initial model of the tunnel operation area into several model areas, and use the model areas as nodes to construct an undirected graph of the operation area;

[0007] Step 4, assigning colors to nodes in the undirected graph of the work area;

[0008] Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter groups corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

[0009] Furthermore, the measured data of the tunnel operation area includes geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data and construction parameter data.

[0010] Furthermore, the method of establishing the initial model of the tunnel operation area includes:

[0011] Classify and sort out the various types of measured data collected, archive them according to the data type, obtain corresponding data packets, perform smoothing interpolation processing on the corresponding data packets, and obtain corresponding smoothed data packets;

[0012] The three-dimensional coordinate data of each position point in the tunnel operation area is obtained, the obtained point cloud data is registered, and merged into a unified point cloud data set, and the point cloud data set is filtered to obtain a standard point cloud data set; the standard point cloud data set is segmented into different point cloud partial data using a segmentation algorithm based on region growing;

[0013] Fitting the corresponding surface model from each segmented point cloud data; meshing the fitted surface model to generate a mesh model composed of triangles or quadrilateral patches; optimizing the generated mesh model; obtaining a preliminary mesh model corresponding to each point cloud data; registering and connecting the preliminary mesh models corresponding to each point cloud data to obtain a preliminary tunnel operation area model;

[0014] The data points in the smoothed data packet are aligned to the space where the preliminary tunnel operation area model is located according to the corresponding coordinate values, that is, the preliminary tunnel operation area model is given initial parameters to obtain the initial model.

[0015] Furthermore, the method of performing smooth interpolation processing on the corresponding data packets includes:

[0016] The coordinates of the starting and ending points of the tunnel line are used as the preliminary reference points of the spatial range of the data points in the data packet, the spatial range of the data packet is determined based on the preliminary reference points, the number of division levels of the spatial range is defined, the entire spatial range is used as the first level, and the first level is divided into a regular three-dimensional grid;

[0017] Define the order of the mapping curve M = ceil(log2(max(n_x,n_y,n_z,α×(-∑(p(i)×log2(p(i))))))); where n_x, n_y and n_z are the number of three-dimensional grids in the x, y and z directions of the three-dimensional space respectively; p(i) is the probability of the i-th attribute value appearing in the data packet; α is the weight coefficient, and ceil() means rounding up; According to the defined order M, calculate the coding range [0, 2 (3×M) -1], that is, the minimum encoding value of the encoding range is 0 and the maximum encoding value is 2 (3×M) -1;

[0018] A virtual mapping curve is constructed based on the order M, and the data points in the data packet are mapped to the mapping curve using the Butz algorithm as a curve point. The value of the curve point corresponding to each data point is obtained, that is, the encoding value on the mapping curve;

[0019] The code value corresponding to the data point is divided into several intervals according to the code range. Each interval corresponds to a data block. The size of the interval Where S is the number of data blocks; according to the calculated interval size, determine the lower and upper bounds of each interval; traverse all data points in the data packet and assign them to corresponding intervals according to their encoding values; each interval corresponds to a first-level data block;

[0020] Obtain the number of data points in the data block. If it is greater than the preset data volume threshold, the corresponding data block is used as a new level and the data is re-divided in the same way as the first level. The division is stopped until the number of data points in the data block is less than or equal to the preset data volume threshold or reaches the preset number of division levels. The data blocks are strategically merged to obtain the final data block.

[0021] The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations are performed in parallel for the distributed data blocks to obtain interpolation results of each computing node; the interpolation results on each computing node are merged to obtain corresponding smoothed data packets.

[0022] Furthermore, the method of strategically merging the data blocks includes:

[0023] For each data block, calculate the centroid coordinates of all data points inside it; the centroid coordinates represent the spatial position center of the data block;

[0024] Construct a spatial index from the centroid coordinates of all data blocks; set an adaptive spatial distance threshold KL = α1×mean+α2×std, where α1 and α2 are two distance adjustment coefficients; mean is the average Euclidean distance between the centroid coordinates of all data points in all data blocks; std is the standard deviation of the Euclidean distance between the centroid coordinates of all data points in all data blocks;

[0025] Preset the search radius, use the spatial index to find the data blocks within the search radius for each data block, and record them as adjacent data blocks; for each pair of adjacent data blocks, calculate the Euclidean distance between their centroid coordinates; for two adjacent data blocks whose Euclidean distance is less than the spatial distance threshold, merge them into a new data block. The merged new data block contains all the data points of the original two data blocks, and recalculate the centroid coordinates of the new data block. Repeat the merging until there are no data block pairs that can be merged; after each merge, update the spatial index and data block.

[0026] Furthermore, the fitting method of the surface model includes:

[0027] Part of the point cloud data is divided into a three-dimensional voxel space, and the voxel space is subjected to median filtering; the filtered voxel space is obtained; the surface voxels are extracted from the filtered voxel space; the extracted surface voxels are subjected to surface fitting to generate a continuous surface model.

[0028] Furthermore, the method of extracting surface voxels from the filtered voxel space includes:

[0029] Define a sliding window, perform sliding window scanning on the filtered voxel space, and calculate the curvature characteristics and density characteristics of the voxels in each sliding window; the curvature characteristics include the principal curvature and Gaussian curvature; the density characteristics include the voxel occupancy rate and the number of voxel neighbors; normalize the curvature characteristics and density characteristics, and perform weighted calculation to obtain the comprehensive eigenvalue; according to the distribution of the comprehensive eigenvalue, divide the filtered voxel space into several voxel regions; for each voxel region, define a seed point number Nse; for each voxel region, randomly select Nse seed points;

[0030] Nse=a1×(Nv) β1 ×exp(γ×σc)×(1+γ1×σd)×(1+δ×σn); where,

[0031] a1, γ, γ1, δ and β1 are all adjustment coefficients; Nv is the number of voxels in the corresponding voxel area, σc is the standard deviation of the comprehensive characteristic value in the corresponding voxel area; σd is the standard deviation of the density value of the voxels in the corresponding voxel area; σn is the standard deviation of the normal vector of the voxels in the corresponding voxel area;

[0032] Take Nse randomly selected seed points as the preliminary area, preset the curvature threshold, and start from the randomly selected seed point to calculate the curvature values ​​of the voxels in its 26 neighborhoods. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current preliminary area and marked as visited. The voxel marked as visited is taken as the new seed point, and the curvature values ​​of the voxels in its 26 neighborhoods are repeatedly calculated. The voxel to be added to the preliminary area next is determined based on the curvature value; until no voxels are added to the preliminary area; the final preliminary area is obtained; the final preliminary area is marked as a surface area; all voxels marked as surface areas are gathered together, which is the extracted surface voxel.

[0033] Furthermore, the method of dividing the initial model of the tunnel operation area into a plurality of model areas includes:

[0034] Analyze the initial model of the tunnel operation area and identify different structural functional components; the structural functional components include the tunnel body, initial support, permanent lining, pipelines and ventilation facilities; initial support, permanent lining, pipelines and ventilation facilities are each considered as a model area;

[0035] Extract the corresponding part of the tunnel body from the initial model and record it as the tunnel grid model. Use the tunnel grid model as the wave source, and evenly sample a set of points on the surface of the tunnel grid model as the wave source point set. Place an initial wave front at each wave source point, define the number of simulation time steps, and update the position and shape of the wave front according to the propagation equation at each simulation time step.

[0036] At each time step, the positions and intersection times of the intersection points between different wavefronts are detected; all the intersection points are arranged in the order of the intersection time; that is, the point set on the center line is obtained; the point set of the center line is smoothed; the point set of the smoothed center line is parameterized to obtain the center line of the tunnel grid model;

[0037] Perform spatial sampling along the extracted center line according to a preset sampling interval; obtain a number of sampling points; each sampling point corresponds to a cross-section position of the tunnel grid model; at each sampling point, extract the cross section of the tunnel grid model, recorded as the cross section; extract the shape descriptor of each cross section; the shape descriptor includes perimeter, area, geometric moment and central moment; perform weighted summation of the perimeter, area, geometric moment and central moment to obtain the shape feature value corresponding to the cross section; define the buried depth value of each cross section as the distance from the center of the cross section to the ground surface;

[0038] Initially, each section is regarded as a model area; a shape threshold and a burial depth threshold are preset. If the absolute value of the difference between the shape feature values ​​of any two sections is less than the shape threshold, and the absolute value of the difference between the burial depth values ​​of the two sections is less than the burial depth threshold, the corresponding two sections are classified as the same model area.

[0039] Furthermore, the method of constructing the undirected graph of the operation area includes:

[0040] Determine a unique identifier for each model region; and calculate the geometric center or centroid of each model region as the representative point of the region; calculate the average value of the Euclidean distance between all pairs of points on the boundaries of any two model regions; record it as the proximity measure between the two model regions;

[0041] A neighboring threshold is preset. If the neighboring measure of two model areas is less than the neighboring threshold, the two model areas are adjacent. Each model area is taken as a node, and the attributes of the node include a unique identifier and the coordinates of the representative point. An empty undirected graph G = (V, E) is constructed, where V is a set of nodes and E is a set of unweighted edges. If two model areas are adjacent, an unweighted edge is added to the undirected graph G to connect the nodes corresponding to the two model areas. The attributes of the unweighted edge include the neighboring measure between the two model areas. The contact area and the intersection volume between the two adjacent model areas are calculated as the weight of the unweighted edge. At this time, the undirected graph G is the undirected graph of the operation area.

[0042] Furthermore, the method of assigning colors to nodes in the undirected graph of the work area includes:

[0043] Calculate the degree of each node in the undirected graph of the work area, define the coloring rule that the node with the largest degree is colored first, and use the greedy strategy to gradually assign colors to each node in the undirected graph of the work area until all nodes are assigned colors; that is, obtain the initial solution for assigning colors to the nodes in the undirected graph of the work area; calculate the number of colorings of the initial solution; count the number of node conflicts in the initial solution, that is, the number of adjacent pairs of nodes assigned the same color;

[0044] Define a set of nested neighborhood structures, each neighborhood structure N_k describes the transition from the current solution to the new solution; k is the index of the neighborhood structure; randomly select a neighborhood structure from a set of nested neighborhood structures, and randomly generate a new solution x' in the selected neighborhood structure; the new solution x' is obtained by locally perturbing the current solution, and the local perturbation method includes exchanging node colors and recoloring some nodes; starting from the new solution x', find the local optimal solution x" in the selected neighborhood structure N_k; calculate the target value of the new solution x', that is, the weighted sum of the number of node conflicts and the number of colorings, and generate all neighborhood solution sets N(x') of the new solution x' according to the definition of the selected neighborhood structure; each neighborhood solution in the neighborhood solution set N(x') is obtained by locally perturbing x';

[0045] For each neighborhood solution x', calculate its target value and record the currently found optimal neighborhood solution and its target value; the optimal neighborhood solution is the one with the minimum weighted sum of the number of node conflicts and the number of colorings; compare the target values ​​of the new solution x' and the optimal neighborhood solution; if the target value of the optimal neighborhood solution is less than the target value of the new solution x', move x' to the optimal neighborhood solution x_best, that is, x'=x_best; otherwise, keep x' unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimal solution x';

[0046] If the target value of the local optimal solution x" is smaller than the target value of the initial solution, the initial solution is updated to the local optimal solution x"; otherwise, the initial solution is kept unchanged; if no solution better than the local optimal solution x" is found in the selected neighborhood structure N_k, the neighborhood structure N_k is replaced with any one of the nested neighborhood structures; a new solution is regenerated and iterated until the local optimal solution is not updated for several consecutive times; the iteration is stopped, and the current optimal solution and its corresponding node coloring scheme are output; according to the node coloring scheme, colors are assigned to the nodes in the undirected graph of the work area.

[0047] Technical effects and advantages of a tunnel operation area dynamic parameter modeling method of the present invention:

[0048] The present invention can greatly improve the accuracy, reliability and applicability of modeling, and provide strong technical support for the design, construction and operation management of tunnel engineering. It makes full use of the comprehensive measured data of the tunnel operation area, including data in various aspects such as geology, mechanics, environment and construction, and provides reliable basic data support for modeling. It cleverly divides the tunnel operation area into multiple areas, and adopts different parameter combinations for modeling in different areas, which can effectively reflect the spatial differences within the tunnel operation area and avoid the errors and deviations caused by the use of unified parameter modeling. At the same time, a smooth transition area is set between adjacent areas to ensure the continuity and consistency of the overall model. The established dynamic parameterized model can not only accurately describe the current state of the tunnel operation area, but more importantly, it can dynamically adjust parameters according to actual conditions, predict and simulate the response of the tunnel operation area under different working conditions, and provide a scientific decision-making basis for the design optimization, construction plan formulation, risk assessment and operation management of tunnel engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 A schematic diagram of a tunnel operation area dynamic parameter modeling method of the present invention;

[0050] Figure 2 It is a schematic diagram of a tunnel operation area dynamic parameterized modeling system of the present invention. DETAILED DESCRIPTION

[0051] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0052] Example 1

[0053] See also Figure 1 As shown, a tunnel operation area dynamic parameterized modeling method described in this embodiment includes:

[0054] Step 1: Collect measured data of the tunnel operation area;

[0055] Step 2: Based on the collected measured data, an initial model of the tunnel operation area is established;

[0056] Step 3: Divide the initial model of the tunnel operation area into several model areas, and use the model areas as nodes to construct an undirected graph of the operation area;

[0057] Step 4, assigning colors to nodes in the undirected graph of the work area;

[0058] Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter groups corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

[0059] The measured data of the tunnel operation area include geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data and construction parameter data;

[0060] Geological exploration data include geological data such as rock and soil layer, lithology, structural surface, permeability, etc. obtained through core drilling and geophysical testing. These data reflect the geological conditions of the tunnel; in-situ test data include rock mechanics parameters, groundwater pressure and other data obtained through in-situ tests such as rock mechanics tests and water pressure tests. These data reflect the actual mechanical properties of the surrounding rock; topographic data include terrain undulations and elevation along the tunnel line. These data reflect the spatial position and burial depth of the tunnel; stress measurement data include in-situ stress test data such as horizontal stress and vertical stress. These data reflect the stress environment of the tunnel; environmental parameter data include environmental parameter test data such as temperature and humidity. These data will affect the mechanical properties of the surrounding rock; construction parameter data include construction parameter data such as excavation method and support measures. These data are important boundary conditions for establishing the initial model.

[0061] The above measured data comprehensively reflect the geological, mechanical, environmental and construction conditions of the tunnel operation area, and are the basic data for establishing the initial model of the tunnel operation area.

[0062] Methods for establishing an initial model of the tunnel working area include:

[0063] The collected measured data of various types are classified and sorted, and archived according to data types (such as geological data, mechanical data, environmental data, etc.), and corresponding data packets are obtained respectively. The corresponding data packets are smoothed and interpolated to obtain the corresponding continuous data distribution; that is, the corresponding smoothed data packets.

[0064] Specifically, the coordinates of the starting and ending points of the tunnel line are used as the preliminary reference points of the spatial range of the data points in the data packet. The spatial range of the data packet is determined based on the preliminary reference points. The spatial range is the range composed of the minimum coordinate value and the maximum coordinate value of the three coordinate systems in the three-dimensional space. The number of division levels of the spatial range is defined (the appropriate number of levels is determined by trial and error), the entire spatial range is taken as the first level, and the first level is divided into a regular three-dimensional grid (the number of three-dimensional grids is determined according to the preset number of data blocks), and the order of the mapping curve is defined as M=ceil(log2(max(n_x,n_y ,n_z,α×(-∑(p(i)×log2(p(i))))))); wherein n_x, n_y and n_z are the number of three-dimensional grids in the x, y and z directions of the three-dimensional space respectively; p(i) is the probability of the ith attribute value appearing in the data packet, which is calculated by dividing the number of occurrences of the attribute value by the number of data points in the data packet; the attribute value is any one of the x, y and z coordinate values ​​of the data point; α is the weight coefficient, which is adjusted according to the actual situation, and ceil() means rounding up; according to the defined order M, the coding range [0, 2 (3×M) -1], that is, the minimum encoding value of the encoding range is 0 and the maximum encoding value is 2 (3×M) -1.

[0065] A virtual mapping curve is constructed based on the order M, and the data points in the data packet are mapped to the mapping curve using the Butz algorithm as a curve point. The value of the curve point corresponding to each data point is obtained, that is, the encoding value on the mapping curve.

[0066] The code value corresponding to the data point is divided into several intervals according to the code range. Each interval corresponds to a data block. The size of the interval Where S is the number of data blocks; according to the calculated interval size, determine the lower and upper bounds of each interval; specifically, the lower bound of the first interval is the minimum coding value, the upper bound is the minimum coding value + the size of the interval - 1, the lower bound of the subsequent interval is the upper bound of the previous interval + 1, and the upper bound is the lower bound + the size of the interval - 1.

[0067] Traverse all data points in the data packet and assign them to corresponding intervals according to their coding values; each interval corresponds to a first-level data block, which contains all data points falling within the interval; obtain the number of data points in the data block, and if it is greater than the preset data volume threshold, use the corresponding data block as a new level and re-perform the same data division as the first level; until the number of data points in the data block is less than or equal to the preset data volume threshold or reaches the preset number of division levels; stop dividing; and strategically merge the data blocks to obtain the final data block.

[0068] Specifically, for each data block, the centroid coordinates of all the data points inside it are calculated (approximately calculated using the mean of the coordinate values ​​of the data points); the centroid coordinates represent the spatial position center of the data block; the centroid coordinates of all data blocks are used to construct a spatial index, such as an R-tree, a quadtree, etc.; an adaptive spatial distance threshold KL=α1×mean+α2×std is set, where α1 and α2 are two distance adjustment coefficients for controlling the size of the threshold; mean is the average Euclidean distance between the centroid coordinates of all data points inside all data blocks; std is the standard deviation of the Euclidean distance between the centroid coordinates of all data points inside all data blocks.

[0069] Preset the search radius, use the spatial index to find the data blocks within the search radius for each data block, and record them as adjacent data blocks; for each pair of adjacent data blocks, calculate the Euclidean distance between their centroid coordinates; for two adjacent data blocks whose Euclidean distance is less than the spatial distance threshold, merge them into a new data block. The merged new data block contains all the data points of the original two data blocks, and recalculate the centroid coordinates of the new data block. Repeat the merging until there are no data block pairs that can be merged; after each merge, update the spatial index and data block.

[0070] The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations are performed in parallel for the distributed data blocks, such as inverse distance weighted interpolation, Kriging interpolation, or radial basis function interpolation. The interpolation results of each computing node are obtained. The interpolation results on each computing node are merged to complete the smooth interpolation processing. The corresponding smoothed data packets are obtained.

[0071] The three-dimensional coordinate data of each position point in the tunnel operation area is obtained, the obtained point cloud data are aligned and merged into a unified point cloud data set, the point cloud data set is filtered to obtain a standard point cloud data set; the standard point cloud data set is segmented into different point cloud partial data using a segmentation algorithm based on region growing.

[0072] The corresponding surface model is fitted from each segmented point cloud data; the fitted surface model is meshed to generate a mesh model composed of triangles or quadrilateral patches; the generated mesh model is optimized by smoothing, simplifying, etc. to improve the quality and computational efficiency of the model; a preliminary mesh model corresponding to each point cloud data is obtained; the preliminary mesh models corresponding to each point cloud data are aligned and connected to obtain a preliminary tunnel operation area model; the data points in the smoothed data packet are aligned to the space where the preliminary tunnel operation area model is located according to the corresponding coordinate values, that is, the preliminary tunnel operation area model is given initial parameters; and the initial model is obtained.

[0073] The surface model fitting methods include:

[0074] The point cloud data is divided into a three-dimensional voxel space, and each voxel records whether the spatial position is occupied by the point cloud data; the voxel space is subjected to median filtering to remove isolated voxels and noise voxels; and the filtered voxel space is obtained.

[0075] Extract surface voxels from the filtered voxel space; specifically, define a sliding window, perform sliding window scanning on the filtered voxel space, and calculate the curvature features and density features of the voxels in each sliding window; the curvature features include the principal curvature and Gaussian curvature; the density features include the voxel occupancy rate (how many voxels are occupied or filled by point cloud data in a given three-dimensional voxel grid) and the number of voxel neighbors (how many occupied neighboring voxels are around each occupied voxel in a given voxel grid); normalize the curvature features and density features , and perform weighted calculation to obtain comprehensive eigenvalues; according to the distribution of comprehensive eigenvalues, the filtered voxel space is divided into several voxel regions; specifically, the K-Means algorithm is used to divide the voxels into K clusters according to the comprehensive eigenvalues; the mean difference of the comprehensive eigenvalues ​​between each two clusters is calculated, and the inter-cluster distance threshold is set. If the mean difference of the two clusters is less than the inter-cluster distance threshold, they are merged into a new cluster, and the merging is repeated until there is no cluster pair that can be merged. The final cluster is the voxel region; each voxel region represents a different combination of curvature features and density features.

[0076] For each voxel region, define a number of seed points Nse; for each voxel region, randomly select Nse seed points;

[0077] Nse=a1×(Nv) β1 ×exp(γ×σc)×(1+γ1×σd)×(1+δ×σn); where,

[0078] a1, γ, γ1, δ and β1 are all adjustment coefficients; Nv is the number of voxels in the corresponding voxel area, σc is the standard deviation of the comprehensive eigenvalue in the corresponding voxel area; σd is the standard deviation of the density value of the voxels in the corresponding voxel area, reflecting the degree of density change; σn is the standard deviation of the normal vector of the voxels in the corresponding voxel area, reflecting the degree of change of the normal vector.

[0079] Take Nse randomly selected seed points as the preliminary area, preset the curvature threshold, and start from the randomly selected seed point to calculate the curvature values ​​of the voxels in its 26 neighborhoods. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current preliminary area and marked as visited. The voxel marked as visited is taken as a new seed point, and the curvature values ​​of the voxels in its 26 neighborhoods are repeatedly calculated. The voxel to be added to the preliminary area next is determined based on the curvature value; until no voxels are added to the preliminary area; the final preliminary area is obtained; the final preliminary area is marked as a surface area; all voxels marked as surface areas are collected together, namely the extracted surface voxels; surface fitting is performed on the extracted surface voxels to generate a continuous surface model; surface fitting methods include moving least squares method, Gaussian kernel regression, etc.

[0080] The initial model of the tunnel operation area can be divided into several model areas:

[0081] The initial model of the tunnel operation area is analyzed to identify different structural functional components; the structural functional components include the tunnel body, initial support, permanent lining, pipelines and ventilation facilities; the initial support, permanent lining, pipelines and ventilation facilities are each regarded as a model area.

[0082] The part corresponding to the tunnel body is extracted from the initial model and recorded as the tunnel grid model. The tunnel grid model is used as the wave source, and a group of points are uniformly sampled on the surface of the tunnel grid model as the wave source point set; an initial wavefront is placed at each wave source point; the initial wavefront can be a sphere, a plane or other simple shapes, and the number of simulation time steps is defined. At each simulation time step, the position and shape of the wavefront are updated according to the propagation equation; the propagation equation is a set equation or a fast subequation.

[0083] At each time step, the positions and intersection times of the intersection points between different wave fronts are detected; all the intersection points are arranged in the order of intersection time; that is, the point set on the center line is obtained; the point set of the center line is smoothed (such as moving average, B-Spli ne curve fitting, etc.); the point set of the smoothed center line is parameterized (such as B-Spli ne curve, piecewise linear curve) to obtain the center line of the tunnel grid model.

[0084] Along the extracted center line, spatial sampling is performed according to a preset sampling interval; a number of sampling points are obtained; each sampling point corresponds to a cross-sectional position of the tunnel grid model; at each sampling point, the cross section of the tunnel grid model is extracted (cross-sectional extraction is performed using methods such as plane cutting and ray projection), which is recorded as the cross section; the shape descriptor of each cross section is extracted; the shape descriptor includes perimeter, area, geometric moment and central moment; the perimeter, area, geometric moment and central moment are weightedly summed to obtain the shape feature value corresponding to the cross section; the burial depth value of each cross section is defined as the distance from the center of the cross section to the ground surface; the geometric moment refers to a set of values ​​obtained by multiplying the coordinate value of each point by the corresponding weight and then integrating it for an area, and the central moment is a set of moment values ​​calculated on the basis of calculating the geometric moment by first translating the coordinate origin of the shape to the centroid position of the shape.

[0085] Initially, each section is regarded as a model area; a shape threshold and a burial depth threshold are preset. If the absolute value of the difference between the shape feature values ​​of any two sections is less than the shape threshold, and the absolute value of the difference between the burial depth values ​​of the two sections is less than the burial depth threshold, the corresponding two sections are classified as the same model area.

[0086] Determine a unique identifier (ID) for each model region; and calculate the geometric center or centroid of each model region as the representative point of the region; calculate the average value of the Euclidean distance between all pairs of points on the boundary of any two model regions; record it as the proximity measure between the two model regions; preset a proximity threshold, if the proximity measure of the two model regions is less than the proximity threshold, the two model regions are adjacent.

[0087] Each model area is regarded as a node, and the attributes of the node include a unique identifier and the coordinates of the representative point; an empty undirected graph G = (V, E) is constructed, where V is a set of nodes and E is a set of unweighted edges; if two model areas are adjacent, an unweighted edge is added to the undirected graph G to connect the nodes corresponding to the two model areas; the attributes of the unweighted edge include the adjacent measure between the two model areas; the contact area and the intersection volume between the two adjacent model areas are calculated as the weight of the unweighted edge; the undirected graph G at this time is the work area undirected graph.

[0088] Ways to assign colors to nodes in the workspace undirected graph include:

[0089] Calculate the degree of each node in the undirected graph of the work area, define the coloring rule that the node with the largest degree is colored first, and use the greedy strategy to gradually assign colors to each node in the undirected graph of the work area until all nodes are assigned colors; that is, obtain the initial solution (node ​​coloring scheme) for assigning colors to the nodes in the undirected graph of the work area; calculate the number of colorings of the initial solution, that is, the number of color types used; count the number of node conflicts in the initial solution, that is, the number of adjacent pairs of nodes assigned the same color.

[0090] Define a set of nested neighborhood structures, each neighborhood structure N_k describes the transition from the current solution to the new solution; k is the index of the neighborhood structure; for example, N_1 is to exchange the colors of two adjacent nodes; N_2 is to recolor a node and all its adjacent nodes; N_3 is to recolor all nodes in a maximum group; N_k_max is to regenerate a completely new coloring scheme.

[0091] A neighborhood structure is randomly selected from a set of nested neighborhood structures, and a new solution x' is randomly generated in the selected neighborhood structure; the new solution x' is obtained by locally perturbing the current solution, and the local perturbation method includes exchanging node colors and recoloring some nodes.

[0092] Starting from the new solution x', find the local optimal solution x" in the selected neighborhood structure N_k; specifically, calculate the target value of the new solution x', that is, the weighted sum of the number of node conflicts and the number of colorings, and generate all neighborhood solution sets N(x') of the new solution x' according to the definition of the selected neighborhood structure; each neighborhood solution in the neighborhood solution set N(x') is obtained by locally perturbing x'.

[0093] For each neighborhood solution x', calculate its target value (the weighted sum of the number of node conflicts and the number of colorings), and record the best neighborhood solution and its target value found so far; the best neighborhood solution is the one with the minimum weighted sum of the number of node conflicts and the number of colorings.

[0094] Compare the target values ​​of the new solution x' and the optimal neighborhood solution; if the target value of the optimal neighborhood solution is less than the target value of the new solution x', move x' to the optimal neighborhood solution x_best, that is, x'=x_best; otherwise, keep x' unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimal solution x".

[0095] If the target value of the local optimal solution x" is smaller than the target value of the initial solution, the initial solution is updated to the local optimal solution x"; otherwise, the initial solution is kept unchanged; if no solution better than the local optimal solution x" (with a smaller target value) is found in the selected neighborhood structure N_k, the neighborhood structure N_k is replaced with any of the nested neighborhood structures; a new solution is regenerated and iterated until the local optimal solution is not updated for several consecutive times; the iteration is stopped and the current optimal solution (with the smallest target value) and its corresponding node coloring scheme are output.

[0096] According to the node coloring scheme, colors are assigned to the nodes in the undirected graph of the work area. In the above process, the neighborhood structure determines the search scope and accuracy, which usually needs to be gradually expanded from a small range to a large range. By systematically exploring different neighborhoods, it is possible to effectively escape from the local optimum and obtain a better node coloring and parameter allocation scheme.

[0097] Based on the color assigned to each model area, several groups of different parameter sets corresponding to different colors are prepared in advance, and each group of parameter sets contains several specific parameters (adjusted based on the measured data of the tunnel operation area) to comprehensively describe the geological, mechanical, environmental and construction conditions of the corresponding area of ​​the group; traverse all model areas, and for the j-th model area, obtain its assigned color c_j; among several groups of different parameter sets, find the parameter set corresponding to color c_j, and determine this parameter set as the parameter set corresponding to the j-th model area.

[0098] An overlapping area is set at the junction of two adjacent model areas. In the overlapping area, the parameter groups of the two adjacent model areas are interpolated to achieve a smooth transition. The interpolation method can be linear interpolation, nonlinear interpolation (such as Gaussian process regression), etc.

[0099] Through the above steps, each model area in the entire tunnel operation area obtains the corresponding parametric model, and realizes the smooth transition of parameters at the junction of the areas, thus completing the dynamic parametric modeling of the entire operation area.

[0100] The initial model is split according to the division of the model area to obtain the initial sub-model corresponding to each model area. For the initial sub-model of each model area, the relevant parameters in the initial sub-model are updated according to the parameter group corresponding to the model area to obtain a parameterized sub-model; for example, the formation model parameters are updated, such as rock and soil strata, lithology, structural surface, etc.; the mechanical parameters are updated, such as rock strength, deformation modulus, etc.; the environmental parameters are updated, such as temperature, humidity, etc.; the construction parameters are updated, such as excavation method, support measures, etc.

[0101] The parameterized sub-model is spliced ​​with the initial model to obtain an overall parameterized model that includes all model areas. The parameterized models of adjacent areas are smoothly interpolated using the overlapping area method to achieve seamless transition of parameters. It should be noted that in the process of updating parameters, the correlation and constraints between parameters need to be fully considered to ensure the rationality and consistency of the parameters. At the same time, the parameters also need to be calibrated and optimized according to actual conditions to improve the accuracy and reliability of modeling.

[0102] This embodiment can greatly improve the accuracy, reliability and applicability of modeling, provide strong technical support for the design, construction and operation management of tunnel engineering, make full use of the comprehensive measured data of the tunnel operation area, including data in various aspects such as geology, mechanics, environment and construction, provide reliable basic data support for modeling, cleverly divide the tunnel operation area into multiple areas, and use different parameter combinations for modeling in different areas, which can effectively reflect the spatial differences within the tunnel operation area and avoid the errors and deviations caused by the use of unified parameter modeling; at the same time, a smooth transition area is set between adjacent areas to ensure the continuity and consistency of the overall model. The established dynamic parameterized model can not only accurately describe the current state of the tunnel operation area, but more importantly, it can dynamically adjust parameters according to actual conditions, predict and simulate the response of the tunnel operation area under different working conditions, and provide a scientific decision-making basis for the design optimization, construction plan formulation, risk assessment and operation management of tunnel engineering.

[0103] Example 2

[0104] See also Figure 2 As shown, the part not described in detail in this embodiment is described in Example 1, which provides a tunnel operation area dynamic parameter modeling system, including:

[0105] Data acquisition module, used to collect measured data in the tunnel operation area;

[0106] The preliminary modeling module builds an initial model of the tunnel operation area based on the collected measured data;

[0107] The area division and composition module is used to divide the initial model of the tunnel operation area into several model areas and construct an undirected graph of the operation area using the model areas as nodes;

[0108] A coloring module, used to assign colors to nodes in the undirected graph of the work area;

[0109] The parameterization module is used to determine the corresponding parameter group for each model area according to its assigned color, and to update the overall parameterized model based on the initial model of the tunnel operation area and the parameter groups corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area; each module is connected by wired and / or wireless means to realize data transmission between modules.

[0110] Example 3

[0111] This embodiment discloses an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the operation mode of the above-mentioned method for dynamic parameterized modeling of a tunnel operation area is implemented.

[0112] Since the electronic device introduced in this embodiment is an electronic device used to implement a method for dynamic parameterized modeling of a tunnel operation area in the embodiment of the present application, based on the method for dynamic parameterized modeling of a tunnel operation area introduced in the embodiment of the present application, a person skilled in the art can understand the specific implementation of the electronic device of the present embodiment and its various variations, so how the electronic device implements the method in the embodiment of the present application is not described in detail here. As long as a person skilled in the art implements the electronic device used in the method for dynamic parameterized modeling of a tunnel operation area in the embodiment of the present application, it belongs to the scope of protection of the present application.

[0113] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters and thresholds in the formula are set by technicians in this field according to actual conditions.

[0114] The above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technical users in this technical field, some improvements and modifications without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A method for dynamic parameterized modeling of a tunnel operation area, characterized in that: include: Step 1: Collect measured data of the tunnel operation area; Step 2: Based on the collected measured data, an initial model of the tunnel operation area is established; Step 3: Divide the initial model of the tunnel operation area into several model areas, and use the model areas as nodes to construct an undirected graph of the operation area; Step 4, assigning colors to nodes in the undirected graph of the work area; Step 5: For each model area, determine the corresponding parameter group according to its assigned color, and update the overall parameterized model based on the initial model of the tunnel operation area and the parameter groups corresponding to each model area; that is, the dynamic parameterized modeling result of the entire tunnel operation area.

2. A tunnel operation area dynamic parameter modeling method according to claim 1, characterized in that: The measured data of the tunnel operation area includes geological exploration data, in-situ test data, topographic data, stress measurement data, environmental parameter data and construction parameter data.

3. A tunnel operation area dynamic parameter modeling method according to claim 2, characterized in that: The method of establishing the initial model of the tunnel operation area includes: Classify and sort out the various types of measured data collected, archive them according to the data type, obtain corresponding data packets, perform smoothing interpolation processing on the corresponding data packets, and obtain corresponding smoothed data packets; The three-dimensional coordinate data of each position point in the tunnel operation area is obtained, the obtained point cloud data is registered, and merged into a unified point cloud data set, and the point cloud data set is filtered to obtain a standard point cloud data set; the standard point cloud data set is segmented into different point cloud partial data using a segmentation algorithm based on region growing; Fitting the corresponding surface model from each segmented point cloud data; meshing the fitted surface model to generate a mesh model composed of triangles or quadrilateral patches; optimizing the generated mesh model; obtaining a preliminary mesh model corresponding to each point cloud data; registering and connecting the preliminary mesh models corresponding to each point cloud data to obtain a preliminary tunnel operation area model; The data points in the smoothed data packet are aligned to the space where the preliminary tunnel operation area model is located according to the corresponding coordinate values, that is, the preliminary tunnel operation area model is given initial parameters to obtain the initial model.

4. A tunnel operation area dynamic parameter modeling method according to claim 3, characterized in that: The method of performing smooth interpolation processing on the corresponding data packets includes: The coordinates of the starting and ending points of the tunnel line are used as the preliminary reference points of the spatial range of the data points in the data packet, the spatial range of the data packet is determined based on the preliminary reference points, the number of division levels of the spatial range is defined, the entire spatial range is used as the first level, and the first level is divided into a regular three-dimensional grid; Define the order of the mapping curve M = ceil(log2(max(n_x,n_y,n_z,α×(-∑(p(i)×log2(p(i))))))); where n_x, n_y and n_z are the number of three-dimensional grids in the x, y and z directions of the three-dimensional space respectively; p(i) is the probability of the i-th attribute value appearing in the data packet; α is the weight coefficient, and ceil() means rounding up; According to the defined order M, calculate the coding range [0, 2 (3×M) -1], that is, the minimum encoding value of the encoding range is 0 and the maximum encoding value is 2 (3×M) -1; A virtual mapping curve is constructed based on the order M, and the data points in the data packet are mapped to the mapping curve using the Butz algorithm as a curve point. The value of the curve point corresponding to each data point is obtained, that is, the encoding value on the mapping curve; The code value corresponding to the data point is divided into several intervals according to the code range. Each interval corresponds to a data block. The size of the interval Where S is the number of data blocks; according to the calculated interval size, determine the lower and upper bounds of each interval; traverse all data points in the data packet and assign them to corresponding intervals according to their encoding values; each interval corresponds to a first-level data block; Obtain the number of data points in the data block. If it is greater than the preset data volume threshold, the corresponding data block is used as a new level and the data is re-divided in the same way as the first level. The division is stopped until the number of data points in the data block is less than or equal to the preset data volume threshold or reaches the preset number of division levels. The data blocks are strategically merged to obtain the final data block. The final data blocks are distributed to different computing nodes. On each computing node, interpolation calculations are performed in parallel for the distributed data blocks to obtain interpolation results of each computing node; the interpolation results on each computing node are merged to obtain corresponding smoothed data packets.

5. A method for dynamic parameterized modeling of a tunnel operation area according to claim 4, characterized in that: The method of strategically merging data blocks includes: For each data block, calculate the centroid coordinates of all data points inside it; the centroid coordinates represent the spatial position center of the data block; Construct a spatial index from the centroid coordinates of all data blocks; set an adaptive spatial distance threshold KL = α1×mean+α2×std, where α1 and α2 are two distance adjustment coefficients; mean is the average Euclidean distance between the centroid coordinates of all data points in all data blocks; std is the standard deviation of the Euclidean distance between the centroid coordinates of all data points in all data blocks; Preset the search radius, use the spatial index to find the data blocks within the search radius for each data block, and record them as adjacent data blocks; for each pair of adjacent data blocks, calculate the Euclidean distance between their centroid coordinates; for two adjacent data blocks whose Euclidean distance is less than the spatial distance threshold, merge them into a new data block. The merged new data block contains all the data points of the original two data blocks, and recalculate the centroid coordinates of the new data block. Repeat the merging until there are no data block pairs that can be merged; after each merge, update the spatial index and data block.

6. A method for dynamic parameterized modeling of a tunnel operation area according to claim 5, characterized in that: The fitting method of the surface model includes: Part of the point cloud data is divided into a three-dimensional voxel space, and the voxel space is subjected to median filtering; the filtered voxel space is obtained; the surface voxels are extracted from the filtered voxel space; the extracted surface voxels are subjected to surface fitting to generate a continuous surface model.

7. A method for dynamic parameterized modeling of a tunnel operation area according to claim 6, characterized in that: The method of extracting surface voxels from the filtered voxel space includes: Define a sliding window, perform sliding window scanning on the filtered voxel space, and calculate the curvature characteristics and density characteristics of the voxels in each sliding window; the curvature characteristics include the principal curvature and Gaussian curvature; the density characteristics include the voxel occupancy rate and the number of voxel neighbors; normalize the curvature characteristics and density characteristics, and perform weighted calculation to obtain the comprehensive eigenvalue; according to the distribution of the comprehensive eigenvalue, divide the filtered voxel space into several voxel regions; for each voxel region, define a seed point number Nse; for each voxel region, randomly select Nse seed points; Nse=a1×(Nv) β1 ×exp(γ×σc)×(1+γ1×σd)×(1+δ×σn); where, a1, γ, γ1, δ and β1 are all adjustment coefficients; Nv is the number of voxels in the corresponding voxel area, σc is the standard deviation of the comprehensive characteristic value in the corresponding voxel area; σd is the standard deviation of the density value of the voxels in the corresponding voxel area; σn is the standard deviation of the normal vector of the voxels in the corresponding voxel area; Take Nse randomly selected seed points as the preliminary area, preset the curvature threshold, and start from the randomly selected seed point to calculate the curvature values ​​of the voxels in its 26 neighborhoods. If the absolute value of the curvature value is greater than the curvature threshold, the corresponding voxel is added to the current preliminary area and marked as visited. The voxel marked as visited is taken as the new seed point, and the curvature values ​​of the voxels in its 26 neighborhoods are repeatedly calculated. The voxel to be added to the preliminary area next is determined based on the curvature value; until no voxels are added to the preliminary area; the final preliminary area is obtained; the final preliminary area is marked as a surface area; all voxels marked as surface areas are gathered together, which is the extracted surface voxel.

8. A method for dynamic parameterized modeling of a tunnel operation area according to claim 7, characterized in that: The method of dividing the initial model of the tunnel operation area into a plurality of model areas includes: Analyze the initial model of the tunnel operation area and identify different structural functional components; the structural functional components include the tunnel body, initial support, permanent lining, pipelines and ventilation facilities; initial support, permanent lining, pipelines and ventilation facilities are each considered as a model area; Extract the corresponding part of the tunnel body from the initial model and record it as the tunnel grid model. Use the tunnel grid model as the wave source, and evenly sample a set of points on the surface of the tunnel grid model as the wave source point set. Place an initial wave front at each wave source point, define the number of simulation time steps, and update the position and shape of the wave front according to the propagation equation at each simulation time step. At each time step, the positions and intersection times of the intersection points between different wavefronts are detected; all the intersection points are arranged in the order of the intersection time; that is, the point set on the center line is obtained; the point set of the center line is smoothed; the point set of the smoothed center line is parameterized to obtain the center line of the tunnel grid model; Perform spatial sampling along the extracted center line according to a preset sampling interval; obtain a number of sampling points; each sampling point corresponds to a cross-section position of the tunnel grid model; at each sampling point, extract the cross section of the tunnel grid model, recorded as the cross section; extract the shape descriptor of each cross section; the shape descriptor includes perimeter, area, geometric moment and central moment; perform weighted summation of the perimeter, area, geometric moment and central moment to obtain the shape feature value corresponding to the cross section; define the buried depth value of each cross section as the distance from the center of the cross section to the ground surface; Initially, each section is regarded as a model area; a shape threshold and a burial depth threshold are preset. If the absolute value of the difference between the shape feature values ​​of any two sections is less than the shape threshold, and the absolute value of the difference between the burial depth values ​​of the two sections is less than the burial depth threshold, the corresponding two sections are classified as the same model area.

9. A method for dynamic parameterized modeling of a tunnel operation area according to claim 8, characterized in that: Ways to construct an undirected graph of the work area include: Determine a unique identifier for each model region; and calculate the geometric center or centroid of each model region as the representative point of the region; calculate the average value of the Euclidean distance between all pairs of points on the boundaries of any two model regions; record it as the proximity measure between the two model regions; A neighboring threshold is preset. If the neighboring measure of two model areas is less than the neighboring threshold, the two model areas are adjacent. Each model area is taken as a node, and the attributes of the node include a unique identifier and the coordinates of a representative point. An empty undirected graph G = (V, E) is constructed, where V is a set of nodes and E is a set of unweighted edges. If two model areas are adjacent, an unweighted edge is added to the undirected graph G to connect the nodes corresponding to the two model areas. The attributes of the unweighted edge include the neighboring measure between the two model areas. The contact area and intersection volume between the two adjacent model areas are calculated as the weight of the unweighted edge. The undirected graph G at this time is the undirected graph of the working area.

10. A tunnel operation area dynamic parameter modeling method according to claim 9, characterized in that: The method of assigning colors to nodes in the undirected graph of the work area includes: Calculate the degree of each node in the undirected graph of the work area, define the coloring rule that the node with the largest degree is colored first, and use the greedy strategy to gradually assign colors to each node in the undirected graph of the work area until all nodes are assigned colors; that is, obtain the initial solution for assigning colors to the nodes in the undirected graph of the work area; calculate the number of colorings of the initial solution; count the number of node conflicts in the initial solution, that is, the number of adjacent pairs of nodes assigned the same color; Define a set of nested neighborhood structures, each neighborhood structure N_k describes the transition from the current solution to the new solution; k is the index of the neighborhood structure; randomly select a neighborhood structure from a set of nested neighborhood structures, and randomly generate a new solution x' in the selected neighborhood structure; the new solution x' is obtained by locally perturbing the current solution, and the local perturbation method includes exchanging node colors and recoloring some nodes; starting from the new solution x', find the local optimal solution x" in the selected neighborhood structure N_k; calculate the target value of the new solution x', that is, the weighted sum of the number of node conflicts and the number of colorings, and generate all neighborhood solution sets N(x') of the new solution x' according to the definition of the selected neighborhood structure; each neighborhood solution in the neighborhood solution set N(x') is obtained by locally perturbing x'; For each neighborhood solution x', calculate its target value and record the currently found optimal neighborhood solution and its target value; the optimal neighborhood solution is the one with the minimum weighted sum of the number of node conflicts and the number of colorings; compare the target values ​​of the new solution x' and the optimal neighborhood solution; if the target value of the optimal neighborhood solution is less than the target value of the new solution x', move x' to the optimal neighborhood solution x_best, that is, x'=x_best; otherwise, keep x' unchanged; until the preset maximum number of iterations is reached, output the current optimal neighborhood solution as the local optimal solution x'; If the target value of the local optimal solution x" is smaller than the target value of the initial solution, the initial solution is updated to the local optimal solution x"; otherwise, the initial solution is kept unchanged; if no solution better than the local optimal solution x" is found in the selected neighborhood structure N_k, the neighborhood structure N_k is replaced with any one of the nested neighborhood structures; a new solution is regenerated and iterated until the local optimal solution is not updated for several consecutive times; the iteration is stopped, and the current optimal solution and its corresponding node coloring scheme are output; according to the node coloring scheme, colors are assigned to the nodes in the undirected graph of the work area.

Images