Mountain high-precision displacement monitoring method based on frequency domain digital image cross correlation

Abstract

The invention discloses a mountain high-precision displacement monitoring method based on frequency domain digital image cross correlation, and relates to the technical field of remote sensing measurement and digital image processing. Comprising the steps of obtaining multi-stage optical images; pre-processing the obtained multi-stage optical image; implementing orthotopography correction by using a digital elevation model (DEM); interference caused by vegetation and shadow is suppressed; constructing a multi-level resolution pyramid and carrying out frequency domain cross-correlation operation; performing sub-pixel peak value positioning to obtain a sub-pixel level displacement value; and generating a continuous displacement field, performing quality control and uncertainty expression at the same time, and finally outputting a single-source displacement field result. According to the method, through DEM orthographic projection, the two-stage projection difference caused by mountain fluctuation is reduced.

Description

Technical Field

[0001] This invention relates to the fields of remote sensing measurement and digital image processing technology, and more specifically to a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation. Background Technology

[0002] Existing digital image correlation (DIC) methods mostly rely on spatial domain grayscale matching. However, in large-format mountain remote sensing images, spatial domain methods are prone to problems such as low peak signal-to-noise ratio, matching drift, and insufficient subpixel accuracy due to the influence of terrain undulation, viewing angle differences, non-rigid movement of vegetation, and changes in shadow and illumination.

[0003] Frequency-domain DIC (FFT-DIC) performs cross-correlation operations in the frequency domain using the convolution theorem, determining pixel-level displacements using phase correlation peaks, and then performing sub-pixel interpolation, which can significantly improve matching robustness. However, without accompanying steps such as DEM orthorectification, vegetation / shading suppression, multi-level resolution pyramids, and quality control, it is still difficult to obtain continuous and reliable displacement fields in mountainous environments.

[0004] The following technical problems exist in the existing technology: 1) Terrain error: There is a projection difference between the two images on undulating terrain. If only image registration is used, the residual system offset is difficult to eliminate. 2) Cross-scale instability: Satellite and UAV imagery have significant differences in resolution and imaging geometry, making direct matching prone to failure or local mismatch; even within a single source, multi-scale texture differences can affect search efficiency and stability. 3) Vegetation / Shadow Interference: Seasonal changes in vegetation and wind-induced movement disrupt the assumption of local rigidity; changes in shadow and illumination bring low-frequency brightness differences, weakening the significance of correlation peaks; 4) Computational cost and search range: Direct full-domain search of high-resolution images involves a large amount of computation, and peak identification is easily affected by noise.

[0005] Therefore, proposing a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation to address the difficulties of existing technologies is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0006] In view of this, the present invention provides a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation, which is used to solve the technical problems existing in the prior art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation includes the following steps: S1. Acquire multiple optical images; S2. Preprocess the acquired multi-phase optical images; S3. Implement orthophoto terrain correction using a digital elevation model (DEM); S4. Suppress the disturbance caused by vegetation and shadows; S5. Construct a multi-level resolution pyramid and perform frequency domain cross-correlation operations; S6. Perform subpixel peak positioning to obtain subpixel-level displacement values; S7. Generate a continuous displacement field, perform quality control and uncertainty expression, and finally output the single-source displacement field result.

[0008] Optionally, the specific preprocessing steps for the acquired multi-stage optical images in S2 are as follows: Image acquisition: Select two or more periods of UAV imagery or satellite multispectral imagery; Unified processing: Unify the optical images from multiple periods to the WGS84 / UTM coordinate system; Radiometric / ensemble processing: Radiometric calibration and atmospheric correction for satellite multispectral imagery; Lens distortion correction and brightness / color equalization for UAV imagery; Coarse registration and resampling: using feature points or imaging geometry models to achieve preliminary alignment of influences in the same region.

[0009] Optionally, the specific details of orthophoto correction using the Digital Elevation Model (DEM) in S3 are as follows: Orthophoto topographic correction employs a standardized digital elevation model (DEM), a unified sensor model, back projection sampling within the same grid, and a quality check process, as follows: First, a digital elevation model (DEM) covering the study area is selected, the horizontal / vertical benchmarks are defined, the geodetic height to orthographic height conversion is performed, and gap filling, outlier smoothing, and micro-depression / peak correction are completed. The digital elevation model (DEM) was then resampled to the same grid resolution as the target orthophoto and aligned to a unified pixel origin, row and column start points, and range to ensure the pixel-level correspondence between the two orthophotos. For satellite multispectral imagery, the RPC imaging geometry provided by the manufacturer is read and RPC offset correction is performed based on stable ground control points or corresponding points between images to improve ground feature positioning accuracy. For UAV imagery, the camera interior and exterior orientation elements obtained from field work or SfM / aerial triangulation are used, and the distortion coefficients are fine-tuned and verified.

[0010] Optionally, the specific details of suppressing the disturbance caused by vegetation and shadows in S4 are as follows: Vegetation masking: Calculated based on near-infrared and red light, using the following formula: ; in, It is in the near-infrared band. It is in the red light band; Vegetation masks are generated by segmenting using adaptive or empirical thresholds; then, the segmented boundaries are expanded or contracted. Shadow masking: Identifies shadow areas using brightness thresholds and local contrast, and normalizes the shadow index; and NDVI Mask merging is used to mask unstable pixels.

[0011] Optionally, the specific details of constructing a multi-level resolution pyramid and performing frequency domain cross-correlation operations in S5 are as follows: Pyramid construction: For the two orthophotos and the mask / weight, first apply a light Gaussian low-pass filter, then downsample at a ratio of ×1 / 2 to form resolution levels of 1 / 8, 1 / 4, 1 / 2, and 1×, with pixels between layers kept aligned; Search shrinkage: The low-resolution layer performs a coarse match first, and the resulting initial displacement value is upsampled and used as the initial search value for the high-resolution layer to narrow the search range, reduce computational cost and the probability of mismatch. Window and step size: A sliding window is used to traverse the grid one by one. If the effective pixels in the window are less than 60-70%, the window is skipped or temporarily enlarged. Frequency domain phase correlation: The cross-power spectrum of the reference window f(x,y) and the target window g(x,y) is calculated and normalized. The phase correlation peak provides the initial value of the pixel-level displacement, as shown in the following formula:

[0012] in, This is the inverse Fourier transform. Fourier spectrum It is a complex conjugate.

[0013] Optionally, the specific details of sub-pixel peak localization to obtain sub-pixel level displacement values ​​in S6 are as follows: After obtaining the integer pixel peak value of the relevant surface C, the sub-pixel offset is estimated in the neighborhood. , The final displacement is obtained by adding the integer displacement to the integer displacement. Logarithmic parabolic interpolation: ; in, The correlation value is 1 pixel to the left of the peak point in the x-direction. The correlation value is 1 pixel to the right of the peak point in the x-direction. The correlation value is located 1 pixel below the peak in the y-direction. The correlation value is located 1 pixel above the peak point in the y-direction. These are the relevant values ​​at the peak points, used as the central values ​​in the x-axis formula. This is the relevant value at the peak point, used as the center value in the y-axis formula; This results in sub-pixel offset. ; Alternatively, a two-dimensional Gaussian surface can be used for fitting: ; in, For the local coordinates of the relevant surface, , The subpixel coordinates of the Gaussian peak. For background / bias term, peak That is, subpixel displacement.

[0014] Optionally, the specific details of generating a continuous displacement field in S7, along with quality control and uncertainty expression, are as follows: Interpolation and Rasterization: Interpolate discrete displacement vectors to a regular grid, outputting the east / north component, magnitude, and direction; Anomaly mitigation: Anomaly vectors are filtered out by combining peak height, PCE, and local RMSE thresholds; 3×3 or 5×5 spatial median / directional consistency filtering is used. Uncertainty raster: Outputs pixel-level quality factor and uncertainty, used for result classification and early warning threshold setting.

[0015] As can be seen from the above technical solution, compared with the prior art, the present invention discloses a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation, the beneficial effects of which are: 1) By orthophotoning the DEM, the difference in projection between the two periods caused by the undulation of the mountain can be reduced; 2) By using NDVI and shadow masks, the influence of non-rigid bodies and lighting interference on peak localization is suppressed; 3) By leveraging a multi-level resolution pyramid, the search range is narrowed, computational costs are reduced, and matching robustness is improved; 4) Subpixel interpolation is used to improve peak positioning accuracy and form a continuous displacement field. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0017] Figure 1 A flowchart of a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation provided by the present invention; Figure 2The images provided in this embodiment of the invention are UAV images of the glacier terminus from June and October 2018, in which the upper image is the main image and the lower image is the secondary image. Figure 3 The images provided in this embodiment of the invention are UAV images of the glacier terminus from October 2018 and October 2019, where the upper image is the main image and the lower image is the secondary image. Figure 4 The images provided in this embodiment of the invention are UAV images of the glacier terminus from October 2019 and September 2020, where the upper image is the main image and the lower image is the secondary image. Figure 5 The images provided in this embodiment of the invention are UAV images of the glacier terminus from September 2020 and November 2020, in which the upper image is the main image and the lower image is the secondary image. Figure 6 The left column provided for the embodiments of the present invention is a two-dimensional displacement amplitude map obtained from the analysis of each master-slave image pair; the right column is a two-dimensional displacement vector direction map obtained from the analysis of each master-slave image pair. Figure 7 The velocity field based on Landsat 8 imagery is provided for embodiments of the present invention. Detailed Implementation

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

[0019] See Figure 1 As shown, this invention discloses a high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation, comprising the following steps: S1. Acquire multiple optical images; S2. Preprocess the acquired multi-phase optical images; S3. Implement orthophoto terrain correction using a digital elevation model (DEM); S4. Suppress the disturbance caused by vegetation and shadows; S5. Construct a multi-level resolution pyramid and perform frequency domain cross-correlation operations; S6. Perform subpixel peak positioning to obtain subpixel-level displacement values; S7. Generate a continuous displacement field, perform quality control and uncertainty expression, and finally output the single-source displacement field result.

[0020] Furthermore, the specific preprocessing steps for the acquired multi-stage optical images in S2 are as follows: Image acquisition: Select two or more periods of UAV imagery or satellite multispectral imagery; Unified processing: Unify the optical images from multiple periods to the WGS84 / UTM coordinate system; Radiometric / ensemble processing: Radiometric calibration and atmospheric correction for satellite multispectral imagery; Lens distortion correction and brightness / color equalization for UAV imagery; Coarse registration and resampling: using feature points or imaging geometry models to achieve preliminary alignment of influences in the same region.

[0021] Specifically, to ensure the robustness and comparability of subsequent frequency domain cross-correlation, this invention first completes unified preprocessing of multiple image periods at the data level: for the study area and monitoring target, two (or more) optical images of the same season / similar solar altitude angle and azimuth angle are selected first; for UAV data, the same model / lens and consistent or correctable flight altitude and swivel are used to ensure that the ground resolution (GSD) difference is ≤10%, and the heading / lateral overlap is maintained as much as possible (e.g., ≥70% / 60%) to facilitate subsequent geometric constraints; for satellite data, products of the same sensor and the same processing level (e.g., L1TP / L2A) are selected first, the cloud cover is controlled within an acceptable threshold (e.g., ≤10%), and geometric metadata such as RPC / attitude and orbit are recorded. After importing the raw data, standardization at the radiometric and geometric levels is first performed: satellite imagery converts DN to surface reflectance using coefficients provided by the sensor and performs atmospheric correction (e.g., using existing physical or empirical models); UAV imagery undergoes radial / tangential distortion correction based on camera intrinsic / distortion parameters or calibration plate results, and exposure / white balance homogenization and brightness-color equalization (mild parameters such as histogram matching, quantile matching, or Wallis contrast normalization) are used to reduce low-frequency brightness differences across multiple periods while avoiding excessive smoothing of textures.

[0022] Coarse registration is then performed: satellite imagery is used as initial values ​​with RPC geometric models or orbital attitudes, and combined with SIFT / SURF / ORB feature point matching and RANSAC robust estimation to solve for full-frame initial alignment of translation / affine or thin-plate splines; if UAV imagery is formed by stitching multiple images into an orthorectified mosaic, an orthorectified base map is first generated after SfM / aerial triangulation and bundle adjustment, and then full-frame preliminary alignment is performed for multiple periods using the same feature matching + RANSAC. To avoid scale effects, a unified pixel grid is introduced: the two periods of imagery are resampled to a consistent GSD (e.g., 10m / px for satellites, 5cm / px for UAVs) and the same projection coordinate system (WGS84 / UTM, area code determined according to the study area), and the same raster origin and pixel row and column starting points are fixed; resampling preferentially uses bilinear or cubic convolution (to avoid mosaic caused by nearest neighbors, but also to avoid introducing excessive smoothing), and a NoData mask is explicitly defined.

[0023] Finally, the image range is cropped using AOI to unify boundaries and extend buffers; collinearity residuals and mismatch rates are checked (e.g., average feature residual ≤ 0.5 pixels, coarse registration mismatch rate ≤ 5%), and the images are saved as GeoTIFF files with complete metadata as DIC input. At this point, multiple image sets are unified in resolution, radiometric scale, geometric framework, and grid reference, providing consistent and traceable input conditions for subsequent DEM orthorectification and frequency domain cross-correlation matching.

[0024] Furthermore, the specific details of orthophoto correction using the Digital Elevation Model (DEM) in S3 are as follows: Orthophoto topographic correction employs a standardized digital elevation model (DEM), a unified sensor model, back projection sampling within the same grid, and a quality check process, as follows: First, select a digital elevation model (DEM) covering the study area (preferably with a resolution ≤1m–5m), define the horizontal / vertical datum (projection such as WGS84 / UTM, vertical datum such as orthographic height / geodetic height), perform geodetic height to orthographic height conversion, and complete gap filling, outlier smoothing, and micro-depression / peak correction (only perform mild filtering on steep slopes to preserve the true terrain undulations). The digital elevation model (DEM) was then resampled to the same grid resolution as the target orthophoto and aligned to a unified pixel origin, row and column start points, and range to ensure the pixel-level correspondence between the two orthophotos. For satellite multispectral imagery, the RPC imaging geometry provided by the manufacturer is read and RPC offset correction is performed based on stable ground control points or corresponding points between images to improve ground feature positioning accuracy. For UAV imagery, the camera interior and exterior orientation elements obtained from field work or SfM / aerial triangulation are used, and the distortion coefficients are fine-tuned and verified.

[0025] Specifically, orthophoto projection uses a map-to-image method: for each ground pixel (X,Y) in the output orthophoto grid, the elevation Z=DEM(X,Y) is obtained from the Digital Elevation Model (DEM). (X,Y,Z) is substituted into the sensor model (RPC for satellites, collinearity equations for UAVs) to obtain the original image coordinates (r,c). Then, bilinear or cubic convolution interpolation is used to sample grayscale / reflectivity and write it into the orthophoto pixel (avoiding jagged edges / mosaic caused by nearest-neighbor interpolation). For pixels with terrain occlusion, visibility is determined based on the line-of-sight direction and the gradient of the neighboring DEM, and an invalid / occluded mask is generated to prevent subsequent DIC from generating false matches in the occluded area. Considering the impact of DEM errors on planar positioning, an error propagation estimate is given: under imaging geometry with an off-axis angle of α, the planar displacement error caused by the elevation error ΔZ is approximately Δ... R≈ΔZ·tanα, corresponding to Δp≈ΔR / GSD in pixels (therefore, a high-precision DEM should be prioritized, or approximately vertical imaging should be used as much as possible); after completing the orthorectification of the two phases of images, sub-pixel correlation verification and offset evaluation are performed on rigid control areas such as stable bare rock / building roofs (statistical mean offset, 95th percentile error and local RMSE). When systematic translation / rotation residuals are found, small-scale rigid / affine fine-tuning is performed (without changing the DEM reference, only correcting the image); finally, the two phases of orthorectified images and the corresponding NODATA / occlusion mask are output. The two strictly share the same raster reference (projection, resolution, origin, and range are consistent), and the source and accuracy of DEM, sensor model / offset correction parameters, interpolation kernel type and quality indicators are recorded in the metadata, providing geometrically consistent and traceable input for subsequent frequency domain cross-correlation.

[0026] Furthermore, the specific details of S4 regarding the suppression of disturbances caused by vegetation and shadows are as follows: Vegetation masking: Calculated based on near-infrared and red light, using the following formula: ; in, It is in the near-infrared band. It is in the red light band; Vegetation masks are generated by segmenting using adaptive or empirical thresholds; then, the segmented boundaries are expanded or contracted. Shadow masking: Identifies shadow areas using brightness thresholds and local contrast, and normalizes the shadow index; and NDVI Mask merging is used to mask unstable pixels.

[0027] Specifically, the mask area is shielded in the cross-correlation calculation to reduce the interference of non-rigid bodies and low-frequency brightness differences on peak positioning.

[0028] To reduce the interference of non-rigid vegetation and illumination / shadow on phase correlation peak localization, this invention constructs a vegetation and shadow dual mask based on orthophotos, and synchronously downsamples and buffers the boundaries according to the pyramid hierarchy to ensure effective texture within the matching window of each layer: First, calculation is performed based on multispectral data. NDVI Sentinel-2 can be set to B8 / B4, and Landsat8 / 9 to B5 / B4. The vegetation binary map is obtained by segmentation using scene-adaptive thresholds or empirical thresholds [0.30, 0.45]. For drone RGB cameras without NIR, VARI=(GR) / (G+RB) or ExG(2G-RB) is used as an approximation, and thresholding is performed in the same way. Subsequently, for shadows, low-light areas are identified by combining global brightness thresholds with local contrast / entropy in the reflectivity domain. Cross-validation is then performed using hillshade thresholds (such as shadow index or shadow probability maps) calculated from the DEM and solar geometry (azimuth angle, solar altitude angle) to suppress pseudo-dark areas caused solely by material differences. The two types of masks are calculated separately for the main image and the secondary image, and the union is taken as the final masking area to cover unstable pixels in either period. The mask boundary is morphologically dilated by 1-2 pixels according to the matching window size or 5% of the window width to avoid aliasing caused by the window crossing the stable / unstable boundary. Small isolated patches (area less than 10% of the typical window area) are subjected to opening operation / connected component denoising to eliminate burrs and holes. Considering the influence of water and snow on the correlation, NDWI / NDSI masks can be optionally superimposed to exclude them. At the same time, the shadow mask is further updated using cloud / cloud shadow quality bands (if available).

[0029] To reduce the impact of hard threshold boundaries on abrupt changes in correlation peaks, soft weighted bands (such as 1-2 pixel Gaussian attenuation) are generated at the final mask boundaries. In the correlation calculations, the masked pixels are assigned zero weights, while the buffer bands participate in the attenuation according to their weights, thus achieving "soft masking." When constructing a multi-level pyramid, the mask and weight map are downsampled at the same magnification and aligned with the pixels of each layer of the image to ensure consistency between layers.

[0030] In terms of quality control, the mask coverage is statistically analyzed (recommended range: 0.1-0.7), and the improvement in peak quality (e.g., peak height) and local RMSE before and after masking is compared. If over-masking leads to insufficient effective texture, the threshold is automatically lowered or the expansion radius is reduced. The final output is a binary mask + soft weight map (two-stage union) of the same size as the image, which is uniformly used in subsequent frequency domain cross-correlation and subpixel fitting. This ensures stable suppression of low-frequency brightness differences and non-rigid body disturbances caused by vegetation and shadows in both UAV and satellite single-source pipelines.

[0031] Furthermore, the specific details of constructing a multi-level resolution pyramid and performing frequency domain cross-correlation operations in S5 are as follows: Pyramid construction: For the two orthophotos and the mask / weight, first apply a light Gaussian low-pass filter, then downsample at a ratio of ×1 / 2 to form resolution levels of 1 / 8, 1 / 4, 1 / 2, and 1×, with pixels between layers kept aligned; Specifically, to reduce aliasing and edge effects, the window is gently windowed (e.g., Hann / Hamming) and normalized to zero mean-unit variance (ZNCC preprocessing) before correlation. In practice: we first perform a light blurring (Gaussian low-pass) on the two images and their corresponding masks to remove overly fine textures. Then, the images are successively reduced by half, generating multiple resolution levels (e.g., original image, 1 / 2, 1 / 4, 1 / 8 size), while maintaining consistent pixel positions. This facilitates subsequent layer-by-layer matching. Before performing correlation calculations, to avoid jagged edges or edge errors, we apply a "smoothing window" (e.g., Hann or Hamming) to each window and unify the pixel values ​​to a standard range of "mean 0, variance 1". This reduces the interference of illumination differences and boundary effects on matching.

[0032] Search shrinkage: The low-resolution layer performs a coarse match first, and the resulting initial displacement value is upsampled and used as the initial search value for the high-resolution layer to narrow the search range, reduce computational cost and the probability of mismatch. Specifically, a coarse matching is first performed on a low-resolution image layer to obtain a rough displacement estimate. This estimate is then upsampled to a higher-resolution layer, and within that layer, a precise match is searched within a small area centered on this estimate. This reduces computational cost and the possibility of matching errors.

[0033] Window and step size: Use a sliding window to traverse the grid (satellite is recommended to be 128×128 pixels, step size 64; drone is recommended to be 64×64 pixels, step size 32; for weak textures, it can be enlarged to 160×160 / 80 (satellite) or 96×96 / 48 (drone)). If the effective pixels in the window are less than 60-70%, skip or temporarily enlarge the window. Specifically, for satellite multispectral imagery, the recommended window size is 128×128 pixels, with each movement being 64 pixels. For drone imagery, the recommended window size is 64×64 pixels, moving 32 pixels at a time. If the texture of the image area is weak (e.g., the texture is very simple), the window can be enlarged appropriately. For example, use 160×160 (step size 80) for satellite imagery and 96×96 (step size 48) for UAV imagery.

[0034] Additionally, if the number of usable pixels (excluding the mask) in a window is less than 60%–70% of the total, it means that the window lacks effective information, so it should be skipped directly, or the window should be temporarily enlarged before calculation.

[0035] Frequency domain phase correlation: The cross-power spectrum of the reference window f(x,y) and the target window g(x,y) is calculated and normalized. The phase correlation peak provides the initial value of the pixel-level displacement, as shown in the following formula:

[0036] in, This is the inverse Fourier transform. Fourier spectrum It is a complex conjugate.

[0037] Specifically, in the frequency domain, a Fourier transform is performed on the reference window and the target window, and their cross-power spectrum is calculated (equivalent to comparing the similarity of the two images in the frequency domain). The result is then normalized. This result yields a clear peak position, which corresponds to the initial pixel-level displacement of the two images.

[0038] For areas with significant lighting differences, Wallis local contrast normalization or bandpass weighting can be enabled within a window to weaken extremely low / high frequency components, thereby improving the signal-to-noise ratio and stability of the correlation peak. In practice: In areas with significant lighting differences, a "local contrast adjustment" (e.g., using the Wallis algorithm) can be performed on the image within each window, or weights can be added to components of different frequencies. The purpose of this is to reduce the impact of excessively dark or bright areas on the matching results, thus making the correlation peak clearer and more stable.

[0039] Furthermore, the specific details of sub-pixel peak localization to obtain sub-pixel level displacement values ​​in S6 are as follows: In digital image correlation (DIC) or phase correlation, the degree of matching between two images is calculated. This matching result is usually not a single number, but a two-dimensional matrix, where each value represents the similarity between the two images at a certain translation amount. This two-dimensional matrix is ​​called the correlation surface, usually denoted as C.

[0040] After obtaining the integer pixel peak value of the relevant surface C, the sub-pixel offset is estimated in the neighborhood. , The final displacement is obtained by adding the integer displacement to the integer displacement. Logarithmic parabolic interpolation: ; in, The correlation value is 1 pixel to the left of the peak point in the x-direction. The correlation value is 1 pixel to the right of the peak point in the x-direction. The correlation value is located 1 pixel below the peak in the y-direction. The correlation value is located 1 pixel above the peak point in the y-direction. These are the relevant values ​​at the peak points, used as the central values ​​in the x-axis formula. This is the relevant value at the peak point, used as the center value in the y-axis formula; This results in sub-pixel offset. ; Alternatively, a two-dimensional Gaussian surface can be used for fitting: ; in, The local coordinates of the relevant surface (in rows and columns / pixels of the grid in the neighborhood of the peak point). , These are the sub-pixel coordinates of the Gaussian peak (continuous positions relative to the neighborhood coordinate system). For background / bias term, peak That is, subpixel displacement.

[0041] Specifically, after finding the integer pixel peak on the relevant surface C, further calculations can be performed in a small area around it (usually a 3×3 or 5×5 neighborhood) to obtain a more accurate displacement value (sub-pixel offset) than a single pixel. The final result is an integer displacement plus sub-pixel correction.

[0042] Two common methods are logarithmic parabolic interpolation (simple and stable) and two-dimensional Gaussian surface fitting (more suitable for high-precision scenarios). In practical applications, both methods can improve accuracy to approximately 0.05 pixels. You can choose one method based on your needs, or perform cross-validation to ensure both robustness and accuracy.

[0043] Furthermore, the specific details of generating a continuous displacement field in S7, along with quality control and uncertainty expression, are as follows: Interpolation and Rasterization: Interpolate discrete displacement vectors to a regular grid, outputting the east / north component, magnitude, and direction; Specifically, these scattered displacement vectors are interpolated onto a regular grid to obtain a complete displacement field. The output typically includes east-west displacement, north-south displacement, displacement magnitude (amplitude), and direction.

[0044] Anomaly mitigation: Anomaly vectors are filtered out by combining peak height, PCE, and local RMSE thresholds; 3×3 or 5×5 spatial median / directional consistency filtering is used. Specifically, the inspection results may contain error vectors. Unreliable data is filtered out using metrics such as peak height, PCE (peak energy ratio), and local RMSE (residual energy error). Then, 3×3 or 5×5 median filtering or directional consistency filtering is used to further remove noise points.

[0045] Uncertainty raster: Outputs pixel-level quality factor and uncertainty, used for result classification and early warning threshold setting.

[0046] Specifically, assigning a quality factor and uncertainty index to each pixel allows us to clearly identify which areas yield reliable results and which require caution. This information can also be used for outcome grading and setting early warning thresholds.

[0047] In a specific embodiment: Case study of surface displacement at the tip of the Hailuogou Glacier tongue based on UAV imagery (June 30, 2018 – November 5, 2020, glacier surface displacement estimation based on frequency domain cross-correlation of multiple UAV images). Figures 2-5 In the middle, the upper part is the main image; the lower part is the secondary image. Figure 6 The left column shows the two-dimensional displacement amplitude, and the right column shows the displacement vector distribution. Figures 2 to 5 The yellow dashed polygon in the image represents the area near the leading edge of a glacier collapse. Figure 6 The black dashed polygon is... Figures 2 to 5 The corresponding areas marked by the yellow dashed lines are used to compare the consistency between displacement amplitude and vector direction.

[0048] Example: Unmanned aerial vehicle (UAV) imagery scene (Hailuogou Glacier in southeastern Qinghai-Tibet Plateau)

[0049] Two drone images from the study area, dated June 30, 2018 and November 5, 2020, were selected and processed according to the procedure of this invention: distortion correction, coordinate unification, and DEM orthorectification were completed; vegetation / shadows were masked according to the sensor band conditions (the drone model was DJI MAVIC series, equipped with an RGB camera, so VARI can be used instead of NDVI here); a multi-level resolution pyramid was constructed and frequency domain phase correlation was performed, followed by logarithmic parabolic (or two-dimensional Gaussian) subpixel fitting in the neighborhood of the correlation peak.

[0050] The results are as follows Figure 2 Wherever Figure 5 As shown: Differential calculation is performed on the master and slave images to obtain a two-dimensional displacement amplitude map (e.g.) Figure 6 (as shown in the left column) and two-dimensional displacement vector direction pattern ( Figure 6 As shown in the right column, this reflects the displacement gradient between the glacier body and its leading edge. The displacement vectors show a consistent direction in the collapse leading edge region, and anomalous vectors are effectively removed after screening. This case demonstrates that this method can stably obtain a continuous displacement field even under single-source UAV imagery conditions.

[0051] Example (Single Satellite Source, Landsat 8): Multiple Landsat 8 OLI images of the study area were selected as input data and processed according to the procedure of this invention: atmospheric and radiative corrections, unified projection coordinates, and DEM orthorectification were performed; after generating vegetation and shadow masks through NDVI thresholding and brightness / contrast identification, frequency domain phase correlation matching was performed on a multi-level resolution pyramid, and sub-pixel fitting was performed using a logarithmic parabola (or two-dimensional Gaussian). The resulting displacement / velocity field is as follows: Figure 7 As shown, the east-west and north-south components are illustrated. The results were used to identify displacement gradients and hotspot regions in the glacier body and its leading edge; anomaly vectors were filtered out using peak quality, local RMSE, and spatial consistency criteria.

[0052] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0053] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A high-precision displacement monitoring method for mountainous areas based on frequency domain digital image cross-correlation, characterized in that, Includes the following steps: S1. Acquire multiple optical images; S2. Preprocess the acquired multi-phase optical images; S3. Implement orthophoto terrain correction using a digital elevation model (DEM); S4. Suppress the disturbance caused by vegetation and shadows; S5. Construct a multi-level resolution pyramid and perform frequency domain cross-correlation operations; S6. Perform subpixel peak positioning to obtain subpixel-level displacement values; S7. Generate a continuous displacement field, perform quality control and uncertainty expression, and finally output the single-source displacement field result.

2. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 1, characterized in that, The specific preprocessing steps for the acquired multi-stage optical images in S2 are as follows: Image acquisition: Select two or more periods of UAV imagery or satellite multispectral imagery; Unified processing: Unify the optical images from multiple periods to the WGS84 / UTM coordinate system; Radiometric / ensemble processing: Radiometric calibration and atmospheric correction for satellite multispectral imagery; Lens distortion correction and brightness / color equalization for UAV imagery; Coarse registration and resampling: using feature points or imaging geometry models to achieve preliminary alignment of influences in the same region.

3. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 2, characterized in that, The specific details of orthophoto correction using the Digital Elevation Model (DEM) in S3 are as follows: Orthophoto topographic correction employs a standardized digital elevation model (DEM), a unified sensor model, back projection sampling within the same grid, and a quality check process, as follows: First, a digital elevation model (DEM) covering the study area is selected, the horizontal / vertical benchmarks are defined, the geodetic height to orthographic height conversion is performed, and gap filling, outlier smoothing, and micro-depression / peak correction are completed. The digital elevation model (DEM) was then resampled to the same grid resolution as the target orthophoto and aligned to a unified pixel origin, row and column start points, and range to ensure the pixel-level correspondence between the two orthophotos. For satellite multispectral imagery, the RPC imaging geometry provided by the manufacturer is read and RPC offset correction is performed based on stable ground control points or corresponding points between images to improve ground feature positioning accuracy. For UAV imagery, the camera interior and exterior orientation elements obtained from field work or SfM / aerial triangulation are used, and the distortion coefficients are fine-tuned and verified.

4. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 1, characterized in that, The specific measures taken in S4 to suppress interference from vegetation and shadows are as follows: Vegetation masking: Calculated based on near-infrared and red light, using the following formula: ； in, It is in the near-infrared band. It is in the red light band; Vegetation masks are generated by segmenting using adaptive or empirical thresholds; then, the segmented boundaries are expanded or contracted. Shadow masking: Identifies shadow areas using brightness thresholds and local contrast, and normalizes the shadow index; and NDVI Mask merging is used to mask unstable pixels.

5. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 1, characterized in that, The specific steps for constructing a multi-level resolution pyramid and performing frequency domain cross-correlation operations in S5 are as follows: Pyramid construction: For the two orthophotos and the mask / weight, first apply a light Gaussian low-pass filter, then downsample at a ratio of ×1 / 2 to form resolution levels of 1 / 8, 1 / 4, 1 / 2, and 1×, with pixels between layers kept aligned; Search shrinkage: The low-resolution layer performs a coarse match first, and the resulting initial displacement value is upsampled and used as the initial search value for the high-resolution layer to narrow the search range, reduce computational cost and the probability of mismatch. Window and step size: A sliding window is used to traverse the grid one by one. If the effective pixels in the window are less than 60-70%, the window is skipped or temporarily enlarged. Frequency domain phase correlation: The cross-power spectrum of the reference window f(x,y) and the target window g(x,y) is calculated and normalized. The phase correlation peak provides the initial value of the pixel-level displacement, as shown in the following formula: in, For inverse Fourier transform, Fourier spectrum It is a complex conjugate.

6. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 1, characterized in that, The specific steps for obtaining sub-pixel level displacement values ​​through sub-pixel peak localization in S6 are as follows: After obtaining the integer pixel peak value of the relevant surface C, the sub-pixel offset is estimated in the neighborhood. , The final displacement is obtained by adding the integer displacement to the integer displacement. Logarithmic parabolic interpolation: ； in, The correlation value is 1 pixel to the left of the peak point in the x-direction. The correlation value is 1 pixel to the right of the peak point in the x-direction. The correlation value is located 1 pixel below the peak in the y-direction. The correlation value is located 1 pixel above the peak point in the y-direction. These are the relevant values ​​at the peak points, used as the central values ​​in the x-axis formula. This is the relevant value at the peak point, used as the center value in the y-axis formula; This results in sub-pixel offset. ; Alternatively, a two-dimensional Gaussian surface can be used for fitting: ； in, For the local coordinates of the relevant surface, , The subpixel coordinates of the Gaussian peak. For background / bias term, peak That is, subpixel displacement.

7. The method for high-precision displacement monitoring in mountainous areas based on frequency domain digital image cross-correlation according to claim 1, characterized in that, The specific details of generating a continuous displacement field in S7, along with quality control and uncertainty expression, are as follows: Interpolation and Rasterization: Interpolate discrete displacement vectors to a regular grid, outputting the east / north component, magnitude, and direction; Anomaly mitigation: Anomaly vectors are filtered out by combining peak height, PCE, and local RMSE thresholds; 3×3 or 5×5 spatial median / directional consistency filtering is used. Uncertainty raster: Outputs pixel-level quality factor and uncertainty, used for result classification and early warning threshold setting.

Images