Aggregate grain shape and gradation identification method and device based on image processing and storage medium

Abstract

The application provides an aggregate grain shape and grading identification method and device based on image processing and a storage medium, a dynamic falling image stream of an aggregate particle group is acquired, a closed contour image of a single aggregate particle is separated from the image stream, a boundary constraint harmonic field is constructed to generate a harmonic field distribution, a stable critical structure is extracted based on the harmonic field distribution, a critical turning point of a field strength value change is identified, the critical turning point is coded into a topological invariant sequence, a shape semantic mapper calibrated by a concrete interface performance is called, semantic mapping processing is performed on the topological invariant sequence, a quantitative representation of a single aggregate particle to an aggregate-paste interface occlusion ability is output, spatial density weighted fusion is performed according to the quantitative representations of the whole batch of aggregate particles, after normalization processing, a structure performance index reflecting the functional characteristics of the aggregate grain shape is generated through weighted fusion. The application can comprehensively reflect the functional characteristics of the whole batch of aggregate, and provides a direct performance evaluation basis for the optimization of the aggregate grading in the concrete mix design.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a method, apparatus and storage medium for identifying aggregate particle shape and gradation based on image processing. Background Technology

[0002] In the field of concrete material performance research and engineering applications, aggregate particle shape and gradation identification are key aspects affecting the workability, strength, and durability of concrete. Accurately obtaining the shape characteristics and gradation distribution of aggregate particles is of great significance for optimizing concrete mix design. Traditional identification methods mostly rely on manual sieving or static image acquisition, extracting geometric parameters such as the perimeter, area, and roundness of particles to describe particle shape. However, these methods are difficult to cope with contour extraction errors caused by particle adhesion and motion blur during dynamic falling. Moreover, relying solely on geometric parameters cannot fully characterize the microscopic morphological features such as particle surface texture and angular distribution, resulting in a lack of direct correlation in the assessment of aggregate-paste interface interlocking ability. At the same time, traditional gradation analysis focuses more on particle size distribution, ignoring the impact of the relative positional relationship between particles on the overall structural performance, making it difficult to directly provide quantitative basis for particle shape functional characteristics for concrete mix optimization. Summary of the Invention

[0003] In view of this, the present invention provides a method, apparatus, and storage medium for identifying aggregate particle shape and gradation based on image processing. The technical solution of the embodiments of the present invention is implemented as follows:

[0004] On one hand, this invention provides an image processing-based method for identifying aggregate particle shape and gradation. The method includes: acquiring a dynamic falling image stream of aggregate particles; separating a closed contour image of a single aggregate particle from the dynamic falling image stream, the closed contour image containing a set of continuous pixel coordinates of the particle boundary; constructing a boundary-constrained harmonic field on the closed contour image, generating a harmonic field distribution that satisfies the contour boundary gradient condition in the particle's internal region, the harmonic field distribution representing the field strength value change at each point inside the particle through a two-dimensional spatial field strength value matrix; and extracting stable critical structures based on the harmonic field distribution to identify the field strength value. The critical inflection point of change is encoded into a topologically invariant sequence by encrypting the spatial coordinate sequence of the critical inflection point. A shape semantic mapper calibrated by concrete interface performance is invoked to perform semantic mapping processing on the topologically invariant sequence, outputting a quantitative characterization of the interlocking ability of a single aggregate particle at the aggregate-paste interface. The quantitative characterization includes the proportion of interface contact area and the distribution parameters of interlocking depth. Based on the quantitative characterization of the entire batch of aggregate particles, spatial density weighted fusion is performed, and a density weight factor is calculated based on the relative positional relationship between particles. The quantitative characterization is then normalized, and a structural performance index reflecting the functional characteristics of aggregate particle shape is generated through weighted fusion.

[0005] On the other hand, the present invention provides an aggregate particle shape and gradation identification device, comprising:

[0006] The image acquisition module is used to acquire a dynamic falling image stream of aggregate particles, and to separate a closed contour image of a single aggregate particle from the dynamic falling image stream. The closed contour image contains a set of continuous pixel coordinates of the particle boundary.

[0007] The harmonic field construction module is used to construct a boundary-constrained harmonic field for the closed contour image, and generate a harmonic field distribution that satisfies the contour boundary gradient condition in the internal region of the particle. The harmonic field distribution is characterized by the field intensity value change of each point inside the particle through a two-dimensional spatial field intensity value matrix.

[0008] The critical encoding module is used to extract stable critical structures based on the harmonic field distribution, identify critical inflection points of field strength value changes, and encode the spatial coordinate sequence of the critical inflection points into a topologically invariant sequence.

[0009] The semantic mapping module is used to call the shape semantic mapper calibrated by the concrete interface performance to perform semantic mapping processing on the topological invariant sequence and output a quantitative characterization of the interlocking ability of a single aggregate particle at the aggregate-paste interface. The quantitative characterization includes the interface contact area ratio and interlocking depth distribution parameters.

[0010] The weighted fusion module is used to perform spatial density weighted fusion based on the quantitative characterization of the entire batch of aggregate particles, calculate the density weight factor based on the relative positional relationship between particles, normalize the quantitative characterization, and generate a structural performance index that reflects the functional characteristics of aggregate particle shape through weighted fusion.

[0011] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described above.

[0012] This invention acquires the dynamic falling image stream of aggregate particles and separates the closed contour images of individual aggregate particles, effectively overcoming the contour extraction errors caused by particle adhesion and motion blur in traditional static image analysis. This provides a high-precision single-particle morphology basis for particle shape analysis. By constructing a boundary-constrained harmonic field on the closed contour images, the external contour features are transformed into internal field strength distribution patterns. Furthermore, the critical inflection points of field strength changes are extracted and encoded with topological invariance. This overcomes the limitations of traditional geometric parameters in describing particle shape, achieving a robust representation of the essential structural features of particle morphology, unaffected by translation, rotation, scaling, etc. The interference of geometric transformations is eliminated; by calling the shape semantic mapper calibrated by concrete interface performance, the correlation between particle topology and aggregate-paste interface interlocking ability is directly established, avoiding the indirectness and error accumulation of performance prediction by traditional empirical formulas; by performing spatial density-weighted fusion through quantitative characterization of the entire batch of aggregate particles, the density weight factor is calculated by introducing the relative positional relationship between particles, realizing the spatial weighted integration of single particle quantitative characterization, which can comprehensively reflect the particle shape functional characteristics of the entire batch of aggregates, and provides a direct performance evaluation basis for the optimization of aggregate gradation in concrete mix design, rather than simply size distribution data. Attached Figure Description

[0013] Figure 1 This is a schematic diagram illustrating the implementation process of an image processing-based method for identifying aggregate particle shape and gradation, provided in an embodiment of the present invention.

[0014] Figure 2 This is a schematic diagram of the composition structure of an aggregate particle shape and gradation identification device provided in an embodiment of the present invention.

[0015] Figure 3 This is a schematic diagram of the hardware entity of a computer system provided in an embodiment of the present invention. Detailed Implementation

[0016] This invention provides an image processing-based method for identifying aggregate particle shape and gradation, which can be executed by a processor in a computer system. The computer system can refer to a device with data processing capabilities, such as a laptop, tablet, desktop computer, or distributed system.

[0017] Figure 1 This is a schematic diagram illustrating the implementation process of an image processing-based aggregate particle shape and gradation identification method provided in an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes:

[0018] Step S100: Obtain the dynamic falling image stream of the aggregate particle group, and separate the closed contour image of a single aggregate particle from the dynamic falling image stream. The closed contour image contains a set of continuous pixel coordinates of the particle boundary.

[0019] A dynamic falling image stream is a sequence of images captured continuously by image acquisition devices such as high-speed cameras during the falling process of aggregate particles. This image stream can capture the different states and positions of the aggregate particles during their fall. The closed contour image of a single aggregate particle is extracted from the dynamic falling image stream and completely depicts the boundary of that individual aggregate particle. The boundary in this image is composed of a continuous set of pixel coordinates, which precisely define the shape of the aggregate particle.

[0020] Optionally, step S100 can be specifically implemented as the following steps S110-S160:

[0021] Step S110: Perform multi-source shadow removal processing on the dynamic falling image stream, and eliminate the high-light reflection area on the particle surface by separating the three primary color channels to generate a shadow-suppressed aggregate image sequence.

[0022] When performing multi-source shadow removal processing, multiple light sources at different angles can be arranged around the aggregate falling channel, such as at the top and sides of the channel. For example, a main light source can be set at the top to provide primary illumination, and two auxiliary light sources can be set on the sides to illuminate the aggregate particles at different angles. While capturing the dynamic falling image stream, images under different light sources are recorded simultaneously. Then, an image fusion algorithm is used to fuse the images from different light sources to eliminate shadows. For the separation of the three primary color channels, the channel separation function in image processing software can be used to separate the image into red, green, and blue channels. For each channel, the distribution of pixel values ​​is analyzed, and areas with excessively high pixel values ​​(i.e., high-brightness reflection areas) are eliminated by adjusting the pixel values. For example, a thresholding method can be used to reduce pixel values ​​above a certain threshold, thereby eliminating high-brightness reflection areas. Finally, a sequence of aggregate images with shadow suppression is generated, which more accurately reflects the true shape of the aggregate particles.

[0023] Step S120: Perform inter-frame difference operation on the aggregate image sequence to extract the motion region mask. Remove noise holes in the mask through morphological opening and closing operations to generate a binary image of the particle motion region. The foreground pixel value in the binary image of the particle motion region represents the area occupied by the particle.

[0024] During inter-frame difference operations, pixel-by-pixel subtraction can be performed on adjacent frames in the aggregate image sequence after shadow suppression. For example, for the nth frame and the (n+1)th frame, the difference between their corresponding pixel values ​​is calculated to obtain a difference image. Then, by setting an appropriate threshold, the difference image is converted into a binary image, i.e., a motion region mask. For morphological opening and closing operations, an opening operation is performed first, using a structuring element of an appropriate size (such as a 3×3 square structuring element) to erode the motion region mask, removing small noise points, followed by a dilation operation to restore the eroded motion region. Next, a closing operation is performed, using structuring elements of the same or different sizes to dilate the image after the opening operation, filling noise holes, followed by an erosion operation to connect adjacent motion regions. Finally, a binary image of the particle motion region is generated, which clearly shows the area occupied by the aggregate particles in the image.

[0025] Step S130: The watershed algorithm is used to segment the binary image of the particle motion region into cohesive particles. A set of foreground marker points is generated based on the distance transform. The particle boundary is separated by the marker-controlled watershed transform to obtain the initially separated single particle region.

[0026] When using the watershed algorithm for segmenting adherent particles, a distance transform is first applied to the binary image of the particle movement region. The Euclidean distance transform algorithm can be used to calculate the Euclidean distance from each pixel to the nearest background pixel, resulting in a distance-transformed image. In the distance-transformed image, regions with larger distance values ​​typically correspond to the center positions of the particles. Then, by setting an appropriate threshold, a set of foreground marker points is extracted from the distance-transformed image. For example, pixels with distance values ​​greater than a certain threshold are used as foreground marker points. Next, a marker-controlled watershed transform is used, employing the foreground marker point set as marker information to segment the binary image of the particle movement region. During the segmentation process, starting from the marker points, the flow of water is simulated to merge or separate adjacent regions, ultimately achieving the separation of particle boundaries and obtaining initially separated individual particle regions.

[0027] Step S140: Perform boundary tracking processing on the separated individual particle regions. Use the neighborhood connected component analysis algorithm to extract the outer contour pixel coordinates of the particle regions and generate an initial boundary coordinate sequence. The initial boundary coordinate sequence is arranged in a clockwise direction.

[0028] When performing boundary tracking, start from a boundary pixel in the granular region and traverse adjacent pixels according to certain rules (such as clockwise) until returning to the starting point. During the traversal, record the coordinates of each boundary pixel. For neighborhood connectivity analysis algorithms, four-neighbor or eight-neighbor analysis methods can be used. Taking four-neighbor analysis as an example, check the top, bottom, left, and right adjacent pixels of each pixel. If adjacent pixels also belong to the granular region, they are considered to belong to the same connected region. In this way, the outer contour pixel coordinates of the granular region are extracted. Then, these outer contour pixel coordinates are arranged clockwise to generate an initial boundary coordinate sequence. For example, a depth-first search algorithm can be used to implement boundary tracking and coordinate sorting, starting from a boundary pixel and recursively visiting its adjacent boundary pixels until the entire boundary is traversed, recording the coordinates of the visited pixels in sequence.

[0029] Step S150: Redundant points in the initial boundary coordinate sequence are removed using a boundary point simplification algorithm, while significant boundary feature points are retained, generating a simplified set of continuous pixel coordinates.

[0030] Redundant points are points on the boundary that contribute little to the description of the boundary shape; removing these points does not significantly change the overall shape of the boundary. Significant boundary feature points are points on the boundary with obvious characteristics, such as corner points or points with large curvature changes. The simplified set of continuous pixel coordinates is the set of coordinates obtained after processing by the boundary point simplification algorithm. This set contains fewer coordinate points but can still accurately represent the boundary shape of the particles.

[0031] In practical applications, the Douglas-Peucker algorithm can be used to simplify boundary points. First, connect the start and end points of the initial boundary coordinate sequence to form a straight line. Then, calculate the perpendicular distance from each point in the sequence to this line. If the distance from a point to the line is greater than a preset threshold, the point is considered a significant boundary feature and is retained; otherwise, the point is considered redundant and removed. Next, the above process is repeated for the subsequences divided by the retained points until all points have been processed.

[0032] Step S160: Construct a closed polygon based on a continuous set of pixel coordinates, process the contour breakage area through a boundary pixel point interpolation completion algorithm, and generate a closed contour image that satisfies topological closure.

[0033] When constructing a closed polygon based on a continuous set of pixel coordinates, the first and last coordinate points in the set can be connected to form a closed polygon. For contour break regions, a linear interpolation algorithm can be used. For example, assuming the pixel coordinates on both sides of the break region are (x1, y1) and (x2, y2), the coordinates of the pixel points to be inserted within the break region can be calculated based on these coordinates and the distance between them. Specifically, for the i-th pixel to be inserted, its horizontal coordinate xi = x1 + i*(x2 - x1) / n, and its vertical coordinate yi = y1 + i*(y2 - y1) / n, where n is the number of pixels to be inserted. In this way, new pixels are inserted into the break region, connecting the broken contour, ultimately generating a closed contour image that satisfies topological closure.

[0034] Step S200: Construct a boundary-constrained harmonic field for the closed contour image, and generate a harmonic field distribution that satisfies the contour boundary gradient condition in the internal region of the particle. The harmonic field distribution is characterized by the field intensity value change of each point inside the particle through a two-dimensional spatial field intensity value matrix.

[0035] Optionally, step S200 can be specifically implemented as the following steps S210-S260:

[0036] Step S210: Convert the set of boundary pixel coordinates of the closed contour image into a continuous curve model, and generate a smooth particle boundary model through curve fitting.

[0037] When converting the set of boundary pixel coordinates of a closed contour image into a continuous curve model, parametric curve representation methods, such as Bézier curves or spline curves, can be used. Taking Bézier curves as an example, based on the characteristics of the boundary pixel coordinate set, an appropriate number and position of control points are selected. Using the calculation formula for Bézier curves, the set of boundary pixel coordinates is converted into a parametric representation of the Bézier curve. For curve fitting, methods such as the least squares method can be used. Assuming a series of discrete points in the boundary pixel coordinate set are known, the least squares method is used to find a curve that minimizes the sum of the squared distances from the curve to these discrete points. For example, for a set containing 100 boundary pixel coordinate points, a cubic Bézier curve is used for fitting. By adjusting the positions of the four control points of the Bézier curve, the curve is made to approximate these discrete points as closely as possible, ultimately generating a smooth granular boundary model.

[0038] Step S220: Calculate the normal vector direction at each point of the boundary based on the particle boundary model, determine the field strength gradient constraint condition at the boundary, and establish the boundary constraint relationship group of the harmonic field. The boundary constraint relationship group includes the field strength value and normal derivative constraint at each point of the boundary.

[0039] When calculating the normal vector direction at each point of the boundary based on a particle boundary model, the curve equation of the particle boundary model can be differentiated to obtain the tangent vector direction at each point of the boundary. Then, the normal vector direction can be obtained through orthogonal transformation. For example, for a particle boundary model represented by a Bézier curve, the parametric equation of the Bézier curve can be differentiated to obtain the parametric representation of the tangent vector. Then, through orthogonal transformation of the vector, the tangent vector can be converted into a normal vector. The field intensity gradient constraint conditions at the boundary can be set according to the physical properties of the particles and actual needs. For example, for certain aggregate particles with specific shapes and materials, the magnitude of the field intensity gradient on the boundary can be set to be proportional to the curvature of the boundary, and the direction can be directed towards the interior of the particle. Through these settings, a set of boundary constraint relationships for the harmonic field is established. This set of relationships can be represented by a set of linear or nonlinear equations for subsequent solution of the harmonic field.

[0040] Optionally, step S220 can be specifically implemented as the following steps S221-S226:

[0041] Step S221: Calculate the first derivative of the curve relationship of the particle boundary model to obtain the tangent vector direction parameters at each point of the boundary. Generate the normal vector direction parameters through orthogonal vector transformation. The normal vector direction points to the internal region of the particle to ensure the internal directionality of the field strength distribution.

[0042] When calculating the first derivative of the curve relationship in a particle boundary model, an appropriate differentiation method should be selected based on the specific form of the particle boundary model. For example, if the particle boundary model is represented by parametric equations, the coordinate components of the tangent vector can be obtained by differentiating the parametric equations with respect to the parameters. For orthogonal vector transformations, the dot product and cross product operations of vectors can be used. Assuming the coordinate components of the tangent vector are known to be (x, y), by exchanging the coordinate components and changing the sign of one of the components, the coordinate components of the normal vector (-y, x) can be obtained, and then normalized so that the magnitude of the normal vector is 1. In this way, the normal vector direction parameter is generated, ensuring that the normal vector direction points towards the interior region of the particle.

[0043] Step S222: Determine the field strength gradient direction based on the normal vector direction parameter, set the reference value of the field strength gradient magnitude at the boundary, and generate the boundary gradient constraint vector. The reference value of the gradient magnitude is positively correlated with the curvature of the particle boundary to reflect the influence of the boundary curvature on the field strength.

[0044] The electric field gradient direction is the direction in which the electric field intensity changes most rapidly at each point on the boundary. It can be determined by the normal vector direction parameter, as the electric field gradient direction is usually consistent with the normal vector direction. The reference value for the electric field gradient magnitude is a reference value for the magnitude of the electric field gradient at each point on the boundary. This reference value is set considering its positive correlation with the curvature of the particle boundary, because the greater the curvature of the boundary, the more drastic the change in electric field intensity tends to be. The boundary gradient constraint vector is a vector containing information on the electric field gradient direction and magnitude, used to constrain the behavior of the harmonic field at the boundary. When determining the electric field gradient direction, since the normal vector direction points to the interior region of the particle and the electric field gradient direction is usually consistent with the normal vector direction, the normal vector direction can be directly used as the electric field gradient direction. To set the reference value for the electric field gradient magnitude at the boundary, the curvature of the particle boundary needs to be calculated first. This can be done by calculating the second derivative of the curve equation of the particle boundary model and combining the result of the first derivative. For example, for a boundary curve represented by a parametric equation, the curvature value at each point on the boundary can be obtained using relevant curvature calculation formulas. Then, based on a preset positive correlation function, the curvature value is mapped to a reference value for the field strength gradient magnitude. For example, a linear function can be set so that the reference value for the field strength gradient magnitude increases linearly with the increase of the curvature value. Finally, the field strength gradient direction and the reference value for magnitude are combined to generate a boundary gradient constraint vector. Based on this, the boundary gradient constraint vector can accurately reflect the characteristics of the field strength gradient at the boundary, demonstrating the influence of the boundary curvature on the field strength.

[0045] Step S223: Divide the internal region of the closed contour image into a finite element calculation mesh, number the mesh elements, and establish a relationship table between nodes and elements. The relationship table records the node index contained in each mesh element to realize mesh information management.

[0046] The finite element computational mesh divides the internal region of a closed contour image into several small units, which can be triangular, quadrilateral, or other shapes. The harmonic field is solved by analyzing and calculating these units.

[0047] When dividing the interior region of a closed contour image into a finite element computational mesh, different meshing algorithms can be used, such as the Delaunay triangulation algorithm or the quadrilateral meshing algorithm. Taking the Delaunay triangulation algorithm as an example, this algorithm can divide a set of points on a plane into a series of non-overlapping triangles, and the circumcircles of these triangles do not contain other points. First, a certain number of points are randomly distributed in the interior region of the closed contour image, and then the Delaunay triangulation algorithm is used to connect these points into triangular mesh elements. Next, each mesh element is numbered, either according to its position or the order in which it was generated. To establish the association table between nodes and elements, each mesh element can be traversed, recording the index of the nodes it contains.

[0048] Step S224: For each boundary grid node, establish a partial differential relationship between the field strength value and the field strength values ​​of adjacent nodes based on the boundary gradient constraint vector, and generate a node field strength constraint relationship. The node field strength constraint relationship includes a linear combination of the node's own field strength value and the field strength values ​​of adjacent nodes.

[0049] When establishing the partial differential relationship between the field strength value of each boundary grid node and the field strength values ​​of its neighboring nodes, this can be achieved using the boundary gradient constraint vector and the finite difference method. The finite difference method is a numerical method that transforms partial differential equations into difference equations. For a boundary grid node, a difference approximation relationship for the field strength value is established based on the positions of its neighboring nodes and the boundary gradient constraint vector. For example, for a boundary grid node, it may have two or three neighboring nodes. Based on the direction and magnitude of the boundary gradient constraint vector, combined with the field strength values ​​of the neighboring nodes, a difference equation for the field strength value of that boundary node is established. By performing this process on each boundary grid node, node field strength constraint relationships are generated. These node field strength constraint relationships can be represented by a set of linear equations, each equation containing a linear combination of the boundary node's own field strength value and the field strength values ​​of its neighboring nodes, reflecting the relationship between the field strength values ​​of the boundary node and its neighboring nodes.

[0050] Step S225: Collect the nodal field strength constraint relationships of all boundary nodes, construct a linear relationship set containing boundary conditions, and use row pivoting elimination to initialize the coefficient matrix of the relationship set to improve the solution stability.

[0051] When collecting the nodal field strength constraints of all boundary nodes, the nodal field strength constraints of each boundary node are arranged in a certain order to form a system of linear equations. The coefficient matrix of this system of equations contains the coefficients of the field strength values ​​of the boundary nodes and the field strength values ​​of adjacent nodes, while the constant terms are determined based on the boundary gradient constraint vector and other known conditions. For initializing the coefficient matrix of the relation group using row pivoting elimination, firstly, a row of the coefficient matrix is ​​selected as the pivot row. Then, through row transformation, the pivot coefficients of the pivot row are changed to 1, and the coefficients of the columns containing the pivots in other rows are eliminated to 0. When selecting the pivot, elements with larger absolute values ​​are preferred to improve the stability of the calculation. Through multiple such row transformations, the coefficient matrix is ​​transformed into a form that is easy to solve, thus completing the initialization of the coefficient matrix of the relation group.

[0052] Step S226: Strengthen the boundary conditions of the linear relationship group by using the penalty function method to transform the gradient constraint conditions into additional relations, thereby generating a complete set of boundary constraint relationships.

[0053] When strengthening boundary conditions for a set of linear relationships, the form of the penalty function is first determined. The penalty function is typically a function related to the gradient constraints; its value increases when the gradient constraints are not met. For example, the penalty function can be set as the sum of squares of the deviations from the gradient constraints. Then, the penalty function is added to the objective function of the linear relationship set, ensuring that solving for the linear relationship not only satisfies the original linear relationships but also minimizes the value of the penalty function. In this way, the gradient constraints are transformed into additional relations, which are then combined with the original linear relationships to generate a complete set of boundary constraint relations.

[0054] Step S230: Construct a two-dimensional orthogonal mesh in the internal region of the particle, map the pixel coordinate system of the closed contour image to the standardized mesh coordinate system, and generate a set of mesh node coordinates. The spacing of the mesh nodes is dynamically adjusted according to the particle size to adapt to particles of different sizes.

[0055] When constructing a two-dimensional orthogonal mesh, the orientation and starting position of the mesh are first determined. Typically, horizontal and vertical directions are chosen as the mesh line directions, and the starting position can be determined based on the position and size of the closed contour image. Then, the spacing of the mesh nodes is dynamically adjusted according to the particle size. The spacing of the mesh nodes can be determined according to preset rules by calculating the dimensions of the closed contour image (e.g., length and width). For example, for larger particles, the mesh node spacing can be set to a certain proportion of the maximum particle size; for smaller particles, this proportion is reduced accordingly. A linear transformation method can be used to map the pixel coordinate system of the closed contour image to the normalized mesh coordinate system. By determining the proportional relationship and offset between the pixel coordinate system and the normalized mesh coordinate system, the pixel coordinates are converted to normalized mesh coordinates. Finally, a set of mesh node coordinates is generated, in which the coordinate points accurately represent the position of each node in the two-dimensional orthogonal mesh, providing a foundation for subsequent harmonic field calculations.

[0056] Step S240: Solve the boundary constraint relationship set using the finite difference method, calculate the field strength distribution that satisfies the Laplace relation on the grid node coordinate set, generate the preliminary harmonic field, and normalize the field strength range of the preliminary harmonic field to unify the dimensions.

[0057] The finite difference method (FDM) is a method for numerically solving partial differential equations by transforming them into difference equations. This method can solve for a set of boundary constraint relations to obtain the field strength values ​​of the harmonic field at the grid nodes. The Laplace relation is a mathematical relationship satisfied by the harmonic field, where the sum of the second derivatives of the field strength is zero. The field strength distribution describes the distribution of the field strength values ​​at each node of the grid node coordinate set. The preliminary harmonic field is the initial result obtained through the finite difference method. Its field strength value range may vary; normalization can unify the field strength values ​​to a preset range, eliminating the influence of dimensions.

[0058] Optionally, step S240 can be specifically implemented as the following steps S241-S246:

[0059] Step S241: Classify the set of grid node coordinates, distinguish between boundary nodes and internal nodes, and establish a node type label array. In the node type label array, boundary nodes are marked with special identifiers to distinguish them from internal nodes.

[0060] When classifying the set of mesh node coordinates, the node type is determined based on the positional relationship between the mesh node's coordinates and the closed contour boundary. Whether a node is a boundary node can be determined by checking if its coordinates satisfy the equations or conditions of the closed contour boundary. For example, for a closed contour represented by a polygon, it can be determined whether the node's coordinates lie on the edge of the polygon. Mesh nodes identified as boundary nodes are marked with a special identifier, such as 1, in the node type label array; internal nodes are marked with 0. By traversing each node in the set of mesh node coordinates, the node type classification and labeling are completed, establishing the node type label array.

[0061] Step S242: For internal nodes, establish the field strength calculation relationship according to the difference scheme of the Laplace relation. The field strength calculation relationship means that the field strength value of an internal node is equal to the average of the field strength values ​​of its four neighboring nodes.

[0062] When establishing the calculation relationship for the field strength values ​​of internal nodes, the difference scheme of the Laplace relation is first clarified. For an internal node in a two-dimensional orthogonal mesh, its four neighboring nodes are typically the four adjacent nodes above, below, left, and right of that node. According to the difference scheme of the Laplace relation, the field strength value of an internal node is expressed as a linear combination of the field strength values ​​of its four neighboring nodes. Specifically, the field strength value of the internal node can be set as an unknown quantity, and the field strength values ​​of its four neighboring nodes can be set as known quantities. An equation is established such that the field strength value of the internal node is equal to the average of the field strength values ​​of its four neighboring nodes. For example, let the field strength value of the internal node be U. i,j The field strength values ​​of its four neighboring nodes are U i-1,j U i+1,j U i,j-1 U i,j+1The relationship for calculating the field strength can be expressed as U i,j =(U i-1,j +U i+1,j +U i,j-1 +U i,j+1 ) / 4.

[0063] Step S243: Merge the boundary constraint relationship group with the internal node field strength value calculation relationship to obtain a combined relationship group containing all nodes.

[0064] When merging the boundary constraint relationship set with the internal node field strength value calculation relationship, the boundary constraint relationship set and the internal node field strength value calculation relationship are arranged in a certain order. Typically, the equations for the boundary nodes are placed first, followed by the equations for the internal nodes. Then, the coefficient matrices and constant terms of these equations are merged to form a new coefficient matrix and constant term vector. The new coefficient matrix contains the coefficients of all node field strength values, and the constant term vector contains the constants determined based on the boundary conditions and the internal node calculation relationship. Through this merging operation, a simultaneous relationship set containing all nodes is obtained. This simultaneous relationship set is a large system of linear equations, which can be solved using appropriate numerical methods to obtain the field strength values ​​for all nodes.

[0065] Step S244: Solve the simultaneous relation set using an iterative method. Initialize the field strength values ​​of the internal nodes to be uniformly distributed. Assign field strength values ​​of the boundary nodes according to the boundary constraints. Iterate and update the field strength values ​​of the nodes to gradually approximate the true solution.

[0066] Iterative methods are numerical methods that continuously update solutions through iteration. For large systems of linear equations, iterative methods can gradually approximate the true solution.

[0067] When using iterative methods to solve simultaneous equations, the first step is to select a suitable iterative algorithm, such as the Jacobi iteration method or the Gauss-Seidel iteration method. Taking the Jacobi iteration method as an example, the linear equation system is solved by continuously updating the nodal field strength values. The internal nodal field strength values ​​are initialized to a uniform distribution; all internal nodal field strength values ​​can be initialized to a preset constant, such as 0. For boundary nodes, their field strength values ​​are determined according to the boundary constraint relationships. In each iteration, the new field strength value for each node is calculated based on the coefficient matrix and constant terms of the simultaneous equations. For internal nodes, the field strength value is updated according to the field strength value calculation relationship (the field strength value of an internal node is equal to the average of the field strength values ​​of its four neighboring nodes); for boundary nodes, their field strength values ​​remain unchanged (because they have already been assigned according to the boundary constraint conditions). This process is repeated continuously, and as the number of iterations increases, the nodal field strength values ​​gradually converge to the true solution. By monitoring the changes in field strength values ​​during the iteration process, it is determined whether the convergence condition has been met. If the overall deviation of the field strength values ​​between two adjacent iterations is less than a certain preset threshold, the iteration is stopped and the converged node field strength values ​​are obtained.

[0068] Step S245: Monitor the convergence of field strength values ​​during the iteration process, calculate the overall deviation of field strength values ​​between two adjacent iterations, stop the iteration when the overall deviation is less than the convergence judgment criterion, and obtain the set of converged node field strength values.

[0069] The convergence of the electric field strength during the iteration process refers to the situation where the nodal electric field strength values ​​gradually approach the true solution as the number of iterations increases during the iterative solution of the simultaneous relation system. The overall deviation of the electric field strength values ​​between two adjacent iterations is the sum of the changes in the electric field strength values ​​of all nodes in two adjacent iterations. By calculating this deviation, it is possible to determine whether the iteration has converged.

[0070] When monitoring the convergence of field strength values ​​during the iteration process, the difference between the field strength values ​​of all nodes in two adjacent iterations is calculated after each iteration. For each node, the absolute value of the difference between its field strength value in the current iteration and the previous iteration is calculated. Then, the differences of all nodes are summed to obtain the overall deviation of the field strength values ​​between two adjacent iterations. This overall deviation is compared with the convergence criterion. If the overall deviation is less than the convergence criterion, the iteration is considered to have converged and the iteration stops; otherwise, the next iteration continues. This process is repeated until the convergence condition is met, obtaining a set of converged node field strength values. This set provides accurate field strength data for the subsequent generation of the initial harmonic field.

[0071] Step S246: Arrange the set of nodal field strength values ​​in order of grid node coordinates to generate a two-dimensional field strength value matrix. The two-dimensional field strength value matrix is ​​determined as the preliminary harmonic field, and each element corresponds to the field strength value of the grid node to fully reflect the spatial distribution of the field strength.

[0072] When arranging the set of nodal field strength values ​​according to the grid node coordinates, the coordinate order of the grid nodes is first determined. This is typically done according to the row and column order of the grid, i.e., first by row, then by column. Then, the field strength values ​​from the set of nodal field strength values ​​are sequentially filled into the corresponding positions in the two-dimensional field strength value matrix according to the coordinate order of the grid nodes. For example, for the node in the i-th row and j-th column of a two-dimensional orthogonal grid, its field strength value is filled into the element position in the i-th row and j-th column of the two-dimensional field strength value matrix. In this way, a two-dimensional field strength value matrix is ​​generated. This matrix is ​​defined as the initial harmonic field, with each element corresponding to the field strength value of a grid node. The matrix form visually displays the spatial distribution of the harmonic field's field strength within a closed contour region.

[0073] Step S250: Perform boundary consistency verification on the preliminary harmonic field, calculate the deviation between the actual value of the field strength gradient at the boundary and the theoretical constraint value, and correct the field strength value of the internal nodes by relaxation iteration method until the deviation meets the accuracy requirements.

[0074] Boundary consistency verification checks the field strength gradient at the boundary of the initial harmonic field to ensure that it meets the boundary constraints. The deviation is the difference between the actual value of the field strength gradient and the theoretical constraint value. By calculating this deviation, it can be determined whether the initial harmonic field meets the requirements at the boundary. The relaxation iteration method is an iterative method used to correct the field strength values ​​at the nodes. By continuously adjusting the field strength values ​​of the internal nodes, the deviation of the field strength gradient at the boundary is gradually reduced until the accuracy requirements are met.

[0075] When performing boundary consistency verification on the preliminary harmonic field, the actual value of the field strength gradient at the boundary is first calculated. This can be achieved by performing a difference calculation on the two-dimensional field strength matrix of the preliminary harmonic field to obtain the gradient value of the field strength at the boundary. Then, the calculated actual value of the field strength gradient is compared with the theoretical constraint value to calculate the deviation between the two. For each boundary node, the absolute value of the difference between its actual field strength gradient and the theoretical constraint value is calculated, and the deviation values ​​of all boundary nodes are added together to obtain the total deviation. If the total deviation is greater than the accuracy requirement, the field strength values ​​of the internal nodes need to be corrected using a relaxation iteration method. In the relaxation iteration method, the field strength values ​​of the internal nodes are adjusted according to the field strength gradient deviation at the boundary. For example, for internal nodes that are close to the boundary and have a significant impact on the boundary field strength gradient, their field strength values ​​are appropriately adjusted according to the magnitude and direction of the deviation. By continuously iterating and updating the field strength values ​​of the internal nodes, the actual value of the field strength gradient and the deviation at the boundary are recalculated until the deviation meets the accuracy requirements.

[0076] Step S260: Convert the corrected grid node field strength values ​​into a two-dimensional spatial field strength value matrix in the original image coordinate system, and generate a harmonic field distribution that satisfies the contour boundary gradient condition. The number of rows and columns of the two-dimensional spatial field strength value matrix is ​​consistent with the size of the closed contour image.

[0077] When converting the corrected grid node field strength values ​​into a two-dimensional spatial field strength matrix in the original image coordinate system, it is first necessary to establish a mapping relationship between the grid node coordinate system and the original image coordinate system. This can be achieved by inverse transformation of the mapping relationship from the pixel coordinate system to the normalized grid coordinate system established in step S230, converting the grid node coordinates to the original image coordinates. Then, based on the field strength values ​​of the grid nodes and the mapping relationship, the field strength values ​​are interpolated to the pixels of the original image. Interpolation can be performed using methods such as bilinear interpolation, calculating the field strength value of each pixel in the original image based on the field strength values ​​of the grid nodes. The interpolated field strength values ​​are arranged in the pixel coordinate order of the original image to form a two-dimensional spatial field strength matrix. The number of rows and columns of this matrix is ​​consistent with the size of the closed contour image, which can completely reflect the distribution of the harmonic field in the original image.

[0078] Step S300: Extract stable critical structures based on harmonic field distribution, identify critical inflection points of field strength value changes, and encode the spatial coordinate sequence of critical inflection points into a topologically invariant sequence.

[0079] Stable critical structure extraction involves identifying structures with specific characteristics from a harmonic field distribution. Critical inflection points in field strength variation are points in the harmonic field where the trend of field strength change alters, such as local maxima or minima. Topologically invariant sequences are sequences obtained by encoding the spatial coordinates of critical inflection points. These sequences possess topological invariance, meaning that regardless of the shape deformation of aggregate particles, the sequence remains unchanged as long as their topological structure remains constant. This facilitates the comparison and analysis of the characteristics of different aggregate particles.

[0080] Optionally, step S300 can be specifically implemented as the following steps S310-S360:

[0081] Step S310: Perform gradient calculation on the two-dimensional spatial field strength value matrix of the harmonic field distribution to generate the field strength gradient magnitude matrix and gradient direction matrix. The gradient magnitude matrix represents the rate of change of the field strength value, and the gradient direction matrix represents the direction of change to reflect the spatial trend of field strength change.

[0082] The electric field gradient is a physical quantity describing the variation of electric field values ​​in space, including the magnitude and direction of the gradient. The electric field gradient magnitude matrix is ​​a matrix composed of the magnitudes of the electric field gradients at each point in the two-dimensional electric field value matrix of the harmonic field distribution. This matrix characterizes the rate of change of the electric field value in space; the larger the magnitude, the faster the change. The gradient direction matrix is ​​a matrix composed of the directions of the electric field gradients at each point. This matrix characterizes the direction of change of the electric field value in space, reflecting the spatial trend of the electric field's variation.

[0083] Optionally, step S310 can be specifically implemented as the following steps S311-S316:

[0084] Step S311: Use the gradient calculation operator to perform convolution operation on the two-dimensional spatial field strength matrix of the harmonic field distribution, and calculate the first-order partial derivative matrices in the horizontal and vertical directions respectively. The horizontal partial derivative matrix represents the rate of change of the horizontal field strength.

[0085] When performing convolution operations on the two-dimensional spatial field intensity matrix of a harmonic field distribution using gradient calculation operators, a suitable gradient calculation operator is typically chosen, such as the Sobel operator or the Prewitt operator. Taking the Sobel operator as an example, this operator contains two convolution kernels, one horizontal and one vertical. For the horizontal convolution kernel, convolution is performed with the two-dimensional spatial field intensity matrix to obtain the first-order partial derivative matrix in the horizontal direction. During the convolution operation, the convolution kernel is multiplied by each local region in the matrix and summed to obtain the horizontal partial derivative at the center point of that local region. By traversing the entire two-dimensional spatial field intensity matrix, the calculation of the first-order partial derivative matrix in the horizontal direction is completed. Each element in this matrix reflects the change in field intensity in the horizontal direction, i.e., the rate of change of the lateral field intensity. Similarly, using the vertical convolution kernel to perform convolution operations with the two-dimensional spatial field intensity matrix yields the first-order partial derivative matrix in the vertical direction, reflecting the rate of change of the field intensity in the vertical direction.

[0086] Step S312: Perform square root operation on the partial derivative matrix in the horizontal direction and the partial derivative matrix in the vertical direction to generate the field strength gradient magnitude matrix. Each element in the gradient magnitude matrix is ​​the magnitude of the gradient vector at the corresponding position to quantify the rate of change.

[0087] The square root operation is a mathematical operation used to calculate the magnitude of a vector. In this step, the magnitude of the field intensity gradient is obtained by performing the square root operation on corresponding elements of the horizontal and vertical partial derivative matrices. The field intensity gradient magnitude matrix is ​​a matrix composed of the magnitudes of the field intensity gradient at each location. Each element in this matrix corresponds to the magnitude of the field intensity gradient at the corresponding location in the two-dimensional field intensity value matrix, quantifying the rate of change of the field intensity value at that location.

[0088] Step S313: Calculate the gradient direction angle using the direction angle calculation method. Convert the ratio of the elements of the horizontal direction partial derivative matrix and the vertical direction partial derivative matrix into angle values ​​to generate the gradient direction matrix. The angle values ​​reflect the direction of field strength change.

[0089] The orientation angle calculation method is used to convert the ratio of vector components into an angle value. In this step, the ratio of the elements of the horizontal and vertical partial derivative matrices is calculated and converted into a gradient orientation angle, thereby generating a gradient orientation matrix. Each element in the gradient orientation matrix corresponds to the orientation angle of the field intensity gradient at the corresponding position in the two-dimensional spatial field intensity matrix. This angle value reflects the direction of change of the field intensity at that position.

[0090] Step S314: Perform non-maximum suppression processing on the gradient direction matrix to retain local maxima points in the gradient direction, enhance edge localization accuracy, and generate a refined gradient direction matrix.

[0091] Non-local maxima suppression (NMS) is a method used for image edge detection. By suppressing non-local maxima points and retaining only local maxima points along the gradient direction, it can refine edges and improve the accuracy of edge localization. Each element in the gradient direction matrix represents the direction of the field intensity gradient. Performing NMS along the gradient direction removes points that are not local maxima, making the edges clearer. The refined gradient direction matrix is ​​the matrix obtained after NMS processing. The elements in this matrix only retain information about local maxima points along the gradient direction, more accurately reflecting the edge positions of field intensity changes.

[0092] Step S315: Smooth the field strength gradient magnitude matrix to generate a smoothed gradient magnitude matrix. The smoothness of the filtering process is dynamically adjusted according to the gradient noise level.

[0093] Smoothing is a method used to reduce image noise, such as Gaussian filtering and mean filtering. By smoothing the field intensity gradient magnitude matrix, noise points in the matrix can be removed, making the changes in gradient magnitudes smoother. The smoothing degree of the filtering process is dynamically adjusted according to the level of gradient noise. The filtering parameters are automatically adjusted based on the degree of noise in the field intensity gradient magnitude matrix. Regions with higher noise levels are treated with smoother filtering, while regions with lower noise levels are treated with relatively weaker filtering, in order to remove noise while preserving the detailed information of the field intensity gradient. The smoothed gradient magnitude matrix is ​​the matrix obtained after smoothing; the changes in gradient magnitudes in this matrix are smoother, which is beneficial for subsequent thresholding and region growing processes.

[0094] Step S316: Align the smoothed gradient magnitude matrix and the refined gradient direction matrix in terms of dimensions to ensure that they have the same number of rows and columns and that the corresponding pixels correspond one-to-one.

[0095] When aligning the smoothed gradient magnitude matrix and the refined gradient direction matrix in terms of dimensions, no adjustment is needed if the number of rows and columns of the two matrices are already the same. If the number of rows and columns differs, cropping or padding operations can be performed depending on the specific situation. For example, if the smoothed gradient magnitude matrix has more rows than the refined gradient direction matrix, the extra rows in the smoothed gradient magnitude matrix can be cropped; if it has fewer rows than the refined gradient direction matrix, zero rows can be padded below the smoothed gradient magnitude matrix. The same applies to the number of columns. Through these operations, it is ensured that the smoothed gradient magnitude matrix and the refined gradient direction matrix have the same number of rows and columns, and that corresponding pixels correspond one-to-one. Based on this, in subsequent processing, such as extracting regions with significant field intensity changes and performing region growing, gradient magnitude and direction information can be conveniently used simultaneously, improving the accuracy and efficiency of the processing.

[0096] Step S320: Set a dynamic threshold in the gradient magnitude matrix, calculate the segmentation threshold of the gradient magnitude using an adaptive threshold determination method, and extract the region where the gradient magnitude is greater than the segmentation threshold as the region of significant field strength change. The region of significant field strength change reflects the location of drastic field strength change.

[0097] When setting a dynamic threshold, an adaptive thresholding method is used to calculate the segmentation threshold based on the gradient magnitude. Various adaptive thresholding methods can be used, such as the Otsu algorithm and local thresholding. Taking the Otsu algorithm as an example, this algorithm determines the optimal segmentation threshold by maximizing the inter-class variance. First, the grayscale distribution of all pixels in the gradient magnitude matrix is ​​statistically analyzed, and the inter-class variance is calculated for different thresholds. The inter-class variance represents the degree of difference between the two classes when the image is divided into two categories (regions with significant field intensity changes and non-significant regions). By traversing all possible thresholds, the threshold that maximizes the inter-class variance is found and used as the segmentation threshold. Then, in the gradient magnitude matrix, regions with gradient magnitudes greater than the segmentation threshold are extracted as regions with significant field intensity changes. This can be achieved by traversing each element in the gradient magnitude matrix, marking pixels with values ​​greater than the segmentation threshold as pixels in regions with significant field intensity changes, and pixels with values ​​less than the segmentation threshold as pixels in regions with non-significant changes.

[0098] Step S330: Perform region growing processing on the region with significant field strength change. Starting from the pixel with the largest gradient direction change rate, merge adjacent pixels according to the gradient direction consistency criterion to generate a connected region with field strength change. The connected region is a continuous region with field strength change.

[0099] When performing region growing on regions with significant field strength changes, the pixel with the largest gradient direction change rate is first identified as the seed point. This can be achieved by calculating the gradient direction change rate of each point in the gradient direction matrix and selecting the point with the largest change rate. Then, starting from the seed point, adjacent pixels are merged according to the gradient direction consistency criterion. For each neighboring pixel of the seed point, the angle between its gradient direction and the seed point's gradient direction is calculated. If the angle is less than a preset threshold, the neighboring pixel is considered to have a similar gradient direction to the seed point and is merged into the current region. Next, the merged region is used as a new seed region, and the process of finding its neighboring pixels is repeated until no neighboring pixels meet the criteria. By continuously merging neighboring pixels, a connected region with varying field strength is generated.

[0100] Step S340: Within each connected region of field strength change, trace the trajectory of field strength value change along the gradient direction, record the local maximum and minimum points of field strength value on the trajectory, and determine the critical turning point of field strength value change. The critical turning point reflects the location of the change in field strength trend.

[0101] The electric field strength (EVS) change trajectory is a path recorded along the gradient direction within each connected region of EVS change. Local maxima and minima are the points on the trajectory where the EVS reaches its local maximum and minimum values. These points reflect the locations where the EVS change trend changes, i.e., the critical turning points. By identifying these critical turning points, the characteristics of EVS change can be accurately grasped.

[0102] Optionally, step S340 can be specifically implemented as the following steps S341-S346:

[0103] Step S341: Within each connected region of field strength change, calculate the gradient magnitude of all pixels, and select the pixel with the largest gradient magnitude as the trajectory starting point. The starting point is the location where the field strength change is most drastic.

[0104] When calculating the gradient magnitude of all pixels within each connected region of field strength change, the field strength gradient magnitude matrix generated in step S312 can be used directly. For each pixel within the connected region of field strength change, its corresponding gradient magnitude is obtained from the field strength gradient magnitude matrix. Then, the gradient magnitudes of all pixels are compared, and the pixel with the largest gradient magnitude is found. This can be achieved by traversing all pixels within the connected region of field strength change, recording the gradient magnitude of each pixel, and using a variable to record the current largest gradient magnitude and the corresponding pixel coordinates. After traversing all pixels, the pixel with the largest gradient magnitude is the starting point of the trajectory. This starting point represents the location where the field strength change is most drastic. Tracing the trajectory of field strength value change from this point can more accurately reflect the situation of field strength change and provide a good starting point for subsequently determining the critical turning point of field strength value change.

[0105] Step S342: Starting from the trajectory starting point, move to the adjacent pixel point along the direction indicated by the gradient direction matrix, record the pixel coordinates and corresponding field strength values ​​on the moving path, and generate the initial field strength change trajectory. The trajectory reflects the change process of the field strength value along the gradient direction.

[0106] When moving along the gradient direction from the starting point of the trajectory, the gradient direction of that point is determined based on the element value of the starting point in the gradient direction matrix. Then, adjacent pixels are found along this gradient direction. For example, if the gradient direction is 0°, the pixel to the right of the current point is moved; if the gradient direction is 90°, the pixel above the current point is moved. The moved pixel is taken as the new current point, and its coordinates and corresponding field strength value are recorded. The field strength value can be obtained from the two-dimensional field strength value matrix. This process is repeated continuously, moving sequentially along the direction indicated by the gradient direction matrix to adjacent pixels until no further movement is possible (e.g., reaching the boundary of a connected region of field strength change). The coordinates and corresponding field strength values ​​of all pixels along the movement path are recorded sequentially to generate the initial field strength change trajectory.

[0107] Step S343: Smooth the initial field strength change trajectory by using a moving average filter to eliminate trajectory jitter and generate a smoothed field strength change trajectory. The window size for smoothing is determined based on the degree of trajectory fluctuation.

[0108] When smoothing the initial field strength variation trajectory, a moving average filtering method is used. First, the window size for smoothing is determined. This can be determined by analyzing the fluctuation of the initial field strength variation trajectory. For example, the standard deviation of the field strength values ​​on the trajectory can be calculated; a larger standard deviation indicates greater trajectory fluctuation, and the window size can be increased accordingly. Then, starting from the first point of the initial field strength variation trajectory, a window of the same size is selected, and the average of all field strength values ​​within the window is calculated. This average is used as the new field strength value at the center of the window. Next, the window is moved one position forward, and the above calculation process is repeated until the entire initial field strength variation trajectory has been traversed. In this way, a smoothed field strength variation trajectory is generated. The field strength value changes in this trajectory are smoother, reducing the impact of jitter and providing more reliable data for subsequently accurately determining the local maxima and minima of the field strength values.

[0109] Step S344: Calculate the change in field intensity value of adjacent pixels along the smoothed field intensity change trajectory. When the change changes from positive to negative, it is marked as a local maximum point, and when it changes from negative to positive, it is marked as a local minimum point.

[0110] The change in field strength is the difference in field strength values ​​between adjacent pixels in the smoothed field strength change trajectory. By calculating the change in field strength, the trend of field strength change can be determined. When the change changes from positive to negative, it indicates that the field strength value has changed from increasing to decreasing, and this point is a local maximum; when the change changes from negative to positive, it indicates that the field strength value has changed from decreasing to increasing, and this point is a local minimum. Marking these local maxima and minima can identify the critical inflection points of field strength value change, providing a basis for subsequent screening and confirmation.

[0111] Step S345: Filter the marked local maxima and minima, remove adjacent extreme points whose distance is less than the distance threshold, and generate a set of candidate critical inflection points.

[0112] When filtering the marked local maxima and minima, a distance threshold is first determined. This threshold can be determined based on the characteristics of the smoothed field strength variation trajectory and actual needs. For example, for cases with small trajectory fluctuations, the distance threshold can be appropriately reduced; for cases with large trajectory fluctuations, the distance threshold can be increased. Then, the marked local maxima and minima are traversed, and for each pair of adjacent extreme points, the distance between them is calculated. The distance can be represented by the difference in pixel indices of the two points on the trajectory. If the distance between adjacent extreme points is less than the distance threshold, one of the extreme points is removed. Extreme points with smaller field strength value changes can be selected for removal. By traversing all adjacent extreme points, the filtering process is completed, generating a set of candidate critical inflection points.

[0113] Step S346: Calculate the second derivative of the field strength value at each point in the candidate critical inflection point set. When the absolute value of the second derivative exceeds the curvature judgment criterion, it is confirmed as a critical inflection point of field strength value change. The second derivative is calculated by the difference method to reflect the curvature characteristics of field strength value change.

[0114] When calculating the second derivative of the electric field strength at each point in the candidate critical inflection point set, a finite difference method is used. For each point in the candidate critical inflection point set, the second derivative is calculated by considering its neighboring points. Taking the central difference method as an example, for a certain candidate critical inflection point, the electric field strength values ​​of its preceding and following points are taken, and the second derivative of that point is approximated by the relationship between the electric field strength values ​​of these three points. Specifically, the difference between the electric field strength value of this point and the preceding point, and the difference between the electric field strength value of this point and the following point are first calculated, and then these two differences are further calculated to obtain an approximate value of the second derivative.

[0115] After calculating the second derivative of the electric field strength at each point, its absolute value is compared with the curvature judgment standard. The curvature judgment standard is pre-set based on actual conditions and experience, used to distinguish between true critical turning points and pseudo-turning points caused by data fluctuations. When the absolute value of the second derivative of the electric field strength at a point exceeds the curvature judgment standard, it indicates that the curvature of the electric field strength change at that point is large, and the trend of electric field strength change has changed significantly. This point is then identified as a critical turning point in the change of electric field strength. Points whose absolute value of the second derivative does not exceed the curvature judgment standard are excluded from being considered critical turning points.

[0116] Step S350: Extract the spatial coordinates of the critical turning points, arrange them in the trajectory tracking order to generate a spatial coordinate sequence of the critical turning points, and the spatial coordinate sequence contains the planar coordinate information of the turning points to reflect their spatial position.

[0117] The spatial coordinates of critical inflection points are the specific location coordinates of these inflection points in the two-dimensional image. By extracting these coordinates and arranging them according to the trajectory tracing order, an ordered sequence can be formed. The trajectory tracing order is the order determined when tracing the trajectory of field strength value changes along the gradient direction in step S340. Arranging the spatial coordinates of critical inflection points in this order can completely reflect the process and path of field strength value changes. Each element in the spatial coordinate sequence represents the planar coordinate information of a critical inflection point. This information can accurately describe the position of the inflection point in the image, which is of great significance for subsequent analysis of aggregate particle characteristics and topological structure encoding.

[0118] Step S360: Perform topological structure encoding on the spatial coordinate sequence, using the relative distance encoding method to convert absolute coordinates into relative coordinates, and generate a topologically invariant sequence.

[0119] When topologically encoding a spatial coordinate sequence, a relative distance encoding method is first employed. For each critical inflection point in the spatial coordinate sequence, its relative distance and relative positional relationship with other inflection points are calculated. Relative distance can be determined by calculating the distance between two points, while relative positional relationship can be described by calculating information such as the angle between two points. For example, for the first critical inflection point in the sequence, its distance and angle with the second inflection point are calculated, then its distance and angle with the third inflection point are calculated, and so on. These relative distances and relative positional relationships are recorded in a specific order to form relative coordinate information. Then, this relative coordinate information is encoded to generate a topologically invariant sequence. This encoding process can use predefined encoding rules to convert relative distances and relative positional relationships into numerical or symbolic sequences. In this way, the spatial coordinate sequence, which originally depended on absolute coordinates, is transformed into a topologically invariant sequence. The generated topologically invariant sequence can be used for subsequent semantic mapping processing. By interacting with a shape semantic mapper calibrated by concrete interface performance, the quantitative characterization of aggregate particles' interlocking ability at the aggregate-paste interface can be further analyzed.

[0120] Step S400: Call the shape semantic mapper calibrated by the concrete interface performance to perform semantic mapping processing on the topological invariant sequence, and output a quantitative characterization of the interlocking ability of a single aggregate particle at the aggregate-paste interface. The quantitative characterization includes the interface contact area ratio and interlocking depth distribution parameters.

[0121] The shape semantic mapper calibrated for concrete interface performance is a trained and calibrated model that maps topologically invariant sequences into semantic information related to concrete interface performance. Semantic mapping involves inputting the topologically invariant sequence into the shape semantic mapper, which then transforms the feature information in the sequence into a specific description of the interlocking ability of aggregate particles at the aggregate-paste interface. The quantitative characterization of the interlocking ability of a single aggregate particle at the aggregate-paste interface is achieved through a series of parameters obtained from the semantic mapping process. These parameters accurately describe the interlocking situation between aggregate particles and the paste. The interface contact area ratio reflects the proportion of the actual contact area between the aggregate and the paste to the total aggregate surface area, while the interlocking depth distribution parameters describe the distribution characteristics of the interlocking depth between the aggregate and the paste, such as central tendency and dispersion.

[0122] Optionally, step S400 can be specifically implemented as the following steps S410-S460:

[0123] Step S410: Input the topologically invariant sequence into the input layer of the shape semantic mapper calibrated by the concrete interface performance, and convert the topologically invariant sequences of different lengths into fixed-dimensional input vectors through sequence length normalization.

[0124] When inputting topologically invariant sequences into the input layer of the shape semantic mapper, the sequence length is first normalized. Various methods can be used for normalization; for example, for shorter topologically invariant sequences, their length can be extended to a fixed dimension by padding with specific values ​​(such as zeros); for longer sequences, their length can be reduced to a fixed dimension by truncation or sampling. This process transforms topologically invariant sequences of different lengths into input vectors with the same dimension. Then, this fixed-dimensional input vector is fed into the input layer of the shape semantic mapper, preparing it for subsequent feature extraction and semantic mapping processing within the mapper.

[0125] Step S420: Based on the hidden layer of the shape semantic mapper, feature extraction is performed on the input vector, and nonlinear transformation is performed through a multi-layer network structure to generate a shape feature vector. The shape feature vector includes the angular features of the particles, surface texture features, and overall morphological features.

[0126] The hidden layers of the shape semantic mapper are network layers within the mapper that perform feature extraction and processing. These layers contain multiple neurons and network structures, enabling in-depth analysis and transformation of the input vector. Feature extraction extracts feature information related to the shape of aggregate particles from the input vector. The multi-layer network structure then combines and transforms the extracted features through a series of nonlinear transformations, generating a more representative shape feature vector. This shape feature vector includes the particle's angular features, surface texture features, and overall morphological features. These features are crucial for describing the shape of aggregate particles and predicting their interlocking ability at the aggregate-slurry interface.

[0127] Optionally, the hidden layer of the shape semantic mapper contains three feature extraction sub-layers, corresponding to the extraction of corner features, surface texture features, and overall shape features, respectively. Each sub-layer uses differentiated network structure parameters to extract different types of features in a targeted manner. Based on this, step S420 can be specifically implemented as the following steps S421-S425:

[0128] Step S421: The input vector is convolved based on the corner feature extraction sub-layer. Local features are extracted from the topologically invariant sequence using convolution kernels of different sizes to generate a corner feature map. The corner feature map reflects the angular change features in the sequence.

[0129] When performing convolution operations on the input vector in the corner feature extraction sublayer, different sizes of convolution kernels are first selected. Different sized kernels can capture local features at different scales, which is very useful for extracting corner features. For example, smaller kernels can capture more subtle angle changes, while larger kernels can capture more macroscopic corner features. Then, the kernels are convolved with the input vector. By sliding the kernel across the input vector, the sum of the products of the kernel and the local regions of the input vector is calculated to obtain the convolution result.

[0130] The convolution results of all convolution kernels are combined to generate an edge feature map. Each element in this feature map represents the local feature response of the input vector at the corresponding position, reflecting the angular variation characteristics in the topologically invariant sequence. Through the edge feature map, the location and feature intensity of potential edges in the topologically invariant sequence can be visually observed, providing important information for further analysis of the edge characteristics of aggregate particles.

[0131] Step S422: Based on the surface texture feature extraction sub-layer, a temporal network structure is used to extract temporal features from the input vector. The texture fluctuation features in the topologically invariant sequence are captured by the temporal memory unit to generate a surface texture feature vector. The surface texture feature vector represents the micro-texture changes on the particle surface.

[0132] When using a temporal network structure to extract temporal features from the input vector in the surface texture feature extraction sublayer, the input vector is sequentially fed into the temporal network structure. The temporal memory unit updates its state based on the current value of the input vector and its previous memory state. For example, for a Long Short-Term Memory (LSTM) network unit, it updates the states of its input gate, forget gate, and output gate based on the input vector, thus determining which information needs to be retained and which needs to be forgotten. By continuously updating the state of the temporal memory unit, the temporal network structure can capture texture fluctuation features in topologically invariant sequences. These features reflect the microscopic texture changes on the particle surface, such as texture roughness and periodicity. The captured texture fluctuation features are then organized and encoded to generate a surface texture feature vector.

[0133] Step S423: The input vector is processed by global pooling based on the overall morphological feature extraction sub-layer. The overall distribution features of the sequence are extracted by combining global max pooling and global average pooling to generate an overall morphological feature vector. The overall morphological feature vector represents the macroscopic shape characteristics of the particles.

[0134] When performing global pooling on the input vector based on the overall morphological feature extraction sublayer, a global max pooling operation is first performed to find the maximum value among all elements of the input vector. These maximum values ​​are then combined into a vector, reflecting the most prominent feature information in the input vector, which is very useful for describing the features of certain key parts of the particle. Next, a global average pooling operation is performed to calculate the average value of all elements in the input vector, resulting in an average value vector. This vector reflects the overall average level of the input vector and can represent the overall size and average characteristics of the particle. The vectors obtained from global max pooling and global average pooling are combined to generate the overall morphological feature vector. This vector integrates the maximum and average value information of the input vector, enabling a more comprehensive characterization of the macroscopic shape properties of the particle.

[0135] Step S424: Concatenate the edge feature map, surface texture feature vector and overall shape feature vector dimensionally to generate a concatenated feature vector.

[0136] When performing dimensional concatenation, the corner feature map, surface texture feature vector, and overall shape feature vector can be arranged in a predetermined order, and then their elements can be combined. For the corner feature map, it may be necessary to convert it into vector form for concatenation with the surface texture and overall shape feature vectors. For example, the corner feature map can be expanded into a vector by rows or columns. Then, the expanded corner feature map vector, surface texture feature vector, and overall shape feature vector are concatenated in sequence to form a new vector, i.e., the concatenated feature vector.

[0137] Step S425: Perform a nonlinear transformation on the spliced ​​feature vector through a fully connected layer, and introduce a nonlinear relationship by using a nonlinear activation function to generate a shape feature vector containing multi-dimensional shape information.

[0138] When performing a nonlinear transformation on the concatenated feature vector using a fully connected layer, the fully connected layer first performs a linear transformation on the concatenated feature vector. Each neuron in the fully connected layer is connected to all elements of the concatenated feature vector. By calculating the weighted sum and adding a bias, the output of that neuron is obtained. The outputs of all neurons are combined to form a new vector. This new vector is the result of the concatenated feature vector after a linear transformation, but it is still linear and may not be able to represent complex shape features well.

[0139] Then, a nonlinear activation function is used to process the output of the fully connected layer. Common nonlinear activation functions include ReLU (Modified Linear Unit) and the Sigmoid function. Taking the ReLU function as an example, it sets values ​​less than 0 to 0 and leaves values ​​greater than 0 unchanged, thus introducing a nonlinear relationship. The output vector of the fully connected layer is input into the nonlinear activation function, and after processing by the function, a new vector is obtained. This vector contains multi-dimensional shape information, which can more accurately describe the comprehensive shape characteristics of aggregate particles, such as the combined performance of features like edges, surface texture, and overall morphology.

[0140] Step S430: Call the concrete interface performance association layer of the shape semantic mapper, match the shape feature vector with the pre-stored concrete interface performance database, and extract the interface interlocking parameters related to the shape features. The concrete interface performance database contains experimental data on interface performance corresponding to different aggregate shapes.

[0141] The concrete interface performance association layer of the shape semantic mapper is specifically designed to match shape feature vectors with a concrete interface performance database. It can find interface performance information corresponding to aggregate shapes similar to the shape feature vectors in the database. The pre-stored concrete interface performance database contains a large amount of experimental data on interface performance corresponding to different aggregate shapes. This data was obtained through actual experimental measurements and reflects the interfacial bonding ability and other properties of different aggregate shapes when in contact with the paste. The interfacial bonding parameters related to the shape features are matched from the database. These parameters describe the bonding between aggregate particles with similar shapes and the paste, and are of great significance for evaluating the bonding ability of individual aggregate particles at the aggregate-paste interface.

[0142] Optionally, the concrete interface performance database stores the shape feature vectors of multiple standard aggregate shape samples and their corresponding interface interlocking parameters. Based on this, step S430 can be specifically implemented as the following steps S431-S435:

[0143] Step S431: Calculate the similarity measure between the shape feature vector and the shape feature vector of each standard sample in the database, and generate a similarity matrix. Each element in the similarity matrix represents the degree of shape similarity between the aggregate to be identified and the standard sample.

[0144] When calculating the similarity measure between the shape feature vector of the aggregate to be identified and the shape feature vectors of each standard sample in the database, a suitable similarity measurement method is first selected, such as Euclidean distance or cosine similarity. The calculated similarity values ​​are then arranged sequentially to form a similarity matrix. For example, for a database with n standard samples, the similarity matrix is ​​an n×1 matrix, where the elements in the i-th row represent the degree of similarity between the shape feature vector of the aggregate to be identified and the shape feature vector of the i-th standard sample. The similarity matrix provides a clear visual representation of the shape similarity between the aggregate to be identified and each standard sample, offering a basis for subsequently selecting a matching sample set.

[0145] Step S432: Select multiple standard samples with the highest similarity in the similarity matrix as the matching sample set. The number of samples is determined according to the sample distribution in the database to ensure the representativeness of the matching.

[0146] When selecting the highest similarity standard samples from the similarity matrix as the matching sample set, the elements in the similarity matrix are first sorted according to their similarity from highest to lowest. Then, the number of samples to be selected is determined based on the distribution of samples in the database. If the sample distribution in the database is relatively uniform, the number of samples can be relatively small; if the sample distribution is uneven, more samples may need to be selected to cover different types of aggregate shapes.

[0147] Step S433: Perform a weighted average on the interface bite parameters of the matching sample set, using the similarity metric as the weight coefficient, to generate the weighted average interface bite parameters.

[0148] When performing a weighted average of the interface engagement parameters of the matching sample set, the similarity metric value of each standard sample is first determined as a weighting coefficient. These similarity metrics, calculated in step S431, reflect the degree of shape similarity between each standard sample and the aggregate to be identified. Then, for the interface engagement parameter of each standard sample in the matching sample set, it is multiplied by the corresponding weighting coefficient. For example, for a certain interface engagement parameter (such as a contact area parameter), the parameter value of the first standard sample is multiplied by its similarity metric value, the parameter value of the second standard sample is multiplied by its similarity metric value, and so on.

[0149] The weighted interface fit parameters of all standard samples are summed, and then divided by the sum of the weight coefficients to obtain the weighted average interface fit parameters. This weighted averaging process allows for a more accurate synthesis of the interface fit parameters of the matching sample set.

[0150] Step S434: Correct the weighted average interface bite parameters using a regression model to eliminate the influence of shape similarity matching error on the parameters and generate corrected interface bite parameters. The regression model is trained based on historical matching data.

[0151] When correcting the weighted average interface bite parameters using a regression model, the regression model trained on historical matching data is first employed. This historical matching data is accumulated during previous matching of shape feature vectors with the database, including shape feature vectors, standard matching samples, and interface bite parameters. The parameters of the regression model are obtained through analysis and training of this historical data. The weighted average interface bite parameters are then used as input to the regression model, and the corrected interface bite parameters are calculated based on these parameters. The regression model can learn the relationship between shape similarity matching error and interface bite parameters from historical data. This relationship is used to adjust the weighted average interface bite parameters, eliminating the influence of errors. For example, the regression model might discover from historical data that when there is a certain error in shape similarity matching, the interface bite parameters will have a corresponding deviation. Through the calculations of the regression model, this deviation can be corrected to obtain more accurate corrected interface bite parameters.

[0152] Step S435: Divide the corrected interface occlusion parameters into preset categories, namely contact area parameters, occlusion depth parameters and interface strength parameters, and extract the contact area parameters and occlusion depth parameters as interface occlusion parameters related to shape features.

[0153] When classifying the corrected interfacial interlocking parameters according to preset categories, they are categorized into contact area parameters, interlocking depth parameters, and interfacial strength parameters based on their definitions and meanings. For example, parameters representing the actual size or proportion of the contact area are categorized as contact area parameters; parameters describing the average value and standard deviation of the interlocking depth are categorized as interlocking depth parameters; and parameters representing the interfacial bonding strength are categorized as interfacial strength parameters. Then, contact area parameters and interlocking depth parameters are extracted from the categorized groups. These parameters are shape-related interfacial interlocking parameters, which can more directly reflect the influence of aggregate particle shape on its interlocking ability with the slurry. These extracted parameters provide crucial input for subsequent weighted fusion of interfacial interlocking parameters, calculation of the proportion of interfacial contact area, and generation of interlocking depth distribution parameters.

[0154] Step S440: The extracted interface interlocking parameters are weighted and fused. The weighting coefficients are adjusted according to the concrete strength requirements, mix proportions and aggregate types to generate a weighted set of interface interlocking parameters.

[0155] Weighted fusion is a method that comprehensively processes multiple parameters. By assigning different weight coefficients to each parameter and then summing them in a weighted manner, a comprehensive parameter set is obtained. In this step, the weight coefficients are adjusted according to concrete strength requirements, mix proportions, and aggregate types because different concrete strength requirements, mix proportions, and aggregate types have varying impacts on the importance of interfacial bonding parameters. The weighted interfacial bonding parameter set is the result of weighted fusion, comprehensively considering the importance of contact area and bonding depth parameters under different conditions.

[0156] Step S450: Calculate the interface contact area ratio based on the weighted set of interface interlocking parameters. The interface contact area ratio is the ratio of the actual contact area between the aggregate and the slurry to the surface area of ​​the aggregate, which is derived from the surface area-related parameters in the shape feature vector.

[0157] When calculating the proportion of interfacial contact area based on a weighted set of interfacial interlocking parameters, parameters related to surface area are first extracted from the shape feature vector. These parameters may include information such as the total surface area of ​​aggregate particles and the surface area of ​​different parts. Then, based on the contact area parameters in the weighted set of interfacial interlocking parameters, combined with surface area-related parameters, the actual contact area between the aggregate and the slurry is calculated. For example, the contact area parameters may contain specific numerical values ​​of the actual contact area or proportional information related to the contact area. Using this information and the surface area parameters, the actual contact area can be accurately calculated. Finally, the calculated actual contact area is divided by the surface area of ​​the aggregate to obtain the proportion of interfacial contact area.

[0158] Step S460: Generate bite depth distribution parameters based on the weighted set of interface bite parameters. The parameters include the central tendency, dispersion and extreme values ​​of bite depth. The spatial distribution reflecting the bite effect is obtained through statistical analysis of the surface features of the particles.

[0159] When generating occlusal depth distribution parameters based on the weighted set of interfacial occlusal parameters, parameters related to occlusal depth are first extracted from the weighted set. These parameters may include occlusal depth measurements at different locations, proportional information about occlusal depth, etc. Then, statistical analysis is performed on these occlusal depth data. For calculating central tendency, the mean can be calculated by summing all occlusal depth measurements and dividing by the number of data points; alternatively, the median can be calculated by arranging the data in ascending order and taking the middle value. For calculating dispersion, the standard deviation can be calculated, reflecting the degree of dispersion of the data relative to the mean; alternatively, the variance can be calculated, which is the square of the standard deviation. For extreme values, the maximum and minimum values ​​in the data are directly identified.

[0160] Simultaneously, the analysis incorporates the surface irregularities of the particles. The unevenness of the particle surface leads to variations in bite depth at different locations. Statistical analysis of these surface irregularities allows for a better understanding of the reasons behind the bite depth distribution. For example, areas with greater surface irregularities may have larger and more dispersed bite depths; conversely, areas with smoother surfaces may have more concentrated bite depths.

[0161] Step S500: Based on the quantitative characterization of the entire batch of aggregate particles, perform spatial density weighted fusion, calculate the density weight factor based on the relative positional relationship between particles, normalize the quantitative characterization, and generate a structural performance index that reflects the functional characteristics of aggregate particle shape through weighted fusion.

[0162] The quantitative characterization of the entire batch of aggregate particles is a quantitative description of the interlocking ability of each aggregate particle at the aggregate-paste interface, obtained through the previous steps. This includes the proportion of the interface contact area and parameters such as the distribution of interlocking depth. Spatial density-weighted fusion is a method that comprehensively considers the spatial distribution and quantitative characterization of aggregate particles. By calculating density weighting factors, different weights are assigned to particles at different locations, thereby more accurately evaluating the performance of the entire batch of aggregate. The density weighting factors are calculated based on the relative positional relationships between particles, reflecting the relative importance of each particle in the entire batch of aggregate. Normalization is a process that unifies the quantitatively characterized data, giving it the same dimensions and range, facilitating weighted fusion. The structural performance index is a comprehensive indicator generated after spatial density-weighted fusion, reflecting the contribution of the functional characteristics of aggregate particle shape to concrete performance.

[0163] Optionally, step S500 can be specifically implemented as the following steps S510-S570:

[0164] Step S510: Obtain the quantitative characterization of the entire batch of aggregate particles and establish a particle quantitative characterization set, which includes the proportion of interface contact area and the interlocking depth distribution parameters of each particle.

[0165] When obtaining the quantitative characterization of the entire batch of aggregate particles, the interfacial contact area ratio and interlocking depth distribution parameters of each particle are extracted from the previous calculation results. The interfacial contact area ratio reflects the ratio of the actual contact area between the aggregate and the paste to the aggregate surface area; the interlocking depth distribution parameters include information such as the central tendency, dispersion, and extreme values ​​of the interlocking depth. These quantitative characterization information for each particle are recorded in a certain order (such as the particle numbering order) to form a particle quantitative characterization set.

[0166] Step S520: Perform spatial position modeling for the entire batch of aggregate particles, determine the relative position coordinates of each particle based on the coordinate information in the dynamic falling image stream, and generate a particle spatial position matrix. The matrix records the planar coordinates of each particle.

[0167] When modeling the spatial position of the entire batch of aggregate particles, the coordinate information of each particle is first extracted from the dynamic falling image stream. Image processing and target tracking algorithms can be used to identify and track the particles in the images, recording the coordinate position of each particle in different frames. Then, the relative position coordinates of each particle are determined based on this coordinate information. Alternatively, a reference particle can be selected, and the position coordinates of other particles relative to that reference particle can be calculated, or the relative position coordinates can be determined by calculating the centroid position of all particles and using the centroid as a reference point.

[0168] Arrange the relative position coordinates of each particle in order to form a particle spatial position matrix. For example, for a batch of aggregate with n aggregate particles, the particle spatial position matrix is ​​an n×2 matrix, where the two elements in the i-th row represent the x-coordinate and y-coordinate of the i-th particle, respectively. This matrix records the planar coordinate information of the entire batch of aggregate particles.

[0169] Step S530: Calculate the relative distance between particles. Calculate the spatial distance between any two particles based on the particle spatial position matrix to generate a particle distance matrix. The particle distance matrix has a symmetric matrix structure.

[0170] When calculating the relative distance between particles, the planar coordinates of each particle are obtained from the particle spatial position matrix. For any two particles i and j, their planar coordinates are (xi, yi) and (xj, yj) respectively, and the spatial distance between them can be calculated using the Euclidean distance formula. The calculated spatial distances between any two particles are then filled into the corresponding positions in the particle distance matrix. For example, the element in the i-th row and j-th column of the matrix represents the spatial distance from particle i to particle j. Due to the symmetry of the distance, the element in the j-th row and i-th column of the matrix is ​​equal to the element in the i-th row and j-th column, so only the elements of the upper or lower triangular part of the matrix need to be calculated, and then the result is symmetrically copied to the other half. The final generated particle distance matrix can intuitively display the relative distance relationships between the aggregate particles in the entire batch.

[0171] Step S540: Calculate the density weight factor based on the particle distance matrix. Particles with smaller distances are assigned a larger weight factor. Use a kernel function to convert the distance values ​​into weight values ​​and generate a particle weight vector. The weight vector represents the relative importance of particles in the population.

[0172] The density weighting factor is used to measure the relative importance of each particle in the entire aggregate batch. Particles with smaller distances have greater mutual influence and are therefore assigned a larger weighting factor. The kernel function is a function that converts distance values ​​into weight values. It maps distance values ​​in the particle distance matrix to weight values, thus generating a particle weight vector. Each element in the particle weight vector corresponds to the weight value of an aggregate particle, reflecting its relative importance within the group.

[0173] Optionally, step S540 can be implemented as the following steps S541-S546:

[0174] Step S541: Perform row normalization on the particle distance matrix by dividing the distance value of each row by the maximum distance value of that row to generate a normalized distance matrix.

[0175] Row normalization is a method of uniformly processing each row of data in a matrix. By dividing the distance value of each row by the maximum distance value of that row, the range of each row of data can be unified to the interval [0,1], which facilitates subsequent calculations and comparisons. The normalized distance matrix is ​​a matrix obtained after row normalization, with its element values ​​between [0,1], reflecting the relative distance relationship between each particle and other particles.

[0176] Step S542: Obtain the bandwidth parameter of the kernel function. The bandwidth parameter is adjusted according to the distribution characteristics of the aggregate particles in the whole batch to adapt to particle groups with different densities.

[0177] When obtaining the bandwidth parameter of the kernel function, first analyze the distribution characteristics of the entire batch of aggregate particles. Distribution characteristics can be evaluated by calculating some statistics in the particle distance matrix, such as the average distance and the standard deviation of the distance. If the average distance and standard deviation are both small, it indicates that the particle group is relatively dense; if the average distance and standard deviation are both large, it indicates that the particle group is relatively sparse. Based on the distribution characteristics obtained from the analysis, adjust the bandwidth parameter of the kernel function. An empirical formula or experimental methods can be used to determine a suitable bandwidth parameter. For example, for the Gaussian kernel function, a suitable σ value can be selected as the bandwidth parameter based on the statistical information of particle distance. By adjusting the bandwidth parameter, the kernel function can better adapt to the distribution of the entire batch of aggregate particles, reasonably converting distance values ​​into weight values, and accurately reflecting the relative importance of particles in the group.

[0178] Step S543: Apply a kernel function to each element in the normalized distance matrix to calculate the similarity weights between particles and generate a similarity weight matrix. The element values ​​in the matrix represent the degree of mutual influence between two particles.

[0179] Similarity weights are weight values ​​calculated using a kernel function based on the distance between particles, reflecting the degree of mutual influence between two particles. The similarity weight matrix is ​​a matrix where each element represents the similarity weight between any two particles. This matrix visually displays the relationships between particles in the entire batch of aggregate. When applying a kernel function to each element in the normalized distance matrix, a suitable kernel function should be chosen, such as the Gaussian kernel function or the exponential kernel function. Taking the Gaussian kernel function as an example, for the element d in the i-th row and j-th column of the normalized distance matrix... ij * The distance value after row normalization is substituted into the Gaussian kernel function to calculate the similarity weight between particle i and particle j. The calculated similarity weights are then filled into the corresponding positions in the matrix to form the similarity weight matrix. The element in the i-th row and j-th column of this matrix represents the degree of mutual influence between particle i and particle j. Due to the symmetry of the distance, the similarity weight matrix is ​​also a symmetric matrix.

[0180] Step S544: Perform row summation on the similarity weight matrix, calculate the total influence weight of each particle, and generate a particle influence weight vector. The element value in the particle influence weight vector is the total influence weight of the corresponding particle.

[0181] By performing row summation on the similarity weight matrix, the sum of the similarity weights between each particle and all other particles can be calculated, which is the total influence weight of each particle. The particle influence weight vector is a vector in which each element corresponds to the total influence weight of an aggregate particle, reflecting the overall influence of that particle in the entire batch of aggregate.

[0182] Step S545: Normalize the particle influence weight vector so that the sum of the weight values ​​of all particles is a uniform value, and generate the density weight factor vector.

[0183] Normalization is a method of unifying the elements in a vector. By normalizing the particle influence weight vector, the sum of the weight values ​​of all particles becomes a uniform value (usually 1). This eliminates differences in the range of weight values ​​among different particles, making weighted fusion operations more convenient. The density weight factor vector is the vector obtained after normalization. Each element corresponds to the density weight factor of an aggregate particle, reflecting the relative importance of that particle in the entire batch of aggregates, and the sum of all elements is 1.

[0184] Step S546: Smooth the particle weight vector to generate the final particle weight vector. The window size for smoothing is determined based on the uniformity of particle distribution to balance local and global weight characteristics.

[0185] When smoothing the particle weight vector, a smoothing method is selected, such as moving average or Gaussian smoothing. Taking the moving average method as an example, the window size for smoothing is first determined. Based on the uniformity of particle distribution, the uniformity of particle distribution is evaluated by analyzing information in the particle distance matrix or particle spatial location matrix. If the particle distribution is relatively uniform, the window size can be set to a smaller value, such as 3 or 5; if the particle distribution is uneven, the window size can be set to a larger value, such as 7 or 9. For each element in the particle weight vector, the average value of elements within the window size is calculated with that element as the center, and the average value is used as the smoothed weight value of the i-th element. The smoothed weight values ​​of all elements are arranged in order to form the final particle weight vector.

[0186] Step S550: The weighted summation of the interface contact area ratios in the particle quantization characterization set is performed, and the overall interface contact area ratio is generated using the particle weight vector as the weight coefficient.

[0187] When performing weighted summation, the interfacial contact area ratio of each particle is sequentially extracted from the particle quantification characterization set, and the corresponding weight coefficient is extracted from the particle weight vector. The interfacial contact area ratio of each particle is multiplied by its corresponding weight coefficient to obtain the particle's contribution to the overall interfacial contact area ratio. Then, the contribution values ​​of all particles are summed to obtain the overall interfacial contact area ratio. This indicator comprehensively considers the situation of each particle in the entire batch of aggregate, taking into account the relative positions and mutual influences between particles. Compared to the interfacial contact area ratio of a single particle, it better reflects the performance of the entire batch of aggregate in practical applications. For example, in concrete, a larger overall interfacial contact area ratio means that the bond between the aggregate and the paste is likely to be tighter, and the concrete performance is likely to be better.

[0188] Step S560: Perform statistical analysis on the occlusal depth distribution parameters in the particle quantification characterization set, calculate the weighted average occlusal depth, weighted occlusal depth dispersion and weighted maximum occlusal depth, normalize the overall interface contact area ratio and weighted occlusal depth parameters, and generate normalized interface contact area index and normalized occlusal depth index.

[0189] In statistical analysis, the interlocking depth distribution parameter for each particle is calculated in conjunction with its weighting coefficient in the particle weight vector. For example, when calculating the weighted average interlocking depth, the average interlocking depth of each particle is multiplied by its corresponding weighting coefficient, and then all results are summed. A similar weighted calculation method is used for calculating the weighted interlocking depth dispersion and the weighted maximum interlocking depth. The overall interface contact area ratio and weighted interlocking depth parameters are normalized because these parameters may have different value ranges; normalization unifies them to a common scale, facilitating subsequent weighted fusion. Through normalization, normalized interface contact area indices and normalized interlocking depth indices are generated. These normalized indices eliminate the influence of different parameter value ranges, more accurately reflecting the relative performance of the entire batch of aggregate in terms of interface contact and interlocking depth.

[0190] Step S570: The normalized interface contact area index and the normalized interlocking depth index are weighted and fused to generate the structural performance index. The structural performance index comprehensively represents the contribution of aggregate particle shape to concrete performance.

[0191] When performing weighted fusion, the weighting coefficients of the normalized interfacial contact area index and the normalized interlocking depth index are first determined based on specific engineering needs and concrete design requirements. For example, for projects with high requirements for the integrity and bond strength of concrete, a higher weight may be assigned to the normalized interfacial contact area index; for projects requiring high shear strength, the normalized interlocking depth index may be given more weight. Then, the normalized interfacial contact area index is multiplied by its corresponding weighting coefficient, and the normalized interlocking depth index is multiplied by its corresponding weighting coefficient. The two results are then added together to obtain the structural performance index. The structural performance index is a comprehensive indicator that integrates the two important factors of interfacial contact area and interlocking depth, and can comprehensively characterize the contribution of aggregate particle shape to concrete performance.

[0192] The various algorithms involved in the above descriptions of the embodiments of this invention, such as the Euclidean distance algorithm, the cosine distance algorithm, the Gaussian kernel function, etc., can all be obtained from relevant content in the prior art, and will not be elaborated on in the embodiments of this invention. Furthermore, those skilled in the art can supplement the details based on common knowledge in the field when implementing the solutions of this invention. For example, based on common knowledge in the field, they can use normalization to eliminate dimensional conflicts before feature fusion, use interpolation to eliminate dimensional differences, reasonably set thresholds based on historical data, experience, or business scenario requirements, train the model based on a general model training method, set the number of layers in the model structure based on actual needs, select activation functions, etc. This invention will not provide redundant descriptions of overly detailed implementation processes here.

[0193] Figure 2 This is a schematic diagram of the composition structure of an aggregate particle shape and gradation identification device provided in an embodiment of the present invention, as shown below. Figure 2 As shown, the aggregate particle shape and gradation identification device 200 includes:

[0194] Image acquisition module 210 is used to acquire a dynamic falling image stream of aggregate particle group, and separate a closed contour image of a single aggregate particle from the dynamic falling image stream. The closed contour image contains a set of continuous pixel coordinates of the particle boundary.

[0195] The harmonic field construction module 220 is used to construct a boundary-constrained harmonic field for the closed contour image, and generate a harmonic field distribution that satisfies the contour boundary gradient condition in the internal region of the particle. The harmonic field distribution is characterized by the field intensity value change of each point inside the particle through a two-dimensional spatial field intensity value matrix.

[0196] The critical encoding module 230 is used to extract stable critical structures based on the harmonic field distribution, identify critical inflection points of field strength value changes, and encode the spatial coordinate sequence of the critical inflection points into a topologically invariant sequence.

[0197] The semantic mapping module 240 is used to call the shape semantic mapper calibrated by the concrete interface performance to perform semantic mapping processing on the topological invariant sequence and output a quantitative characterization of the interlocking ability of a single aggregate particle at the aggregate-paste interface. The quantitative characterization includes the interface contact area ratio and interlocking depth distribution parameters.

[0198] The weighted fusion module 250 is used to perform spatial density weighted fusion based on the quantitative characterization of the whole batch of aggregate particles, calculate the density weight factor based on the relative positional relationship between particles, normalize the quantitative characterization, and generate a structural performance index that reflects the functional characteristics of aggregate particle shape through weighted fusion.

[0199] The descriptions of the apparatus embodiments above are similar to those of the method embodiments above, and have similar beneficial effects. In some embodiments, the functions or modules included in the apparatus provided by the present invention can be used to perform the methods described in the method embodiments above. For technical details not disclosed in the apparatus embodiments of the present invention, please refer to the descriptions of the method embodiments of the present invention for understanding.

[0200] Figure 3 A hardware entity diagram of a computer system provided as an embodiment of the present invention, such as... Figure 3 As shown, the hardware entity of the computer system 1000 includes a processor 1001 and a memory 1002, wherein the memory 1002 stores a computer program that can run on the processor 1001, and the processor 1001 executes the program to implement the steps in the method of any of the above embodiments.

[0201] It should be noted that the descriptions of the storage medium and device embodiments above are similar to those of the method embodiments above, and have similar beneficial effects. For technical details not disclosed in the storage medium and device embodiments of the present invention, please refer to the descriptions of the method embodiments of the present invention for understanding.

Claims

1. An image processing-based aggregate grain shape and gradation identification method, characterized by, The method comprises: acquiring a dynamic falling image stream of the aggregate particle group, separating a closed contour image of a single aggregate particle from the dynamic falling image stream, the closed contour image containing a continuous pixel coordinate set of a particle boundary; boundary-constrained harmonic field construction is performed on the closed contour image to generate a harmonic field distribution satisfying the gradient condition of the contour boundary in the internal region of the particle, the harmonic field distribution represents the change of field strength value of each point in the internal region of the particle through a two-dimensional space field strength value matrix; specifically, the boundary pixel coordinate set of the closed contour image is converted into a continuous curve model, a smooth particle boundary model is generated through curve fitting; the normal vector direction of each point on the boundary is calculated based on the particle boundary model, the field strength gradient constraint condition at the boundary is determined, a boundary constraint relationship group of the harmonic field is established, the boundary constraint relationship group contains the field strength value and the normal derivative constraint of each point on the boundary; a two-dimensional orthogonal grid is constructed in the internal region of the particle, the pixel coordinate system of the closed contour image is mapped to a standardized grid coordinate system to generate a grid node coordinate set, the spacing of the grid nodes is dynamically adjusted according to the size of the particle to adapt to particles of different sizes; the boundary constraint relationship group is solved by using the finite difference method, the field strength value distribution satisfying the Laplace relationship is calculated on the grid node coordinate set to generate a preliminary harmonic field, the field strength value range of the preliminary harmonic field is processed by normalization to unify the dimension; boundary consistency verification is performed on the preliminary harmonic field, the deviation amount of the actual value of the field strength gradient at the boundary from the theoretical constraint value is calculated, the field strength value of the internal node is corrected by the relaxation iteration method until the deviation amount meets the accuracy requirement; the corrected grid node field strength value is converted into a two-dimensional space field strength value matrix in the original image coordinate system to generate a harmonic field distribution satisfying the gradient condition of the contour boundary, the number of rows and columns of the two-dimensional space field strength value matrix is consistent with the size of the closed contour image; stable critical structure extraction is performed based on the harmonic field distribution, a critical turning point of the field strength value change is identified, and a spatial coordinate sequence of the critical turning point is encoded into a topologically invariant sequence; a shape semantic mapper calibrated by the concrete interface performance is called to perform semantic mapping processing on the topologically invariant sequence, and a quantitative representation of the occlusion ability of a single aggregate particle to the aggregate-paste interface is output, the quantitative representation contains an interface contact area ratio and a occlusion depth distribution parameter; spatial density weighted fusion is performed on the quantitative representation of the whole batch of aggregate particles, a density weight factor is calculated based on the relative positional relationship between the particles, the quantitative representation is normalized, and a structure performance index reflecting the functional characteristics of the aggregate particle shape is generated through weighted fusion.

2. The method of claim 1, wherein, The acquiring a dynamic falling image stream of the aggregate particle group, separating a closed contour image of a single aggregate particle from the dynamic falling image stream, the closed contour image containing a continuous pixel coordinate set of a particle boundary, comprises: multi-light-source shadow elimination processing is performed on the dynamic falling image stream, a particle surface highlight reflection area is separated and eliminated through three primary color channels to generate a shadow-suppressed aggregate image sequence; The granular motion region binary image is generated by performing inter-frame difference operation on the aggregate image sequence, extracting a motion region mask, and removing noise holes in the mask through morphological opening and closing operation, wherein a foreground pixel value in the granular motion region binary image represents a granular occupied region; The granular motion region binary image is subjected to connected granule segmentation processing by using a watershed algorithm, a foreground marker point set is generated based on distance transformation, and granule boundary separation is achieved through marker-controlled watershed transformation, thereby obtaining an initially separated single granule region; Boundary tracking processing is performed on the separated single granule region, a neighborhood connected domain analysis algorithm is used to extract an outer contour pixel coordinate of the granule region, an initial boundary coordinate sequence is generated, and the initial boundary coordinate sequence is arranged in a clockwise direction; Redundant point elimination processing is performed on the initial boundary coordinate sequence, and boundary feature salient points are retained, thereby generating a simplified continuous pixel coordinate set; A closed polygon is constructed based on the continuous pixel coordinate set, a boundary pixel point interpolation completion algorithm is used to process a contour broken region, and a closed contour image satisfying topological closure is generated.

3. The method of claim 1, wherein, The normal vector direction of each point on the granule boundary model is calculated based on the granule boundary model, a field strength gradient constraint condition at the boundary is determined, and a boundary constraint relationship group of the harmonic field is established, including: First derivative calculation is performed on the curve relationship of the granule boundary model, a tangent vector direction parameter of each point on the boundary is obtained, and the normal vector direction parameter is generated through orthogonal vector transformation, wherein the normal vector direction points to the interior region of the granule to ensure the interior directivity of the field strength distribution; The field strength gradient direction is determined according to the normal vector direction parameter, a field strength gradient amplitude reference value at the boundary is set, a boundary gradient constraint vector is generated, and the gradient amplitude reference value is positively correlated with the curvature of the granule boundary to reflect the influence of the boundary bending degree on the field strength; The interior region of the closed contour image is divided into a finite element calculation grid, the grid elements are numbered, and a node and element association table is established, which records the node index contained by each grid element; For each boundary grid node, a partial differential relationship between the field strength value and the field strength values of adjacent nodes is established according to the boundary gradient constraint vector, a node field strength constraint relationship is generated, and the node field strength constraint relationship includes a linear combination of the node's own field strength value and the field strength values of adjacent nodes; The node field strength constraint relationships of all boundary nodes are collected to construct a linear relationship group containing boundary conditions, and a row pivot elimination method is used to initialize the relationship group coefficient matrix; The boundary condition is strengthened by performing boundary condition strengthening processing on the linear relationship group, the gradient constraint condition is converted into an additional relationship through a penalty function method, and a complete boundary constraint relationship group is generated.

4. The method of claim 3, wherein, The boundary constraint relationship group is solved by using a finite difference method, the field strength value distribution satisfying the Laplace relationship is calculated on the grid node coordinate set, and a preliminary harmonic field is generated, including: The grid node coordinate set is subjected to classification processing to distinguish boundary nodes from internal nodes, and a node type marker array is established, wherein the boundary nodes are marked as special identifiers to distinguish from internal nodes. For the internal nodes, a field strength value calculation relationship is established according to a difference format of Laplace relationship, and the field strength value calculation relationship represents that the field strength value of the internal node is equal to the average value of the field strength values of four neighboring nodes of the internal node; The boundary constraint relationship group and the internal node field strength value calculation relationship are combined to obtain a simultaneous relationship group containing all nodes; An iterative method is used to solve the simultaneous relationship group, the field strength values of the internal nodes are initialized as a uniform distribution, the field strength values of the boundary nodes are assigned according to the boundary constraint condition, and the field strength values of the nodes are updated through iteration to gradually approach the true solution; The convergence of the field strength values in the iteration process is monitored, the overall deviation of the field strength values of two adjacent iterations is calculated, the iteration is stopped when the overall deviation is less than a convergence judgment standard, and a converged node field strength value set is obtained; The node field strength value set is arranged in order of the grid node coordinates to generate a two-dimensional field strength value matrix, and the two-dimensional field strength value matrix is determined as a preliminary harmonic field, and each element of the two-dimensional field strength value matrix corresponds to the field strength value of a grid node.

5. The method of claim 1, wherein, The stable critical structure is extracted based on the harmonic field distribution, and a critical turning point of the field strength value change is identified, and a spatial coordinate sequence of the critical turning point is encoded as a topological invariant sequence, including: Gradient calculation is performed on the two-dimensional spatial field strength value matrix of the harmonic field distribution to generate a field strength gradient modulus value matrix and a gradient direction matrix, the field strength gradient modulus value matrix represents a field strength value change rate, and the gradient direction matrix represents a change direction to reflect a field strength spatial change trend; A dynamic threshold is set in the gradient modulus value matrix, a gradient modulus value segmentation threshold is calculated by an adaptive threshold determination method, a region with a gradient modulus value greater than the segmentation threshold is extracted as a field strength change significant region, and the field strength change significant region reflects a position of a field strength change. Region growing processing is performed on the field strength change significant region, starting from a pixel point with the maximum gradient direction change rate, adjacent pixels are merged according to a gradient direction consistency criterion to generate a field strength change connected region, and the connected region is a continuous region of the field strength change. In each field strength change connected region, a field strength value change trajectory is tracked along the gradient direction, local maximum points and minimum points of the field strength value on the trajectory are recorded, a critical turning point of the field strength value change is determined, and the critical turning point reflects a change position of the field strength change trend. Spatial coordinates of the critical turning points are extracted, and a spatial coordinate sequence of the critical turning points is generated in the order of the trajectory tracking, and the spatial coordinate sequence contains plane coordinate information of the turning points to reflect spatial positions of the turning points. The spatial coordinate sequence is topologically encoded, and absolute coordinates are converted into relative coordinates to generate a topological invariant sequence.

6. The method of claim 5, wherein, The two-dimensional spatial field strength value matrix of the harmonic field distribution is subjected to gradient calculation to generate a field strength gradient modulus value matrix and a gradient direction matrix, the field strength gradient modulus value matrix represents a field strength value change rate, and the gradient direction matrix represents a change direction, including: Convolution operation is performed on the two-dimensional spatial field strength value matrix of the harmonic field distribution to calculate a first-order partial derivative matrix in a horizontal direction and a first-order partial derivative matrix in a vertical direction, and the first-order partial derivative matrix in the horizontal direction represents a horizontal field strength change rate. The horizontal partial derivative matrix and the vertical partial derivative matrix are subjected to square sum and square root operations to generate a field strength gradient modulus value matrix, each element of the gradient modulus value matrix being the modulus length of the gradient vector at the corresponding position to quantify the change rate; The gradient direction angle is calculated by a direction angle calculation method, the element ratio of the horizontal partial derivative matrix and the vertical partial derivative matrix is converted into an angle value to generate a gradient direction matrix, and the angle value reflects the direction of the change in the field strength; The gradient direction matrix is subjected to non-maximum suppression processing to retain local maximum points in the gradient direction, thereby enhancing the edge positioning accuracy and generating a refined gradient direction matrix; The field strength gradient modulus value matrix is subjected to smoothing processing to generate a smoothed gradient modulus value matrix, and the smoothing degree is dynamically adjusted according to the gradient noise level; The smoothed gradient modulus value matrix and the refined gradient direction matrix are subjected to dimension alignment to ensure that the number of rows and columns of the two matrices is the same, and the corresponding pixel points are one-to-one corresponding; The method further comprises the following steps: In each of the field strength change connected regions, the gradient modulus values of all the pixel points are calculated, and the pixel point with the maximum gradient modulus value is selected as a starting point of the trajectory, and the starting point is the position with the most severe change in the field strength; Starting from the starting point of the trajectory, the trajectory is moved to an adjacent pixel point in the direction indicated by the gradient direction matrix, and the coordinates of the pixel points on the moving path and the corresponding field strength values are recorded to generate an initial field strength change trajectory, and the trajectory reflects the change process of the field strength value along the gradient direction; The initial field strength change trajectory is subjected to smoothing processing to eliminate trajectory jitter, and a smoothed field strength change trajectory is generated, and the window size of the smoothing processing is determined according to the fluctuation degree of the trajectory; Along the smoothed field strength change trajectory, the change amount of the field strength value of the adjacent pixel points is calculated, and when the change amount changes from positive to negative, a local maximum point is marked, and when the change amount changes from negative to positive, a local minimum point is marked; The marked local maximum points and minimum points are screened to remove adjacent extreme points with a distance less than a distance threshold to generate a candidate critical turning point set; The second-order derivative of the field strength value of each point in the candidate critical turning point set is calculated, and when the absolute value of the second-order derivative exceeds a curvature judgment standard, the point is confirmed as a critical turning point of the change in the field strength value, and the second-order derivative is calculated by a difference method.

7. The method of claim 1, wherein, The shape semantic mapper calibrated by the concrete interface performance is called to perform semantic mapping processing on the topologically invariant sequence, and a quantitative representation of the occlusion ability of a single aggregate particle to the aggregate-paste interface is output, and the quantitative representation includes an interface contact area ratio and an occlusion depth distribution parameter, and the method comprises the following steps: The topologically invariant sequence is input into the input layer of the shape semantic mapper calibrated by the concrete interface performance, and different lengths of the topologically invariant sequence are converted into fixed-dimension input vectors through sequence length normalization processing. The shape semantic mapper extracts features of the input vector based on a hidden layer and generates a shape feature vector through a multi-layer network structure, the shape feature vector containing corner features, surface texture features and overall morphology features of the particle; The concrete interface performance correlation layer of the shape semantic mapper is called to match the shape feature vector with a pre-stored concrete interface performance database to extract interface engagement parameters related to the shape features, the concrete interface performance database containing interface performance experimental data corresponding to different aggregate shapes; The extracted interface engagement parameters are weighted and fused, and weight coefficients are adjusted according to concrete strength requirements, mixing ratio parameters and aggregate types to generate a weighted interface engagement parameter set; The interface contact area ratio is calculated based on the weighted interface engagement parameter set, the interface contact area ratio being a ratio of an actual contact area of the aggregate and the paste to a surface area of the aggregate, and being derived from surface area-related parameters in the shape feature vector; The engagement depth distribution parameter is generated according to the weighted interface engagement parameter set, the parameter containing a central tendency, a dispersion degree and an extreme value of the engagement depth, and being obtained through statistical analysis of the concave-convex features of the particle surface to reflect the spatial distribution of the engagement effect.

8. An apparatus for identifying the shape and gradation of aggregate particles, characterized by The image acquisition module is configured to acquire a dynamic falling image stream of the aggregate particle group, separate a closed contour image of a single aggregate particle from the dynamic falling image stream, and the closed contour image contains a continuous pixel coordinate set of a particle boundary. ​ The harmonic field construction module is configured to perform boundary-constrained harmonic field construction on the closed contour image, generate a harmonic field distribution satisfying the gradient condition of the contour boundary in the internal region of the particle, and represent the change of field strength value of each point in the internal region of the particle by a two-dimensional spatial field strength value matrix; specifically, the boundary pixel coordinate set of the closed contour image is converted into a continuous curve model, and a smooth particle boundary model is generated by curve fitting; the normal vector direction of each point on the boundary is calculated based on the particle boundary model, the field strength gradient constraint condition at the boundary is determined, the boundary constraint relation set of the harmonic field is established, and the boundary constraint relation set includes the field strength value and the normal derivative constraint of each point on the boundary; a two-dimensional orthogonal grid is constructed in the internal region of the particle, the pixel coordinate system of the closed contour image is mapped to a standardized grid coordinate system, a grid node coordinate set is generated, and the spacing of the grid nodes is dynamically adjusted according to the size of the particle to adapt to particles of different sizes; the boundary constraint relation set is solved by using the finite difference method, the field strength value distribution satisfying the Laplace relation is calculated on the grid node coordinate set, a preliminary harmonic field is generated, and the field strength value range of the preliminary harmonic field is processed by normalization to unify the dimension; the boundary consistency of the preliminary harmonic field is verified, the deviation amount of the actual value of the field strength gradient at the boundary from the theoretical constraint value is calculated, the field strength value of the internal node is corrected by the relaxation iteration method until the deviation amount meets the accuracy requirement; the corrected grid node field strength value is converted into a two-dimensional spatial field strength value matrix in the original image coordinate system, and the harmonic field distribution satisfying the contour boundary gradient condition is generated, and the number of rows and columns of the two-dimensional spatial field strength value matrix is consistent with the size of the closed contour image. The critical coding module is configured to perform stable critical structure extraction based on the harmonic field distribution, identify a critical turning point of the field strength value change, and encode the spatial coordinate sequence of the critical turning point into a topologically invariant sequence. The semantic mapping module is configured to call a shape semantic mapper calibrated by a concrete interface performance, perform semantic mapping processing on the topologically invariant sequence, and output a quantitative representation of the occlusion ability of a single aggregate particle to the aggregate-paste interface, wherein the quantitative representation includes an interface contact area proportion and an occlusion depth distribution parameter. The weighted fusion module is configured to perform spatial density weighted fusion according to the quantitative representation of the entire batch of aggregate particles, calculate a density weight factor based on the relative positional relationship between the particles, perform normalization processing on the quantitative representation, and generate a structure performance index reflecting the functional characteristics of the aggregate particle shape by weighted fusion.

9. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by a processor to implement the steps in the method of any one of claims 1 to 7.

Images