Non-contact long-distance monitoring method and system for structural settlement

Abstract

The present invention proposes a non-contact long-distance monitoring method and system for structural settlement, which includes obtaining an original image set and performing preprocessing to obtain a preprocessed image set; performing feature point recognition and feature point tracking to obtain a feature point vertical trajectory set; constructing a topological relationship between feature points to obtain a time-varying feature graph; calculating a graph structure change index and a local vertical deformation field, constructing and solving a global settlement field equation to obtain a settlement field function; calculating the predicted settlement field at a future moment; quantifying the uncertainty of the settlement estimate to obtain an uncertainty field; performing multi-scale error decomposition to obtain an error scale coefficient set; constructing and applying an adaptive Kalman filter to obtain an optimized settlement field function; and using a sub-pixel interpolation algorithm to calculate a high-precision settlement field function. The present invention not only realizes non-contact, large-scale, and long-distance structural settlement monitoring, but also improves the resolution and accuracy of the settlement monitoring results.

Description

Technical Field

[0001] The invention belongs to the field of engineering structure settlement measurement, in particular to a non-contact long-distance monitoring method and system for structure settlement. Background Art

[0002] There is a great demand for non-contact, long-distance, high-precision structural settlement monitoring (monitoring) technology. For example, in the construction of bridges across waters (road bridges, high-speed railway bridges, etc.), there will be a "window period" for pier settlement monitoring before the upper bridge is built. During this "window period", close-range observation / monitoring schemes lack a supporting platform and are inconvenient to implement; if traditional methods such as GPS are used, the accuracy has great limitations, and there is an urgent need to implement long-distance settlement observation technology in a targeted manner; in addition, non-contact, long-distance settlement monitoring methods have good adaptability and convenience in more application scenarios.

[0003] Currently, Digital Image Correlation (DIC) technology has been widely used in non-contact monitoring and detection. Its basic principle is to obtain deformation information by comparing two digital images of an object before and after deformation using correlation calculation. However, traditional DIC technology applications are mostly short-distance monitoring at laboratory scale. How to achieve larger-scale and long-distance DIC technology applications is a technical problem that needs to be solved urgently. Summary of the invention

[0004] The purpose of the invention is to propose a non-contact long-distance monitoring method and system for structural settlement to solve the above-mentioned problems existing in the prior art.

[0005] The technical solution, a non-contact remote monitoring method for structural settlement, includes the following steps:

[0006] S1. Shooting the heat source matrix based on a multi-view synchronous acquisition method to obtain an original image set, and preprocessing the original image set to obtain a preprocessed image set;

[0007] S2. Based on the preprocessed image set, an improved fast corner detection algorithm is used to identify feature points to obtain a feature point set; based on the feature point set, an improved optical flow tracking algorithm is used in combination with a time consistency constraint to track feature points to obtain a feature point vertical trajectory set; based on the feature point vertical trajectory set, an adaptive triangulation algorithm is applied to construct a topological relationship between feature points to obtain a time-varying feature map;

[0008] S3. Based on the time-varying feature graph, a spectral graph analysis algorithm is used to calculate the graph structure change index and obtain the graph topology change index set; based on the feature point vertical trajectory set and the time-varying feature graph, a local rigidity-preserving deformation algorithm is used to estimate the local vertical deformation of the feature point and obtain the local vertical deformation field; based on the local vertical deformation field and the graph topology change index set, a global settlement field equation is constructed and solved to obtain the settlement field function; based on the settlement field function, a dynamic mode decomposition algorithm is used to calculate the dynamic mode set and the corresponding growth rate; based on the dynamic mode set and the growth rate, a nonlinear operator prediction method is used to calculate the predicted settlement field at the future moment;

[0009] S4. Based on the settlement field function and the predicted settlement field, the polynomial chaos expansion method is used to quantify the uncertainty of the settlement estimation and obtain the uncertainty field; based on the uncertainty field, the wavelet packet transform algorithm is used to perform multi-scale error decomposition to obtain a set of error scale coefficients; based on the error scale coefficient set, an adaptive Kalman filter is constructed and applied to optimize the settlement field function to obtain the optimized settlement field function; based on the optimized settlement field function, a sub-pixel interpolation algorithm of surface fitting is used to calculate a high-precision settlement field function.

[0010] Non-contact remote monitoring system for structural settlement, including:

[0011] at least one processor; and,

[0012] a memory communicatively connected to at least one of the processors; wherein,

[0013] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the non-contact remote monitoring method of structural settlement.

[0014] Beneficial effect: the present invention not only realizes non-contact, large-scale, and long-distance structural settlement monitoring, but also improves the resolution and accuracy of settlement monitoring results. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a flow chart of the present invention.

[0016] Figure 2 This is a flow chart of step S1 of the present invention.

[0017] Figure 3 This is a flow chart of step S2 of the present invention.

[0018] Figure 4 This is a flow chart of step S3 of the present invention.

[0019] Figure 5 This is a flow chart of step S4 of the present invention. DETAILED DESCRIPTION

[0020] like Figure 1 As shown, the present application proposes a non-contact remote monitoring method for structural settlement, comprising the following steps:

[0021] S1. Shooting the heat source matrix based on a multi-view synchronous acquisition method to obtain an original image set, and preprocessing the original image set to obtain a preprocessed image set;

[0022] S2. Based on the preprocessed image set, an improved fast corner detection algorithm is used to identify feature points to obtain a feature point set; based on the feature point set, an improved optical flow tracking algorithm is used in combination with a time consistency constraint to track feature points to obtain a feature point vertical trajectory set; based on the feature point vertical trajectory set, an adaptive triangulation algorithm is applied to construct a topological relationship between feature points to obtain a time-varying feature map;

[0023] S3. Based on the time-varying feature graph, a spectral graph analysis algorithm is used to calculate the graph structure change index and obtain the graph topology change index set; based on the feature point vertical trajectory set and the time-varying feature graph, a local rigidity-preserving deformation algorithm is used to estimate the local vertical deformation of the feature point and obtain the local vertical deformation field; based on the local vertical deformation field and the graph topology change index set, a global settlement field equation is constructed and solved to obtain the settlement field function; based on the settlement field function, a dynamic mode decomposition algorithm is used to calculate the dynamic mode set and the corresponding growth rate; based on the dynamic mode set and the growth rate, a nonlinear operator prediction method is used to calculate the predicted settlement field at the future moment;

[0024] S4. Based on the settlement field function and the predicted settlement field, the polynomial chaos expansion method is used to quantify the uncertainty of the settlement estimation and obtain the uncertainty field; based on the uncertainty field, the wavelet packet transform algorithm is used to perform multi-scale error decomposition to obtain a set of error scale coefficients; based on the error scale coefficient set, an adaptive Kalman filter is constructed and applied to optimize the settlement field function to obtain the optimized settlement field function; based on the optimized settlement field function, a sub-pixel interpolation algorithm of surface fitting is used to calculate a high-precision settlement field function.

[0025] like Figure 2 As shown, according to one aspect of the present application, step S1 is further:

[0026] S11, acquiring real-time image stream data from a plurality of pre-arranged infrared optical observation devices, numbering and time-stamping the real-time image stream data, and forming an original image set;

[0027] S12, based on the original image set, using a histogram equalization adaptive exposure algorithm to obtain a contrast-optimized image set;

[0028] S13, based on the contrast-optimized image set, performing an inverse operation using an atmospheric turbulence correction algorithm based on Zernike polynomials to obtain a corrected image set;

[0029] S14, based on the corrected image set, using the wavelet transform multi-resolution analysis method to perform multi-scale image fusion to obtain a fused high-quality image;

[0030] S15. Based on the fused high-quality image and the pre-stored reference image, a feature matching algorithm based on local feature descriptors is applied to identify corresponding point pairs; based on the identified corresponding point pairs, a RANSAC algorithm is used to estimate the transformation matrix; based on the transformation matrix, an affine transformation is performed on the fused high-quality image to obtain a preprocessed image set.

[0031] This embodiment improves the image quality and accuracy of long-distance structural settlement monitoring. Multi-view synchronous acquisition effectively reduces the occlusion and distortion problems caused by a single view, expands the monitoring range, and improves the integrity and reliability of the data. The histogram equalization adaptive exposure algorithm improves the contrast and detail performance of the image under different lighting conditions by dynamically adjusting the pixel value mapping relationship, making it easier to detect tiny structural deformations. The Zernike polynomial atmospheric turbulence correction algorithm effectively compensates for the image distortion caused by atmospheric turbulence in long-distance observations, improves the clarity and stability of the image, and lays the foundation for subsequent precise analysis. The wavelet transform multi-resolution analysis method achieves high-quality image fusion by selecting the optimal image features at different scales, which not only retains large-scale structural information, but also maintains small-scale details, enhancing the information content and analysis value of the image. This embodiment can achieve high-precision and high-reliability long-distance structural settlement monitoring in complex environments, and provides strong technical support for the safety monitoring of important infrastructure such as large buildings, bridges, and tunnels.

[0032] According to one aspect of the present application, step S12 is further:

[0033] S121, based on the original image set, calculating the grayscale histogram of each original image; performing cumulative summation on the grayscale histogram to obtain a cumulative distribution function; normalizing the cumulative distribution function to obtain a grayscale mapping function;

[0034] S122, based on the grayscale mapping function, processing each pixel in the original image to obtain a new pixel value; combining the new pixel values ​​to form a balanced image; based on the balanced image, calculating an average brightness value; based on the average brightness value and a preset target brightness value, calculating a brightness adjustment factor;

[0035] S123, performing linear adjustment based on the brightness adjustment factor and the equalized image to obtain adjusted pixel values; combining the adjusted pixel values ​​to form a brightness adjusted image; dividing the brightness adjusted image into a predetermined number of sub-blocks, and calculating the local contrast of each sub-block;

[0036] S124, based on the local contrast, the preset contrast threshold and the brightness adjusted image, perform adaptive gamma correction on each sub-block to obtain a gamma-corrected sub-block; based on the gamma-corrected sub-block, use a bilinear interpolation method to perform a smooth transition to obtain a smoothed sub-block; recombine the smoothed sub-blocks to form a final adaptive exposure adjusted image; combine all the adaptive exposure adjusted images to form a contrast optimized image set.

[0037] In one embodiment of the present application, the original image set {I i (t)}, calculate the grayscale histogram H(k) of the image, where k is the grayscale level (0-255). Cumulative summation of the histogram H(k) is performed to obtain the cumulative distribution function CDF(k). Normalize the cumulative distribution function CDF(k) to the range of [0, 255] to obtain the grayscale mapping function M(k).

[0038] Read the grayscale mapping function M(k) and the original image I, apply the mapping function M to each pixel p(x, y) in the image I, and obtain the new pixel value p'(x, y) = M(p(x, y)). Combine all the processed pixel values ​​to form the equalized image I eq . Get the equalized image I eq , calculate the average brightness value avg of the image brightness . Read the preset target brightness value target from the system configuration file brightness , calculate the brightness adjustment factor α=target brightness / avg brightness .

[0039] Read the brightness adjustment factor α and the equalized image I eq , for image I eq Each pixel value p in eq (x, y) is linearly adjusted to obtain the adjusted pixel value p adj (x, y) = min(255, α*p eq (x, y)). All adjusted pixel values ​​are combined to form the brightness-adjusted image I adj . Get the brightness adjusted image I adj , divide the image into m×n sub-blocks. For each sub-block B ij Calculate the local contrast Cij = (max(B ij ) - min(B ij )) / (max(B ij ) + min(B ij )) Read the contrast threshold T from the system configuration file contrast Based on the local contrast C ij , contrast threshold T contrast and the brightness adjusted image I adj For each sub-block B ij , if C ij < T contrast , then adaptive gamma correction is applied: γ ij = 1 + (T contrast - C ij ) / T contrast For each pixel p in the sub-block adj (x, y) Apply gamma correction: p gamma (x, y) = 255 * (p adj (x, y) / 255) 1 / γij .

[0040] Get all the gamma-corrected sub-blocks and use the bilinear interpolation method to smoothly transition the boundary pixels of adjacent sub-blocks. Reassemble all the processed sub-blocks to form the final adaptive exposure adjusted image I'. Add I' to the new image set {I' i (t)}. The processed original image set {I i (t)}, and finally output the contrast-optimized image set {I' i (t)}, where each image undergoes histogram equalization and adaptive exposure adjustment.

[0041] This embodiment improves the quality and availability of long-distance structural settlement monitoring images through a histogram equalization adaptive exposure algorithm. By calculating and normalizing the cumulative distribution function, the global grayscale is equalized, the overall contrast of the image is effectively enhanced, and the structural details that were originally difficult to distinguish become clearly visible. By introducing an adaptive brightness adjustment mechanism, the image is dynamically adjusted according to the preset target brightness, which effectively solves the common underexposure or overexposure problems in long-distance observations and ensures that the overall brightness of the image is within the optimal observation range. By dividing the image into multiple sub-blocks and applying local contrast adjustment, the problem of uneven lighting in complex scenes is successfully addressed, so that the details of the shadow area and the highlight area are well preserved. In particular, by introducing adaptive gamma correction, the contrast can be dynamically adjusted according to the characteristics of each sub-block, which not only improves the local detail performance of the image, but also maintains the global visual consistency. The smooth transition of the sub-block boundaries is achieved by the bilinear interpolation method, which effectively eliminates the artificial traces that may be caused by block processing. This embodiment can obtain high-quality structural images under various complex lighting conditions and long-distance observation scenarios, providing a reliable data basis for subsequent precise settlement analysis and improving the applicability and accuracy of long-distance non-contact structural settlement monitoring.

[0042] According to one aspect of the present application, step S13 is further:

[0043] S131, based on the contrast-optimized image set, performing a fast Fourier transform on each image to obtain a frequency domain image; calculating a power spectrum of the frequency domain image; performing a logarithmic transform on the power spectrum to obtain a logarithmic power spectrum; performing radial averaging on the logarithmic power spectrum to obtain a one-dimensional power spectrum density function;

[0044] S132. Based on the one-dimensional power spectral density function, the Kolmogorov turbulence model is fitted using the least squares method, and the best fitting coefficient is calculated to obtain an estimated turbulence intensity parameter; based on the turbulence intensity parameter, the atmospheric coherence length is calculated; based on the atmospheric coherence length, the order of the Zernike polynomial is calculated, and the corresponding Zernike polynomial is generated;

[0045] S133, based on the frequency domain image, calculating the coefficient of each Zernike polynomial; based on the coefficient of each Zernike polynomial, constructing the wavefront phase of atmospheric turbulence, and calculating the atmospheric transfer function; based on the atmospheric transfer function and the frequency domain image, sequentially performing a deconvolution operation and an inverse Fourier transform to obtain a corrected spatial domain image;

[0046] S134, based on the corrected spatial domain image, edge sharpening processing is performed to obtain an edge-sharpened image; based on the edge-sharpened image, the image gradient is calculated using the Sobel operator to obtain a gradient amplitude; the gradient amplitude is added to the corrected spatial domain image to obtain a final corrected image; all the final corrected images are combined to form a corrected image set.

[0047] In one embodiment of the present application, a contrast-optimized image set {I' i (t)}. Perform a fast Fourier transform (FFT) on image I' to obtain a frequency domain image F(u, v). Calculate the power spectrum of F(u, v) P(u, v) = |F(u, v)| 2 Based on the power spectrum P(u, v), logarithmic transformation is applied to obtain the logarithmic power spectrum L(u, v) = log(P(u, v)). The logarithmic power spectrum L(u, v) is radially averaged to obtain the one-dimensional power spectrum density function PSD(f), where f is the spatial frequency.

[0048] Based on the power spectral density function PSD(f), the Kolmogorov turbulence model PSD is fitted using the least squares method model (f) = A * f -5 / 3 , where A is the coefficient to be determined. Calculate the best fitting coefficient A best , get the estimated turbulence intensity parameter C n 2 Based on the turbulence intensity parameter C n 2 and the optical parameters (wavelength λ, aperture D) stored in the system configuration file, calculate the atmospheric coherence length r0 = (0.423 * k 2 * C n 2 * L) -3 / 5 , where k = 2π / λ and L is the optical path length. Calculate the order of the Zernike polynomial J = (D / r0) 2 Generate the J-order Zernike polynomial Z j (ρ, θ), where j = 1, 2, ..., J, ρ and θ are polar coordinates.

[0049] Based on Zernike polynomial Z j (ρ, θ) and the frequency domain image F(u, v), calculate the coefficient a of each Zernike polynomial j = <F(u,v), Z j (ρ, θ)>, where <·, ·> represents the inner product operation. Based on the Zernike coefficient a j, construct the wavefront phase of atmospheric turbulence φ(x, y) =Σ(j=1toJ) a j * Z j (ρ, θ). Calculate the atmospheric transfer function H(u, v) = exp(-i *φ(x, y)).

[0050] Based on the atmospheric transfer function H(u, v) and the frequency domain image F(u, v), a deconvolution operation F is performed corrected (u, v) = F (u, v) / H (u, v). For image F corrected (u, v) performs an inverse Fourier transform (IFFT) to obtain the corrected spatial domain image I corrected .

[0051] For the corrected image I corrected , perform edge sharpening. Use the Sobel operator to calculate the image gradient G x and G y , we get the gradient magnitude G = sqrt(G x 2 + G y 2 ). Add the gradient magnitude G to the corrected image I corrected , and the final corrected image I'' is obtained. The image set {I' i (t)}, output the corrected image set {I'' i (t)}, where each image is corrected for atmospheric turbulence using Zernike polynomials.

[0052] This embodiment effectively solves the image distortion problem caused by atmospheric turbulence in long-distance optical observations through the Zernike polynomial atmospheric turbulence correction algorithm, and improves the accuracy and reliability of structural settlement monitoring. By performing fast Fourier transform and power spectrum analysis on the image, the intensity parameters of atmospheric turbulence can be accurately estimated, providing an accurate reference basis for subsequent correction. By calculating the atmospheric coherence length and generating Zernike polynomials of appropriate order, a mathematical model that can accurately describe the wavefront distortion caused by atmospheric turbulence is successfully constructed. In particular, by calculating the Zernike polynomial coefficients and constructing the atmospheric transfer function, the influence of atmospheric turbulence on optical imaging can be accurately simulated, providing a theoretical basis for image correction. By performing deconvolution operations in the frequency domain, the blur and distortion caused by atmospheric turbulence are effectively removed, and the clarity and geometric accuracy of the image are improved. By introducing edge sharpening processing, not only the influence of atmospheric turbulence is corrected, but also the detail performance of the image is enhanced, so that the slight deformation of the structure is easier to detect. This embodiment can obtain high-quality structural images under long-distance and complex atmospheric conditions, expanding the application scope and accuracy of non-contact settlement monitoring. By effectively eliminating the influence of atmospheric turbulence, a reliable image basis is provided for the precise quantification of long-distance structural settlement, thus improving the accuracy and credibility of the monitoring results.

[0053] According to one aspect of the present application, step S14 is further:

[0054] S141, based on the corrected image set, using the wavelet basis function to perform a two-dimensional discrete wavelet transform on each image to obtain a first-layer decomposition result, including an approximate coefficient matrix and a detail coefficient matrix; based on the first-layer decomposition result, perform a two-dimensional discrete wavelet transform again to obtain a second-layer decomposition result; repeat until a preset layer decomposition result is obtained; combine all decomposition results to obtain a multi-scale decomposition coefficient;

[0055] S142, calculating the energy distribution of the detail coefficient matrix based on the multi-scale decomposition coefficients; determining the significance scale set based on the energy distribution; performing soft threshold processing based on the significance scale set and the multi-scale decomposition coefficients to obtain processed coefficients;

[0056] S143, based on the processed coefficients, using inverse wavelet transform to reconstruct the image to obtain an enhanced image; calculating the structural similarity index between the enhanced image and the images in the corrected image set; if the structural similarity index is lower than a preset threshold, adjusting the parameters of the soft threshold processing, repeating steps S142 to S143 until the maximum number of iterations is reached; otherwise, no processing is performed; outputting the final enhanced image;

[0057] S144, performing weighted fusion on the final enhanced image and the images in the corrected image set to obtain a fused high-quality image.

[0058] In one embodiment of the present application, the corrected image set {I'' i (t)}. Select a suitable wavelet basis function ψ(x), such as Daubechies wavelet. Perform a two-dimensional discrete wavelet transform (2D-DWT) on image I'' to obtain the approximate coefficient matrix A1 and the detail coefficient matrix {H1, V1, D1}. Among them, A1 is the low-frequency component, and H1, V1, and D1 are the high-frequency components in the horizontal, vertical, and diagonal directions, respectively. Read the coefficient matrix A1 and perform 2D-DWT again to obtain the second-level decomposition results: the approximate coefficient matrix A2 and the detail coefficient matrix {H2, V2, D2}. Repeat this process until the preset decomposition layer number L is reached to obtain the multi-scale decomposition coefficient {A L , {H l , V l , D l}, l=1, 2, ..., L}.

[0059] Based on the multi-scale decomposition coefficients, the detail coefficients {H l , V l , D l}Calculate the energy E l =Σ(H l 2 +V l 2 +D l 2 ). According to the energy distribution, the saliency scale set S = {l | E l >T E}, where T E is the preset energy threshold. Read the significant scale set S and multi-scale decomposition coefficients. Perform soft threshold processing on the scale coefficients belonging to the set S: C l ' = sign(C l ) * max(0, |C l | -λ l ), where C l Represents H l 、V l or D l ,λ l is the adaptive threshold.

[0060] Get the processed coefficients {A L , {H l ', V l ', D l '}, l∈S}. Perform inverse wavelet transform (IDWT) to reconstruct the image and obtain the enhanced image I enhanced. Calculate the structural similarity index (SSIM) between the enhanced image and the original image. If the SSIM is lower than the preset threshold T SSIM , then adjust the parameter λ of the soft threshold processing l , repeat the soft threshold processing and image reconstruction process until the SSIM meets the requirement or the maximum number of iterations is reached.

[0061] Get the final enhanced image I enhanced and the original image I''. Calculate the weighted average I final =w*I enhanced +(1-w)*I'', where the weight w is dynamically adjusted according to the SSIM value. final Add to new image collection {I f (t)}. The processed and rectified image set {I'' i (t)}, and output the fused high-quality image set {I f (t)}, where each image is subjected to wavelet transform multi-resolution analysis and enhancement.

[0062] Where ψ(x) is the wavelet basis function; A l is the approximate coefficient matrix of the l-th layer decomposition; H l 、V l , D l are the horizontal, vertical and diagonal detail coefficient matrices of the lth layer decomposition respectively; L represents the total number of layers of wavelet decomposition; E l is the energy of the detail coefficient of the lth layer; T E represents the energy threshold; S is the set of significance scales; C l Represents H l 、V l or D l The general symbol for l is the adaptive threshold for soft threshold processing; T SSIM is the SSIM threshold; w represents the weight of the weighted average; I f (t) represents the final fused high-quality image.

[0063] This embodiment achieves high-quality image fusion and enhancement through the wavelet transform multi-resolution analysis method, and improves the information content and analysis value of the structural settlement monitoring image. By selecting suitable wavelet basis functions and multi-level two-dimensional discrete wavelet transforms, the low-frequency approximate components and high-frequency detail components of the image can be effectively separated, laying the foundation for subsequent multi-scale analysis. By calculating the energy distribution of each scale and determining the significance scale, the frequency components containing key structural information are successfully identified, effectively improving the pertinence and efficiency of subsequent processing. In particular, by introducing adaptive soft threshold processing, the threshold can be dynamically adjusted according to the characteristics of different scales, and useful information can be retained to the maximum extent while suppressing noise. By iteratively optimizing the threshold parameters and monitoring the structural similarity index, it is ensured that the enhanced image not only maintains the structural characteristics of the original image, but also improves the clarity of the details. By weighted averaging the fusion of the original image and the enhanced image, the enhancement effect is maximized on the basis of retaining the original information. This embodiment can extract and integrate the most valuable information from images collected from multiple angles to generate a high-quality fused image. This not only improves the spatial resolution and signal-to-noise ratio of settlement monitoring, but also enhances the detection capability of tiny structural deformation. By providing clearer, more detailed and more informative images, it creates favorable conditions for subsequent accurate settlement analysis and improves the overall performance and reliability of long-distance non-contact structural settlement monitoring.

[0064] According to another aspect of the present application, in step S1, the heat source matrix is ​​photographed based on a multi-view synchronous acquisition method to obtain a set of original images as follows:

[0065] S1a, arrange active heat source matrix on the observation target, and form a fixed temperature gradient field on the surface of the target structure by precisely controlling the heat source output, so that the target and the ambient temperature form a stable contrast;

[0066] S1b. Acquire real-time image streams of multiple infrared optical observation devices, and use a multi-view synchronous acquisition method to shoot the temperature gradient field formed by the heat source matrix to obtain a set of original images.

[0067] In one embodiment of the present application, a plurality of controllable heat sources (such as infrared heaters) are used to form a matrix and arranged on the observation target surface; the output power of each heat source is precisely adjusted by the control system to form a fixed temperature gradient field on the target surface. By adjusting the output of the heat source, the temperature of the target surface is stably contrasted with the ambient temperature. A temperature sensor is used to monitor the temperature distribution of the target surface in real time to ensure the stability of the temperature gradient field. Multiple infrared cameras are configured to ensure that they cover the heat source matrix from different angles; a synchronous control system is used to ensure that all infrared cameras start shooting at the same time. The real-time image stream of the infrared camera is collected to record the temperature gradient field formed by the heat source matrix. The collected original image set is preprocessed, such as denoising, correction, etc. An image registration algorithm is used to align images from different perspectives to generate a complete temperature gradient field image.

[0068] This embodiment improves the image quality and accuracy of long-distance structural settlement monitoring through active heat source matrix control technology and multi-view synchronous acquisition method. The application of active heat source matrix solves the difficulty of obtaining stable images at a long distance. By forming a fixed temperature gradient field on the surface of the target structure, a stable infrared imaging condition is created, which improves the recognition and positioning accuracy of feature points. It effectively overcomes the influence of ambient temperature changes and atmospheric interference on infrared imaging, and provides high-quality raw data for subsequent image processing and analysis. Multi-view synchronous acquisition technology further enhances the integrity and reliability of the data, and effectively reduces the occlusion and distortion problems caused by a single perspective. Combined with the histogram equalization adaptive exposure algorithm, the contrast and detail performance of the image under different lighting conditions are improved. This embodiment can achieve high-precision and high-reliability long-distance structural settlement monitoring in complex environments, and provides strong technical support for the safety monitoring of important infrastructure such as large buildings, bridges, and tunnels.

[0069] like Figure 3 As shown, according to one aspect of the present application, step S2 is further:

[0070] S21. Based on the preprocessed image set, an improved fast corner point detection algorithm is used to identify the significant corner points in the heat source matrix; non-maximum suppression and sub-pixel precision optimization are performed on the significant corner points, and finally a feature point set is obtained;

[0071] S22. Based on the feature point set and the preprocessed image set, an enhanced binary descriptor generation algorithm is used to generate a descriptor for each feature point to obtain a descriptor set; based on the feature point set and descriptor set of the current frame, and the pre-stored feature point set and descriptor set of the previous frame, an improved optical flow tracking algorithm is applied to estimate the vertical displacement of the feature points between adjacent frames; based on the estimated vertical displacement, a time consistency constraint is introduced to filter out the vertical displacement that does not conform to the expected vertical trajectory; and finally a vertical trajectory set of the feature points is obtained;

[0072] S23, based on the feature point vertical trajectory set, calculating the multidimensional feature vector of each vertical trajectory, including the vertical trajectory length, average speed and acceleration change; based on the multidimensional feature vector, calculating the Mahalanobis distance of each vertical trajectory to its k nearest neighbors in the feature space; based on the Mahalanobis distance, using an adaptive threshold method to identify outliers, marking the outliers as abnormal vertical trajectories; removing the abnormal vertical trajectories from the feature point vertical trajectory set, and finally obtaining a filtered vertical trajectory set;

[0073] S24. Based on the filtered vertical trajectory set, Delaunay triangulation is performed on the plane, and shape constraints of triangles in the Delaunay triangulation are dynamically adjusted to obtain a triangulation result; the triangulation result is converted into a graph structure to obtain a time-varying feature graph.

[0074] In one embodiment of the present application, the registered image I is read r (t), an improved fast corner detection algorithm is applied to identify the significant corners in the heat source matrix by sliding a small window on the image and comparing the intensity difference between the pixels in the window and the central pixel. Then, the detected corners are subjected to non-maximum suppression and sub-pixel precision optimization, and finally the feature point set P(t) = {p1(t), p2(t), ..., p n (t)}, where p i (t) represents the exact coordinates of the i-th feature point at time t.

[0075] Get the feature point set P(t) and the registered image I r (t), apply the enhanced binary descriptor generation algorithm to each feature point: define a set of pixel pairs around the feature point, compare the intensity values ​​of each pair of pixels, generate a binary string, and dynamically adjust the spatial distribution and contrast threshold of the pixel pairs to improve the recognition of the descriptor. Finally, the descriptor set D(t)={d1(t), d2(t), ..., d n (t)}, where d i (t) is the corresponding p i (t) is a binary descriptor of

[0076] Read the feature point set P(t) and descriptor set D(t) of the current frame, as well as the corresponding data P(t-1) and D(t-1) of the previous frame obtained from the system memory. Apply the improved optical flow tracking algorithm: combine descriptor matching and local image block similarity measurement to estimate the vertical displacement of feature points between adjacent frames. At the same time, introduce temporal consistency constraints, and filter out matching results that do not meet the expected vertical trajectory by analyzing the motion patterns of feature points in multiple consecutive frames. Finally, the feature point vertical trajectory set T(t) = {T1( t), T2(t), ..., T m (t)}, where T i (t) represents the complete vertical trajectory of the i-th feature point from the initial frame to the current frame t.

[0077] Get the feature point vertical trajectory set T(t), and apply the robust statistical distance anomaly detection algorithm to filter the vertical trajectories: calculate the multidimensional feature vector of each vertical trajectory, including vertical trajectory length, average speed, acceleration change, etc.; calculate the Mahalanobis distance of each vertical trajectory to its k nearest neighbors in the feature space, and use the adaptive threshold method to identify outliers. The vertical trajectories marked as abnormal will be removed from the set, and finally the filtered vertical trajectory set T'(t) is obtained.

[0078] Read the filtered vertical trajectory set T'(t), and apply the adaptive triangulation algorithm to construct the topological relationship between the feature points: project all feature points on the plane, and then perform Delaunay triangulation. In order to adapt to the possible deformation of the heat source matrix, the shape constraints of the triangles are dynamically adjusted to allow irregular triangles to be formed in high deformation areas; the triangulation result is converted into a graph structure to obtain the feature graph G(t) = (V(t), E(t)), where V(t) is the vertex set (corresponding to the feature points) and E(t) is the edge set (representing the topological connection between the feature points).

[0079] This embodiment achieves efficient and accurate feature extraction and tracking through algorithms such as improved fast corner detection, enhanced binary descriptor generation, improved optical flow tracking, and robust statistical distance anomaly detection. The improved fast corner detection algorithm improves the speed and accuracy of feature point recognition by optimizing the window sliding strategy and non-maximum suppression method, and can quickly locate the key points in the structure, providing a stable reference point for settlement analysis. The enhanced binary descriptor generation algorithm adopts an adaptive sampling mode and a dynamic threshold strategy to improve the recognition and robustness of the descriptor, making the feature matching at different times and different perspectives more accurate and reliable. The improved optical flow tracking algorithm combines time consistency constraints and covariance propagation methods to improve the accuracy and stability of feature point tracking, and can effectively deal with small structural deformations and environmental interference. The robust statistical distance anomaly detection algorithm successfully identifies and eliminates inconsistent vertical trajectories through multidimensional feature analysis and density accessibility principles, thereby improving the reliability of tracking results. This embodiment can accurately capture small deformations of the structure in a complex actual monitoring environment, provide high-quality data support for subsequent settlement analysis, and improve the accuracy and reliability of long-distance non-contact structural settlement monitoring.

[0080] In another embodiment of the present application, a feature point set P(t) = {p1(t), p2(t), ..., p n (t)} and the fused high-quality image I f (t). For each feature point p i (t), determine its neighborhood radius r, extract i (t) is the center of the circular image block B with a radius of r i Based on image block B i , apply Gaussian filter G(σ) to B i Perform smoothing to obtain the filtered image block B i '. Calculate B i 'The gradient magnitude M(x, y) and direction θ(x, y).

[0081] Get the gradient information M(x, y) and θ(x, y) at the feature point p i (t) generates a sampling pattern S = {(x k , y k ) |k=1, 2, ..., K}, where K is the number of sampling point pairs. The positions of the sampling point pairs are determined according to the logarithmic polar coordinate distribution. Based on the sampling mode S and the gradient information, for each sampling point pair (x k , y k ), calculate the gradient difference value d k =M(x k , y k )-M(-x k , -y k ). Generate the initial binary string b={b k | b k = 1 if d k > 0, else 0; k = 1, 2, ..., K}.

[0082] Get the initial binary string b and calculate the variance var(b) of b. If the variance var(b) is less than the preset threshold T var , then adjust the spatial distribution of the point pairs in the sampling pattern S, and repeat until var(b)≥T var Or the maximum number of iterations is reached. Based on the final binary string b and the gradient direction θ(x, y), calculate the main direction θ main = mode(θ(x, y)), where mode represents the mode. main Perform a cyclic vertical shift on b to obtain a rotationally invariant binary string b rot .

[0083] Based on the binary string b rot , calculate the Hamming weight wH (b rot ). If the Hamming weight w H (b rot ) is not within the preset range [w min , w max ], then for the binary string b rot Perform a bit flip operation to obtain the adjusted binary string b adj , so that w min ≤w H (b adj )≤w max . The adjusted binary string b adj Convert to decimal number d i . Calculate the feature point p i The spatial coordinates (x) of (t) i ,y i ). Combination d i and (x i , y i ) forms the final descriptor d i (t) = (d i , x i , y i ). i (t) is added to the descriptor set D(t). Process all feature points in the feature point set P(t) and output the final descriptor set D(t) = {d1(t), d2(t), ..., d n (t)}, where each descriptor is obtained through an enhanced binary descriptor generation algorithm.

[0084] Where P(t) represents the set of feature points; p i (t) represents the i-th feature point; I f (t) represents the fused high-quality image; r represents the neighborhood radius of the feature point; B i represents a circular image block centered on the feature point; G(σ) represents a Gaussian filter, σ is the standard deviation; M(x, y) represents the gradient amplitude; θ(x, y) represents the gradient direction; S represents the sampling mode; K represents the logarithm of the sampling points; d k represents the gradient difference value; b represents the binary string; T var represents the variance threshold; θ main Indicates the main direction; b rot represents a rotationally invariant binary string; w H represents the Hamming weight; w min , w max Indicates the preset range of Hamming weight; b adj Represents the adjusted binary string; d iRepresents the binary representation of the descriptor (converted to decimal); (x i , y i ) represents the spatial coordinates of the feature points; d i (t) represents the final descriptor; D(t) represents the descriptor set.

[0085] This embodiment improves the accuracy and robustness of feature point description through an enhanced binary descriptor generation algorithm, and provides a reliable basis for feature matching and tracking in long-distance structural settlement monitoring. By adopting a dynamically adjusted sampling mode around the feature points, it can adapt to local structural features of different scales and directions, and improve the distinguishing ability of the descriptor. By calculating the gradient difference between the sampling point pairs and generating the initial binary string, the complex local structure of the image is successfully encoded into a compact binary representation, which not only improves the computational efficiency but also retains the key structural information. In particular, by introducing the variance threshold and iterative optimization mechanism, the sampling mode can be dynamically adjusted to ensure that the generated descriptor has sufficient information entropy, effectively improving the stability of the descriptor under different viewing angles and lighting conditions. By calculating the main direction and performing cyclic vertical displacement, the rotation invariance of the descriptor is achieved, so that feature matching can effectively cope with the rotational deformation of the structure. By adjusting the Hamming weight of the binary string, the statistical characteristics of the descriptor are optimized, and its performance in feature matching is further enhanced. This embodiment can accurately and efficiently describe and match structural feature points in a complex long-distance monitoring environment, and provide high-quality corresponding point data for subsequent settlement analysis. The accuracy and real-time performance of long-distance non-contact structural settlement monitoring are improved by generating highly discriminative, rotation-invariant and computationally efficient binary descriptors.

[0086] In another embodiment of the present application, the current frame descriptor set D(t) and the previous frame descriptor set D(t-1) read from the system memory are obtained. For each descriptor d in D(t) i (t), search for the most similar descriptor in D(t-1), using the Hamming distance H(d i (t), d j (t-1)) as the similarity measure. Establish the initial matching pair set M={(d i (t), d j (t-1))}. Read the initial matching pair set M and the current frame image I f (t) and the previous frame image I f (t-1), for each matching pair (d i (t), d j (t-1)), extract the image block B centered on the corresponding feature point i (t) and B j (t-1). Calculate the brightness constraint equation ΔI=I for each pair of image blocksx * u + I y * v + I t , where I x ,I y is the spatial gradient, I t is the time gradient, (u, v) is the optical flow vector to be calculated.

[0087] Based on the brightness constraint equations, construct the augmented matrix A = [I x , I y ] and vector b = -I t The initial optical flow estimate w is obtained by solving the equation Aw = b using weighted least squares, where the weights are based on the similarity of the descriptors. init = (u init , v init ). Apply the covariance propagation method to estimate the uncertainty of optical flow and calculate the covariance matrix Σ w =(A T *W*A) -1 , where W is the weight matrix. According to the covariance matrix Σ w Compute confidence intervals for optical flow estimates.

[0088] Based on the covariance matrix Σ w , the optical flow estimation is optimized using iterative reweighted least squares (IRLS). In each iteration, the weights are updated according to the residual and the optical flow vector w is recalculated. k and the covariance matrix Σ wk . Get the optimized optical flow vector w opt and the covariance matrix Σ wopt ,Apply the temporal consistency constraint, compare the current estimate with the optical flow estimates of the previous frames, and use Kalman filtering to correct it if the current estimate deviates too much from the historical vertical trajectory.

[0089] Read the corrected optical flow vector w final . Update the position of feature points: p i (t) = p j (t-1) + w final . Calculate the tracking quality index Q i , based on the confidence and temporal consistency of optical flow estimation. i (t), Q i ) is added to the vertical trajectory set T(t). Process all feature point matching pairs and output the updated vertical trajectory set T(t)={(p1(t), Q1), (p2(t), Q2), ..., (p m (t), Q m)}, where each element contains the updated feature point position and the corresponding tracking quality indicator.

[0090] This embodiment realizes high-precision and high-reliability feature point tracking through an improved optical flow tracking algorithm, providing key technical support for continuous monitoring of structural settlement. By combining descriptor matching and brightness constraint equations, it is possible to maintain stable tracking performance under large vertical displacement and illumination changes, effectively solving the common feature point loss problem in long-distance monitoring. By introducing the weighted least squares method to solve the optical flow equation, the reliability differences in different regions are successfully considered, and the accuracy of optical flow estimation is improved. In particular, by applying the covariance propagation method to estimate the uncertainty of optical flow, the credibility of the tracking results can be quantified, providing an important basis for subsequent data fusion and decision-making. By optimizing the optical flow estimation through iterative reweighted least squares method, the robustness to outliers is effectively improved, ensuring the tracking stability in complex environments. By introducing time consistency constraints and Kalman filtering, the influence of short-term noise is successfully suppressed, and the smooth continuity of the vertical trajectory of the feature points is achieved. This embodiment can maintain high-precision feature point tracking in long-distance and long-term structural monitoring, providing continuous and reliable vertical displacement data for settlement analysis. By accurately capturing the tiny deformations and movements of the structure, the temporal resolution and sensitivity of non-contact settlement monitoring are improved, providing technical support for the timely detection and assessment of potential structural risks.

[0091] In another embodiment of the present application, a vertical trajectory set T(t) is obtained. i Calculate its eigenvector f i = [Δx,Δy, v x , v y , a x , a y , Q], where Δx and Δy are vertical displacements, v x and v y is the speed, a x and a y is the acceleration, and Q is the average tracking mass.

[0092] Based on the feature vector set {f i}. Calculate the mean μ and covariance matrix Σ of the eigenvector. Use the Minimum Covariance Determinant (MCD) algorithm to estimate the robust μ rob and Σ rob , to reduce the impact of outliers. For each eigenvector f i Calculate the Mahalanobis distance d i =sqrt((f i -μrob ) T *Σ rob -1 *(f i -μ rob )). All d i Stored in the distance set D. The adaptive threshold T is calculated using the median absolute deviation (MAD) method adap =median(D)+k*MAD(D), where k is an adjustable parameter. adap The vertical traces corresponding to the elements of are marked as potential anomalies.

[0093] Based on the set of potential abnormal vertical trajectories, for each potential abnormal vertical trajectory, check the status of its k nearest neighbors. If more than half of the nearest neighbors are also marked as potential abnormalities, the vertical trajectory is confirmed to be abnormal, otherwise it is reclassified as normal. Based on the set of confirmed abnormal vertical trajectories, calculate the local density ρ of these vertical trajectories i and the relative distance Δ i According to ρ i and Δ i The product of identifies the core anomaly in the abnormal vertical trajectory, that is, the vertical trajectory with the largest product value.

[0094] Based on the core abnormal vertical trajectory, with the core anomaly as the center, other abnormal vertical trajectories are clustered by the density accessibility principle to form abnormal vertical trajectory clusters. The representative features of each cluster are calculated, such as cluster center, range, and duration. The vertical trajectories belonging to the abnormal vertical trajectory cluster are removed from the original vertical trajectory set T(t). The mean shift algorithm is applied to the remaining vertical trajectories for smoothing to obtain the filtered vertical trajectory set T'(t). The filtered vertical trajectory set T'(t) and the abnormal vertical trajectory cluster information are output for subsequent analysis and visualization to assist in identifying potential structural anomalies or environmental interference.

[0095] This embodiment effectively identifies and eliminates inconsistent feature point vertical trajectories through a robust statistical distance anomaly detection algorithm, thereby improving the reliability and accuracy of settlement monitoring data. By calculating multidimensional feature vectors, the dynamic characteristics of vertical trajectories such as vertical displacement, velocity, and acceleration are fully considered, providing a rich information basis for anomaly detection. By applying the MCD algorithm to estimate the robust mean and covariance, the influence of outliers on statistical characteristic estimation is successfully reduced, and the accuracy of subsequent anomaly detection is improved. In particular, by calculating the Mahalanobis distance and setting the adaptive threshold using the median absolute deviation (MAD) method, it is possible to flexibly respond to data distribution characteristics under different monitoring scenarios, and achieve highly adaptive anomaly detection. By considering the k-nearest neighbor state and local density, not only isolated outliers are identified, but also clustered abnormal vertical trajectories are successfully detected, improving the comprehensiveness and reliability of anomaly detection. Clustering abnormal vertical trajectories through the density accessibility principle achieves refined classification of different types of anomalies, providing an important basis for subsequent cause analysis and processing. This embodiment can effectively filter out erroneous vertical trajectories caused by equipment errors, environmental interference or local structural anomalies in a complex monitoring environment, and retain valid data that truly reflects the state of structural settlement. By providing high-quality and reliable vertical trajectory data, a solid data foundation is laid for accurately evaluating the structural settlement status and trend, and the overall credibility and practicality of long-distance non-contact structural settlement monitoring is improved.

[0096] like Figure 4 As shown, according to one aspect of the present application, step S3 is further:

[0097] S31. Based on the time-varying feature graph, calculate the Laplace matrix at adjacent moments; perform eigenvalue decomposition on the Laplace matrix to obtain an eigenvalue sequence; based on the eigenvalue sequence, calculate a graph topology change index set;

[0098] S32, based on the vertical trajectory set of feature points, dividing the time-varying feature map into overlapping local areas; based on each local area, constructing an energy function, including a data term and a regularization term; minimizing the energy function through an iterative optimization method to obtain an optimal deformation vector for each feature point; based on the optimal deformation vector of each feature point, forming a local vertical deformation field;

[0099] S33. Based on the local vertical deformation field and the graph topology change index set, a global settlement field equation is constructed, including a local vertical deformation constraint term and a global consistency constraint term; the global settlement field equation is discretized into a large-scale sparse linear system using a variational method, and solved by a conjugate gradient method to obtain a settlement field function;

[0100] S34, based on the settlement field function and the historical settlement field function set pre-stored in the system database, discretize to obtain a high-dimensional data matrix; based on the high-dimensional data matrix, construct a time offset matrix pair; based on the time offset matrix pair, calculate the characteristic decomposition of the dynamic matrix through singular value decomposition and low-rank approximation, and finally obtain a dynamic mode set and a corresponding growth rate;

[0101] S35. Project the settlement field function to the dynamic pattern space to obtain a pattern coefficient vector; predict the pattern coefficient at a future time based on the pattern coefficient vector and the growth rate; linearly combine the pattern coefficient at a future time with the dynamic pattern set to obtain a predicted settlement field at a future time.

[0102] In one embodiment of the present application, the current feature graph G(t) and the previous feature graph G(t-1) obtained from the system memory are read, and the Laplace matrices L(t) and L(t-1) of the two graphs are calculated. The two Laplace matrices are subjected to eigenvalue decomposition to obtain an eigenvalue sequence. By comparing the changes in the corresponding eigenvalues, the graph topology change index Δλ(t) = {Δλ1(t), Δλ2(t), ..., Δλ k (t)}, where Δλ i (t) represents the relative rate of change of the i-th eigenvalue.

[0103] Based on the feature map G(t) and the filtered vertical trajectory set T'(t), the local rigidity preserving deformation algorithm is applied to estimate the local deformation of the feature points: the feature map is divided into overlapping local regions. For each local region, an energy function is constructed, which includes a data term (maintaining vertical trajectory consistency) and a regularization term (maintaining local rigidity). The energy function is minimized by an iterative optimization method to obtain the optimal deformation vector of each feature point. The final output is the local vertical deformation field D(t)={d1(t), d2(t), ..., d m (t)}, where d i (t) represents the local deformation vector of the i-th feature point.

[0104] The local vertical deformation field D(t) and the graph topology change index Δλ(t) are read to construct the global settlement field equation. The equation contains two main terms: local deformation constraint term and global consistency constraint term. The local deformation constraint term ensures the accuracy of local deformation based on D(t), and the global consistency constraint term uses Δλ(t) to ensure the rationality of the overall topological structure. The equation is discretized into a large-scale sparse linear system using the variational method, and solved by the conjugate gradient method to obtain the settlement field function S(x, y, t), which represents the settlement at the coordinate (x, y) at time t.

[0105] The current settlement field function S(x, y, t) and the historical settlement field function set {S(x, y, τ)|τ∈[0, t)} stored in the system database are obtained, and the dynamic mode decomposition algorithm is applied to extract the spatiotemporal dynamic characteristics of settlement: the settlement field function is discretized into a high-dimensional data matrix X, in which each column represents the settlement field state at a time point; the time offset matrix pair X and X' is constructed, where X' is obtained by removing the first column of X and adding the current state at the end; the eigendecomposition of the dynamic matrix A is calculated through singular value decomposition and low-rank approximation, and finally the dynamic mode set Φ={φ1, φ2, ..., φ r} and the corresponding growth rate Λ={λ1,λ2,...,λ r}, where φ i represents the i-th dynamic mode, λ i Indicates its corresponding growth rate.

[0106] The dynamic pattern set Φ, growth rate Λ and current settlement field function S(x, y, t) are read, and the nonlinear prediction method based on Koopman operator theory is used to predict settlement: the current settlement field S(x, y, t) is projected into the dynamic pattern space to obtain the pattern coefficient vector a(t); the growth rate Λ is used to predict the pattern coefficient a(t+Δt) = diag(exp(Λ·Δt))·a(t) at future times; the settlement field at future times is reconstructed by linearly combining the predicted pattern coefficients and the dynamic pattern to obtain the predicted settlement field function S'(x, y, t+Δt), which represents the settlement prediction result after Δt time.

[0107] This embodiment achieves high-precision vertical displacement estimation and settlement analysis through methods such as spectral analysis, local rigidity-preserving deformation, variational method to solve the global settlement field equation, dynamic mode decomposition and nonlinear prediction based on Koopman operator theory. The spectral analysis algorithm sensitively captures the slight changes in the structural topology by calculating the changes in the eigenvalues ​​of the Laplace matrix of the graph structure, providing a global perspective for settlement analysis. The local rigidity-preserving deformation algorithm accurately estimates the local deformation of the feature points by balancing the local rigidity and observation data, and effectively handles the non-uniform settlement problem of the structure. The variational method to solve the global settlement field equation combines local deformation constraints and global consistency constraints, realizes accurate reconstruction from discrete feature points to continuous settlement fields, and provides an important basis for comprehensively evaluating the settlement state of the structure. The dynamic mode decomposition algorithm deeply reveals the inherent laws of the settlement process by extracting the spatiotemporal dynamic characteristics of the settlement, and provides a scientific basis for predicting future settlement trends. The nonlinear prediction method based on Koopman operator theory overcomes the limitations of traditional linear prediction methods, can accurately capture the nonlinear dynamic characteristics of the settlement process, and improves the accuracy and reliability of prediction. This embodiment not only realizes the accurate estimation of the current settlement state, but also can reliably predict the future settlement trend, provides strong decision-making support for structural safety assessment and preventive maintenance, and improves the foresight and scientific nature of structural settlement monitoring.

[0108] According to one aspect of the present application, the graph topology change index set calculated based on the eigenvalue sequence in step S31 is further:

[0109] S311, based on the eigenvalue sequence, calculate the relative change rate of each pair of adjacent eigenvalues; store all relative change rates into the eigenvalue change rate array, calculate the mean and standard deviation of the eigenvalue change rate array; based on the mean and standard deviation, use the three times standard deviation principle to identify the eigenvalues ​​with significant changes, and generate a set of significant change eigenvalue indexes;

[0110] S312, extracting corresponding eigenvectors based on the significant change eigenvalue index set; calculating the angle of the corresponding eigenvectors; calculating the mean and standard deviation of the angle; identifying the significantly rotated eigenvectors based on the mean and standard deviation of the angle, and generating a significant change set;

[0111] S313. Calculate the comprehensive change index based on the significant change set, the eigenvalue sequence and the eigenvalue change rate array; normalize the comprehensive change index to obtain the final graph topology change index, and form a graph topology change index set.

[0112] In one embodiment of the present application, the feature graph G(t) = (V(t), E(t)) at the current moment and the feature graph G(t-1) = (V(t-1), E(t-1) at the previous moment read from the system memory are obtained. For the feature graphs G(t) and G(t-1), the degree matrices D(t) and D(t-1) are calculated respectively, where D ii Denotes the degree of vertex i. Based on the degree matrices D(t) and D(t-1), and the original adjacency matrices A(t) and A(t-1), calculate the Laplacian matrices L(t) = D(t) - A(t) and L(t-1) = D(t-1) -A(t-1). Use the QR algorithm to perform eigenvalue decomposition on the Laplacian matrices L(t) and L(t-1) to obtain the eigenvalue sequence λ(t) = {λ1(t), λ2(t), ..., λ n (t)} and λ(t-1) = {λ1(t-1), λ2(t-1), ..., λ n (t-1)}.

[0113] Based on the eigenvalue sequences λ(t) and λ(t-1), calculate the relative change rate Δλ of each pair of corresponding eigenvalues i =(λ i (t)-λ i (t-1)) / λ i (t-1). The calculated Δλ i Stored in the array Δλ. Calculate the mean μ of the array Δλ Δλ and standard deviation σ Δλ The triple standard deviation (3σ) principle is used to identify the eigenvalues ​​with significant changes, i.e. |Δλ i -μ Δλ |>3σ Δλ The index i of .

[0114] Read the index set of significantly changed eigenvalues. For each significantly changed eigenvalue, extract its corresponding eigenvector v i (t) and v i (t-1). Calculate the angle θ between these eigenvectors i =arccos((v i (t)·v i (t-1)) / (||v i (t)||*||v i (t-1)||)). Based on the eigenvector angle θ i . Calculate the mean value μ of the angle θ and standard deviation σ θ , identify the significantly rotated eigenvectors, i.e., θ i >μ θ +2σ θThe indices i are stored in the significant change set S.

[0115] Read the significant change set S, as well as the eigenvalue sequence and eigenvalue change rate array. For each index i in the significant change set S, calculate the comprehensive change index C i = w1* |Δλ i | + w2 *θ i , where w1 and w2 are predefined weights. i} is normalized to obtain the final graph topology change index Δλ(t) = {Δλ1(t), Δλ2(t), ..., Δλ k (t)}, where k is the size of the significant change set S. The output graph topology change index Δλ(t) is used for subsequent analysis.

[0116] This embodiment achieves high-sensitivity detection of structural topological changes through the spectrum analysis algorithm, and provides strong technical support for comprehensively evaluating the impact of settlement on the integrity of the structure. By constructing and analyzing the Laplace matrix, the complex structural topological relationship is successfully encoded into a mathematical expression, laying the foundation for subsequent fine analysis. By comparing the eigenvalue sequence before and after the moment, it is possible to sensitively capture the slight changes in the structural topology and effectively identify early signs of settlement that may be ignored by traditional methods. In particular, by calculating the relative rate of change of the eigenvalue and the angle of the eigenvector, not only the intensity of the topological change is quantified, but also the direction and pattern of the change are revealed, providing important clues for a deep understanding of the settlement mechanism. By introducing the 3σ principle and the adaptive threshold strategy, the detection sensitivity and reliability are successfully balanced, effectively reducing false positives and false negatives. By comprehensively considering the eigenvalue changes and eigenvector rotations, a comprehensive graph topological change index is generated, providing a multi-dimensional evaluation basis for subsequent settlement analysis. This embodiment can accurately identify and quantify the local and global topological changes caused by settlement in a complex structural system, providing a scientific basis for evaluating the overall stability of the structure and predicting potential risks. By capturing subtle changes in structural topology, the sensitivity and foresight of long-distance non-contact structural settlement monitoring are significantly improved, creating conditions for timely detection and prevention of serious structural problems.

[0117] In another embodiment of the present application, a feature map G(t) = (V(t), E(t)) and a filtered vertical track set T'(t) are obtained. The feature map G(t) is divided into overlapping local regions {R1, R2, ..., R m}, each region contains a set of adjacent vertices and edges.

[0118] For the local area set {R i Each region R in i, extract the vertex set V it contains i and edge set E i . Construct the local rigid energy function E rigid (R i ) =Σ j，k∈Ei w jk ||(v j -v k ) - (v j 0 -v k 0 )|| 2 , where v j and v k is the position after deformation, v j 0 and v k 0 is the initial position. Based on the local rigid energy function E rigid (R i ) and the vertical trajectory set T'(t), construct the data item energy function E data (R i ) =Σ j∈Vi w j ||v j -p j || 2 , where p j is the observation position in the corresponding vertical trajectory.

[0119] Based on the function E rigid (R i ) and E data (R i ), construct the total energy function E total (R i )=α*E rigid (R i )+(1-α)*E data (R i ), where α is the equilibrium parameter. Use gradient descent to minimize E total (R i ), and obtain the local optimal deformation v i *. Based on local optimal deformation v i *, calculate the affine transformation matrix A before and after deformation i , so that v i *≈A i *v i 0 . For the affine transformation matrix A i Perform singular value decomposition A i =U i *Σ i *V iT , extract the rotation component R i =U i *V i T and the scaling component S i =V i *Σ i *V i T .

[0120] Based on the rotation component R i and the scaling component S i , calculate the stiffness measure ρ of the local deformation i = ||R i || F / ||S i || F , where ||·|| F represents the Frobenius norm. i With the preset threshold T ρ Compare, if ρ i < T ρ , then mark the area as a non-rigid deformation area. For the area marked as non-rigid, increase its rigid weight w rigid , repeat the above steps until the stiffness measure ρ of all regions is i ≥T ρ Or the maximum number of iterations is reached, and the final local deformation result is obtained. The local deformation result is interpolated to the entire feature map using the moving least squares method to obtain the global vertical deformation field D(t) = {d1(t), d2(t), ..., d m (t)}, where d i (t) represents the local deformation vector of the i-th feature point. Based on the global vertical deformation field D(t), the divergence div(D) and curl(D) of the vertical deformation field are calculated to characterize the expansion and rotation characteristics of the deformation. The global vertical deformation field D(t) and its feature quantities div(D) and curl(D) are output for subsequent analysis.

[0121] This embodiment achieves accurate estimation of local deformation of the structure through a local rigidity-preserving deformation algorithm, and provides key technical support for accurately evaluating the impact of non-uniform settlement on the structure. By dividing the structure into overlapping local regions, the algorithm successfully balances global consistency and local flexibility, and can accurately capture differentiated settlement patterns in complex structures. By constructing local rigid energy functions and data item energy functions, the prior knowledge and observation data are effectively combined to improve the reliability and accuracy of deformation estimation. In particular, by introducing equilibrium parameters to dynamically adjust the weights of rigid constraints and data fitting, it is possible to flexibly respond to settlement characteristics in different regions, ensuring the local adaptability of the estimation results. By calculating the affine transformation matrix before and after deformation and performing singular value decomposition, the rotation and scaling components are successfully separated, providing an important basis for in-depth understanding of the structural deformation mechanism caused by settlement. By calculating the rigidity metric of local deformation and applying an iterative optimization strategy, the degrees of freedom of deformation and the overall stability of the structure are effectively balanced, ensuring the physical rationality of the estimation results. This embodiment can accurately reconstruct the local vertical deformation field of the structure in a complex non-uniform settlement scenario, providing detailed and reliable data support for evaluating the impact of settlement on structural integrity. By accurately quantifying local deformation, the spatial resolution and accuracy of long-distance non-contact structural settlement monitoring are improved, providing a scientific basis for formulating targeted maintenance and reinforcement strategies.

[0122] In another embodiment of the present application, the global vertical deformation field D(t) = {d1(t), d2(t), ..., d m (t)} and the graph topology change index Δλ(t) = {Δλ1(t), Δλ2(t), ..., Δλ k (t)}. Construct the energy functional E[S] = E data [S] +α* E smooth [S] +β* E topo [S], where S(x, y, t) is the required sedimentation field function, α and β represent weight coefficients; based on the energy functional E[S], define the data item E data [S] = Σ i w i * ||S(x i , y i , t) - d i (t)|| 2 , where (x i , y i ) is the coordinate of the feature point, w i is the weight. Calculate the weight w of each feature point i , based on its reliability in vertical deformation fields.

[0123] Based on data item E data[S], construct the smoothing term E smooth [S] =∫∫(|▽S| 2 + |▽ 2 S| 2 ) dxdy, where ▽S and ▽ 2 S represents the gradient and Hessian matrix of S respectively. The smooth term is discretized using the finite difference method. Based on the smooth term E smooth [S] and the graph topology change index Δλ(t). Define the topology constraint term E topo [S] = Σ j γ j *|∫∫φ j (x, y)*S(x, y, t) dxdy-Δλ j (t)| 2 , where φ j (x, y) corresponds to Δλ j (t) is the characteristic function of

[0124] Based on the topological constraint E topo [S], discretize the energy functional E[S] into a large-scale sparse linear system Ax = b, where x represents the discretized sedimentation field value, A is the coefficient matrix, and b is the right-hand side vector. Use the conjugate gradient method to solve the linear system. Read the discrete sedimentation field value x obtained by the solution, and use bicubic spline interpolation to reconstruct the discrete value into a continuous function S(x, y, t). Calculate the gradient field ▽S and Laplacian field ▽ of the continuous function S(x, y, t) 2 S is used to characterize the direction and acceleration of the settlement. Critical points in the settlement field are identified, including local extreme points and saddle points. The positions and values ​​of these critical points are calculated for subsequent analysis.

[0125] Based on the critical point information and the settlement field function S(x, y, t), the set of contour lines of the settlement field is calculated. i | S(x,y, t) = h i}, where h i is the preset settlement threshold. Analyze the shape and distribution of contour lines to identify abnormal settlement areas. Comprehensive settlement field function S(x, y, t), gradient field ▽S, Laplacian field ▽ 2 S, critical point and contour information, generate the settlement field feature descriptor F S Output the sedimentation field function S(x, y, t) and feature descriptor F S For subsequent analysis.

[0126] This embodiment solves the global settlement field equation by the variational method, realizes the high-precision reconstruction from discrete observation points to the continuous settlement field, and provides a powerful analysis tool for comprehensively evaluating the settlement state of the structure. By constructing an energy functional containing data terms, smoothing terms and topological constraints, the discrete observation data, physical prior knowledge and global topological features are successfully integrated into a unified mathematical framework, laying a theoretical foundation for accurately reconstructing the settlement field. By introducing weight coefficients to dynamically adjust the reliability of different observation points, the influence of outliers and noise is effectively reduced, and the robustness of the reconstruction results is improved. In particular, by constructing and solving large-scale sparse linear systems, the continuous problem is successfully discretized into a computable form, achieving efficient and accurate numerical solutions. By calculating the gradient and Laplacian of the settlement field, not only the settlement amount is reconstructed, but also the spatial distribution information of the settlement rate and acceleration is provided, creating conditions for in-depth analysis of the dynamic characteristics of the settlement. By identifying critical points and calculating contour lines, the key features of the settlement field are successfully extracted, providing strong support for the rapid positioning of high-risk areas. This embodiment can accurately reconstruct a continuous global settlement field from limited observation data, providing a detailed and reliable data basis for comprehensively evaluating the settlement state of the structure and predicting future trends. By generating high-resolution and high-precision settlement field functions, the analysis depth and prediction capabilities of long-distance non-contact structural settlement monitoring are improved, providing strong decision-making support for structural health management and risk assessment.

[0127] In another embodiment of the present application, the settlement field function S(x, y, t) and the historical settlement field function set {S(x, y, τ) |τ∈[0, t)} stored in the system database are obtained. The continuous function is discretized into a matrix X = [x1, x2, ..., x n ], where x i Indicates time point t i The state vector of the sedimentation field. Construct the time offset matrix X1 = [x1, x2, ..., x( n-1 )] and X2 = [x2, x3, ..., x n ]. Calculate the SVD decomposition of the time offset matrix X1 and X2: X1 = U * Σ *V T . Select the first r singular values ​​and the corresponding singular vectors to construct the truncated matrix U r ,Σ r and V r , calculate the projection matrix P = U r T .

[0128] Read the projection matrix P and the constructed time offset matrix pair X1 and X2. Calculate the matrix A tilde =P*X2*V r *Σr -1 . For the matrix A tilde Perform eigenvalue decomposition: A tilde *W =W*Λ, where Λ is the eigenvalue diagonal matrix and W is the eigenvector matrix. Calculate the dynamic mode matrix Φ = X2*V r *Σ r -1 *W. Normalize each column of the dynamic pattern matrix Φ to obtain the standardized dynamic pattern set Φ norm = {φ1, φ2, ..., φ r}. Calculate each dynamic mode φ i The growth rate λ i = log(|Λ ii |) / Δt, where Δt is the time step. Store the growth rate in the vector Λ growth middle.

[0129] Get the growth rate vector Λ growth and dynamic pattern set Φ norm Calculate the energy contribution E of each dynamic mode i = ||φ i || 2 *|Λ ii | 2 . Sort the dynamic modes according to their energy contribution and select the top k main modes. Based on the top k main dynamic modes and the corresponding growth rates, for each main mode φ i , calculate its spatial correlation C i (x, y) = corr(S(x, y,:), φ i ), where corr represents the time series correlation coefficient, S(x, y, :) represents the value of the sedimentation field function at the position (x, y) at all time points, and : represents all possible values. Generate a spatial correlation atlas {C i (x, y)}.

[0130] Get Spatial Correlation Atlas {C i (x, y)} and the main dynamic mode set. The dynamic mode, growth rate and spatial correlation information are integrated to generate the dynamic feature descriptor F D The output dynamic pattern set Φ = {φ1, φ2, ..., φ k}, the corresponding growth rate Λ = {λ1, λ2, ..., λ k} and feature descriptor F D For subsequent analysis.

[0131] This embodiment successfully extracts the spatiotemporal dynamic characteristics of sedimentation through the dynamic mode decomposition algorithm, providing an advanced analytical tool for in-depth understanding of the sedimentation mechanism and predicting future trends. By discretizing the continuous sedimentation field function into a high-dimensional data matrix, the complex spatiotemporal dynamics are successfully compressed into an analyzable numerical form, laying the foundation for subsequent mode decomposition. By constructing a time offset matrix pair and performing singular value decomposition, the main dynamic modes in the sedimentation process are effectively captured, reducing the dimension of the data while retaining key information. In particular, by calculating the characteristic decomposition of the dynamic matrix, the main dynamic modes and their evolutionary characteristics in the sedimentation process are successfully identified, providing important clues for understanding the internal mechanism of sedimentation. By calculating the growth rate of each dynamic mode, not only the dominant sedimentation mode is identified, but also the development trend of each mode is quantified, providing a scientific basis for predicting future sedimentation behavior. By analyzing the spatial correlation of the dynamic mode, the dynamic characteristics of the time domain are successfully mapped to the spatial domain, realizing the precise positioning of the sedimentation hotspot area. This embodiment can extract dynamic modes with physical significance from long-term series sedimentation data, providing a new perspective for in-depth understanding of the complex dynamic behavior of the sedimentation process. By revealing the main modes and evolution laws of settlement, the analysis depth and prediction ability of long-distance non-contact structural settlement monitoring are improved, providing important theoretical support for the formulation of long-term structural maintenance strategies and risk management plans.

[0132] In another embodiment of the present application, a dynamic pattern set Φ={φ1, φ2, ..., φ k}, growth rate Λ={λ1, λ2,..., λ k} and the current settlement field function S(x, y, t). Define the observation function set g = {g1, g2, ..., g m}, including polynomial functions and radial basis functions.

[0133] Based on the observation function set g and the current settlement field function S(x, y, t), calculate the observation vector y(t) = [g1(S), g2(S), ..., g m (S)] T . Use the least squares method to project y(t) into the dynamic pattern space and obtain the pattern coefficient vector a(t). Based on the pattern coefficient vector a(t) and the growth rate Λ, construct the finite-dimensional approximation matrix K of the Koopman operator, whose diagonal elements are exp(λ i *Δt). Calculate the mode coefficient a(t+Δt) = K * a(t) at the future time.

[0134] Based on the future model coefficient a(t+Δt) and the dynamic model set Φ, the future sedimentation field is reconstructed by linear combination: S'(x, y, t+Δt) = Σ i ai (t+Δt) * φ i (x, y). Calculate the gradient ▽S' and Laplacian ▽ of the predicted sedimentation field 2 S'. Calculate the difference between the predicted sedimentation field and the current sedimentation field ΔS = S' - S. Identify the significant change area in the difference ΔS, that is, |ΔS| > μ ΔS + 2σ ΔS Apply local polynomial regression to the significant change areas to improve the prediction accuracy in these areas. Update the predicted settlement field to obtain the optimized S''(x, y, t+Δt).

[0135] Based on the optimized predicted settlement field S''(x, y, t+Δt) and the current settlement field S(x, y, t). Calculate the predicted confidence interval CI(x, y) = [S''(x, y, t+Δt) -ε(x, y), S''(x, y, t+Δt) +ε(x, y)], where ε(x, y) is estimated based on the local prediction error. Based on the predicted confidence interval CI(x, y) and the optimized predicted settlement field S''(x, y, t+Δt), generate a set of contour lines of the predicted settlement field {C' i | S''(x, y, t+Δt) = h i}. Compare the predicted contours with the current contours to identify potential areas of abnormal development.

[0136] Obtain potential abnormal development area information and predicted settlement field S''(x, y, t+Δt). Comprehensively predict settlement field, gradient field, Laplacian field, confidence interval and abnormal development area information to generate prediction feature descriptor F P Output predicted settlement field function S''(x, y, t+Δt), confidence interval CI(x, y) and feature descriptor F P For subsequent analysis and decision support.

[0137] This embodiment achieves high-precision prediction of complex settlement dynamics through a nonlinear prediction method based on Koopman operator theory, providing strong technical support for early assessment and prevention of potential risks. By defining an observation function set containing polynomial functions and radial basis functions, the nonlinear settlement dynamics is successfully converted into a high-dimensional linear system, laying a theoretical foundation for subsequent prediction analysis. By projecting the current settlement field into the dynamic mode space, the dimension of the problem is effectively reduced, while retaining key dynamic information, improving computational efficiency and prediction accuracy. In particular, by constructing a finite-dimensional approximate matrix of the Koopman operator, the nonlinear evolution characteristics of the system can be accurately captured, overcoming the limitations of traditional linear prediction methods. By applying exponential mapping to predict future mode coefficients, long-term prediction of settlement dynamics is successfully achieved, providing an important basis for evaluating the long-term stability of the structure. By introducing local polynomial regression and prediction confidence intervals, not only the local accuracy of the prediction is improved, but also the uncertainty of the prediction results is quantified, providing comprehensive information support for risk assessment. This embodiment can achieve high-precision, long-term scale predictions in complex nonlinear settlement systems, creating conditions for timely discovery of potential structural risks. By accurately predicting future settlement trends and patterns, the foresight and preventive nature of long-distance non-contact structural settlement monitoring is improved, providing a scientific basis for formulating proactive maintenance strategies and emergency plans, and improving the safety management level of important infrastructure.

[0138] like Figure 5 As shown, according to one aspect of the present application, step S4 is further:

[0139] S41, based on the settlement field function and the predicted settlement field, a linear combination of orthogonal polynomial basis functions is expanded; based on the linear combination, the expansion coefficient is calculated by using the Galerkin projection method, and finally an uncertainty field is obtained;

[0140] S42, based on the uncertainty field and the preprocessed image set, performing wavelet packet decomposition to obtain a predetermined number of sub-bands; calculating the energy and the cross-correlation coefficient of each sub-band, and selecting a significant sub-band through a thresholding method based on the energy and the cross-correlation coefficient of each sub-band; based on the significant sub-bands, reconstructing error components of multiple scales, and finally forming an error scale coefficient set;

[0141] S43, discretizing the settlement field function into a state vector, constructing a state transfer equation and an observation equation based on the state vector; dynamically adjusting the covariance matrix of the noise in the state transfer equation and the observation equation based on the error scale coefficient set to obtain the adjusted state transfer equation and the observation equation; performing Kalman filtering based on the adjusted state transfer equation and the observation equation to iteratively optimize the state estimation; converting the state estimation into a function form to obtain an optimized settlement field function;

[0142] S44. Based on the optimized sedimentation field function, a high-resolution sub-grid is constructed; based on the high-resolution sub-grid, the sedimentation data of each sub-grid area is fitted to obtain a local quadratic polynomial fitting function; based on the local quadratic polynomial fitting function, the polynomial coefficients are estimated using the least squares method; based on the polynomial coefficients, the sedimentation values ​​of the sub-grid points are calculated; based on the grayscale information of the preprocessed image set, the sedimentation values ​​are fine-tuned to obtain the fine-tuned sedimentation values; based on the fine-tuned sedimentation values, a high-precision sedimentation field function is generated.

[0143] In one embodiment of the present application, the settlement field function S(x, y, t) and the predicted settlement field S'(x, y, t+Δt) are read, and the uncertainty of the settlement estimation is quantified by applying the polynomial chaos expansion method: the uncertainty of the input parameters is expressed as a random variable set ξ={ξ1,ξ2,...,ξn}; the settlement field function is expanded into a linear combination of orthogonal polynomial basis functions Ψi(ξ): S(x, y, t, ξ)≈Σ i=0 to P αi(x, y, t)Ψi(ξ), where αi is the expansion coefficient. The expansion coefficient is calculated by the Galerkin projection method, and the uncertainty field U(x, y, t) is finally obtained, which represents the variance or confidence interval of the settlement estimate at each position (x, y) at time t.

[0144] Get the uncertainty field U(x, y, t) and the registered image I r (t), the wavelet packet transform algorithm is used to perform multi-scale error decomposition: for the uncertainty field U(x, y, t) and Figure I r (t) Perform wavelet packet decomposition to obtain multiple subbands. Calculate the energy and cross-correlation coefficient of each subband to evaluate its contribution to the overall error. Then, select the significant subbands through the thresholding method and reconstruct the error components of multiple scales. The final output error scale coefficient set {ε1, ε2, ..., ε L}, where ε i Represents the error contribution weight of the i-th scale.

[0145] Read the sedimentation field function S(x, y, t) and the error scale coefficient set {ε1, ε2, ..., ε L}, construct and apply adaptive Kalman filter to optimize the settlement field; discretize the settlement field function into state vector x(t). Construct state transfer equation x(t+1) = F(t)x(t) + w(t) and observation equation z(t) = H(t)x(t) + v(t), where F(t) is the state transfer matrix, H(t) is the observation matrix, w(t) and v(t) are process noise and observation noise, respectively. According to the error scaling coefficients {ε1, ε2, ..., ε LDynamically adjust the noise covariance matrix Q(t) and R(t). Perform the prediction and update steps of the Kalman filter and iteratively optimize the state estimation. Finally, the optimized settlement field function S*(x, y, t) is obtained.

[0146] Get the optimized sedimentation field function S*(x, y, t) and the registered image I r (t), the spatial resolution of the sedimentation estimate is improved using a sub-pixel interpolation algorithm based on surface fitting: a high-resolution sub-grid is defined around the original grid points. For each sub-grid area, a local quadratic polynomial function f(x, y) = ax 2 + by 2 + cxy + dx + ey + f fits the sedimentation data. The polynomial coefficients {a, b, c, d, e, f} are estimated by the least squares method. The sedimentation values ​​of the subgrid points are calculated using the fitted polynomial function. At the same time, the interpolation results are fine-tuned in combination with the grayscale information of the image Ir(t). Finally, a high-precision sedimentation field function S**(x, y, t) is output, which has higher spatial resolution and accuracy.

[0147] This embodiment realizes error analysis and precision improvement of settlement estimation through algorithms such as polynomial chaos expansion, wavelet packet transform, adaptive Kalman filtering and sub-pixel interpolation of surface fitting. The polynomial chaos expansion method effectively quantifies the uncertainty of settlement estimation by expressing the uncertainty of input parameters as a linear combination of orthogonal polynomial basis functions, providing a reliable basis for risk assessment. The wavelet packet transform algorithm accurately identifies the error contribution of different frequency components through multi-scale error decomposition, pointing out the direction for targeted improvement of monitoring accuracy. The adaptive Kalman filter combines the error scale coefficient and dynamically adjusts the covariance matrix of process noise and observation noise, improves the adaptability and robustness of filtering, and effectively suppresses the influence of various interference factors. The sub-pixel interpolation algorithm of surface fitting successfully improves the spatial resolution of the settlement field to the sub-pixel level through high-order polynomial fitting and image gradient consistency constraints, and enhances the detection capability of micro-sedimentation. This embodiment not only improves the accuracy and reliability of settlement estimation, but also provides comprehensive error analysis and uncertainty quantification, laying a solid foundation for the scientific interpretation and application of monitoring results. Through this series of error analysis and accuracy improvement measures, highly reliable settlement monitoring results can be provided in complex actual monitoring environments, providing strong technical support for the safe management of important engineering structures.

[0148] In one embodiment of the present application, the settlement field function S(x, y, t) and the predicted settlement field S''(x, y, t+Δt) are obtained. A set of random variables ξ = {ξ1, ξ2, ..., ξn} is defined to represent the uncertainty of the input parameters. M sample points {ξ(i)} are sampled from a predefined probability distribution. i=1 M . Choose an orthogonal polynomial basis function set {Ψj(ξ)} j=0 P , where P is the truncation order. Calculate the value of the basis function at the sampling point Ψj(ξ(i)) and construct the basis function matrix Ψ.

[0149] Based on the basis function matrix Ψ and the settlement field function S(x, y, t). For each spatial point (x, y), calculate the function S(x, y, t, ξ(i)), i = 1, ..., M. Construct the response matrix R, where R ij = S(x i , y i , t, ξ(j)). Use the least squares method to solve the expansion coefficient: α= (Ψ T Ψ) -1 Ψ T R. Get the expansion coefficient vector α(x, y) = {α0(x, y), α1(x, y), ..., α P (x, y)}.

[0150] Based on the expansion coefficient α(x, y), calculate the mean μ of the sedimentation field S (x, y) = α0(x, y) and variance σ S 2 (x, y) = Σ j=1 P α j 2 (x, y) E[Ψ j 2 ], where E[Ψ j 2 ] is the second-order moment of the basis function. Construct the 95% confidence interval of the sedimentation field CI(x, y) = [μ S (x, y) -1.96σ S (x, y), μ S (x, y) + 1.96σ S (x, y)]. Calculate the width of the confidence interval W(x, y) = 3.92σ S (x, y).

[0151] Based on the confidence interval width W(x, y), calculate the global average width Wavg =∫∫W(x, y)dxdy / A, where A is the area of ​​the study area. Identify W(x, y) > W avg The area with + 2σW is regarded as the high uncertainty area H. Based on the expansion coefficient α(x, y), for each point in the identified high uncertainty area H, the Sobol sensitivity index Si = Var(E[S|ξi]) / Var(S) is calculated to identify the main contributing sources of uncertainty.

[0152] Get Sobol sensitivity index and mean field μ S (x, y), variance field σ S 2 (x, y), generate uncertainty field descriptor F U , including mean, variance, high uncertainty area and main uncertainty sources. Output uncertainty field U(x, y, t) = {μ S , σ S 2 , CI, H, Si} and descriptor F U . Where ξ represents a random variable; Ψ j represents the orthogonal polynomial basis function; α represents the expansion coefficient; μ S represents the mean value of the sedimentation field; σ S 2 represents the variance of the sedimentation field; CI represents the confidence interval; W represents the width of the confidence interval; H represents the high uncertainty area; Si represents the Sobol sensitivity index; F U Represents an uncertainty field descriptor.

[0153] In this embodiment, the polynomial chaos expansion method is used to comprehensively quantify the uncertainty of settlement estimation, providing a strong theoretical support for reliability analysis and risk assessment. By representing the uncertainty of input parameters as a set of random variables, the complex uncertainty sources are successfully modeled, laying a foundation for subsequent probabilistic analysis. By selecting an appropriate set of orthogonal polynomial basis functions, the random process is effectively expanded into a linear combination of deterministic polynomials, simplifying the computational complexity. In particular, by using the least squares method to solve the expansion coefficients, an efficient and accurate uncertainty propagation analysis is achieved, overcoming the problem of high computational cost of the traditional Monte Carlo method. By calculating the mean and variance of the settlement field, not only the expected value of settlement prediction is provided, but also the degree of dispersion of the prediction results is quantified, providing an important index for comprehensively evaluating the reliability of monitoring results. By constructing confidence intervals and identifying high-uncertainty regions, the risk regions that need to be focused on are successfully located, providing a scientific basis for optimizing the monitoring strategy. This embodiment can accurately evaluate the reliability and risk level of settlement prediction results considering multiple uncertainty sources, providing comprehensive and detailed uncertainty information for decision-makers. By systematically quantifying and analyzing the uncertainty of prediction, the credibility and practicality of long-distance non-contact structural settlement monitoring are improved, providing a solid theoretical foundation for formulating risk-oriented monitoring and maintenance strategies.

[0154] In one embodiment of the present application, the uncertainty field U(x, y, t) and the registered image I r (t) are obtained. An appropriate wavelet basis function ψ(x), such as the Daubechies wavelet, is selected. The maximum decomposition level L and the integrity parameter p of the wavelet packet tree are determined. An initial two-dimensional discrete wavelet transform is performed on the uncertainty field U(x, y, t) to obtain the approximation coefficient A1 and the detail coefficients {H1, V1, D1}. An initial wavelet packet tree node set T1 = {A1, H1, V1, D1} is constructed. For each node in the initial wavelet packet tree node set T1, its energy E(node)=Σ|coef| 2 and entropy S(node)=-Σ|coef| 2 log(|coef| 2 ) are calculated. It is decided whether to further decompose the node according to the energy and entropy. If E(node) > p*E(parent) or S(node)<p*S(parent), the node is further decomposed. This process is repeated until the maximum level L is reached or no node needs to be further decomposed. Here, coef represents the coefficient set of the current node, E(parent) represents the energy of the parent node, and S(parent) represents the entropy of the parent node.

[0155] The wavelet packet tree after completion of decomposition is obtained. The normalized energy E of each leaf node is calculatednorm (node) = E(node) / E total , where E total is the total energy. Select leaf nodes whose energy ratio exceeds the threshold τ to form a significant node set S. For each node in the significant node set S, calculate its support region R(node) in the original space. Construct a feature graph F(x, y), where F(x, y) = Σ( node∈S ) E norm (node) * χ R (node)(x, y), χ is the characteristic function. Calculate the feature map F(x, y) and the registered image I r The correlation coefficient of (t)ρ(x, y) = Cov(F, I r ) / (σF *σI r ). Identify the region where ρ(x, y) > μρ+ 2σρ as the high correlation region C. For each point (x, y) in the high correlation region C, extract the corresponding significant wavelet packet coefficient {coef i (x, y)}. Using these coefficients, the local signal is reconstructed to obtain the fine structure field D(x, y).

[0156] Based on the fine structure field D(x, y) and the feature map F(x, y), calculate the multi-scale error index ME(l) = || D - F l || 2 / ||F l || 2 , where F l is the approximation of F at the lth layer. The output error scale coefficient set {ε l = ME(l) / Σ k ME(k), l = 1, ..., L} and the fine structure field D(x, y).

[0157] Where ψ(x) represents the wavelet basis function; L represents the maximum decomposition level; p represents the wavelet packet tree integrity parameter; A1, H1, V1, D1 all represent wavelet transform coefficients; E(node) represents node energy; S(node) represents node entropy; τ represents energy threshold; ρ represents the mutual correlation coefficient; ME represents the multi-scale error index; ε l represents the error scale coefficient.

[0158] This embodiment realizes multi-scale decomposition and refined analysis of the sedimentation field error through the wavelet packet transform algorithm, which provides important technical support for improving monitoring accuracy and optimizing monitoring strategies. By selecting appropriate wavelet basis functions and determining the optimal decomposition layer number, the complex error field is successfully decomposed into components of different frequencies and spatial scales, creating conditions for a comprehensive understanding of the error structure. By calculating the energy and entropy of each wavelet packet node, the significant components in the error are effectively identified, which improves the pertinence and efficiency of subsequent analysis. In particular, by adaptively selecting significant nodes, the amount of data can be greatly reduced while retaining key information, and efficient extraction of error features is achieved. By constructing a feature map and calculating the mutual correlation coefficient with the original image, the error features are successfully associated with the actual observation data, providing important clues for understanding the source of the error. By reconstructing local signals and calculating multi-scale error indicators, not only the quantitative description of the error is realized, but also the distribution characteristics of the error at different scales are revealed, which points out the direction for targeted improvement of monitoring accuracy. This embodiment can extract error patterns and scale characteristics with physical significance from complex error fields, providing a scientific basis for optimizing monitoring systems and improving data processing algorithms. By deeply analyzing the multi-scale structure of errors, the accuracy assessment and adaptive optimization capabilities of long-distance non-contact structural settlement monitoring are improved, creating conditions for continuously improving the performance and reliability of the monitoring system.

[0159] In one embodiment of the present application, the settlement field function S(x, y, t) and the error scale coefficient set {ε l}. Discretize the continuous settlement field function into a state vector x(t) = [S(x1, y1, t), S(x2, y2, t), ..., S(xN, yN, t)] T , where (xi, yi) are discrete grid points. Construct the state transfer equation x(t+1) = F(t)x(t) + w(t), where F(t) is the state transfer matrix, initialized to the identity matrix I. w(t) is the process noise, which is assumed to obey a Gaussian distribution with a mean of 0 and a covariance of Q(t). Based on the state transfer equation and the fused high-quality image I f (t), construct the observation equation z(t) = H(t)x(t) + v(t), where H(t) is the observation matrix and v(t) is the observation noise, which is assumed to obey a Gaussian distribution with a mean of 0 and a covariance of R(t).

[0160] Read the observation equation and error scale coefficient set {ε l Initialization process noise covariance matrix Q(0) = diag({ε l}) and the observation noise covariance matrix R(0) =σ 2 *I, where σ 2is the estimated observation error variance.

[0161] Based on the initialized covariance matrices Q(0) and R(0), perform the prediction step of the Kalman filter: x'(t+1|t)=F(t)x'(t|t), P(t+1|t) = F(t)P(t|t)F(t) T + Q(t), where x' is the state estimate and P is the estimation error covariance matrix. Calculate the Kalman gain K(t+1) = P(t+1|t)H(t+1) T [H(t+1)P(t+1|t)H(t+1) T + R(t+1)] -1 .

[0162] Get the Kalman gain K(t+1) and the observed data z(t+1). Perform the update step of the Kalman filter: x'(t+1|t+1) = x'(t+1|t) + K(t+1)[z(t+1) - H(t+1)x'(t+1|t)], P(t+1|t+1) = [I - K(t+1)H(t+1)]P(t+1|t).

[0163] Based on the update results x'(t+1|t+1) and P(t+1|t+1), calculate the normalized innovation sequence d(t+1) = [z(t+1) - H(t+1)x'(t+1|t)] / sqrt(H(t+1)P(t+1|t)H(t+1) T + R(t+1)). Check the white noise characteristics of the innovation sequence d(t+1). If the innovation sequence d(t+1) does not meet the white noise assumption, the matrices Q(t+1) and R(t+1) are adaptively adjusted: Q(t+1) = (1-α)Q(t) + α[K(t+1)d(t+1)d(t+1) T K(t+1) T ], R(t+1) = (1-β)R(t) + β[d(t+1)d(t+1) T ], where α and β are adaptive factors. Repeat until the filter converges. Output the optimized settlement field function S*(x, y, t) = reshape(x'(t+1|t+1), [sqrt(N), sqrt(N)]) and the estimated error covariance P(t+1|t+1). Where Q(t) is the process noise covariance matrix; R(t) is the observation noise covariance matrix.

[0164] This embodiment realizes high-precision estimation and dynamic optimization of the settlement field through the construction and application of adaptive Kalman filter, and provides key algorithmic support for providing stable and reliable monitoring results. By discretizing the continuous settlement field function into a state vector, the complex spatial distribution problem is successfully converted into a processable state estimation problem, laying the foundation for the application of Kalman filtering technology. By constructing the state transfer equation and observation equation containing process noise, the dynamic characteristics and measurement uncertainty of the settlement process are effectively simulated, and the authenticity and applicability of the model are improved. In particular, by introducing the error scale coefficient to dynamically adjust the noise covariance matrix, adaptive processing of errors of different scales and sources is achieved, and the robustness and accuracy of the filter are improved. By calculating the Kalman gain and performing state updates, the model prediction and actual observation are successfully integrated to achieve the optimization of settlement field estimation. By analyzing the normalized innovation sequence and adaptively adjusting the process noise and observation noise covariance, it can dynamically adapt to changes in system characteristics and ensure the stability and reliability of long-term monitoring. This embodiment can provide high-precision, low-noise settlement field estimation in a complex and changeable monitoring environment, providing reliable data support for structural health assessment. By adaptively optimizing the filtering process, the anti-interference ability and long-term stability of long-distance non-contact structural settlement monitoring are improved, creating conditions for achieving continuous and high-quality structural monitoring.

[0165] In one embodiment of the present application, the optimized sedimentation field function S*(x, y, t) and the registered image I are obtained. r (t). Define a high-resolution subgrid with a grid spacing of 1 / k of the original grid, where k is the super-resolution factor. Initialize the high-resolution settlement field S**(x, y, t). For each original grid cell, extract the settlement values ​​of (2m+1)×(2m+1) adjacent grid points to form a local settlement data set D = {(xi, yi, S*(xi, yi, t))}. Construct a quadratic polynomial fitting function f(x, y) = ax 2 + by 2 + cxy + dx + ey + f. Use the least squares method to estimate the coefficient vector θ = [a, b, c, d, e, f] T , solve the equation (X T X)θ = X T S, where X is the design matrix of coordinates in D and S is the corresponding sedimentation value vector.

[0166] For each point (x', y') in the high-resolution subgrid, calculate its fitted sedimentation value S fit (x', y', t) = f(x', y') = [x' 2 , y' 2, x'y', x', y', 1]θ. The fitted sedimentation value S fit (x', y', t) is stored in the high-resolution sedimentation field S**(x, y, t).

[0167] Obtain high-resolution fitting sedimentation field S**(x, y, t) and registered image I r (t). For the registered image I r (t) Perform bicubic interpolation to obtain a high-resolution image I rhr (t). Calculate the local image gradient ▽I rhr (x', y', t) = [ΨI / Ψx, ΨI / Ψy], where Ψ represents the partial derivative. Construct the consistency constraint of the sedimentation field gradient: E grad (x', y') = ||▽S**(x', y', t) - λ▽I rhr (x', y', t)|| 2 , where λ is the scaling factor.

[0168] Based on the gradient consistency constraint E grad (x', y'). Define the overall energy function E total (x', y') = E fit (x', y')+γE grad (x', y'), where E fit (x', y') = (S**(x', y', t) - S* fit (x', y', t)) 2 , γ is a trade-off parameter. Use gradient descent to minimize E total (x', y'), update S**(x', y', t). Read the optimized high-resolution sedimentation field S**(x, y, t), calculate the local curvature κ(x', y') = |▽ 2 S**(x', y', t)| / (1 + ||▽S**(x', y', t)|| 2 ) 3 / 2 Identify the region where κ(x', y') > μκ + 2σκ as the high curvature region H κ .

[0169] Get the identified high curvature area H κ And the optimized high-resolution sedimentation field S**(x, y, t). In the high curvature area H κAnisotropic diffusion filtering is applied to the image to reduce noise while maintaining edge sharpness: ΨS** / Ψt = div(c(||▽S**||)▽S**), where c(·) is the diffusion coefficient function. Iterate until convergence or the maximum number of iterations is reached.

[0170] Read the final high-resolution sedimentation field S**(x, y, t). Calculate the root mean square error RMSE and structural similarity index SSIM with the original sedimentation field S*(x, y, t). Generate a high-precision sedimentation field descriptor F H , including RMSE, SSIM, high curvature area H κ Output high-precision sedimentation field function S**(x, y, t) and descriptor F H .

[0171] Where S*(x, y, t) represents the optimized sedimentation field function; S**(x, y, t) represents the high-precision sedimentation field function; k represents the super-resolution factor; m represents the local fitting window radius; f(x, y) represents the quadratic polynomial fitting function; θ represents the polynomial coefficient vector; X represents the design matrix; S* fit Represents the fitted sedimentation value; I rhr represents a high-resolution image; ▽ represents a gradient operator; E grad represents the gradient consistency constraint; E fit represents the fitting error; E total represents the overall energy function; γ represents the trade-off parameter; κ represents the local curvature; H κ represents the high curvature region; div represents the divergence operator; c(·) represents the diffusion coefficient function; RMSE represents the root mean square error; SSIM represents the structural similarity index; F H Represents a high-precision sedimentation field descriptor.

[0172] This embodiment achieves the improvement of the spatial resolution of the sedimentation field through the sub-pixel interpolation algorithm of surface fitting, and provides important technical support for the refined analysis of structural sedimentation. By defining high-resolution subgrids and extracting local sedimentation data sets, discrete sedimentation data are successfully converted into continuous high-precision sedimentation fields, creating conditions for capturing small sedimentation changes. By constructing a quadratic polynomial fitting function and using the least squares method to estimate the coefficients, the noise in the original data is effectively smoothed, while retaining important structural features and improving the reliability of the interpolation results. In particular, by introducing image gradient information and constructing gradient consistency constraints, the detail information of high-resolution images is successfully integrated into the sedimentation field reconstruction process, improving the authenticity and accuracy of the interpolation results. By defining the overall energy function and applying the gradient descent method to optimize, the best balance between fitting accuracy and gradient consistency is achieved, ensuring the physical rationality of the reconstruction results. By identifying high curvature areas and applying anisotropic diffusion filtering, important structural features in the sedimentation field are successfully retained, while noise and artifacts are effectively suppressed. This embodiment can reconstruct a high-resolution, high-precision sedimentation field from raw data of limited resolution, providing a solid data foundation for the refined analysis of structural sedimentation. By improving the spatial resolution of the settlement field, the analytical capability and application value of long-distance non-contact structural settlement monitoring are enhanced, creating conditions for early detection of minor settlement anomalies and formulation of precise maintenance strategies.

[0173] According to another aspect of the present application, a non-contact remote monitoring method for structural settlement comprises the following steps:

[0174] S0. Method Overview. By establishing a data model in advance, accurate output of structural settlement is achieved. That is, a big data model of "image displacement" and "actual settlement" is established: F (y1, y2, U), where y1 is the image displacement difference; y2 is the actual settlement; U is the observation parameter group (distance, inclination, etc.).

[0175] S1. Experimental design. A series of experiments are designed using a special experimental design method. The confidence interval of y2 is 0~100mm. The number of experimental groups is not less than 100 (n≥100).

[0176] S2. Conduct the test. Select a suitable test site (with the required observation and working condition change conditions), conduct the test in sequence according to the working conditions of the test design, collect image displacement, and fill in the table.

[0177] S3. Establish a model. With [observation parameters, actual settlement] as input and [image displacement] as output, establish a data model: the model reflects the mapping relationship between input data and output data.

[0178] S4. Accuracy test. Input a set of observation parameters and a settlement distance y2 test (0~100mm) into the model at random to obtain a calculated image displacement Vmeter. Then, based on the observation parameters and the settlement distance y2 test, conduct field tests to obtain the displacement difference Vtrue between the two images. Accuracy requirement: (Vmeter-Vtrue) / Vtrue ≤5%; if the accuracy does not meet the requirements, the number of test groups should be increased, or the inversion accuracy of the model should be improved.

[0179] S5. Real application: Input “observation parameters” and “image displacement difference” into the data model to obtain real settlement data.

[0180] According to another aspect of the present application, a non-contact remote monitoring method for structural settlement comprises the following steps:

[0181] Step 1: Deploy the infrared optical observation system.

[0182] The infrared optical observation system consists of three parts: reinforced concrete pier; lifting control and rotating platform; infrared optical observation equipment.

[0183] When selecting a site for the observation platform, the following should be noted: there should be no physical obstructions between the observation platform and the observation target (single target or multiple targets); in order to ensure the observation accuracy of this application, the straight-line distance of observation should be controlled within 1500m; the location of the observation platform should be convenient for personnel to reach; the foundation should be in a basically stable state.

[0184] System component 1: Concrete pier. The height shall not be less than 1.5m (to the ground where personnel are located), and the cross-sectional area shall not be less than 0.3m×0.3m; a forced centering plate shall be installed on the top of the concrete pier; the concrete pier foundation shall be constructed to ensure the stability of the concrete pier, and the construction shall be carried out in accordance with the relevant standards of the "temporary leveling point". It is recommended to adopt relevant measures such as deep burial, pile foundation, and strong tamping to ensure the stability of the observation platform; protective measures shall be set up for the observation platform, and a movable colored steel house can be built as a protective measure, and the front end can be freely opened and closed to ensure that the observation line of sight is not affected.

[0185] System component 2: lifting control and rotating table. The top of the concrete pier is forced to be aligned with the centering plate, and the Z-axis lifting control console and the R-axis rotating platform are connected to the optical observation instrument. The moving accuracy of the lifting control console and the rotating table is required to be less than 0.02mm.

[0186] System component 3: Infrared optical observation instrument. It is used to observe the target object. The technical parameters of the infrared optical observation instrument should meet the following requirements: when ① the temperature difference between the observed object and the ambient temperature is 50 degrees Celsius; ② the distance from the observation platform is 1200m; ③ the observed object is a 10cm×10cm heat source, the image resolution is not less than 1024×768.

[0187] Step 2: Lay out the heat source matrix system.

[0188] On the observed target object, a heat source matrix system is deployed, which mainly includes three parts:

[0189] System component 1: heat source matrix. Silicone heating sheets are pasted on the insulation board in an array. The heating sheets can be heated to a specified temperature after being powered on. The heating sheets are square, and the size of each heating sheet should not be less than 10cm×10cm, and the interval should not be less than 5cm. The size and interval of the heating sheets can be increased according to specific observation needs. The thickness of the insulation board is greater than 1cm, and the edge of the insulation board should be greater than 5cm away from the heating sheet (at least 5cm is reserved). The heat source matrix is ​​hung on the observation target, and the hanging height is located in the middle of the structure. It should ensure that the line of sight with the observation platform is unobstructed, there is no physical obstruction, and the straight-line observation distance is less than 1500m.

[0190] System component 2: Intelligent thermostat. The main output line of the heat source matrix is ​​connected to the thermostat; the temperature sensor is connected to one of the heating plates; the maximum setting temperature of the thermostat should be greater than or equal to 100 degrees Celsius; in order to ensure the stability of observation, it is necessary to monitor the ambient temperature (this temperature refers to the ambient temperature near the observation object during the observation time window):

[0191] When the ambient temperature is between -20 and -10 degrees Celsius, the thermostat is set to 30 degrees Celsius;

[0192] When the ambient temperature is between -10 and 0 degrees Celsius, the thermostat is set to 40 degrees Celsius;

[0193] When the ambient temperature is between 0 and 10 degrees Celsius, the thermostat is set to 50 degrees Celsius;

[0194] When the ambient temperature is between 10 and 20 degrees Celsius, the thermostat is set to 60 degrees Celsius;

[0195] When the ambient temperature is between 20 and 30 degrees Celsius, the thermostat is set to 70 degrees Celsius;

[0196] When the ambient temperature is between 30 and 40 degrees Celsius, the thermostat is set to 80 degrees Celsius;

[0197] When the ambient temperature is greater than 40 degrees Celsius, the thermostat is directly set to 100 degrees Celsius.

[0198] System component 3: Solar remote control power supply system. Use solar panels as power source, the power should be greater than or equal to 60 watts; equip the power supply system with a remote communication remote switch, which can be used to remotely control the power supply on the mobile phone; set up a waterproof box for the thermostat and communication switch, and hang it near the equipment.

[0199] Step 3: Conduct observations and perform image analysis.

[0200] Observation method: Use mobile phone applications to control infrared optical observation equipment and watch the monitoring screen on the mobile phone screen in real time; take photos and videos on the mobile phone or other mobile device port (direct operation on the instrument is prohibited to avoid affecting the observation accuracy); observe at least once every 24 hours; when observing multiple targets, the total time for each observation should be controlled within 30 minutes. When the observation time is long, a higher-power solar power supply should be equipped; after the observation is completed, remotely control the switch equipment and turn off the power supply.

[0201] The observation method of this embodiment has great adaptability to weather conditions and brightness of day and night, and the specific observation time window can be customized; it is recommended to conduct observations in combination with the ambient temperature to maintain the consistency and stability of monitoring sampling; when observing multiple objects, the displacement readings of the lifting console and the rotating platform should be accurately recorded and controlled, and accurately returned to the original position when repeated observations are made.

[0202] Image processing method: The observation data (i.e., the acquired infrared digital image) is transferred to the computer via the mobile phone for storage, and the target object number and shooting time are marked for the image file; image analysis is performed on the computer, and the digital image correlation (DIC) analysis algorithm is used for image comparison and analysis, so as to conduct regular monitoring of structural settlement.

[0203] In one embodiment of the present application, in a certain intercity high-speed railway construction project, the high-speed railway line crosses a river of about 2,400 m, and there are 45 high-speed railway bridge piers in the construction area. A temporary construction traffic bridge is built on the construction site, and the bridge piers will be removed after the construction is completed. This embodiment intends to install observation platforms on both sides of the river to monitor the settlement of the bridge piers from a distance.

[0204] Step 1: Build a shore optical observation platform;

[0205] Taking the right bank upstream observation platform as an example, a site was selected on the right bank 100m upstream from the shore pier to build the observation platform. The farthest observation distance in this implementation case is about 1200m.

[0206] Connect the adapter plate to the forced centering disk, first connect the Z-axis lifting control console; then directly connect the R-axis rotating platform through the table; and then connect the infrared optical observation instrument through the adapter plate.

[0207] The infrared optical observation instrument adopts the existing high-precision infrared sight, model PARD-SA61-45L (45mm lens); the technical parameters of the infrared sight are as follows: when the distance is 1200m, it can observe a 10cm×10cm heat source, and the image resolution is not less than 1024×768.

[0208] In order to verify the accuracy of the equipment observation, this embodiment first carried out a feasibility test between the pedestrian bridge deck and the lakeshore of XX Lake Park in a certain city (about 1350m measured on the map). The instrument was fixed on the pedestrian bridge deck through a tripod, and the heat source matrix point (10cm×10cm heat source) was vertically fixed on the ground on the opposite bank. The observation results verified the accuracy of the instrument equipment.

[0209] Step 2: Lay out the pier heat source matrix system;

[0210] The heating material is made of adjustable temperature silicone rubber material and is customized for DC 24V power supply; this embodiment is equipped with 9 heating sheets, each of which is 10cm×10cm in size and 5cm in equal intervals; the power consumption of each heating sheet is 6W; there is 3M glue on the back of the heating sheet, and the heating sheet is directly pasted on the insulation board in an array, and the size of the insulation board is 50cm×50cm×1cm.

[0211] The main output line of the heat source matrix is ​​connected to the thermostat, the temperature of the thermostat is set to 70 degrees Celsius, and the output line of the thermostat is connected to the electrical end of the remote switch; the output line of the solar power supply is connected to the power end of the remote switch, the specification is DC 12V, and an inverter needs to be connected to provide DC 24V power.

[0212] After connecting the equipment, check that the power of the battery behind the solar panel is turned on, and check that the double-sided indicator lights of the remote control switch are normal; open the mobile cloud smart APP, connect the switch equipment, and control the power supply on the mobile phone. This operation can be used at any distance (only 4G network coverage is required).

[0213] Use four expansion screws to hang the heat source matrix on the right bank side of the pier, slightly upstream, at the middle of the pier. Install the solar panels on the surface of the pier (fixed with expansion screws) to ensure that the solar panels can effectively obtain light. Provide a waterproof box for the thermostat and remote control switch. The size of the waterproof box is 30cm×20cm×30cm and it is hung near the equipment.

[0214] In this embodiment, a heat source matrix system is deployed for the three bridge piers that are relatively farthest from the shore, and the observation distance is about 1120m~1200m.

[0215] Step 3: Conduct regular observations and image analysis;

[0216] Use the mobile phone application "PARD" to control the infrared optical observation equipment, and you can see the monitoring image on the mobile phone screen in real time; take photos and videos directly on the mobile phone; observe once at 8 o'clock every morning; for multi-object observation, turn the knob of the rotating platform to aim the instrument at the next pier target, and ensure that the displacement readings of the lifting platform and the rotating platform are consistent.

[0217] The observed infrared digital images are exported to the computer for storage and marked with the pier number and shooting time. Based on the DIC calculation module (graphic deformation measurement program) in MATLAB, the computer performs image calculation and analysis at different times to obtain the settlement data of the piers, thereby carrying out regular monitoring of the pier settlement.

[0218] According to one aspect of the present application, a non-contact remote monitoring system for structural settlement includes:

[0219] at least one processor; and,

[0220] a memory communicatively connected to at least one of the processors; wherein,

[0221] The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the non-contact long-distance monitoring method for structural settlement described in any one of the above embodiments.

[0222] The preferred embodiments of the present invention are described in detail above; however, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all belong to the protection scope of the present invention.

Claims

1. A non-contact remote monitoring method for structural settlement, characterized in that: The steps include: S1. Shooting the heat source matrix based on a multi-view synchronous acquisition method to obtain an original image set, and preprocessing the original image set to obtain a preprocessed image set; S2, based on the preprocessed image set, using an improved fast corner point detection algorithm to identify feature points and obtain a feature point set; Based on the feature point set, an improved optical flow tracking algorithm combined with time consistency constraints is used to track the feature points to obtain a set of vertical trajectories of the feature points. Based on the set of vertical trajectories of the feature points, an adaptive triangulation algorithm is applied to construct the topological relationship between the feature points to obtain a time-varying feature map. S3. Based on the time-varying feature graph, a spectral graph analysis algorithm is used to calculate the graph structure change index and obtain the graph topology change index set; based on the feature point vertical trajectory set and the time-varying feature graph, a local rigidity-preserving deformation algorithm is used to estimate the local vertical deformation of the feature point and obtain the local vertical deformation field; based on the local vertical deformation field and the graph topology change index set, a global settlement field equation is constructed and solved to obtain the settlement field function; based on the settlement field function, a dynamic mode decomposition algorithm is used to calculate the dynamic mode set and the corresponding growth rate; based on the dynamic mode set and the growth rate, a nonlinear operator prediction method is used to calculate the predicted settlement field at the future moment; S4. Based on the settlement field function and the predicted settlement field, the polynomial chaos expansion method is used to quantify the uncertainty of the settlement estimation and obtain the uncertainty field; based on the uncertainty field, the wavelet packet transform algorithm is used to perform multi-scale error decomposition and obtain the error scale coefficient set; based on the error scale coefficient set, an adaptive Kalman filter is constructed and applied to optimize the settlement field function and obtain the optimized settlement field function; based on the optimized settlement field function, a sub-pixel interpolation algorithm of surface fitting is used to calculate a high-precision settlement field function; Step S1 is further as follows: S11, acquiring real-time image stream data from a plurality of pre-arranged infrared optical observation devices, numbering and time-stamping the real-time image stream data, and forming an original image set; S12, based on the original image set, using a histogram equalization adaptive exposure algorithm to obtain a contrast-optimized image set; S13, based on the contrast-optimized image set, performing an inverse operation using an atmospheric turbulence correction algorithm based on Zernike polynomials to obtain a corrected image set; S14, based on the corrected image set, using the wavelet transform multi-resolution analysis method to perform multi-scale image fusion to obtain a fused high-quality image; S15. Based on the fused high-quality image and the pre-stored reference image, a feature matching algorithm based on a local feature descriptor is applied to identify corresponding point pairs; based on the identified corresponding point pairs, a RANSAC algorithm is used to estimate a transformation matrix; based on the transformation matrix, an affine transformation is performed on the fused high-quality image to obtain a pre-processed image set; Step S2 is further as follows: S21. Based on the preprocessed image set, an improved fast corner point detection algorithm is used to identify the significant corner points in the heat source matrix; non-maximum suppression and sub-pixel precision optimization are performed on the significant corner points, and finally a feature point set is obtained; S22. Based on the feature point set and the preprocessed image set, an enhanced binary descriptor generation algorithm is used to generate a descriptor for each feature point to obtain a descriptor set; based on the feature point set and descriptor set of the current frame, and the pre-stored feature point set and descriptor set of the previous frame, an improved optical flow tracking algorithm is applied to estimate the vertical displacement of the feature points between adjacent frames; based on the estimated vertical displacement, a time consistency constraint is introduced to filter out the vertical displacement that does not conform to the expected vertical trajectory; and finally a vertical trajectory set of the feature points is obtained; S23, based on the vertical trajectory set of feature points, calculating a multi-dimensional feature vector of each vertical trajectory, including vertical trajectory length, average velocity, and acceleration change; Based on the multidimensional feature vector, the Mahalanobis distance from each vertical trajectory to its k nearest neighbors is calculated in the feature space; based on the Mahalanobis distance, an adaptive threshold method is used to identify outliers, and the outliers are marked as abnormal vertical trajectories; the abnormal vertical trajectories are removed from the feature point vertical trajectory set, and finally a filtered vertical trajectory set is obtained; S24, based on the filtered vertical trajectory set, performing Delaunay triangulation on the plane, and dynamically adjusting the shape constraints of triangles in the Delaunay triangulation to obtain a triangulation result; converting the triangulation result into a graph structure to obtain a time-varying feature graph; Step S3 is further as follows: S31. Based on the time-varying feature graph, calculate the Laplace matrix at adjacent moments; perform eigenvalue decomposition on the Laplace matrix to obtain an eigenvalue sequence; based on the eigenvalue sequence, calculate a graph topology change index set; S32, based on the vertical trajectory set of feature points, dividing the time-varying feature map into overlapping local areas; based on each local area, constructing an energy function, including a data term and a regularization term; minimizing the energy function through an iterative optimization method to obtain an optimal deformation vector for each feature point; based on the optimal deformation vector of each feature point, forming a local vertical deformation field; S33, based on the local vertical deformation field and the graph topology change index set, construct the global settlement field equation, including the local vertical deformation constraint term and the global consistency constraint term; The global sedimentation field equation is discretized into a large-scale sparse linear system using the variational method, and solved by the conjugate gradient method to obtain the sedimentation field function. S34, discretizing based on the settlement field function and the historical settlement field function set pre-stored in the system database to obtain a high-dimensional data matrix; Based on the high-dimensional data matrix, a time-shift matrix pair is constructed; based on the time-shift matrix pair, the eigendecomposition of the dynamic matrix is ​​calculated through singular value decomposition and low-rank approximation, and finally the dynamic mode set and the corresponding growth rate are obtained; S35, projecting the settlement field function to the dynamic pattern space to obtain a pattern coefficient vector; predicting the pattern coefficient at a future time based on the pattern coefficient vector and the growth rate; linearly combining the pattern coefficient at a future time with the dynamic pattern set to obtain a predicted settlement field at a future time; Step S4 is further as follows: S41, based on the settlement field function and the predicted settlement field, a linear combination of orthogonal polynomial basis functions is expanded; based on the linear combination, the expansion coefficient is calculated by using the Galerkin projection method, and finally an uncertainty field is obtained; S42, performing wavelet packet decomposition based on the uncertainty field and the preprocessed image set to obtain a predetermined number of subbands; Calculate the energy and cross-correlation coefficient of each sub-band, and select the significant sub-band through thresholding method based on the energy and cross-correlation coefficient of each sub-band; Based on the significant subbands, error components of multiple scales are reconstructed, and finally a set of error scale coefficients is formed; S43, discretizing the settlement field function into a state vector, and constructing a state transfer equation and an observation equation based on the state vector; Based on the error scale coefficient set, the covariance matrix of the noise in the state transfer equation and the observation equation is dynamically adjusted to obtain the adjusted state transfer equation and the observation equation; Based on the adjusted state transfer equation and observation equation, Kalman filtering is performed to iteratively optimize the state estimation; The state estimation is converted into a function form to obtain an optimized settlement field function; S44, constructing a high-resolution subgrid based on the optimized settlement field function; fitting the settlement data of each subgrid area based on the high-resolution subgrid to obtain a local quadratic polynomial fitting function; estimating polynomial coefficients based on the local quadratic polynomial fitting function using the least squares method; Based on the polynomial coefficients, the settlement values ​​of the sub-grid points are calculated; based on the grayscale information of the preprocessed image set, the settlement values ​​are fine-tuned to obtain the fine-tuned settlement values; based on the fine-tuned settlement values, a high-precision settlement field function is generated.

2. The non-contact long-distance monitoring method for structural settlement according to claim 1, characterized in that: Step S12 is further as follows: S121, based on the original image set, calculating the grayscale histogram of each original image; performing cumulative summation on the grayscale histogram to obtain a cumulative distribution function; normalizing the cumulative distribution function to obtain a grayscale mapping function; S122, processing each pixel in the original image based on the grayscale mapping function to obtain a new pixel value; The new pixel values ​​are combined to form an equalized image; Based on the equalized image, calculate the average brightness value; Calculate a brightness adjustment factor based on the average brightness value and a preset target brightness value; S123, performing linear adjustment based on the brightness adjustment factor and the equalized image to obtain an adjusted pixel value; Combining the adjusted pixel values ​​to form a brightness-adjusted image; dividing the brightness-adjusted image into a predetermined number of sub-blocks, and calculating the local contrast of each sub-block; S124, based on the local contrast, the preset contrast threshold and the brightness adjusted image, perform adaptive gamma correction on each sub-block to obtain a gamma-corrected sub-block; based on the gamma-corrected sub-block, use a bilinear interpolation method to perform a smooth transition to obtain a smoothed sub-block; recombine the smoothed sub-blocks to form a final adaptive exposure adjusted image; combine all the adaptive exposure adjusted images to form a contrast optimized image set.

3. The non-contact long-distance monitoring method for structural settlement according to claim 1 is characterized in that: Step S13 is further as follows: S131, based on the contrast-optimized image set, performing a fast Fourier transform on each image to obtain a frequency domain image; calculating a power spectrum of the frequency domain image; performing a logarithmic transform on the power spectrum to obtain a logarithmic power spectrum; performing radial averaging on the logarithmic power spectrum to obtain a one-dimensional power spectrum density function; S132. Based on the one-dimensional power spectral density function, the Kolmogorov turbulence model is fitted using the least squares method, and the best fitting coefficient is calculated to obtain an estimated turbulence intensity parameter; based on the turbulence intensity parameter, the atmospheric coherence length is calculated; based on the atmospheric coherence length, the order of the Zernike polynomial is calculated, and the corresponding Zernike polynomial is generated; S133, based on the frequency domain image, calculating the coefficient of each Zernike polynomial; based on the coefficient of each Zernike polynomial, constructing the wavefront phase of atmospheric turbulence, and calculating the atmospheric transfer function; based on the atmospheric transfer function and the frequency domain image, sequentially performing a deconvolution operation and an inverse Fourier transform to obtain a corrected spatial domain image; S134, performing edge sharpening processing based on the corrected spatial domain image to obtain an edge-sharpened image; based on the edge-sharpened image, using the Sobel operator to calculate the image gradient to obtain a gradient amplitude; The gradient magnitude is added to the corrected spatial domain image to obtain a final corrected image; and all final corrected images are combined to form a corrected image set.

4. The non-contact long-distance monitoring method for structural settlement according to claim 1 is characterized in that: Step S14 is further as follows: S141, based on the corrected image set, using the wavelet basis function to perform a two-dimensional discrete wavelet transform on each image to obtain a first-layer decomposition result, including an approximate coefficient matrix and a detail coefficient matrix; Based on the first-layer decomposition result, two-dimensional discrete wavelet transform is performed again to obtain the second-layer decomposition result: repeat until the preset layer decomposition result is obtained; All decomposition results are combined to obtain multi-scale decomposition coefficients; S142, calculating energy distribution of detail coefficient matrix based on multi-scale decomposition coefficients; Based on the energy distribution, determine the set of significant scales; Based on the significant scale set and multi-scale decomposition coefficients, soft threshold processing is performed to obtain the processed coefficients; S143, based on the processed coefficients, using inverse wavelet transform to reconstruct the image to obtain an enhanced image; calculating a structural similarity index between the enhanced image and the images in the corrected image set; If the structural similarity index is lower than the preset threshold, the parameters of the soft threshold processing are adjusted, and steps S142 to S143 are repeated until the maximum number of iterations is reached; otherwise, no processing is performed; Output the final enhanced image; S144, performing weighted fusion on the final enhanced image and the images in the corrected image set to obtain a fused high-quality image.

5. The non-contact long-distance monitoring method for structural settlement according to claim 1 is characterized in that: The graph topology change index set calculated based on the eigenvalue sequence in step S31 is further: S311, based on the eigenvalue sequence, calculate the relative change rate of each pair of adjacent eigenvalues; store all relative change rates into the eigenvalue change rate array, calculate the mean and standard deviation of the eigenvalue change rate array; based on the mean and standard deviation, use the three times standard deviation principle to identify the eigenvalues ​​with significant changes, and generate a set of significant change eigenvalue indexes; S312, extracting corresponding eigenvectors based on the significant change eigenvalue index set; calculating the angle of the corresponding eigenvectors; calculating the mean and standard deviation of the angle; identifying the significantly rotated eigenvectors based on the mean and standard deviation of the angle, and generating a significant change set; S313. Calculate the comprehensive change index based on the significant change set, the eigenvalue sequence and the eigenvalue change rate array; normalize the comprehensive change index to obtain the final graph topology change index, and form a graph topology change index set.

6. Non-contact long-distance monitoring system for structural settlement, characterized in that: include: at least one processor; as well as, a memory communicatively connected to at least one of the processors; wherein, The memory stores instructions that can be executed by the processor, and the instructions are used to be executed by the processor to implement the non-contact long-distance monitoring method for structural settlement according to any one of claims 1 to 5.

Images