CN114594503A - Shallow sea terrain inversion method, computer equipment and storage medium - Google Patents
Shallow sea terrain inversion method, computer equipment and storage medium Download PDFInfo
- Publication number
- CN114594503A CN114594503A CN202210202515.4A CN202210202515A CN114594503A CN 114594503 A CN114594503 A CN 114594503A CN 202210202515 A CN202210202515 A CN 202210202515A CN 114594503 A CN114594503 A CN 114594503A
- Authority
- CN
- China
- Prior art keywords
- photon
- sea
- photons
- data
- height
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 59
- 238000001914 filtration Methods 0.000 claims abstract description 18
- 238000012876 topography Methods 0.000 claims abstract description 12
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 230000006870 function Effects 0.000 claims description 55
- 230000002902 bimodal effect Effects 0.000 claims description 22
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 21
- 238000004590 computer program Methods 0.000 claims description 20
- 238000012937 correction Methods 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 7
- 230000011218 segmentation Effects 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 6
- 238000004140 cleaning Methods 0.000 claims description 4
- 239000013535 sea water Substances 0.000 claims description 3
- 238000005259 measurement Methods 0.000 abstract description 8
- 238000012887 quadratic function Methods 0.000 abstract description 6
- 230000003287 optical effect Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 238000012216 screening Methods 0.000 description 4
- 238000000605 extraction Methods 0.000 description 3
- 238000013507 mapping Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000005406 washing Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
The invention discloses a shallow sea topography inversion method, computer equipment and a storage medium, wherein a sliding window algorithm is adopted to segment zonal photon data of a height measurement satellite, the segmented photon data is finely processed, noise is removed, high-precision sea surface photons and seabed photons are obtained, then a relation is established between the depth of the seabed photons and the blue-green waveband logarithmic ratio of hyperspectral image pixels in a corresponding time period by a quadratic function model, inversion parameters are calculated, and finally the whole hyperspectral image is inverted by the inversion parameters to obtain the high-precision shallow sea topography. The invention adopts a mode of 'double-peak Gaussian fitting + mean filtering' to stably and accurately distinguish photons above the sea surface, photons below the sea surface and photons above the sea surface, and then adopts a mode of 'single-peak Gaussian fitting + median filtering' to extract the submarine photons, thereby more completely extracting the submarine photons, ensuring the submarine photon data quantity and further ensuring the inversion precision.
Description
Technical Field
The invention belongs to the field of ocean remote sensing and shallow sea terrain, and particularly relates to a shallow sea terrain inversion method based on altimetry satellite data (ICESat-2) and a hyperspectral remote sensing image (Sentinel-2), computer equipment and a storage medium.
Background
The submarine topography, particularly the shallow sea topography is one of the main targets of ocean remote sensing mapping, and a water depth topographic map drawn according to the shallow sea topography can be applied to the fields of navigation, shipping, fishery fishing, ocean environment management, ocean engineering, ocean resource development and the like.
In recent years, with the vigorous development of satellite height measurement technology, the resolution and the accuracy of the satellite height measurement technology are greatly improved, the precise inversion can be performed on the ground level surface, the gravity anomaly and the shallow sea terrain by using the satellite height measurement data, and the height measurement satellite data and the optical image become one of the main research methods for the shallow sea terrain inversion.
For example, Yue Ma and others propose "Satellite-driven bathymetry using the ocean at-2lidar and Sentinel-2 imaging data" (published: isps Journal of photonic and Remote Sensing, 2020), in which a method for extracting sea-surface sea-floor photons by filtering noise by using a DBSCAN algorithm is disclosed, the method processes the photon data of day and night by adjusting a threshold value, but in the original photon data, a considerable part of the photon density of sea floor is greatly different, so that a part of the sea-floor photon data with a larger difference between the density and the set threshold value is lost, the sea-floor photon data is insufficient, and the inverted parameters are inaccurate.
For example, Hsiao-Journal Hsu et al propose "A semi-empirical scheme for basal mapping in shallow water by ICESat-2and Sentinel-2: A case study in the South China Sea" (published: ISPRS Journal of photon mapping and Remote Sensing, 2021), wherein a method for inverting shallow Sea terrain by establishing a relationship between the ratio of the blue-green band pair (R value) of the optical image and the depth measurement of the laser is disclosed, which employs a method for taking the average depth of the corresponding photons when dealing with the condition that one pixel corresponds to multiple photons, and which loses some high-precision photon information during the processing, so that the number of actually participating photons is reduced, and the inverted parameters are not accurate enough to some extent.
Because photon data of a laser altimetry satellite (ICESat-2) generally has more noise, the data with higher confidence coefficient can be obtained only by manual denoising treatment, and compared with an optical image, the data quantity of the altimetry data is less, taking photon data of the ICESat-2 altimetry satellite as an example (ICESat-2 is the highest satellite in the current laser altimetry satellite), the ICESat-2 satellite measures the distance from the ICESat-2 satellite to the earth surface along the satellite orbit direction by using 3 pairs of green lasers (each pair of lasers consists of two strong and weak lasers which are 90 meters away from each other and 3.3 kilometers away from each pair of lasers), and calculates the earth surface height of each photon, and the data are relatively sparse. The resolution of a hyperspectral image (such as a Sentinel-2 image) can reach 10m, the contained information is wider, and the key of how to obtain high-precision height measurement satellite data and how to fully utilize high-precision photon data is shallow sea topography inversion.
Disclosure of Invention
The invention aims to provide a shallow sea terrain inversion method, computer equipment and a storage medium, so as to solve the problem that parameters obtained by inversion in the prior art are inaccurate. According to the invention, through sliding window segmentation, sea surface photons and photons below the sea surface are extracted by adopting a mode of double-peak Gaussian fitting and mean filtering, and then seabed photons in the photon data below the sea surface are extracted by adopting a mode of single-peak Gaussian fitting and median filtering, so that the seabed photons can be extracted more completely even if the density of the seabed photons is different, and the seabed photon data quantity is ensured; a corresponding R value is found for each photon by a bilinear interpolation method, and all high-precision photon data are adopted to participate in water depth inversion, so that the high-precision photon data are not lost, and the accuracy of inversion parameters is ensured.
The invention solves the technical problems through the following technical scheme: a shallow sea terrain inversion method comprises the following steps:
acquiring a hyperspectral image, and calculating the R value of each pixel in the hyperspectral image;
acquiring photon data, cleaning the photon data, and removing photon data, non-sea area data and noise data which are not in a research area;
carrying out sectional processing on the cleaned photon data to obtain a plurality of sections of photon data;
calculating a histogram of each section of photon data, wherein each section of photon data corresponds to one histogram;
fitting each histogram by adopting a bimodal Gaussian function to obtain an optimal bimodal Gaussian function corresponding to each histogram;
filtering the photon data above the sea surface from the corresponding histogram by adopting a Gaussian mean filtering method according to the parameters of the optimal double-peak Gaussian function, and extracting the photon data above the sea surface and the photon data below the sea surface;
fitting the photon data below the sea surface by adopting a unimodal Gaussian function to obtain an optimal unimodal Gaussian function;
extracting submarine photon data from photon data below the sea surface by adopting a median filtering method according to parameters of the optimal unimodal Gaussian function;
performing refraction correction on the submarine photon data to obtain submarine photon data after refraction correction;
aligning the refraction-corrected submarine photon data with the R value of each pixel in the hyperspectral image according to the longitude and latitude;
calculating the R value of each submarine photon according to the R values of the adjacent pixels of each submarine photon;
and fitting the R values and the water depths of all the submarine photons to determine the coefficients of a fitting function, wherein the coefficients of the fitting function are inversion parameters.
Preferably, the washed photon data is processed in a segmented manner by adopting a sliding window algorithm to obtain a plurality of segments of photon data.
Further, the specific implementation process of calculating the histogram of each segment of photon data includes:
calculating the height of the vertical subsection according to the height range and the number of photons in each section of photon data;
calculating the number of photons in each vertical subsection height in the section of photon data in the height direction of the section of photon data;
and forming a histogram of the photon data according to the vertical segmentation heights and the number of photons in each vertical segmentation height.
Further, the calculation expression of the vertical segment height is as follows:
wherein h isiIs the vertical segment height, H, of the i-th segment of photon dataiIs the height range of all photons in the ith segment of photon data, niIs the number of photons in the ith segment of photon data.
Further, the expression of the bimodal gaussian function is:
wherein, f (x)1Representing a bimodal Gaussian function, a1、a2All represent the amplitude u1、u2All indicate expectation, u1＞u2，σ1、σ2All represent standard deviations; when x represents the photon height, the optimal bimodal Gaussian function f (x)1Amplitude a of1Approximately as the maximum number of sea surface photons, amplitude a, falling within the vertical segment height in the corresponding segment photon data2Approximately the maximum number of seafloor photons, expected u, of the corresponding segment photon data that fall within the vertical segment height1For approximate sea surface height in the corresponding segment photon data, u is expected2To approximate sea floor height, σ, in the corresponding segment of photon data1Standard deviation, σ, of the peak of a photon at sea2The standard deviation of the ocean bottom photon peak.
Preferably, when the photon height x > u1+3σ1Photons above the sea surface;
when u is1-3σ1The photon height x is less than or equal to u1+3σ1The time is sea surface photons;
when the photon height x < u1-3σ1Photons below the sea surface; wherein u is1Expectation of sea surface photon peak, σ, as an optimal bimodal gaussian function1The standard deviation of the sea surface photon peak of the optimal bimodal gaussian function.
Further, the expression of the unimodal gaussian function is:
wherein, f (x)2Representing a unimodal Gaussian function, a3Represents the amplitude u3Indicates expectation, σ3Represents the standard deviation; when x represents the photon height, the optimal unimodal Gaussian function f (x)2Amplitude a of3Approximately the maximum number of seafloor photons, expected u, of the corresponding segment photon data that fall within the vertical segment height3To approximate sea floor height, σ, in the corresponding segment of photon data3The standard deviation of the ocean bottom photon peak.
Preferably, when mj-3σ3The height x of the photons is less than or equal to mj+3σ3Time is a photon of the sea bottom, where3Standard deviation of the optimal single and double peak Gaussian function, mjIs the median value in the j-th median filtering.
Further, a specific expression of the refraction correction is as follows:
wherein z is the depth of the refraction-corrected seabed photon water, z0Depth of ocean bottom photon water before refraction correction, naIs the refractive index of the laser beam in air, nwIs the refractive index of the laser beam in seawater.
Further, calculating the R value of each submarine photon by using a bilinear interpolation method, wherein the specific calculation expression is as follows:
wherein (x)p,yp) Is the latitude and longitude (x) of the photon P at the sea bottom1,y1) The latitude and longitude of the adjacent pixel at the lower left corner of the photon P at the sea floor, (x)1,y2) Latitude and longitude of the adjacent pixel at the upper left corner of the photon P at sea floor, (x)2,y1) Latitude and longitude of the adjacent pixel at the bottom right corner of the photon P at sea floor, (x)2,y2) Latitude and longitude R of the adjacent pixel at the upper right corner of the photon P on the sea bottom11Is a pixel (x)1,y1) R value of (A), R21Is a pixel (x)2,y1) R value of (A), R12Is a pixel (x)1,y2) R value of (A), R22Is a pixel (x)2,y2) R value of (A), RPThe R value of the ocean bottom photon P.
Further, the method further comprises: and calculating the water depth of each pixel according to the inversion parameters and the R value of each pixel in the hyperspectral image, and acquiring the shallow-sea topography according to the water depth of each pixel and the sea surface elevation.
The present invention also provides a computer apparatus comprising: a memory for storing a computer program; a processor for implementing the steps of the shallow sea terrain inversion method as described above when executing the computer program.
The invention also provides a storage medium having stored thereon a computer program which, when executed by a processor, carries out the steps of the shallow sea terrain inversion method as described above.
Advantageous effects
Compared with the prior art, the invention has the advantages that:
the shallow sea topography inversion method, the computer equipment and the storage medium provided by the invention are characterized in that a sliding window algorithm is adopted to segment zonal photon data of a height measurement satellite, the segmented photon data is finely processed to remove noise to obtain high-precision sea surface photons and seabed photons, then a relation is established between the depth of the seabed photons and the blue-green waveband logarithmic ratio (namely R value) of hyperspectral image pixels in a corresponding time period by a quadratic function model, inversion parameters are calculated, and finally the whole hyperspectral image is inverted by the inversion parameters to obtain the high-precision shallow sea topography;
the invention adopts a mode of 'double-peak Gauss fitting + mean filtering' to steadily and accurately distinguish photons (noise data) above the sea surface, photons on the sea surface and photons (including submarine photons and noise data) below the sea surface, and then adopts a mode of 'single-peak Gauss fitting + median filtering' to extract the submarine photons, so that the submarine photons can be more completely extracted even if the density of the submarine photons is different, and the submarine photon data volume is ensured; aiming at the condition that one pixel corresponds to a plurality of photon data, a bilinear interpolation method is adopted to replace an averaging method to obtain all photon information, so that high-precision information of all photon on the sea bottom is fully utilized to invert high-precision shallow-sea topography.
The method aims at shallow sea terrain inversion, and has higher inversion accuracy compared with other existing inversion methods.
Drawings
In order to more clearly illustrate the technical solution of the present invention, the drawings needed to be used in the description of the embodiments are briefly introduced below, and it is obvious that the drawings in the following description are only one embodiment of the present invention, and it is obvious for those skilled in the art to obtain other drawings based on the drawings without creative efforts.
FIG. 1 is a flow chart of a shallow sea terrain inversion method in an embodiment of the present invention;
FIG. 2 is a hyperspectral image corrected by atmosphere after cloud cover screening in an embodiment of the invention;
FIG. 3 is a schematic view of a sliding window segment and a vertical segment in an embodiment of the invention;
FIG. 4 is a schematic diagram of a fitting of a bimodal Gaussian function according to an embodiment of the present invention;
FIG. 5 is a diagram illustrating bilinear interpolation in an embodiment of the present invention;
FIG. 6 is a schematic diagram of fitting inversion parameters by different methods in an embodiment of the present invention, where (a) is inversion by a mean method, the number of photons involved in inversion is 769, and the fitting degree is 0.85; the figure (b) adopts bilinear interpolation method to carry out inversion, the number of photons participating in the inversion is 5771, and the fitting degree is 0.86;
FIG. 7 is a shallow sea depth map obtained using inversion parameters in an embodiment of the present invention.
Detailed Description
The technical solutions in the present invention are clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be obtained by a person skilled in the art without inventive efforts based on the embodiments of the present invention, shall fall within the scope of protection of the present invention.
The technical solution of the present application will be described in detail below with specific examples. The following several specific embodiments may be combined with each other, and details of the same or similar concepts or processes may not be repeated in some embodiments.
As shown in fig. 1, the shallow sea terrain inversion method provided in this embodiment includes the following steps:
1. hyperspectral image processing
(1.1) hyperspectral image acquisition: the hyperspectral image (image of Sentinel-2) can be obtained from published image data (https:// scihub. copernius. eu/dhus/#/home).
(1.2) hyperspectral image screening: the shallow sea terrain inversion needs the reflectivity of pixels in a shallow sea area of an optical image, in order to ensure that the shallow sea area in the image is not shielded by a cloud layer and the image is clear, the cloud amount is less than 10%, particularly the shallow sea area is cloud-free or cloud-less and the image with sufficient illumination is selected, and the image at early morning or late night is avoided as much as possible because the sun altitude is smaller, sunlight obliquely irradiates the sea surface and the energy reflected to a satellite sensor is less.
The hyperspectral image screening ensures that the shallow sea area of the selected image is not shielded by a cloud layer and is sufficiently illuminated, so that the shallow sea area can be completely imaged, and the inversion calculation can be conveniently carried out by utilizing the wave band reflection value of the pixel in the hyperspectral image shallow sea area.
(1.3) atmospheric correction: because the published Sentinel-2 data is not corrected by the atmosphere, the atmosphere correction needs to be performed by using official Sentinel-2L2A product generation and formatting processing software Sen2Cor (http:// step. esa. int/main/snap-supported-plugs/Sen 2 Cor) to eliminate the influence of the atmosphere on the transmission process of the remote sensing signal, and fig. 2 shows a hyperspectral image corrected by the atmosphere after cloud amount screening.
(1.4) calculating the R value of each pixel in the hyperspectral image: and calculating the R value of each pixel according to the blue-green wave band value of each pixel in the hyperspectral image.
2. Optical sub-data processing
(2.1) photon data reading: photon Data is extracted from a photon file downloaded from NSIDC (National Snow & Ice Data Center, National Snow Data Center) (https:// NSIDC. org/Data/atl03), and the photon Data comprises key information required by photons, such as longitude and latitude, height and distance along a track, and is grouped into a multi-dimensional list.
(2.2) photon data washing: and (4) performing photon data cleaning according to the longitude and latitude and the height, and removing photon data, non-sea area data with too high height and too low height and noise data which are not in the research area. Photon data which are not in the research area are useless in the inversion process and are to be rejected, while noise data with too high or too low photon height influence the photon extraction process and are also to be rejected, and the photon extraction efficiency and the photon extraction precision are improved by data cleaning.
(2.3) photon data segmentation: and performing segmented processing on the cleaned photon data by adopting a sliding window algorithm to obtain a plurality of segments of photon data.
As shown in fig. 3, the strip-shaped photon data is segmented in the horizontal direction by using the sliding window to obtain a plurality of segments of photon data, and each segment of photon data corresponds to a sliding window segment. The segmented photon data can be processed more finely, so that noise removal is facilitated, and high-precision sea surface photons and seabed photons can be obtained.
(2.4) calculating a histogram of each segment of photon data, wherein the specific implementation process comprises the following steps:
(2.41) calculating the height of the vertical segment according to the height range and the number of photons of all photons in each segment of photon data, wherein the specific calculation formula is as follows:
wherein h isiIs the vertical segment height, H, of the i-th segment of photon dataiIs the height range of all photons in the ith segment of photon data, niIs the number of photons in the ith segment of photon data.
(2.42) calculating the number of photons in each vertical segment height in the segment of photon data in the height direction of the segment of photon data. As shown in fig. 3, for a certain segment of sliding window segment (i.e. a certain segment of photon data), vertical segmentation is performed in the height direction according to the height of the vertical segment, and the number of photons in each vertical segment is calculated.
(2.43) forming a histogram of the segment of photon data according to the vertical segment heights and the number of photons within each vertical segment height, i.e. forming a histogram of the segment of photon data according to the number of photons within each vertical segment in the segment of photon data, as shown in fig. 4.
(2.5) for each histogram, fitting with a bimodal Gaussian function: performing least square fitting on each histogram by adopting a bimodal Gaussian function, selecting the bimodal Gaussian function with the minimum residual square sum and parameters thereof as the optimal bimodal Gaussian function, namely the optimal fitting function, of the corresponding histogram, wherein the expression of the bimodal Gaussian function is as follows:
wherein the content of the first and second substances,f(x)1representing a bimodal Gaussian function, a1、a2All represent the amplitude u1、u2All represent expectation, σ1、σ2All represent standard deviations; when x represents the photon height, the optimal bimodal Gaussian function f (x)1Amplitude a of1Approximately as the maximum number of sea surface photons, amplitude a, falling within the vertical segment height in the corresponding segment photon data2Approximately the maximum number of submarine photons, expected u, falling within the height of the vertical segment in the corresponding segment photon data1For approximate sea surface height in the corresponding segment photon data, u is expected2To approximate sea floor height, σ, in the corresponding segment of photon data1Standard deviation of sea surface photon peak (related to the vertical range of sea surface photons), σ2Is the standard deviation of the peak of the ocean bottom photon (related to the vertical extent of the ocean bottom photon). As shown in fig. 4, a peak with a large expected value is a sea-surface photon peak, such as the right peak in fig. 4, and a peak with a small expected value is a sea-bottom photon peak, such as the left peak in fig. 4.
(2.6) adopting a Gaussian mean filtering method to extract photons on the sea surface and photons below the sea surface
When the photon height x > u1+3σ1Photons above the sea surface; when u is1-3σ1The photon height x is less than or equal to u1+3σ1The time is sea surface photons; when the photon height x < u1-3σ1Photons below the sea surface are present.
After the sea surface photons are extracted, the remaining photons below the sea surface comprise seabed photons and noise, and then the unimodal Gaussian function is adopted to extract accurate seabed photons.
(2.7) fitting with a unimodal Gaussian function: the method comprises the following steps of performing least square fitting on the sub-sea surface photon data by adopting a unimodal Gaussian function, selecting the unimodal Gaussian function with the minimum residual square sum and parameters thereof as an optimal unimodal Gaussian function, namely an optimal fitting function, wherein the expression of the unimodal Gaussian function is as follows:
wherein, f (x)2Representing a single-peak Gaussian function, a3Represents the amplitude u3Indicates expectation, σ3Represents the standard deviation; when x represents the photon height, the optimal unimodal Gaussian function f (x)2Amplitude a of3Approximately the maximum number of seafloor photons, expected u, of the corresponding segment photon data that fall within the vertical segment height3To approximate sea floor height, σ, in the corresponding segment of photon data3Is the standard deviation of the peak of the ocean bottom photon (related to the vertical extent of the ocean bottom photon). The seafloor typically contains more noise and the fitting is more accurate using unimodal gaussians.
(2.8) extracting the submarine photons by adopting a cubic median filtering method
When m isj-3σ3The height x of the photons is less than or equal to mj+3σ3Time is a photon of the sea bottom, where mjIn this embodiment, j is 3, and a good denoising effect can be achieved only by a small number of times of median filtering, so as to obtain submarine photons with high precision and high reliability. .
(2.9) refraction correction of the submarine photons, wherein the specific calculation formula is as follows:
wherein z is the refraction corrected depth of the photon at the sea bottom, z0Depth of ocean bottom photon water before refraction correction, naIs the refractive index of the laser beam in air, nwIs the refractive index of the laser beam in seawater. The photon data comprises the longitude and latitude, the height, the distance along the track and the like of photons, and the depth of the submarine photon is equal to the height of the photons on the sea surface minus the height of the photons on the sea bottom.
3. Data alignment: aligning the refraction-corrected submarine photon data with the R value of each pixel in the hyperspectral image according to the longitude and latitude.
4. Calculating the R value of the submarine photons: and calculating the R value of each submarine photon by a bilinear interpolation method according to the R values of the adjacent pixels of each submarine photon.
The R value of each pixel in the hyperspectral image is a grid structure, the submarine photon is a three-dimensional point, when the data are aligned, the longitude and latitude of the submarine photon is not completely consistent with the longitude and latitude of the grid center (namely the pixel center of the hyperspectral image), so that the submarine photon is scattered near the grid center after falling into each grid structure.
As shown in FIG. 5, let the position of a certain photon P in the latitude and longitude coordinate system be (x)p,yp) The four grid centers or pixel centers adjacent to the ocean bottom photon P are respectively (x)1,y1)、(x1,y2)、(x2,y1)、(x2,y2) R values corresponding to the four grid centers are R respectively11、R12、R21、R22And carrying out bilinear interpolation on the position of the submarine photon P to obtain an R value corresponding to the position of the submarine photon P.
(4.1) first, point M (x) is pointed in the x directionp,y1) Point N (x)p,y2) Linear interpolation is performed (M, N two points are respectively a straight line x ═ xpWith line y ═ y1、y＝y2Intersection of):
wherein R isMR is the value of R corresponding to point MNThe value of R corresponding to point N.
(4.2) the point P (x) is pointed again in the y directionp,yp) Linear interpolation is performed, and the formula is:
wherein R isPThe R value of the ocean bottom photon P (i.e., the abscissa R' in fig. 6).
FIG. 6 is a schematic diagram of fitting inversion parameters by different methods, wherein (a) is an inversion performed by a mean method, the number of photons participating in the inversion is 769, and the fitting degree is 0.85; and (b) the inversion is carried out by adopting a bilinear interpolation method, the number of photons participating in the inversion is 5771, and the fitting degree is 0.86. It can be seen from fig. 6 that the number of photons participating in the inversion in the graph (b) is larger, the parameters of the inversion are more accurate, and the reliability is higher.
5. Parameter inversion: and performing quadratic function fitting on the R values and the water depths of all the submarine photons to determine a primary term coefficient, a secondary term coefficient and a constant term of the fitting function, wherein the primary term coefficient, the secondary term coefficient and the constant term are inversion parameters.
The quadratic function model is:
zP＝a(RP)2+bRP+c (8)
wherein z isPFor the water depth of the photon P at the sea bottom (i.e. the ordinate Z in fig. 6), a, b, and c are the quadratic term coefficient, the first-order term coefficient, and the constant term of the quadratic function model, respectively, i.e. a, b, and c are the inversion parameters to be solved. Illustratively, as in fig. 6(a), a is 6386.188, b is 12553.308, c is 6169.736; in FIG. 6(b), a is 6742.714, b is 13281.673, and c is 6541.305.
6. Shallow sea elevation model: and calculating the water depth of each pixel according to the inversion parameters and the R value of each pixel in the hyperspectral image, acquiring the shallow sea topography according to the water depth and the sea surface elevation of each pixel, and acquiring the shallow sea depth map according to the water depth of each pixel.
The inversion parameters can be obtained in step 5, the coefficients of the quadratic functions in the formula (8) are known, the water depth of each pixel can be calculated according to the known R value of each pixel, and the digital elevation model and the shallow sea depth map of the whole research area can be obtained by subtracting the sea surface elevation from the water depth of the pixel, as shown in fig. 7.
Embodiments of the present invention further provide a computer device, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor implements the steps of the shallow sea terrain inversion method as described above when executing the computer program.
Illustratively, the computer program may be partitioned into one or more modules/units that are stored in the memory and executed by the processor to implement the invention. The one or more modules/units may be a series of computer program instruction segments capable of performing specific functions, which are used to describe the execution of the computer program in the computer device.
The device may be a desktop computer, a notebook, a palm top computer, a cloud server, or other computing device. The apparatus may include, but is not limited to, a processor, a memory. Those skilled in the art will appreciate that the computer apparatus of the present invention is merely exemplary of an apparatus and is not intended to be limiting and may include more or less components than the apparatus, or some of the components may be combined, or different components, such as input output devices, network access devices, buses, etc.
The Processor may be a Central Processing Unit (CPU), other general purpose Processor, a Digital Signal Processor (DSP), an Application Specific Integrated Circuit (ASIC), an off-the-shelf Programmable Gate Array (FPGA) or other Programmable logic device, discrete Gate or transistor logic device, discrete hardware component, etc. A general purpose processor may be a microprocessor or the processor may be any conventional processor or the like.
The memory may be used to store the computer program and/or module, and the processor may implement each step of the shallow sea terrain inversion method by running or executing the computer program and/or module stored in the memory, and invoking data stored in the memory. The memory may mainly include a storage program area and a storage data area, wherein the storage program area may store an operating system, an application program required by at least one function (such as a sound playing function, an image playing function, etc.), and the like; the storage data area may store data (such as audio data, a phonebook, etc.) created according to the use of the cellular phone, and the like. In addition, the memory may include high speed random access memory, and may also include non-volatile memory, such as a hard disk, a memory, a plug-in hard disk, a Smart Media Card (SMC), a Secure Digital (SD) Card, a Flash memory Card (Flash Card), at least one magnetic disk storage device, a Flash memory device, or other volatile solid state storage device.
The computer program, when executed by a processor, implements the steps of the shallow sea terrain inversion method.
The computer device integrated modules/units, if implemented in the form of software functional units and sold or used as separate products, may be stored in a computer readable storage medium. Based on such understanding, all or part of the flow of the method according to the embodiments of the present invention may also be implemented by a computer program, which may be stored in a computer-readable storage medium, and when the computer program is executed by a processor, the steps of the method embodiments described above may be implemented. Wherein the computer program comprises computer program code, which may be in the form of source code, object code, an executable file or some intermediate form, etc. The computer-readable medium may include: any entity or device capable of carrying the computer program code, recording medium, usb disk, removable hard disk, magnetic disk, optical disk, computer Memory, Read-Only Memory (ROM), Random Access Memory (RAM), electrical carrier wave signals, telecommunications signals, software distribution medium, and the like. It should be noted that the computer readable medium may contain content that is subject to appropriate increase or decrease as required by legislation and patent practice in jurisdictions, for example, in some jurisdictions, computer readable media does not include electrical carrier signals and telecommunications signals as is required by legislation and patent practice.
The above disclosure is only for the specific embodiments of the present invention, but the scope of the present invention is not limited thereto, and any changes or modifications within the technical scope of the present disclosure may be easily conceived by those skilled in the art and shall be covered by the scope of the present invention.
Claims (10)
1. A shallow sea terrain inversion method is characterized by comprising the following steps:
acquiring a hyperspectral image, and calculating the R value of each pixel in the hyperspectral image;
acquiring photon data, cleaning the photon data, and removing the photon data, non-sea area data and noise data which are not in a research area;
carrying out sectional processing on the cleaned photon data to obtain a plurality of sections of photon data;
calculating a histogram of each section of photon data, wherein each section of photon data corresponds to one histogram;
fitting each histogram by adopting a bimodal Gaussian function to obtain an optimal bimodal Gaussian function corresponding to each histogram;
filtering the photon data above the sea surface from the corresponding histogram by adopting a Gaussian mean filtering method according to the parameters of the optimal double-peak Gaussian function, and extracting the photon data above the sea surface and the photon data below the sea surface;
fitting the photon data below the sea surface by adopting a unimodal Gaussian function to obtain an optimal unimodal Gaussian function;
extracting submarine photon data from photon data below the sea surface by adopting a median filtering method according to parameters of the optimal unimodal Gaussian function;
performing refraction correction on the submarine photon data to obtain submarine photon data after refraction correction;
aligning the refraction-corrected submarine photon data with the R value of each pixel in the hyperspectral image according to the longitude and latitude;
calculating the R value of each ocean bottom photon according to the R values of the adjacent pixels of each ocean bottom photon;
fitting the R values and the water depths of all the submarine photons to determine the coefficients of a fitting function, wherein the coefficients of the fitting function are inversion parameters;
preferably, the washed photon data is processed in a segmented manner by a sliding window algorithm to obtain a plurality of segments of photon data.
2. The shallow sea terrain inversion method of claim 1, wherein the specific implementation process for calculating the histogram of each segment of photon data is as follows:
calculating the height of the vertical subsection according to the height range and the number of photons in each section of photon data;
calculating the number of photons in each vertical subsection height in the section of photon data in the height direction of the section of photon data;
and forming a histogram of the photon data according to the vertical segmentation heights and the number of photons in each vertical segmentation height.
3. The shallow sea terrain inversion method of claim 2, wherein the vertical section height is calculated by the expression:
wherein h isiIs the vertical segment height, H, of the i-th segment of photon dataiIs the height range of all photons in the ith segment of photon data, niIs the number of photons in the ith segment of photon data.
4. The shallow sea terrain inversion method of claim 1, wherein the expression of the bimodal gaussian function is:
wherein, f (x)1Representing a bimodal Gaussian function, a1、a2All represent the amplitude u1、u2All indicate expectation, u1＞u2，σ1、σ2All represent standard deviations; when x represents the photon height, the optimal bimodal Gaussian function f (x)1Amplitude a of1Approximately as the maximum number of sea surface photons, amplitude a, falling within the vertical segment height in the corresponding segment photon data2Approximately the maximum number of seafloor photons, expected u, of the corresponding segment photon data that fall within the vertical segment height1For approximate sea surface height in the corresponding segment photon data, u is expected2To approximate sea floor height, σ, in the corresponding segment of photon data1Standard deviation, σ, of the peak of a photon at sea2Standard deviation of the sea-bottom photon peak;
preferably, when the photon height x > u1+3σ1Photons above the sea surface;
when u is1-3σ1The photon height x is less than or equal to u1+3σ1The time is sea surface photons;
when the photon height x < u1-3σ1Photons below the sea surface; wherein u is1Expectation of sea surface photon peak, σ, as an optimal bimodal gaussian function1The standard deviation of the sea surface photon peak of the optimal bimodal gaussian function.
5. The shallow sea terrain inversion method of claim 1, wherein the unimodal gaussian function has the expression:
wherein, f (x)2Representing a single-peak Gaussian function, a3Represents the amplitude u3Indicates expectation, σ3Represents the standard deviation; when x represents the photon height, the optimal unimodal Gaussian function f (x)2Amplitude a of3Is approximated as corresponding segment lightMaximum number of seafloor photons, expected u, of the subdata falling within the vertical segment height3To approximate sea floor height, σ, in the corresponding segment of photon data3Standard deviation of the sea-bottom photon peak;
preferably, when mj-3σ3The height x of the photons is less than or equal to mj+3σ3Time is a photon of the sea bottom, of which σ3Standard deviation of the optimal single and double peak Gaussian function, mjIs the median value in the j-th median filtering.
6. The shallow sea terrain inversion method of claim 1, wherein the refraction correction is specifically expressed as:
wherein z is the depth of the refraction-corrected seabed photon water, z0Depth of ocean bottom photon water before refraction correction, naIs the refractive index of the laser beam in air, nwIs the refractive index of the laser beam in seawater.
7. The shallow sea terrain inversion method of any one of claims 1-6, characterized in that a bilinear interpolation method is used to calculate the R value of each seafloor photon, and the specific calculation expression is:
wherein (x)p,yp) Is the latitude and longitude (x) of the photon P at the sea bottom1,y1) The latitude and longitude of the adjacent pixel at the lower left corner of the photon P at the sea floor, (x)1,y2) Latitude and longitude of the adjacent pixel at the upper left corner of the photon P at sea floor, (x)2,y1) Latitude and longitude of the adjacent pixel at the bottom right corner of the photon P at sea floor, (x)2,y2) Latitude and longitude R of the adjacent pixel at the upper right corner of the photon P on the sea bottom11Is a pixel (x)1,y1) R value of (A), R21Is a pixel (x)2,y1) R value of (A), R12Is a pixel (x)1,y2) R value of (A), R22Is a pixel (x)2,y2) R value of (A), RPThe R value of the ocean bottom photon P.
8. The shallow sea terrain inversion method of claim 1, further comprising: and calculating the water depth of each pixel according to the inversion parameters and the R value of each pixel in the hyperspectral image, and acquiring the shallow sea topography according to the water depth and the sea surface elevation of each pixel.
9. A computer device, comprising: a memory for storing a computer program; a processor for implementing the steps of the shallow sea terrain inversion method of any of claims 1-8 when executing the computer program.
10. A storage medium having a computer program stored thereon, wherein the computer program when executed by a processor implements the steps of the shallow sea terrain inversion method of any of claims 1-8.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210202515.4A CN114594503B (en) | 2022-03-02 | Shallow sea topography inversion method, computer equipment and storage medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210202515.4A CN114594503B (en) | 2022-03-02 | Shallow sea topography inversion method, computer equipment and storage medium |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114594503A true CN114594503A (en) | 2022-06-07 |
| CN114594503B CN114594503B (en) | 2024-06-25 |
Family
ID=
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117237430A (en) * | 2023-11-10 | 2023-12-15 | 中国地质大学（武汉） | High-precision multi-time-sequence water depth inversion method, computing equipment and storage medium |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109059796A (en) * | 2018-07-20 | 2018-12-21 | 国家***第三海洋研究所 | The multispectral satellite remote sensing inversion method of shallow water depth without depth of water control point region |
| CN110032939A (en) * | 2019-03-13 | 2019-07-19 | 浙江工业大学 | A kind of remote sensing time series data approximating method based on gauss hybrid models |
| CN111781146A (en) * | 2020-06-30 | 2020-10-16 | 自然资源部第一海洋研究所 | Wave parameter inversion method using high-resolution satellite optical image |
| KR102342662B1 (en) * | 2021-10-15 | 2021-12-22 | 강원대학교산학협력단 | Subsea geospatial information construction device and construction method |
| KR102355001B1 (en) * | 2021-08-24 | 2022-01-24 | 서울대학교 산학협력단 | Estimation of spatial distribution of suspended sediment concentration from hyperspectral images using machine learning regression models and probabilistic clustering method in rivers |
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109059796A (en) * | 2018-07-20 | 2018-12-21 | 国家***第三海洋研究所 | The multispectral satellite remote sensing inversion method of shallow water depth without depth of water control point region |
| CN110032939A (en) * | 2019-03-13 | 2019-07-19 | 浙江工业大学 | A kind of remote sensing time series data approximating method based on gauss hybrid models |
| CN111781146A (en) * | 2020-06-30 | 2020-10-16 | 自然资源部第一海洋研究所 | Wave parameter inversion method using high-resolution satellite optical image |
| KR102355001B1 (en) * | 2021-08-24 | 2022-01-24 | 서울대학교 산학협력단 | Estimation of spatial distribution of suspended sediment concentration from hyperspectral images using machine learning regression models and probabilistic clustering method in rivers |
| KR102342662B1 (en) * | 2021-10-15 | 2021-12-22 | 강원대학교산학협력단 | Subsea geospatial information construction device and construction method |
Non-Patent Citations (3)
| Title |
|---|
| 张云生等: "Extracting shallow water depth from the fusion of multi-temporal ICESat-2 data and multi-spectral imageries", 《2023INTERNATIONAL CONFERENCE ON FRONTIERS OF OCEAN SCIENCEAND TECHNOLOGY》, 15 March 2024 (2024-03-15) * |
| 杨振等: "黑龙江佳木斯地块南端及邻区深部三维电性结构及其构造意义", 《地球物理学报》, 4 November 2022 (2022-11-04) * |
| 陈本清;杨燕明;罗凯;: "基于高分一号卫星多光谱数据的岛礁周边浅海水深遥感反演", 热带海洋学报, no. 02, 15 March 2017 (2017-03-15) * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117237430A (en) * | 2023-11-10 | 2023-12-15 | 中国地质大学（武汉） | High-precision multi-time-sequence water depth inversion method, computing equipment and storage medium |
| CN117237430B (en) * | 2023-11-10 | 2024-03-08 | 中国地质大学（武汉） | High-precision multi-time-sequence water depth inversion method, computing equipment and storage medium |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Huang et al. | Discharge estimation in high-mountain regions with improved methods using multisource remote sensing: A case study of the Upper Brahmaputra River | |
| Chen et al. | A multiresolution hierarchical classification algorithm for filtering airborne LiDAR data | |
| Mongus et al. | Parameter-free ground filtering of LiDAR data for automatic DTM generation | |
| Cobby et al. | Image processing of airborne scanning laser altimetry data for improved river flood modelling | |
| Xu et al. | Deriving highly accurate shallow water bathymetry from Sentinel-2 and ICESat-2 datasets by a multitemporal stacking method | |
| Yang et al. | Filtering of airborne LiDAR bathymetry based on bidirectional cloth simulation | |
| CN112433227B (en) | Water capacity change monitoring method and system, terminal equipment and storage medium | |
| Shamsoddini et al. | Improving lidar-based forest structure mapping with crown-level pit removal | |
| Zhao et al. | A comparison of LiDAR filtering algorithms in vegetated mountain areas | |
| Chen et al. | A mathematical morphology-based multi-level filter of LiDAR data for generating DTMs | |
| CN115422981B (en) | Land and water classification method and system for single-frequency airborne laser sounding data and application | |
| Kulawiak et al. | Processing of LiDAR and multibeam sonar point cloud data for 3D surface and object shape reconstruction | |
| WO2024036739A1 (en) | Reservoir water reserve inversion method and apparatus | |
| CN112669333A (en) | Single tree information extraction method | |
| CN116758049A (en) | Urban flood three-dimensional monitoring method based on active and passive satellite remote sensing | |
| Zhang et al. | Automatic extraction of DTM from low resolution DSM by twosteps semi-global filtering | |
| Zieger et al. | Mapping reef features from multibeam sonar data using multiscale morphometric analysis | |
| Shao et al. | Automated searching of ground points from airborne lidar data using a climbing and sliding method | |
| CN113516764B (en) | Lake and reservoir underwater three-dimensional terrain simulation method and device based on digital elevation model | |
| CN114594503A (en) | Shallow sea terrain inversion method, computer equipment and storage medium | |
| CN114594503B (en) | Shallow sea topography inversion method, computer equipment and storage medium | |
| CN113432549B (en) | Tidal trench three-dimensional parameter extraction method and system based on unmanned aerial vehicle photogrammetry | |
| Hua et al. | The research of artificial shoreline extraction based on airborne LIDAR data | |
| Abdeldayem | Automatic weighted splines filter (AWSF): A new algorithm for extracting terrain measurements from raw LiDAR point clouds | |
| Youssef et al. | Morphometric analysis of hillslope evolution in the Kadisha River Basin based on archived aerial photographs |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |