A method for calculating three-dimensional height of sand dunes by coupling multi-source remote sensing and physical model

Abstract

This invention relates to a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model. The method includes: acquiring several initial Sentinel-2 optical images of the same target area at different time phases; sequentially performing orthorectification and radiometric normalization on the initial Sentinel-2 optical images to obtain preprocessed Sentinel-2 optical images; obtaining the migration rate based on the relative displacement between two preprocessed Sentinel-2 optical images; calculating the saturated sand flux in the same time period as the migration rate based on the Martin & Kok sand flux model; and obtaining the dune height through inversion based on the Exner equation in geomorphology, according to the migration rate and the saturated sand flux. This invention enables low-cost and high-efficiency inversion of the dynamic elevation distribution of large-scale desert areas without the need for ground contact measurements, providing accurate scientific basis for wind and sand disaster prevention and environmental evolution research.

Description

Technical Field

[0001] This invention relates to the field of sand dune measurement technology, specifically to a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model. Background Technology

[0002] The deformation and migration of dune landforms not only directly reflect the evolution of regional aeolian environments but also serve as a crucial basis for assessing desertification processes and the effectiveness of sand control projects. Dune migration is essentially a spatial redistribution of sediment mass driven by wind. From the perspective of aeolian physics, a strict mass conservation and energy coupling mechanism exists between the vertical height of dunes, their migration rate, and the local wind transport potential. Therefore, accurately obtaining dune height parameters is essential for revealing the laws governing aeolian movement.

[0003] However, in the practice of long-term remote sensing monitoring of large areas in the heart of deserts, existing mainstream technologies all face insurmountable bottlenecks. Direct measurement technologies, such as airborne or ground-based lidar (LiDAR), while capable of acquiring sub-meter precision three-dimensional point clouds of the Earth's surface by emitting high-frequency laser pulses, are high-cost, discrete observation methods. Their data acquisition heavily relies on flight platforms, facing not only expensive flight costs and stringent low-altitude airspace control, but also limited single-scan swath width and exponentially increasing data processing volume. This makes it difficult to meet the demands of high-frequency, full-coverage dynamic monitoring of tens of thousands of square kilometers of uninhabited areas, resulting in extremely challenging long-term data acquisition.

[0004] Furthermore, while the most widely used optical image cross-correlation technique currently achieves sub-pixel-level horizontal migration rate extraction by tracking the texture features of preceding and following images, it is essentially a two-dimensional kinematic observation based on geometric matching. This method has two significant drawbacks: firstly, it can only sense displacement vectors on a plane, completely missing the vertical height information crucial for landform evolution; secondly, the algorithm is highly dependent on the texture contrast of the images, and in the spectrally homogeneous and texture-poor desert interior or in shadow areas with varying lighting conditions, it is prone to matching noise or spurious displacements due to sand grain aperture issues, leading to a decrease in the physical interpretability of the monitoring results.

[0005] Therefore, establishing a comprehensive monitoring and inversion framework that integrates multi-source remote sensing data, computer vision algorithms, and aeolian physical mechanisms is key to achieving accurate dune height inversion. Summary of the Invention

[0006] To address the aforementioned problems in existing technologies, this invention provides a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model. The technical problem to be solved by this invention is achieved through the following technical solution: This invention provides a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model, comprising: Acquire initial Sentinel-2 optical images of the same target area at several different time phases; The initial Sentinel-2 optical image is subjected to orthorectification and radiometric normalization in sequence to obtain the preprocessed Sentinel-2 optical image; The migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images. Based on the Martin & Kok sand flux model, the saturated sand flux in the same time period as the migration rate is calculated. Based on the Exner equation in geomorphology, the dune height is obtained by inversion according to the migration rate and the saturated sand flux.

[0007] In one embodiment of the invention, all the initial Sentinel-2 optical images are in the near-infrared band.

[0008] In one embodiment of the present invention, the initial Sentinel-2 optical image is sequentially subjected to orthorectification and radiometric normalization to obtain a preprocessed Sentinel-2 optical image, including: Based on the DEM data of the target area, the initial Sentinel-2 optical image is orthorectified to obtain the corrected Sentinel-2 optical image; The corrected Sentinel-2 optical image is subjected to radiometric normalization to obtain the preprocessed Sentinel-2 optical image. In one embodiment of the present invention, the corrected Sentinel-2 optical image is subjected to radiometric normalization to obtain the preprocessed Sentinel-2 optical image, including: Using the grayscale histogram of any one of the corrected Sentinel-2 optical images as a radiometric reference, the histogram matching method is used to perform radiometric normalization on the grayscale histograms of the remaining corrected Sentinel-2 optical images to obtain the preprocessed Sentinel-2 optical images. In one embodiment of the present invention, the migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images, including: The relative displacement is obtained by processing any two preprocessed Sentinel-2 optical images using COSI-Corr. The migration rate is obtained based on the relative displacement and time interval.

[0009] In one embodiment of the present invention, the migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images, including: The relative displacement is obtained by processing any two preprocessed Sentinel-2 optical images using the dense optical flow method. The migration rate is obtained based on the relative displacement and time interval. In one embodiment of the present invention, the saturated sand flux for the same time period as the migration rate is calculated based on the Martin & Kok sand flux model, including: Obtain the 10-meter height wind speed and surface sediment transport dynamic roughness from the ERA5-Land meteorological reanalysis data of the target area; Based on the wind speed at a height of 10 meters and the surface sand transport dynamic roughness, the surface friction wind speed is calculated by inversely using a preset formula corresponding to the logarithmic wind speed profile law. The preset formula is expressed as follows:

[0010] in, Indicates the height above the ground as Wind speed at that time Denotes the Kármán constant. Represents surface sediment transport dynamic roughness. Indicates the wind speed due to surface friction; The critical friction wind speed is obtained based on the critical friction wind speed calculation formula, which is expressed as follows:

[0011] in, This represents the critical frictional wind speed. Represents gravitational acceleration. Indicates the particle size of sand. Indicates the density of sand grains. Indicates air density; Based on the Martin & Kok sand flux model, the saturated sand flux in the same time period as the migration rate is obtained according to the surface friction wind speed and the critical friction wind speed. In one embodiment of the present invention, the saturated sand flux is expressed as:

[0012] in, Indicates saturated sand flux. This represents a dimensionless empirical coefficient. In one embodiment of the present invention, the dune height is obtained by inversion based on the Exner equation in geomorphology, according to the migration rate and the saturated sand flux, including: Based on the Exner equation, a first relationship is constructed between the spatial divergence of the rate of change of surface elevation and sediment flux, which is expressed as:

[0013] in, Indicates the surface elevation. Indicates time, Indicates the porosity of sand dunes. Represents sediment flux Spatial divergence; The second relationship is derived based on the dune migration rate, and is expressed as follows:

[0014] in, Indicates the migration rate of sand dunes. A spatial gradient representing the height of sand dunes; Based on the first and second relationships, by integrating a one-dimensional profile along the wind direction, the inverse coupling relationship between dune height, dune migration rate, and saturated sand flux at the dune crest is derived, resulting in a theoretical formula for calculating dune height. This theoretical formula is expressed as follows:

[0015] in, Indicates the height of the sand dunes. Indicates saturated sand flux. This indicates the overall density of sand dunes in the target area. Represents a dimensionless constant; A shape correction factor is introduced to modify the theoretical dune height calculation formula into an engineering-usable dune height calculation formula. The relative height of the dune is then calculated using this formula, which is expressed as follows:

[0016] in, Represents dimensionless empirical coefficients. This represents the critical frictional wind speed. Indicates wind speed due to surface friction. Represents gravitational acceleration. Indicates air density, This is the shape correction factor. In one embodiment of the present invention, after obtaining the dune height by inversion based on the migration rate and the saturated sand flux, the method further includes: By using flat Gobi desert or sparse ICESat-2 laser altimetry points at the edge of the target area as absolute elevation control points, the relative height field of the dunes is corrected to a reference surface, generating an absolute digital elevation model with geographic coordinates to reconstruct the three-dimensional dune topography. Compared with the prior art, the beneficial effects of the present invention are as follows: This invention proposes a method for calculating the three-dimensional elevation of sand dunes, which incorporates a sediment transport dynamics model to achieve a physical inversion from a two-dimensional motion field to the vertical height of the dunes. For optical data, it utilizes the rich texture features of the near-infrared band to achieve multi-scale and multi-modal data complementarity. Secondly, it introduces a large-displacement dense optical flow algorithm, unlike sparse feature point tracking. This algorithm achieves a "coarse-to-fine" stepwise matching by constructing an image pyramid, solving the problem of large-displacement sand dune tracking. It also introduces smoothing and data terms into the energy function, using a variational framework to minimize the global energy, thus enabling displacement information to be transmitted through neighborhood constraints even in weakly textured regions, generating a physically continuous and dense pixel-level displacement field. Finally, it constructs a physical inversion model for dune height based on the conservation of sediment transport flux. This invention overcomes the limitation of only monitoring planar displacement, establishing a bridge between two-dimensional migration rate, sediment transport flux, and dune height by utilizing the law of conservation of sediment mass. Specifically, by introducing the Martin and Kok model and using meteorological reanalysis wind speed data, the regional sediment transport flux considering the sand-lifting threshold is calculated. Based on the physical mechanism that the dune migration rate is inversely proportional to its height, the vertical height of the dunes is quantitatively inferred using high-precision migration rates retrieved from remote sensing. This step enables the calculation of dune height based solely on planar monitoring data, without the support of DEM data.

[0017] This invention achieves accurate extraction of pixel-level deformation fields. Compared with the fixed window and sparse grid of the traditional POT method, the dense optical flow method used in this invention can calculate the motion vector for each pixel in the image. Combined with variational regularization constraints, it can more accurately capture the differentiated motion features of details such as dune ridges and wings, avoiding the errors caused by window averaging.

[0018] This invention realizes the construction of a physical mechanism model for inverting dune height from two-dimensional dune motion. It goes beyond surface geometric observations, deeply integrating aeolian physics mechanisms and utilizing the principle of mass conservation to couple kinematic parameters from remote sensing observations with sand transport dynamics parameters, successfully inverting dune height. This method provides a novel, low-cost, and large-scale approach to obtaining dune vertical height parameters in areas lacking expensive LiDAR or high-precision DEM data.

[0019] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model, provided by the present invention. Figure 2 This is a schematic diagram of the annual average migration rate calculated from Sentinel-2 optical images using the dense optical flow method provided by the present invention. Figure 3 This is a schematic diagram showing the results of dune height inversion at several research points in the target area, provided by the present invention. Detailed Implementation

[0021] The present invention will be further described in detail below with reference to specific embodiments, but the implementation of the present invention is not limited thereto.

[0022] Example 1 Currently, in the field of dynamic monitoring and height inversion of sand dune landforms, there are two main implementation schemes: simple estimation models based on empirical rules and traditional optical image cross-correlation techniques.

[0023] For estimation models based on empirical rules, their core logic largely follows Bagnold's early classical observations, namely, that under ideal conditions, the migration rate of crescent-shaped dunes is strictly inversely proportional to their height. In practice, these schemes typically measure dune movement distance manually or semi-automatically, and directly calculate relative height using the inverse relationship, assuming the regional sediment transport flux is a spatiotemporal constant. However, this linearization reveals significant physical flaws when facing complex natural environments. It not only ignores the acceleration or shading effects of undulating terrain on local wind fields, resulting in a lack of spatial heterogeneity in flux parameter values, but also frequently employs the traditional Bagnold sediment transport formula. Modern aeolian physics evidence suggests that this formula tends to significantly overestimate actual sediment transport in high-wind-speed ranges, thus introducing systematic inversion errors. More importantly, in engineering practice, this method is often limited to the morphological analysis of individual standard dunes, making it difficult to support the automated, pixel-level continuous inversion requirements for large-scale, complex, parallel dune chains.

[0024] Another mainstream technology widely used in remote sensing monitoring is the conventional sliding window-based image cross-correlation technique. Mechanistically, this technique is a typical region matching method, its core logic based on the statistical similarity assumption of image texture: by setting a search window of fixed size on two consecutive time-series images, the statistical correlation of pixel grayscale distribution within the window is calculated to find the optimal matching position and determine the displacement vector at the window center. Although this technique has formed a standardized workflow for planar displacement extraction, its inherent window smoothing effect and dimensionality loss constitute insurmountable technical barriers when facing the complex dune height inversion task addressed in this invention. Specifically, from the perspective of spatial solution accuracy, this type of algorithm forcibly assumes that the surface units within the search window have a uniform motion state. This assumption is acceptable in flat areas, but it completely fails on the key geomorphic boundary of dune ridges, because the windward and leeward slopes on both sides of the ridge often exhibit drastically different or even opposite directions of material transport. Fixed-window integration inevitably involves spatially weighted averaging of these high-frequency, differentiated motion signals, leading to the neglect of displacement gradients at the most dynamically significant ridges, thus failing to accurately capture the true motion characteristics of micro-topography. From an observational perspective, this technology essentially remains at the level of two-dimensional kinematic geometry, only acquiring positional changes on the horizontal plane, completely lacking the crucial vertical height information for topographic evolution. This deficiency in the vertical dimension prevents the establishment of a complete height measurement model.

[0025] Therefore, existing technologies have significant shortcomings in dune monitoring and inversion, specifically in terms of limited observation accuracy, simplified physical models, and missing dimensions.

[0026] Based on this, the present invention provides a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model. Please refer to [link to relevant documentation]. Figure 1 , Figure 1 This is a flowchart illustrating a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model, provided by the present invention. The present invention provides a method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model, which includes: Step 1: Acquire initial Sentinel-2 optical images of the same target area at several different time phases.

[0027] In addition, in order to compare with the Sentinel-2 optical image selected in this embodiment, several initial Sentinel-1 SAR data that belong to the same target area and have the same time phase as the initial Sentinel-2 optical image were also selected.

[0028] Furthermore, since the initial Sentinel-2 optical image in the near-infrared band is most sensitive to the differences in humidity and particle texture on the sandy surface and can provide the best image contrast, this embodiment selects the near-infrared band from the initial Sentinel-2 optical image of multiple bands, specifically Band 8 with a resolution of 10m.

[0029] Step 2: Perform orthorectification and radiometric normalization on the initial Sentinel-2 optical image in sequence to obtain the preprocessed Sentinel-2 optical image.

[0030] Step 2.1: Based on the DEM (Digital Elevation Model) data of the target area, perform orthorectification on the initial Sentinel-2 optical image to obtain the corrected Sentinel-2 optical image.

[0031] Specifically, to ensure the accuracy of optical flow calculation, it is necessary to eliminate the geometric distortion and radiation differences of the data source itself. First, DEM data belonging to the same target area as the initial Sentinel-2 optical image is acquired, and the DEM data is imported into SNAP software. In SNAP software, the DEM data is used to perform orthorectification on the Sentinel-2 image. Specifically, based on the sensor geometric model of the DEM data and the initial Sentinel-2 optical image, the initial Sentinel-2 optical image is projected onto the specified map coordinate system through resampling to obtain the corrected Sentinel-2 optical image. The registration error between the corrected Sentinel-2 optical image and the ground control points is controlled within 0.1 pixels to eliminate terrain distortion.

[0032] Step 2.2: Perform radiometric normalization on the corrected Sentinel-2 optical image to obtain the preprocessed Sentinel-2 optical image.

[0033] Specifically, using the grayscale histogram of any corrected Sentinel-2 optical image as a radiometric reference, the histogram matching method is used to perform radiometric normalization on the grayscale histograms of the remaining corrected Sentinel-2 optical images to obtain the preprocessed Sentinel-2 optical images.

[0034] In other words, because the solar altitude angle varies in different seasons, the change in the shadow length of the sand dunes can interfere with optical flow calculations. Therefore, this embodiment selects a corrected Sentinel-2 optical image as a radiation reference image. The grayscale histogram of the near-infrared band of this radiation reference image is extracted. Similarly, the grayscale histogram of the corresponding near-infrared band of each of the remaining corrected Sentinel-2 optical images is extracted. Then, a histogram matching algorithm is used to adjust the grayscale histogram of each of the remaining corrected Sentinel-2 optical images to be statistically consistent with the grayscale histogram of the radiation reference image, thus obtaining several pre-processed Sentinel-2 optical images. This embodiment uses histogram matching to perform local contrast enhancement and radiation consistency processing on the images to satisfy the assumption of constant brightness in the optical flow method.

[0035] For comparison, this embodiment also preprocessed the Sentinel-1 SAR data. First, the Sentinel-1 SAR data was preprocessed by steps including orbit refinement to obtain the SLC standard format. Then, in order to remove the speckle noise inherent in the SAR data, nonlocal mean filtering was used and multi-look operation was performed to preserve the edge structure information of the dune ridge to the greatest extent. This obvious texture information will provide a high-quality base map for subsequent intensity flow tracking.

[0036] Step 3: Obtain the migration rate based on the relative displacement between the two preprocessed Sentinel-2 optical images.

[0037] Step 3.1: Process any two preprocessed Sentinel-2 optical images to obtain the relative displacement.

[0038] This embodiment uses two methods to determine the relative displacement between two preprocessed Sentinel-2 optical images.

[0039] The first method uses COSI-Corr (Co-registration of Optically Sensed Images and Correlation) to process any two preprocessed Sentinel-2 optical images to obtain the relative displacement.

[0040] Specifically, the principle of COSI-Corr is based on the Fourier translation theorem, which transforms the displacement problem in the spatial domain into a phase difference problem in the frequency domain. The two preprocessed Sentinel-2 optical images are denoted as I1 and I2, respectively. Therefore, the relationship between I1 and I2 can be obtained from the Fourier translation theorem:

[0041] in, The frequency domain of the preprocessed Sentinel-2 optical image I2. The frequency domain of the preprocessed Sentinel-2 optical image I1. The frequencies are the horizontal and vertical axes, respectively. , These represent the displacements of I2 relative to I1 in the horizontal and vertical directions, respectively. The normalized cross-power spectral density between I1 and I2 can be calculated from the above formula and expressed as:

[0042] in, To normalize the cross-power spectral density, This is the complex conjugate spectrum of I2.

[0043] To calculate the accurate displacement from the normalized cross-power spectral density, COSI-Corr constructs an objective function containing a weighting matrix. Its purpose is to measure the difference between the observed phase difference and the phase difference produced by the theoretical displacement, as shown in the following formula:

[0044] in, This is the frequency weighting function.

[0045] Therefore, when the objective function When the minimum value is reached, the corresponding offset is considered to be the relative displacement between the two effects.

[0046] The second method uses dense optical flow to process any two preprocessed Sentinel-2 optical images to obtain the relative displacement.

[0047] Specifically, the core theory of dense optical flow follows a variational framework of global energy minimization. It achieves a "coarse-to-fine" hierarchical estimation mechanism by constructing a multi-scale image pyramid, significantly improving the accuracy of displacement field calculation and enabling precise modeling of non-rigid dune deformation. The specific solution process involves constructing a multi-resolution pyramid structure for two preprocessed Sentinel-2 optical images, I3 and I4: first, an initial coarse-scale optical flow estimation is performed at the top layer of the pyramid (low-resolution image) to capture large-scale displacement trends; then, the coarse-scale displacement field is upsampled and passed to the next higher-resolution image layer, and this displacement field is used to perform image deformation operations. This "coarse-to-fine" stepwise optimization strategy effectively overcomes the problems of traditional optical flow methods easily getting trapped in local extrema and the ambiguity of weakly textured surface matching when dealing with large dune displacements. The energy function of dense optical flow can be expressed as follows:

[0048]

[0049] in, For data items, For smoothing terms, , All of these are regularization parameters. It is a convex function. The number of image channels. For local window functions, For the first i The horizontal gradient structure tensor of the channel. For the first i The cross-gradient structure tensor of the channel, For displacement field Spatial gradient, For displacement field Spatial gradient.

[0050] The sum of the data term and the smoothing term in the above formula is the total energy function. The data term, based on the assumptions of brightness and gradient conservation, constrains the feature consistency of the same dune surface unit in two images; the smoothing term introduces spatial neighborhood constraints, promoting smooth motion transitions between adjacent pixels. The core objective of the algorithm is to iteratively adjust the motion vector of each pixel based on the previously generated quasi-dense displacement field. , This allows for the minimization of the total energy function, enabling accurate calculation of the true displacement field and thus obtaining the relative displacement.

[0051] In addition, for comparison, this embodiment also uses the dense optical flow method to process the preprocessed Sentinel-1 SAR data to obtain the relative displacement.

[0052] Step 3.2: Obtain the migration rate based on the relative displacement and time interval.

[0053] Specifically, the migration rate is obtained by utilizing the time interval between two preprocessed Sentinel-2 optical images to measure the relative displacement between them.

[0054] Step 4: Based on the Martin & Kok sand flux model, calculate the saturated sand flux in the same time period as the migration rate.

[0055] In this embodiment, accurately calculating the sand flux driving dune movement is the physical basis for successful height inversion. This invention abandons the traditional Bagnold cubic model because it significantly overestimates sand transport at high wind speeds by neglecting the negative feedback effect of particles on the wind field. This invention employs the Martin & Kok sand flux model, which is more consistent with the physical mechanisms of aeolian sand transport and is a physically corrected model describing the abrupt shift in sand flux at high wind speeds.

[0056] Step 4.1: Obtain the 10-meter height wind speed and surface sediment transport dynamic roughness from the ERA5-Land meteorological reanalysis data of the target area.

[0057] Specifically, the ERA5-Land meteorological reanalysis data is a high-resolution land meteorological reanalysis dataset released by the European Centre for Medium-Range Weather Forecasts (ECMWF). Therefore, the 10-meter height wind speed of the target area can be extracted from the ERA5-Land meteorological reanalysis data, that is, the horizontal wind speed at 10 meters above the ground. The surface sand transport dynamic roughness is a parameter that reflects the influence of surface roughness on wind and sand movement.

[0058] Step 4.2: Based on the wind speed at a height of 10 meters and the surface sand transport dynamic roughness, the surface friction wind speed is derived by inversely using the preset formula corresponding to the logarithmic wind speed profile law. The preset formula is expressed as:

[0059] in, Indicates the height above the ground as Wind speed at that time Denotes the Kármán constant. Represents surface sediment transport dynamic roughness. This indicates the wind speed due to surface friction.

[0060] Therefore, by substituting the wind speed at a height of 10 meters and the surface sand transport dynamic roughness into the preset formula, the surface friction wind speed can be obtained. .

[0061] Step 4.3: Obtain the critical friction wind speed based on the critical friction wind speed calculation formula, which is expressed as:

[0062] in, This represents the critical frictional wind speed, which is the threshold at which sand grains are blown by the wind. Represents gravitational acceleration. Indicates the particle size of sand. Indicates the density of sand grains. This indicates air density.

[0063] Step 4.4: Based on the Martin & Kok sand flux model, obtain the saturated sand flux in the same time period as the migration rate, according to the surface friction wind speed and the critical friction wind speed.

[0064] Specifically, the surface frictional wind speed and critical frictional wind speed are substituted into the Martin & Kok sand flux model. This model takes into account that as wind speed increases, the jumping particles extract momentum from the wind field, leading to a decrease in near-surface wind speed, and the particle velocity eventually tends to saturate (i.e., it does not increase linearly with shear velocity). Therefore, the saturated sand flux... Wind speed due to friction with the ground The relationship between them is quadratic, not cubic; the saturated sand flux is expressed as:

[0065] in, Indicates saturated sand flux. This represents a dimensionless empirical coefficient.

[0066] By integrating wind data hourly over the study period, the modulus of the synthetic sediment transport potential (RDP) can be calculated as the input value for the annual average sediment flux. This improved model significantly enhances the physical accuracy of sediment flux estimation under high-energy wind conditions and eliminates the highly inverted systematic bias caused by traditional models.

[0067] Step 5: Based on the Exner equation in geomorphology, the dune height is obtained by inversion according to the migration rate and the saturated sand flux.

[0068] Step 5.1: Based on the Exner equation, construct the first relationship between the rate of change of surface elevation and the spatial divergence of sediment flux.

[0069] Specifically, the Exner equation is the Exner sediment mass conservation equation in geomorphology. The Exner equation describes the relationship between the rate of change of surface elevation and the spatial divergence of sediment flux. This relationship is the first relation, which is expressed as:

[0070] in, Indicates the surface elevation. Indicates time, Indicates the porosity of sand dunes. Represents sediment flux Spatial divergence.

[0071] Step 5.2: Obtain the second relationship based on the migration rate of the sand dunes.

[0072] Specifically, for a morphology that remains relatively stable and at a migration rate For sand dunes that migrate as a whole (i.e., satisfy the solitary wave hypothesis), the time derivative of their surface elevation can be transformed into the space derivative, which yields the second relation, expressed as:

[0073] in, Indicates the migration rate of sand dunes. The spatial gradient representing the height of sand dunes.

[0074] Step 5.3: Based on the first and second relationships, by integrating the one-dimensional profile along the wind direction, the inverse coupling relationship between dune height, dune migration rate and saturated sand flux at the top of the dune is derived to obtain the theoretical calculation formula for dune height.

[0075] Specifically, by substituting the second relation into the first relation and assuming that the sand flux at the dune base is zero, and integrating over a one-dimensional profile along the wind direction, the inverse coupling relationship between dune height, migration rate, and saturated sand flux at the dune crest can be derived, thus determining the theoretical formula for calculating dune height. This theoretical formula for calculating dune height is expressed as follows:

[0076] in, Indicates the height of the sand dunes. This indicates the overall density of sand dunes in the target area. This represents a dimensionless constant.

[0077] Step 5.4: Introduce a shape correction factor to modify the theoretical calculation formula for dune height into an engineering-usable dune height calculation formula, and calculate the relative height of the dune using the dune height calculation formula.

[0078] Specifically, to adapt to the non-ideal shape of sand dunes in actual desert environments (such as flux loss caused by flow field separation on the leeward slope), this embodiment introduces a shape correction factor to modify the above theoretical formula into an engineering-usable inversion model. This inversion model is the sand dune height calculation formula, which is expressed as:

[0079] in, This is the shape correction factor.

[0080] During implementation, the migration rate obtained in step 3 will be used. and the vector sand flux calculated in step 4 By substituting the values ​​into the dune height calculation formula, the relative height of the dunes can be calculated pixel by pixel.

[0081] Step 6: Using the flat Gobi Desert or sparse ICESat-2 laser altimeter points at the edge of the target area as absolute elevation control points, the relative height field of the dunes is corrected to generate an absolute digital elevation model with geographic coordinates, so as to reconstruct the three-dimensional dune landform.

[0082] Specifically, a three-dimensional relative height field of sand dunes is established based on the relative height of sand dunes at various points in the target area. Flat Gobi desert (with a mobility of 0 and height as the reference surface) or sparse ICESat-2 laser altimeter points at the edge of the target area are used as absolute elevation control points to correct the reference surface of the relative height field, thereby generating an absolute digital elevation model (DEM) with geographic coordinates. This process enables the reconstruction of three-dimensional terrain based solely on planar observation data, effectively solving the technical challenge of traditional remote sensing technology: "accurately measuring the plane but not the height."

[0083] This invention introduces a dense optical flow algorithm from the field of computer vision, abandoning the limitation of fixed window matching. This algorithm establishes a global energy functional containing data terms and smoothing terms, and uses a variational optimization framework to minimize this global energy, achieving physically continuous and dense pixel-by-pixel displacement estimation. This mechanism based on neighborhood smoothing constraints not only accurately captures high-frequency motion details of dune ridges, solving the smoothing blurring problem of traditional correlation methods, but also effectively transmits displacement information in weakly textured regions such as the desert interior. More importantly, this invention innovates at the physical model level, abandoning the simple Bagnold model and instead introducing the Martin & Kok improved sand flux model based on particle momentum conservation. This model fully considers the particle velocity saturation effect, significantly improving the physical accuracy of sand flux calculation. Based on this, this invention, using the Exner sediment mass conservation equation, couples a high-precision mobility field measured by remote sensing with a vector sand flux field derived from meteorology to establish a physical inversion mechanism for solving the absolute height distribution of dunes, thus achieving a leap from two-dimensional planar observation to three-dimensional terrain reconstruction. This invention enables the low-cost and high-efficiency inversion of the dynamic elevation distribution of large-scale desert areas without the need for ground-contact measurements, providing accurate scientific basis for wind and sand disaster prevention and environmental evolution research.

[0084] In this embodiment, during the experiments in steps 1, 2, and 3, the feather-shaped dune region of the Kumtag Desert was selected as the study area. The dune migration results of Sentinel-1 SAR data and Sentinel-2 optical imagery were calculated and obtained respectively. Figure 2 The results of annual average migration rate calculated from Sentinel-2 optical images using the dense optical flow method are presented. Figure 2 Figure a shows the migration rate results in the east-west direction, and Figure b shows the migration rate results in the north-south direction. Based on the results obtained in steps 1, 2, and 3, and following the physical derivation process in steps 4, 5, and 6, dune height inversion was performed at several study points in this region. The specific results are as follows: Figure 3 As shown, the result is in good agreement with the dune height calculated from the high-precision DEM data provided by the ZY-3 satellite, confirming the feasibility of the three-dimensional dune elevation calculation method provided in this invention.

[0085] This invention provides a method for calculating dune height based on the coupling of physical equations. The specific technical means are as follows: external meteorological wind speed data is substituted as initial conditions into the Martin & Kok sand flux physical model to calculate the near-surface sand flux; subsequently, the two-dimensional horizontal displacement field of the dune obtained through remote sensing imagery and the calculated sand flux value are input into the Exner geomorphic mass conservation partial differential equation to calculate the vertical height variable of the dune. In the step of solving for the sand flux, this invention explicitly uses the Martin & Kok quadratic model, which includes a negative feedback mechanism for the momentum of aeolian sand flow-transferred particles, to replace the conventional Bagnold cubic model. Through the correction of the physical mechanism, the overestimation bias of the sand flux calculation results by traditional empirical formulas under high wind speed conditions is eliminated. This invention does not require LiDAR or other measured elevation data; it can complete the numerical conversion to three-dimensional vertical height solely based on the displacement in a two-dimensional plane.

[0086] This invention protects a multi-method fusion technique for extracting ground displacement fields from multi-source data. The technique employs both traditional frequency-domain COSI-Corr and dense optical flow methods. It utilizes the energy functional optimization direction constrained by frequency-domain correlation results to prevent divergence in large displacement regions, and leverages the pixel-by-pixel computational capability of dense optical flow to fill in high-frequency details. This addresses the window smoothing and matching noise issues encountered by traditional frequency-domain methods in low-texture areas and dune ridges.

[0087] In the description of this invention, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0088] Although the invention has been described herein in conjunction with various embodiments, those skilled in the art will understand and implement other variations of the disclosed embodiments by reviewing the accompanying drawings, disclosure, and appended claims in carrying out the claimed invention. In the claims, the word "comprising" does not exclude other components or steps, and "a" or "an" does not exclude a plurality. A single processor or other unit can implement several functions listed in the claims. While different dependent claims may recite certain measures, this does not mean that these measures cannot be combined to produce good results.

[0089] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, any modifications made without departing from the inventive concept should be considered within the scope of protection of the present invention.

Claims

1. A method for calculating the three-dimensional elevation of sand dunes by coupling multi-source remote sensing with a physical model, characterized in that, include: Acquire initial Sentinel-2 optical images of the same target area at several different time phases; The initial Sentinel-2 optical image is subjected to orthorectification and radiometric normalization in sequence to obtain the preprocessed Sentinel-2 optical image; The migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images. Based on the Martin & Kok sand flux model, the saturated sand flux in the same time period as the migration rate is calculated; Based on the Exner equation in geomorphology, the dune height is obtained by inversion according to the migration rate and the saturated sand flux.

2. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, All the initial Sentinel-2 optical images are in the near-infrared band.

3. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, The initial Sentinel-2 optical image is sequentially subjected to orthorectification and radiometric normalization to obtain a preprocessed Sentinel-2 optical image, including: Based on the DEM data of the target area, the initial Sentinel-2 optical image is orthorectified to obtain the corrected Sentinel-2 optical image; The corrected Sentinel-2 optical image is subjected to radiometric normalization to obtain the preprocessed Sentinel-2 optical image.

4. The method for calculating the three-dimensional elevation of sand dunes according to claim 3, characterized in that, The corrected Sentinel-2 optical image is subjected to radiometric normalization to obtain the preprocessed Sentinel-2 optical image, including: Using the grayscale histogram of any one of the corrected Sentinel-2 optical images as a radiometric reference, the histogram matching method is used to perform radiometric normalization on the grayscale histograms of the remaining corrected Sentinel-2 optical images to obtain the preprocessed Sentinel-2 optical images.

5. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, The migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images, including: The relative displacement is obtained by processing any two preprocessed Sentinel-2 optical images using COSI-Corr. The migration rate is obtained based on the relative displacement and time interval.

6. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, The migration rate is obtained based on the relative displacement between the two preprocessed Sentinel-2 optical images, including: The relative displacement is obtained by processing any two preprocessed Sentinel-2 optical images using the dense optical flow method. The migration rate is obtained based on the relative displacement and time interval.

7. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, Based on the Martin & Kok sand flux model, the saturated sand flux for the same time period as the migration rate is calculated, including: Obtain the 10-meter height wind speed and surface sediment transport dynamic roughness from the ERA5-Land meteorological reanalysis data of the target area; Based on the wind speed at a height of 10 meters and the surface sand transport dynamic roughness, the surface friction wind speed is calculated by inversely using a preset formula corresponding to the logarithmic wind speed profile law. The preset formula is expressed as follows: in, Indicates the height above the ground as Wind speed at that time Denotes the Kármán constant. Represents surface sediment transport dynamic roughness. Indicates the wind speed due to surface friction; The critical friction wind speed is obtained based on the critical friction wind speed calculation formula, which is expressed as follows: in, This represents the critical frictional wind speed. Represents gravitational acceleration. Indicates the particle size of sand. Indicates the density of sand grains. Indicates air density; Based on the Martin & Kok sand flux model, the saturated sand flux in the same time period as the migration rate is obtained according to the surface friction wind speed and the critical friction wind speed.

8. The method for calculating the three-dimensional elevation of sand dunes according to claim 7, characterized in that, The saturated sand flux is expressed as: in, Indicates saturated sand flux. This represents a dimensionless empirical coefficient.

9. The method for calculating the three-dimensional elevation of sand dunes according to claim 1, characterized in that, Based on the Exner equation in geomorphology, the dune height is obtained by inversion according to the migration rate and the saturated sand flux, including: Based on the Exner equation, a first relationship is constructed between the spatial divergence of the rate of change of surface elevation and sediment flux, which is expressed as: in, Indicates the surface elevation. Indicates time, Indicates the porosity of sand dunes. Represents sediment flux Spatial divergence; The second relationship is derived based on the dune migration rate, and is expressed as follows: in, Indicates the migration rate of sand dunes. A spatial gradient representing the height of sand dunes; Based on the first and second relationships, by integrating a one-dimensional profile along the wind direction, the inverse coupling relationship between dune height, dune migration rate, and saturated sand flux at the dune crest is derived, resulting in a theoretical formula for calculating dune height. This theoretical formula is expressed as follows: in, Indicates the height of the sand dunes. Indicates saturated sand flux. This indicates the overall density of sand dunes in the target area. Represents a dimensionless constant; A shape correction factor is introduced to modify the theoretical dune height calculation formula into an engineering-usable dune height calculation formula. The relative height of the dune is then calculated using this formula, which is expressed as follows: in, Represents dimensionless empirical coefficients. This represents the critical frictional wind speed. Indicates wind speed due to surface friction. Represents gravitational acceleration. Indicates air density, This is the shape correction factor.

10. The method for calculating the three-dimensional elevation of sand dunes according to claim 9, characterized in that, After obtaining the dune height through inversion based on the migration rate and the saturated sand flux, the process further includes: By using flat Gobi desert or sparse ICESat-2 laser altimetry points at the edge of the target area as absolute elevation control points, the relative height field of the dunes is corrected to a reference surface, generating an absolute digital elevation model with geographic coordinates to reconstruct the three-dimensional dune topography.

Images