A water depth inversion method and device coupling wave dynamics mechanism and chart constraints

Abstract

The present application provides a kind of coupling wave dynamics mechanism and chart constraint's water depth inversion method and device, it is related to ocean remote sensing and coastal hydrology technical field.The method includes: obtaining the multispectral remote sensing image of target water area and corresponding chart sounding data respectively extract pixel gray intensity, chart priori water depth field to construct multidimensional input matrix;Pixel gray intensity is carried out spatial domain frequency transformation, and constructs regional wave dynamic characteristic field;In combination with regional wave dynamic characteristic field, construct the mapping model about water depth, with deep water reference area as benchmark to determine the initial inversion water depth field, and to the initial inversion water depth field is carried out abnormal data rejection, obtains candidate inversion water depth field;Candidate inversion water depth field is carried out multi-time phase integration, and in combination with chart priori water depth field is carried out dynamic weight fusion, obtains fusion water depth field, and in combination with measured water depth value is carried out deviation correction, obtains final inversion water depth field.The present application improves the high-precision estimation of water depth in turbid water area.

Description

Technical Field

[0001] This invention relates to the fields of marine remote sensing and coastal hydrology, and in particular to a method and apparatus for water depth inversion that couples wave dynamics mechanisms with nautical chart constraints. Background Technology

[0002] Water depth information in turbid shallow waters is crucial for water resource management, ecological environment monitoring, and waterway planning. However, traditional optical bathymetry methods are highly dependent on water transparency. In turbid waters or medium-deep waters (approximately 10–30 m), the light signal attenuates rapidly, leading to a significant decrease in inversion accuracy or even failure. While wave theory-based water depth inversion methods are theoretically unaffected by water quality, in practical applications, they often suffer from inaccurate extraction of swell characteristics, missing spatial topological relationships, and insufficient physical constraints, resulting in unsatisfactory water depth inversion accuracy. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention proposes a water depth inversion method and apparatus that couples wave dynamics mechanisms with nautical chart constraints. This method integrates multi-temporal satellite imagery with wave dynamics analysis to achieve high-precision estimation and spatial distribution prediction of water depth in turbid waters.

[0004] This invention provides a method for water depth inversion that couples wave dynamics mechanisms with chart constraints, characterized in that it comprises: Acquire multispectral remote sensing images of the target water area and corresponding nautical chart bathymetry data; Pixel gray intensity is extracted from the multispectral remote sensing image, and the prior water depth field of the nautical chart is extracted from the nautical chart bathymetry data. A multidimensional input matrix is ​​constructed by combining the pixel gray intensity and the prior water depth field of the nautical chart. Spatial domain frequency transformation is performed on the pixel grayscale intensity to construct a regional wave dynamic feature field; Based on the wave dynamic characteristic field of the region, a mapping model for water depth is constructed. The global angular frequency is obtained by inversion with the deep water reference area as the reference. The initial inversion water depth field is determined by combining the global angular frequency. Abnormal data is removed from the initial inversion water depth field to obtain the candidate inversion water depth field. The candidate inverted water depth fields are integrated across multiple time phases and dynamically weighted in conjunction with the prior water depth fields from the nautical charts to obtain a fused water depth field. The fused water depth field is then corrected for deviations by incorporating measured water depth values ​​to obtain the final inverted water depth field.

[0005] In one embodiment of the present invention, the steps of extracting pixel grayscale intensity from the multispectral remote sensing image, extracting prior water depth field from the nautical chart bathymetry data, and constructing a multidimensional input matrix by combining the pixel grayscale intensity and the prior water depth field include: Atmospheric correction is performed on the multispectral remote sensing image to obtain the reflectance data of the lower atmosphere; The preset blue and red light bands are extracted from the reflectance data of the lower atmosphere as spectral reflectance data sources; The nautical chart bathymetry data is resampled so that the resampled nautical chart bathymetry data has the same spatial resolution as the multispectral remote sensing image. Spatial coordinate alignment is performed between the spectral reflectance data source and the resampled nautical chart bathymetry data; The spatially aligned spectral reflectance data source is linearly converted into the pixel grayscale intensity. Extract the prior water depth field of the target area from the spatially aligned chart bathymetry data; By integrating the pixel grayscale intensity, prior water depth field of the nautical chart, and image spatial resolution, the multidimensional input matrix is ​​constructed, wherein the multidimensional input matrix is ​​represented as follows: , in, For a multidimensional input matrix, For geospatial pixel coordinate indexing, The pixel grayscale intensity is used to characterize the geometric shape of sea surface ripples; Pre-examination of water depth field for nautical charts This refers to the spatial resolution of the image.

[0006] In one embodiment of the present invention, the step of performing spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field includes: Based on Radon transform, linear integral projection is performed on each pixel of the pixel gray intensity along multiple azimuth angles to map the two-dimensional spatial ripple texture into a one-dimensional projection vector, and the one-dimensional projection vectors from different angles are integrated to construct a sine curve. A cross-spectrum matrix is ​​constructed from the main image and the secondary image of the corresponding band of the sine curve. The phase difference gradient of each frequency component is calculated. The amplitude of the cross-spectrum is weighted and reconstructed using the phase difference gradient to eliminate 180° directional ambiguity and determine the direction of the wave motion vector. Using the direction of the wave motion vector as a constraint, the pixel gray intensity is processed by a two-dimensional discrete Fourier transform within a sliding window to calculate the wave energy spectrum. Identify the spectral peaks of the wave energy spectrum, and calculate the combiner number and main wave direction of the main wave group by combining the frequency domain coordinates of the spectral peaks with the spatial resolution parameters of the image, and calculate the wavelength from the combiner number; By traversing the target water area through a step-by-step sliding window, a regional wavelength distribution field and a main wave direction distribution field are generated, thus constructing a wave dynamic characteristic field.

[0007] In one embodiment of the present invention, the steps of identifying the spectral peaks of the wave energy spectrum, calculating the combiner number and main wave direction of the main wave group by combining the frequency domain coordinates of the spectral peaks with the spatial resolution parameters of the image, and calculating the wavelength from the combiner number include: The peaks in the wave energy spectrum were identified using an extreme value search algorithm; Based on the frequency domain coordinates of the spectral peaks and the spatial resolution of the image, the physical wavenumber components in the horizontal and vertical directions of the main wave group are calculated using the following formula: , in, The physical wavenumber component is in the horizontal direction. The physical wavenumber component is in the vertical direction. The frequency domain coordinates of the spectral peak, For image spatial resolution, This represents the total number of pixels contained in the sliding window in the horizontal direction. This represents the total number of pixels contained in the sliding window in the vertical direction. The combined wavenumber of the main wave group is calculated based on the physical wavenumber components in the horizontal and vertical directions. The calculation formula is: , Based on the physical wavenumber components in the horizontal and vertical directions, the main wave direction of the main wave group is calculated. The calculation formula is: , The wavelength is calculated based on the combined wave number, using the following formula: , in, λ is the wavelength.

[0008] In one embodiment of the present invention, the wave dynamic characteristic field is represented as: , in, This is the discrete coordinate index of the sliding window center in the geospatial grid. For wave dynamic characteristic field, This represents the wavelength value at that coordinate. This represents the dominant wave direction at that coordinate. The peak confidence parameter is defined as the ratio of the peak value of the energy spectrum to the mean value of the background.

[0009] In one embodiment of the present invention, before performing spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field, the method further includes: Adaptive texture enhancement is performed on the pixel grayscale intensity in the multidimensional input matrix to improve the grayscale contrast between peaks and troughs; A bandpass filtering algorithm is used to remove ultra-long period background signals and high-frequency noise, while retaining effective wave signals modulated by terrain, to obtain the pixel gray intensity after denoising and enhancement.

[0010] In one embodiment of the present invention, the steps of constructing a mapping model about water depth by combining the regional wave dynamic characteristic field, obtaining the global angular frequency by inversion based on the deep-water reference area, and determining the initial inverted water depth field by combining the global angular frequency include: A mapping model between angular frequency, resultant wave number, and water depth is established based on linear wave theory, and is expressed as follows: , in, For global angular frequency, It is the acceleration due to gravity. The combination number of the main wave group. This represents the water depth. Identify deep-water reference regions that meet the deep-water conditions in the multi-dimensional input matrix, and extract the deep-water wavenumber of the deep-water reference regions. Based on the deep-water wave number, and under deep-water conditions Based on the constraint characteristics, the global angular frequency is obtained by inversely solving the mapping model, and the calculation formula is as follows: , in, For deep-water wave number; The Newton-Raphson iterative algorithm is used to solve the transcendental equation of the linear wave dispersion relation. After the iteration converges, the initial inverted water depth field is obtained. The calculation formula is as follows: , Wherein, objective function derivative term , Let be a hyperbolic secant function, and the iteration termination condition is: Less than the preset threshold; and The first Next and first The water depth value in the first iteration.

[0011] In one embodiment of the present invention, the step of removing outlier data from the initial inverted water depth field to obtain a candidate inverted water depth field includes: For each pixel in the initial inverted depth field, the significance of the residual is quantified using the Wald statistical test, and the calculation formula is as follows; , in: The Wald statistic is used to assess the degree to which the inversion results deviate from the physical mechanism or prior constraints. This value follows a set of 1 degree of freedom. distributed; This is for retrieving water depth values; This serves as a reference value for water depth. The variance of the inverted water depth values; The inverted water depth values ​​of the pixels that pass the physical consistency check are retained. For the removed abnormal pixels, spatial interpolation is performed using the inverted water depth values ​​of the neighboring valid pixels to fill the space, thus obtaining the candidate inverted water depth field.

[0012] In one embodiment of the present invention, the steps of performing multi-temporal integration on the candidate inverted depth field, dynamically weighting and fusing it with the prior nautical chart depth field to obtain a fused depth field, and correcting the deviation of the fused depth field with measured depth values ​​to obtain the final inverted depth field include: The candidate inversion water depth fields of the target water area in multiple time phases are obtained, and the multi-time phase integrated inversion water depth field is generated by median synthesis technique; By introducing a dynamic weighting function, the multi-temporal integrated inversion water depth field and the prior water depth field of the nautical chart are fused to obtain the fused water depth field; Using local residual correction techniques, and combining the measured water depth values ​​of known control points, the fused water depth field is checked for deviations to obtain the final inverted water depth field. The calculation formula is as follows: , in: This is the final inverted water depth field after verification; To integrate the deep water field; The local residual correction term is defined as the spatial interpolation residual field between the measured water depth and the inverted water depth at the known control point. It is a geospatial cell coordinate index used to ensure that the correction operation corresponds to each cell across the entire target water area.

[0013] In another aspect, the present invention provides a water depth inversion device that couples wave dynamics mechanisms with chart constraints, employing the aforementioned water depth inversion method that couples wave dynamics mechanisms with chart constraints. The device comprises: The acquisition module is used to acquire multispectral remote sensing images of the target water area and the corresponding nautical chart bathymetry data; The dataset construction module is used to extract pixel gray intensity from the multispectral remote sensing image, extract the prior water depth field of the nautical chart from the nautical chart bathymetry data, and construct a multidimensional input matrix by combining the pixel gray intensity and the prior water depth field of the nautical chart. The wave kinematics feature extraction module is used to perform spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field. The inversion modeling module is used to construct a mapping model about water depth by combining the wave dynamic characteristic field of the region, inverting the global angular frequency based on the deep water reference area, determining the initial inversion water depth field by combining the global angular frequency, and removing abnormal data from the initial inversion water depth field to obtain the candidate inversion water depth field. An optimization module is used to perform multi-temporal integration of the candidate inverted water depth field, and to perform dynamic weighted fusion with the prior water depth field of the nautical chart to obtain a fused water depth field. The fused water depth field is then corrected for deviations by combining the measured water depth values ​​to obtain the final inverted water depth field.

[0014] As can be seen from the above solutions, the advantages of the present invention are: The water depth inversion method coupled with wave dynamics mechanism and nautical chart constraints provided by this invention extracts pixel gray intensity from the multispectral remote sensing image and extracts wave dynamic feature field through spatial domain frequency transformation of pixel gray intensity. It breaks through the dependence of traditional remote sensing inversion on the synchronization of multi-source remote sensing data. It can extract complete wave dynamic field features using only a single image without the need for synchronous auxiliary observation data, which greatly reduces the data cost and application threshold of inversion. Combined with the wave dispersion mechanism, a water depth mapping model is constructed, which solves the bottleneck problem of weak optical spectral signals and limited inversion depth in turbid waters, and extends the effective inversion depth to the medium and deep water area. Simultaneously, a priori depth field is extracted from the chart bathymetry data. A multi-dimensional input matrix is ​​constructed by combining the pixel grayscale intensity and the priori depth field. This priori depth field serves as the physical constraint for the entire inversion process, rather than merely as a later verification benchmark. By constraining the range of wave parameters using the topographic framework provided by the chart, and serving as the initial values ​​for iterative solutions and a reference benchmark for the fusion stage, the model's predictive stability in mid-deep water is significantly improved, avoiding the divergence problem of pure data-driven inversion in deep water. Furthermore, the candidate inversion depth fields are integrated across multiple time phases, and dynamic weight fusion is performed in conjunction with the priori depth field. This combination of multi-time phase integration and dynamic weight fusion optimizes the inversion results, utilizing wave motion laws to compensate for the attenuation of optical signals in turbid water. It achieves stable and high-precision depth measurement, and still has good applicability in turbid nearshore waters with high suspended matter. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the overall process of a water depth inversion method that combines wave dynamics mechanism and chart constraints according to an embodiment of the present invention. Figure 2 A schematic diagram of the selected target water area; Figure 3 for Figure 1 A detailed flowchart of step S2 is shown below; Figure 4The parameters for wave depth inversion are set; where (a) is the spatial wave field (BO4), and (b) is the two-dimensional energy spectrum in polar coordinates obtained after converting to the frequency domain by two-dimensional fast Fourier transform (2D-FFT); Figure 5 A schematic diagram illustrating the wave texture enhancement and noise reduction effects; Figure 6 for Figure 1 A detailed flowchart of step S3 is shown below; Figure 7 This is a schematic diagram illustrating the transformation effect of Bin Laden. Figure 8 This is a schematic diagram illustrating the results of cross-spectral analysis. Figure 9 The diagram shows the effect of spatial discrete Fourier transform; where (a) is a wave modulated by terrain refraction / diffraction, (b) is a regular striped wave, (c) is a weak wave with low signal-to-noise ratio, and (d) is a non-uniform shallow water wave. Figure 10 This is a schematic diagram of the polar coordinate spectral space. Figure 11 for Figure 1 A detailed flowchart of step S5 is shown below; Figure 12 This is a schematic diagram illustrating the final inversion of the water depth field. Figure 13 This is a sampling map of spatial reference points for the final inversion of water depth field data; Figure 14 This is a spatial reference point sampling map for multibeam bathymetry data. Figure 15 A comparison chart of inverted water depth values ​​and multibeam measured water depth values; Figure 16 A general structural block diagram of a water depth inversion device provided in another embodiment of the present invention is shown.

[0016] The attached figures are labeled as follows: 510: Acquisition module; 520: Dataset building module; 530: Wave kinematics feature extraction module; 540: Inversion Modeling Module; 550: Optimization module. Detailed Implementation

[0017] It should be noted that, in this invention, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus.

[0018] In the absence of further restrictions, an element defined by the phrase "comprising a..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0019] refer to Figure 1 As shown, Figure 1 The diagram shows the overall flow of the water depth inversion method that combines wave dynamics mechanism and chart constraints provided in one embodiment of the present invention.

[0020] A water depth inversion method that couples wave dynamics mechanisms with chart constraints includes the following steps: Step S1: Acquire multispectral remote sensing images of the target water area and the corresponding nautical chart bathymetry data.

[0021] In one embodiment, Sentinel-2 multispectral remote sensing imagery covering a target water area is acquired. This target water area can be selected from nearshore sea areas, such as... Figure 2 As shown, Figure 2 A schematic diagram of the selected target water area is shown. High spatial resolution near-infrared or visible light bands of multispectral remote sensing imagery are used to capture the texture features of sea surface waves. Simultaneously, nautical chart bathymetry data of the target water area is collected. This data can be obtained from historical nautical charts or global bathymetry data (GEBCO chart data), serving as the physical benchmark and prior constraint for water depth inversion.

[0022] Step S2: Extract pixel gray intensity from the multispectral remote sensing image, extract the prior water depth field from the nautical chart bathymetry data, and construct a multidimensional input matrix by combining the pixel gray intensity and the prior water depth field.

[0023] In one embodiment, reference Figure 3 As shown, Figure 3 A schematic diagram of the specific process for step S2 is shown.

[0024] Step S21: Perform atmospheric correction on the multispectral remote sensing image to obtain bottom atmospheric reflectance data. In one specific implementation, select an observation period with cloud cover below 10%, and use the Sen2Cor plugin to perform atmospheric correction on the multispectral remote sensing image to obtain bottom atmospheric reflectance data (BOA).

[0025] Step S22: Extract the preset blue and red light bands from the bottom atmospheric reflectance data as spectral reflectance data sources. In a specific implementation, the blue light band (B2 band) and red light band (B4 band) with a spatial resolution of 10m are selected as the basic data sources; the high penetration of the B2 band is used to capture topographic fluctuations in the mid-deep water area, and combined with the high contrast characteristics of the B4 band in turbid water, the wave train characteristics modulated by the complex underwater topography are fully covered.

[0026] Step S23: Resample the nautical chart bathymetry data to ensure that the resampled nautical chart bathymetry data has the same spatial resolution as the multispectral remote sensing image. In one specific implementation, a bilinear interpolation algorithm is used to uniformly resample the low-resolution nautical chart bathymetry data to a spatial resolution of 10m, ensuring that the spatial scale is completely aligned with the pixels of the Sentinel-2 multispectral remote sensing image.

[0027] Step S24: Spatial coordinate alignment is performed between the spectral reflectance data source and the resampled nautical chart bathymetry data to unify all data into the WGS84 / UTM coordinate system and ensure that the pixel geographic location error is controlled within 0.5 pixels.

[0028] Step S25: Linearly convert the spatially aligned spectral reflectance data source into the pixel grayscale intensity. It is used to characterize the geometric shape of sea surface ripples.

[0029] Step S26: Extract the prior water depth field of the target area from the spatially aligned nautical chart bathymetry data.

[0030] Step S27: Integrate the pixel grayscale intensity, prior water depth field of the nautical chart, and image spatial resolution to construct the multidimensional input matrix, wherein the multidimensional input matrix is ​​represented as: , in, For a multidimensional input matrix, For geospatial pixel coordinate indexing, The pixel grayscale intensity is used to characterize the geometric shape of sea surface ripples; Pre-examination of water depth field for nautical charts The prior water depth values ​​include the coordinates of all pixels in the target area. A two-dimensional matrix; This refers to the spatial resolution of the image, such as 10m.

[0031] The multidimensional input matrix contains spectral texture and topographic prior features, serving as the core data operator for subsequent wave dynamics evolution and water depth calculations. The relevant parameter settings are as follows: Figure 4 As shown, Figure 4 The wave depth inversion parameter settings are shown.

[0032] Step S3: Perform spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field.

[0033] Spatial domain frequency transformation is performed on the pixel gray intensity to transform the static pixel gray intensity information in the multidimensional input matrix into dynamic features with physical meaning, including key kinematic features such as wavelength, combined wave number and main wave direction. These dynamic features are swell parameters modulated by seabed topography.

[0034] Furthermore, in one embodiment, before constructing the regional wave dynamic feature field by performing spatial domain frequency transformation on the pixel grayscale intensity, wave texture enhancement and denoising processing are performed on the pixel grayscale intensity in the multidimensional input matrix. Specifically, for non-wave noise generated by high suspended matter in turbid waters, an adaptive enhancement algorithm is used to improve the contrast between wave crests and troughs; then, a bandpass filtering algorithm is used to remove ultra-long ocean current background and high-frequency sensor noise, retaining the swell signal modulated by topography, to obtain the denoised and enhanced pixel grayscale intensity, such as... Figure 5 As shown, Figure 5 The diagram illustrates the wave texture enhancement and noise reduction effects.

[0035] Then, spatial domain frequency transformation is performed on the denoised and enhanced pixel grayscale intensity to construct a regional wave dynamic feature field. Through spatial domain frequency transformation, the static pixel grayscale intensity information in the multidimensional input matrix is ​​transformed into dynamic features, including wavelength, wave combination number, and dominant wave direction. This process does not rely on water transparency; its core lies in identifying the spatial evolution characteristics of wave trains in the image.

[0036] Specifically, in one embodiment, reference Figure 6 As shown, Figure 6 A schematic diagram of the specific process for step S3 is shown.

[0037] Step S31: Based on the Radon transform, linear integral projection is performed on each pixel of the pixel gray intensity along multiple azimuth angles to map the two-dimensional spatial ripple texture into a one-dimensional projection vector, and the one-dimensional projection vectors from different angles are integrated to construct a sine map. In a specific embodiment, to extract coherent wave fields from nearshore backgrounds with high spatial non-stationarity, the Radon transform is introduced to perform geometric feature enhancement. By performing linear integral projection on each pixel covered by the circular mask of pixel gray intensity along the azimuth angle θ∈(0,π), the ripple texture in two-dimensional space is mapped into a one-dimensional projection vector. This process utilizes the physical coherent superposition effect along the wavefront direction to effectively suppress bottom heterogeneity and random radiation noise, enabling the wave signal to achieve directional enhancement in complex background textures. Based on the Radon transform, a sine map is constructed by integrating one-dimensional projection vectors from different scanning angles. Figure 7 The bin Laden transform effect for B2 and B4 bands of Sentinel-2 is shown.

[0038] Step S32: Construct a cross-spectral matrix for the main and secondary images corresponding to the sine wave band, calculate the phase difference gradient of each frequency component, and use the phase difference gradient to reconstruct the amplitude of the cross-spectrum using weighted recalculation to eliminate 180° directional ambiguity and determine the direction of wave motion vector. In this embodiment, cross-spectral analysis is introduced to further determine the physical propagation direction of the wave and eliminate 180° directional ambiguity. By constructing a cross-spectral matrix between the main and secondary images, the phase difference gradient of each frequency component is calculated. Using the phase difference gradient to reconstruct the amplitude of the cross-spectrum using weighted recalculation can accurately lock the true direction of wave motion vector, and combined with polynomial fitting, a sub-pixel angular spectral resolution better than 1° is achieved, ensuring the physical robustness of direction extraction. The cross-spectral analysis effect is as follows: Figure 8 As shown.

[0039] Step S33: Using the wave motion vector direction as a constraint, process the pixel grayscale intensity within the sliding window using a two-dimensional discrete Fourier transform to calculate the wave energy spectrum. The specific calculation formula is as follows: , Where f(x,y) is the pixel brightness value within the window, and u,v are the spatial frequency components; This represents the total number of pixels contained in the sliding window in the horizontal direction. This represents the total number of pixels contained in the sliding window in the vertical direction. Figure 9 The spatial discrete Fourier transform effect is shown.

[0040] Step S34: Identify the spectral peaks of the wave energy spectrum, calculate the combination number and main wave direction of the main wave group by combining the frequency domain coordinates of the spectral peaks with the spatial resolution parameters of the image, and calculate the wavelength from the combination number.

[0041] In one specific implementation, an extreme value search algorithm is used to identify spectral peaks in the wave energy spectrum. These peaks represent the points of highest energy concentration in the wave energy spectrum, and their frequency domain coordinates are obtained. This coordinate represents the discrete position of the main wave group in the frequency domain.

[0042] Then, based on the frequency domain coordinates of the spectral peak With the image spatial resolution The physical wavenumber components in the horizontal and vertical directions of the main wave group are calculated using the following formula: , in, The physical wavenumber component is in the horizontal direction. The physical wavenumber component is in the vertical direction. The frequency domain coordinates of the spectral peak, For image spatial resolution, This represents the total number of pixels contained in the sliding window in the horizontal direction. This represents the total number of pixels contained in the sliding window in the vertical direction.

[0043] Based on the physical wavenumber components in the horizontal direction Vertical physical wavenumber components Calculate the combination number of the main wave group. The calculation formula is: , Based on the physical wavenumber components in the horizontal direction Vertical physical wavenumber components Calculate the main wave direction of the main wave group Main wave direction The formula represents the propagation direction of the surge in the geographic coordinate system, and is as follows: , The wavelength is calculated based on the combined wave number, using the following formula: , in, λ is the wavelength.

[0044] Through the above calculations, the texture features of the image are transformed into a core parameter field describing the dynamic evolution of waves, including wavelength, wave number, and dominant wave direction, such as... Figure 10 As shown, Figure 10 A schematic diagram of the polar coordinate spectral space is shown.

[0045] Step S35: Traverse the entire target water area through a step-sliding window to generate a regional wavelength distribution field and a main wave direction distribution field, and construct a wave dynamic characteristic field.

[0046] The matrix representation of the wave dynamic characteristic field is as follows: , in, This is the discrete coordinate index of the sliding window center in the geospatial grid. For wave dynamic characteristic field; The wavelength value at that coordinate is the wavelength distribution field, which is a two-dimensional matrix of wavelength values ​​corresponding to all pixel coordinate points in the entire target area. The dominant wave direction at this coordinate is a two-dimensional matrix of the dominant wave directions corresponding to all pixel coordinate points in the entire target area, which is used to help identify the evolution of the propagation direction caused by wave refraction. The peak confidence parameter is defined as the ratio of the peak value of the energy spectrum to the mean value of the background, and serves as the dynamic weight for the reliability of the weighted fusion in the subsequent step S5.

[0047] Step S4: Combine the wave dynamic characteristic field of the region to construct a mapping model about the water depth, invert the global angular frequency with the deep water reference area as the reference, determine the initial inverted water depth field with the global angular frequency, and remove abnormal data from the initial inverted water depth field to obtain the candidate inverted water depth field.

[0048] In one embodiment, a mapping model between angular frequency, combined wave number, and water depth is first established based on linear wave theory (Airy theory), expressed as: , in, For global angular frequency, It is the acceleration due to gravity. The combination number of the main wave group. This represents the water depth. Because it is difficult to directly obtain the period from single-phase images The water depth value needs to be solved iteratively by combining prior information or wave velocity relationships in deep water areas. .

[0049] Identify deep-water reference regions that meet the deep-water conditions in the multidimensional input matrix, and extract the deep-water wavenumbers of these reference regions. In one specific implementation method, the deep-water conditions are set as follows: .

[0050] Then, based on the deep-water wave number, under deep-water conditions... Based on the constraint characteristics, the global angular frequency is obtained by inversely solving the mapping model, and the calculation formula is as follows: , in, This represents the deep-water wave number.

[0051] Based on the global angular frequency, the transcendental equation of the linear wave dispersion relation is solved using the Newton-Raphson iterative algorithm. After the iteration converges, the initial inverted depth field is obtained, which includes information about the deep-water reference area. The iterative process is based on the objective function. The zero-point optimization is achieved, and the calculation formula is as follows: , , Wherein, objective function The physical essence is the residual between the local wave dynamic parameters and the theoretical dispersion relationship. ; derivative term , , It is a hyperbolic secant function. and The first Next and first The water depth value in one iteration, and the iteration termination condition is: If the depth is less than the preset threshold, the initial value can be determined from the prior water depth value on the nautical chart during the solution process.

[0052] Furthermore, for each pixel in the initial inverted depth field, the Wald statistical test is used to quantify the significance of the residuals. The Wald statistical test quantifies the residuals of the inverted depth values ​​deviating from the linear dispersion curve, and the calculation formula is as follows: , in: The Wald statistic is used to assess the degree to which the inversion results deviate from the physical mechanism or prior constraints. This value follows a set of 1 degree of freedom. distributed; This is for retrieving water depth values; The reference value for water depth can be determined from prior water depth values ​​on nautical charts or from theoretical water depth values ​​calculated using mapping equations based on linear wave theory. The variance of the inverted water depth value is used to characterize the uncertainty in the inversion process. It is obtained by comprehensively evaluating the signal-to-noise ratio of the spectral analysis within the sub-window and the convergence residual of the numerical iteration.

[0053] Remove significance level Abnormal pixel points are identified to ensure that all retained data points conform to the basic laws of fluid mechanics.

[0054] The inverted water depth values ​​of the pixels that pass the physical consistency check are retained. For the removed abnormal pixels, spatial interpolation is performed using the inverted water depth values ​​of the neighboring valid pixels to fill the gap, thereby obtaining the spatially continuous candidate inverted water depth field.

[0055] In this embodiment, based on the linear wave dispersion relation, the extracted wave dynamic features are converted into quantitative water depth values; the Wald test is introduced to verify physical consistency, quantify the residuals of pixels deviating from the linear dispersion curve, eliminate inversion anomalies that have no physical meaning, and ensure that the inversion process conforms to the basic laws of fluid mechanics.

[0056] Step S5: Perform multi-temporal integration on the candidate inverted water depth field, and perform dynamic weight fusion with the prior water depth field of the nautical chart to obtain the fused water depth field. Then, perform deviation correction on the fused water depth field with the measured water depth value to obtain the final inverted water depth field.

[0057] In one embodiment, reference Figure 11 As shown, Figure 11 A schematic diagram of the specific process for step S5 is shown.

[0058] Step S51: Obtain candidate inversion water depth fields for multiple time phases of the target water area, and generate multi-time phase integrated inversion water depth fields using median synthesis technology.

[0059] , in, For multi-temporal integrated inversion of water depth field, Candidate inversion depth fields corresponding to different time phases, candidate inversion depth fields It is a two-dimensional matrix containing the inverted water depth values ​​h corresponding to the coordinates of all pixels in the target area.

[0060] The median synthesis technique is used to generate a multi-temporal integrated inversion depth field, which effectively resists abnormal fluctuations caused by random white waves, instantaneous cloud shadows or extreme sea conditions.

[0061] Step S52: Introduce a dynamic weighting function to fuse the multi-temporal integrated inversion water depth field with the prior water depth field of the nautical chart to obtain the fused water depth field.

[0062] , in, To integrate the deep water field, Pre-examination of water depth field for nautical charts. This is a dynamic cosine weighting function, with a range of values ​​of... Its value is modulated by the distance from the shore or the terrain slope, and is used to balance the reliability of the model's estimated value with historical prior values. In a specific implementation, in near-shore or shallow water areas where the water depth is less than a first water depth threshold, the first water depth threshold is set to, for example, 20m. In this case, the weight... Towards Offset to preserve terrain details; in deep water areas where the water depth is greater than a second water depth threshold (which is greater than a first water depth threshold), for example, if the second water depth threshold is set to 30m, then the weights... Towards Offset is applied to ensure background stability. In regions where the water depth is greater than or equal to a first water depth threshold and less than or equal to a second water depth threshold, the weight is adjusted. It can be set to 0.5.

[0063] Step S53: Using local residual correction technology, and combining the measured water depth values ​​of known control points, the fused water depth field is checked for deviation to obtain the final inverted water depth field. The calculation formula is as follows: , in: This is the final inverted water depth field after verification; To integrate the deep water field; The local residual correction term is defined as the spatial interpolation residual field between the measured water depth and the inverted water depth at the known control point. It is a geospatial cell coordinate index used to ensure that the correction operation corresponds to each cell across the entire target water area.

[0064] The final inverted water depth field output is a spatially continuous, physically interpretable high-resolution water depth grid product that meets the standard geographic format (GeoTIFF), directly serving turbid waterway monitoring, nearshore engineering construction, and coastal zone disaster prevention and mitigation decision-making. Figure 12 The final inverted water depth field effect diagram is shown.

[0065] In this embodiment, a multi-temporal data integration strategy is adopted, and median synthesis technology is used to suppress random noise. A cosine weighting function is introduced to realize the dynamic fusion of remote sensing inversion solution and nautical chart prior field. Topographic micro-features are preserved in the nearshore high-energy area, and nautical chart data is used to ensure smooth transition and background stability in the deep water area.

[0066] The final inverted depth field was further compared and analyzed with multibeam bathymetric depth data to evaluate its effectiveness. One hundred representative spatial control points were selected within the study area to establish a one-to-one correspondence between the inverted depth values ​​and the measured depth values, such as... Figure 13 , Figure 14 As shown, where Figure 13 This shows a spatial comparison point sampling map of the final inverted water depth field data. Figure 14 This is a spatial reference point sampling map of the multibeam bathymetry data for the corresponding region. Based on this, a quantile-quantile plot is used to verify the consistency of the statistical distribution of the two sets of sampling data. The inversion accuracy is comprehensively evaluated from both the overall trend and distribution characteristics perspectives. Figure 15 As shown, Figure 15A comparison chart of the retrieved water depth values ​​and the multibeam bathymetric measured water depth values ​​is shown. All data points are closely distributed along the diagonal (1:1 line), indicating a high degree of matching between the statistical distributions of the retrieved and measured water depths, and stable overall retrieval accuracy. Quantization performance meets standards: root mean square error (RMSE) is 2.57m, and the coefficient of determination is [missing information]. The value of 0.79 indicates that the inversion results can explain 79% of the measured water depth variation, and the inversion error is controlled within a reasonable range, meeting the accuracy requirements for nearshore water depth remote sensing inversion. Furthermore, there are no significant systematic biases across the entire water depth range of 0–35 m, and the data points are more concentrated in shallow water areas, indicating higher reliability of the inversion in nearshore shallow water areas.

[0067] In summary, the water depth inversion method coupled with wave dynamics mechanism and chart constraints provided by this invention extracts pixel gray intensity from the multispectral remote sensing image and extracts wave dynamic feature field through spatial domain frequency transformation of pixel gray intensity. This breaks through the dependence of traditional remote sensing inversion on the synchronization of multi-source remote sensing data. It can extract complete wave dynamic field features using only a single image without the need for synchronous auxiliary observation data, which greatly reduces the data cost and application threshold of inversion. Furthermore, by constructing a water depth mapping model in conjunction with wave dispersion mechanism, it solves the bottleneck problem of weak optical spectral signals and limited inversion depth in turbid waters, extending the effective inversion depth to the medium and deep water area. Simultaneously, a priori depth field is extracted from the chart bathymetry data. A multi-dimensional input matrix is ​​constructed by combining the pixel grayscale intensity and the priori depth field. This priori depth field serves as the physical constraint for the entire inversion process, rather than merely as a later verification benchmark. By constraining the range of wave parameters through the topographic framework provided by the chart, and using it as the initial value for iterative solutions and a reference benchmark for the fusion stage, the model's predictive stability in mid-deep water is significantly improved, avoiding the divergence problem of pure data-driven inversion in deep water. Furthermore, multi-temporal integration of the candidate inversion depth fields is performed, combined with dynamic weight fusion of the priori depth field. This multi-temporal integration and dynamic weight fusion optimize the inversion results. Wave motion laws are used to compensate for the attenuation of optical signals in turbid water, achieving stable high-precision bathymetry. The root mean square error (RMSE) between the inverted depth and the measured depth is approximately 2.57m, ensuring stable bathymetry. It can remain within a reasonable range even under complex water conditions, and has high suspension levels. It remains well-suited for use in turbid nearshore waters.

[0068] In another embodiment of the present invention, a water depth inversion device is also provided. This embodiment is a device embodiment corresponding to the above embodiments, such as... Figure 16 As shown, Figure 16A general structural block diagram of a depth inversion device according to an embodiment of the present invention is shown. This device embodiment can be implemented in conjunction with the above embodiments. The relevant technical details mentioned in the above embodiments remain valid in this device embodiment, and will not be repeated here to avoid repetition.

[0069] A water depth inversion device 500 that couples wave dynamics mechanisms with chart constraints, the device comprising: The acquisition module 510 is used to acquire multispectral remote sensing images of the target water area and the corresponding nautical chart bathymetry data.

[0070] The dataset construction module 520 is used to extract pixel gray intensity from the multispectral remote sensing image, extract the prior water depth field of the nautical chart from the nautical chart bathymetry data, and construct a multidimensional input matrix by combining the pixel gray intensity and the prior water depth field of the nautical chart.

[0071] The wave kinematics feature extraction module 530 is used to perform spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field.

[0072] The inversion modeling module 540 is used to construct a mapping model about water depth by combining the wave dynamic characteristic field of the region, inverting the global angular frequency based on the deep water reference area, determining the initial inversion water depth field by combining the global angular frequency, and removing abnormal data from the initial inversion water depth field to obtain the candidate inversion water depth field.

[0073] The optimization module 550 is used to perform multi-temporal integration on the candidate inverted water depth field, and to perform dynamic weight fusion with the prior water depth field of the nautical chart to obtain the fused water depth field. The fused water depth field is then corrected for deviations by combining the measured water depth values ​​to obtain the final inverted water depth field.

[0074] In addition, embodiments of the present invention also provide a readable storage medium storing a program or instructions, which, when executed by a processor, implements the steps of the above-described method for water depth inversion that couples wave dynamics mechanisms with chart constraints, and achieves the same technical effect.

[0075] It should be noted that the scope of the methods in the embodiments of the present invention is not limited to performing functions in the order shown or discussed, but may also include performing functions substantially simultaneously or in the reverse order, depending on the functions involved. For example, the described methods may be performed in a different order than described, and various steps may be applied, omitted, or combined. In addition, features described with reference to certain examples may be combined in other examples.

[0076] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the present invention without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of the present invention.

Claims

1. A method for water depth inversion that couples wave dynamics mechanisms with chart constraints, characterized in that, Include: Acquire multispectral remote sensing images of the target water area and corresponding nautical chart bathymetry data; Pixel gray intensity is extracted from the multispectral remote sensing image, and the prior water depth field of the nautical chart is extracted from the nautical chart bathymetry data. A multidimensional input matrix is ​​constructed by combining the pixel gray intensity and the prior water depth field of the nautical chart. The pixel grayscale intensity is subjected to spatial domain frequency transformation to construct a regional wave dynamic feature field; the steps include: Based on Radon transform, linear integral projection is performed on each pixel of the pixel gray intensity along multiple azimuth angles to map the two-dimensional spatial ripple texture into a one-dimensional projection vector, and the one-dimensional projection vectors from different angles are integrated to construct a sine curve. A cross-spectrum matrix is ​​constructed from the main image and the secondary image of the corresponding band of the sine curve. The phase difference gradient of each frequency component is calculated. The amplitude of the cross-spectrum is weighted and reconstructed using the phase difference gradient to eliminate 180° directional ambiguity and determine the direction of the wave motion vector. Using the direction of the wave motion vector as a constraint, the pixel gray intensity is processed by a two-dimensional discrete Fourier transform within a sliding window to calculate the wave energy spectrum. The peaks of the wave energy spectrum are identified, and the number of combined waves and the direction of the main wave group are calculated by combining the frequency domain coordinates of the peaks with the spatial resolution parameters of the image. The wavelength is then calculated from the number of combined waves. By traversing the target water area through a step-sliding window, a regional wavelength distribution field and a main wave direction distribution field are generated to construct a wave dynamic characteristic field. Based on the wave dynamic characteristic field of the region, a mapping model for water depth is constructed. The global angular frequency is obtained by inversion with the deep water reference area as the reference. The initial inversion water depth field is determined by combining the global angular frequency. Abnormal data is removed from the initial inversion water depth field to obtain the candidate inversion water depth field. The candidate inverted water depth fields are integrated across multiple time phases and dynamically weighted in conjunction with the prior water depth fields from the nautical charts to obtain a fused water depth field. The fused water depth field is then corrected for deviations by incorporating measured water depth values ​​to obtain the final inverted water depth field.

2. The method according to claim 1, characterized in that, The steps of extracting pixel grayscale intensity from the multispectral remote sensing image, extracting the prior water depth field from the nautical chart bathymetry data, and constructing a multidimensional input matrix by combining the pixel grayscale intensity and the prior water depth field include: Atmospheric correction is performed on the multispectral remote sensing image to obtain the reflectance data of the lower atmosphere; The preset blue and red light bands are extracted from the reflectance data of the lower atmosphere as spectral reflectance data sources; The nautical chart bathymetry data is resampled so that the resampled nautical chart bathymetry data has the same spatial resolution as the multispectral remote sensing image. The spatial coordinates of the spectral reflectance data source and the resampled nautical chart bathymetry data are aligned. The spatially aligned spectral reflectance data source is linearly converted into the pixel grayscale intensity. Extract the prior water depth field of the target area from the spatially aligned chart bathymetry data; By integrating the pixel grayscale intensity, prior water depth field of the nautical chart, and image spatial resolution, the multidimensional input matrix is ​​constructed, wherein the multidimensional input matrix is ​​represented as: ， in, For a multidimensional input matrix, For geospatial pixel coordinate indexing, The pixel grayscale intensity is used to characterize the geometric shape of sea surface ripples; Pre-examination of water depth field for nautical charts This refers to the spatial resolution of the image.

3. The method according to claim 1, characterized in that, The steps of identifying the spectral peaks of the wave energy spectrum, calculating the combiner number and main wave direction of the main wave group by combining the frequency domain coordinates of the spectral peaks with the spatial resolution parameters of the image, and calculating the wavelength from the combiner number include: The peaks in the wave energy spectrum were identified using an extreme value search algorithm; Based on the frequency domain coordinates of the spectral peaks and the spatial resolution of the image, the physical wavenumber components in the horizontal and vertical directions of the main wave group are calculated using the following formula: ， in, The physical wavenumber component is in the horizontal direction. The physical wavenumber component is in the vertical direction. The frequency domain coordinates of the spectral peak, For image spatial resolution, This represents the total number of pixels contained in the sliding window in the horizontal direction. This represents the total number of pixels contained in the sliding window in the vertical direction. The combined wavenumber of the main wave group is calculated based on the physical wavenumber components in the horizontal and vertical directions. The calculation formula is: ， Based on the physical wavenumber components in the horizontal and vertical directions, the main wave direction of the main wave group is calculated. The calculation formula is: ， The wavelength is calculated based on the combined wave number, using the following formula: ， in, λ is the wavelength.

4. The method according to claim 3, characterized in that, The wave dynamic characteristic field is represented as follows: ， in, This is the discrete coordinate index of the sliding window center in the geospatial grid. For wave dynamic characteristic field, This represents the wavelength value at that coordinate. This represents the dominant wave direction at that coordinate. The peak confidence parameter is defined as the ratio of the peak value of the energy spectrum to the mean value of the background.

5. The method according to claim 1, characterized in that, Before constructing the regional wave dynamic feature field by performing spatial domain frequency transformation on the pixel grayscale intensity, the method further includes: Adaptive texture enhancement is performed on the pixel grayscale intensity in the multidimensional input matrix to improve the grayscale contrast between peaks and troughs; A bandpass filtering algorithm is used to remove ultra-long period background signals and high-frequency noise, while retaining effective wave signals modulated by terrain, to obtain the pixel gray intensity after denoising and enhancement.

6. The method according to claim 1, characterized in that, The steps of constructing a mapping model for water depth based on the wave dynamic characteristic field of the region, obtaining the global angular frequency by inversion using the deep-water reference area as a reference, and determining the initial inverted water depth field based on the global angular frequency include: A mapping model between angular frequency, resultant wave number, and water depth is established based on linear wave theory, and is expressed as follows: ， in, For global angular frequency, It is the acceleration due to gravity. The combination number of the main wave group. This represents the water depth. Identify deep-water reference regions that meet the deep-water conditions in the multi-dimensional input matrix, and extract the deep-water wavenumber of the deep-water reference regions. Based on the deep-water wave number, and under deep-water conditions Based on the constraint characteristics, the global angular frequency is obtained by inversely solving the mapping model, and the calculation formula is as follows: ， in, For deep-water wave number; The Newton-Raphson iterative algorithm is used to solve the transcendental equation of the linear wave dispersion relation. After the iteration converges, the initial inverted water depth field is obtained. The calculation formula is as follows: ， Wherein, objective function derivative term , Let be a hyperbolic secant function, and the iteration termination condition is: Less than the preset threshold; and The first Next and first The water depth value in the first iteration.

7. The method according to claim 6, characterized in that, The step of removing outlier data from the initial inverted depth field to obtain candidate inverted depth fields includes: For each pixel in the initial inverted depth field, the significance of the residuals is quantified using the Wald statistical test, calculated as follows: ， in: The Wald statistic is used to assess the degree to which the inversion results deviate from the physical mechanism or prior constraints. This value follows a set of 1 degree of freedom. distributed; To retrieve the water depth value; This serves as a reference value for water depth. The variance of the inverted water depth value; The inverted water depth values ​​of the pixels that pass the physical consistency check are retained. For the removed abnormal pixels, spatial interpolation is performed using the inverted water depth values ​​of the neighboring valid pixels to fill the space, thus obtaining the candidate inverted water depth field.

8. The method according to claim 1, characterized in that, The steps of performing multi-temporal integration on the candidate inverted depth fields, dynamically weighting and fusing them with the prior nautical chart depth fields to obtain a fused depth field, and then correcting the bias of the fused depth field with measured depth values ​​to obtain the final inverted depth field include: The candidate inversion water depth fields of the target water area in multiple time phases are obtained, and the multi-time phase integrated inversion water depth field is generated by median synthesis technique; By introducing a dynamic weighting function, the multi-temporal integrated inversion water depth field and the prior water depth field of the nautical chart are fused to obtain the fused water depth field; Using local residual correction techniques, and combining the measured water depth values ​​of known control points, the fused water depth field is checked for deviations to obtain the final inverted water depth field. The calculation formula is as follows: ， in: This is the final inverted water depth field after verification; To integrate the deep water field; The local residual correction term is defined as the spatial interpolation residual field between the measured water depth and the inverted water depth at the known control point. It is a geospatial cell coordinate index used to ensure that the correction operation corresponds to each cell across the entire target water area.

9. A water depth inversion device that couples wave dynamics mechanism with chart constraints, characterized in that, The apparatus employing the water depth inversion method based on the coupled wave dynamics mechanism and chart constraints as described in any one of claims 1-8 comprises: The acquisition module is used to acquire multispectral remote sensing images of the target water area and the corresponding nautical chart bathymetry data; The dataset construction module is used to extract pixel gray intensity from the multispectral remote sensing image, extract the prior water depth field of the nautical chart from the nautical chart bathymetry data, and construct a multidimensional input matrix by combining the pixel gray intensity and the prior water depth field of the nautical chart. The wave kinematics feature extraction module is used to perform spatial domain frequency transformation on the pixel grayscale intensity to construct a regional wave dynamic feature field. The inversion modeling module is used to construct a mapping model about water depth by combining the wave dynamic characteristic field of the region, inverting the global angular frequency based on the deep water reference area, determining the initial inversion water depth field by combining the global angular frequency, and removing abnormal data from the initial inversion water depth field to obtain the candidate inversion water depth field. An optimization module is used to perform multi-temporal integration of the candidate inverted water depth field, and to perform dynamic weighted fusion with the prior water depth field of the nautical chart to obtain a fused water depth field. The fused water depth field is then corrected for deviations by combining the measured water depth values ​​to obtain the final inverted water depth field.

Images