Multispectral remote sensing random optimal ratio logarithmic water depth inversion method

Abstract

As water body signals are generally less than 5%, this paper proposes a multispectral remote sensing random optimal ratio logarithmic water depth inversion method. It comprehensively considers the influence of solar spectrum, total upward scattering rate, solar zenith angle, global total transmittance, vertical irradiance and total spherical albedo, constructs an atmospheric coupled radiation calibration correction model to correct the atmospheric and water surface reflectivity, and uses a quotient interpolation tide height correction model to correct the tide height. The FNDWI normalized difference water index accurately extracts water bodies, and the random optimal model effectively extracts sub-datasets. The extreme gradient boosting method optimizes model parameters to obtain the optimal ratio logarithmic model for water depth inversion. Compared with the water depth inversion of the Stumpf model, the R of the water depth inversion of this model is better than that of the Stumpf model. 2 The accuracy is improved by 0.592 and the RMSE is improved by 0.841m.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing technology, and in particular to a multi-spectral remote sensing random optimal ratio logarithmic water depth inversion construction method. Background Art

[0002] Bathymetric data provides crucial data support for maritime activities such as navigation, underwater topography, and coastal ecosystem monitoring, playing a crucial role in oceanographic research. Traditional bathymetric measurement methods primarily rely on shipborne single-beam and multi-beam sonar and airborne lidar. While highly accurate, these methods are also plagued by high costs, a major challenge for modern oceanographic researchers. Consequently, the bathymetric data obtained by these methods is extremely limited, restricting the acquisition of bathymetric data over large areas and around multiple islands and reefs. In recent years, bathymetric inversion based on multispectral remote sensing imagery has become an effective method for rapidly and efficiently acquiring large-scale, high-resolution bathymetric data. Multispectral remote sensing bathymetric inversion uses measured water depths as known quantities and combines them with sea surface spectral information to construct and train an inversion model. Currently, obtaining measured water depth data primarily relies on traditional methods such as shipborne sonar and airborne lidar. However, these traditional methods only capture a limited range of bathymetric data. Therefore, the use of multispectral remote sensing technology for satellite bathymetry still faces many challenges, and its accuracy and reliability still need further research.

[0003] The international patent "Shallow sea area waterdepth inversion method based on Sentinel-2A image" authorized by Li Haibin and others from the 92859th Unit of the Chinese People's Liberation Army, with publication number 114201732, proposes a water depth inversion method for shallow waters based on Sentinel-2A imagery. By preprocessing images, determining inversion factors, and constructing and verifying multiple models, the method selects the optimal model for water depth inversion using the highly correlated blue, green, and red bands combined with Sentinel-2A imagery. This method has the characteristics of wide coverage, fast update, and low cost, effectively making up for the shortcomings of traditional ocean measurements.

[0004] The international patent "Method for selecting hyperspectral water depth inverse waveband" authorized by Chen Jin and others from the Beijing Remote Sensing Information Research Institute, with publication number 108240806, proposes a hyperspectral water depth inverse band selection method that automatically selects bands using spectral band correlation and signal-to-noise ratio information, including band correlation analysis, grouping, signal-to-noise ratio analysis, water depth correlation analysis, and output of bands according to the requirements of the water depth inversion model.

[0005] The invention patent applied for by Liu Zhendong and others from Shandong University of Science and Technology, "A method for water depth inversion based on satellite-borne lidar point cloud", has the patent number CN202410912661.5. It proposes using underwater terrain photons to replace in-situ bathymetric point training and constructing a satellite-borne lidar active and passive fusion water depth inversion model, realizing the refined extraction of underwater terrain photons and improving the accuracy of water depth inversion.

[0006] The invention patent applied for by Xu Qing and others from Ocean University of China, "A method for constructing dual-polarization SAR shallow water depth inversion based on physical constraint neural network", has the patent number CN202410739992.3. It introduces physical auxiliary information of the polarization ratio of SAR remote sensing images and the wind field vector at the imaging moment, solving the problem of large errors in depth inversion and terrain inversion in high turbidity waters.

[0007] The invention patent applied for by Yang Zhixiang and others from China Railway Water Resources and Hydropower Planning and Design Group Co., Ltd. is "A method for water depth inversion using satellite remote sensing images". The patent number is CN202311666078.2. It is based on the water area representativeness index. By selecting the optimal combination of measurement points in sub-areas that do not meet the standards, a water depth inversion model is established according to the reflectivity, so that high-precision water depth data can still be obtained when the water area changes.

[0008] The invention patent applied for by Yang Fanlin and others from Shandong University of Science and Technology is "A method for measuring water depth inversion based on satellite-borne active and passive remote sensing information at high and low tides". The patent number is CN202311666078.2. The water depth inversion model constructed by substituting multi-spectral data at high tide and low tide is used to obtain water depth inversion results in deep water areas. It provides support for water depth inversion of Class II water bodies in coastal zones and water depth data learning in areas without measured water depth data.

[0009] The invention patent applied for by Le Yuan and others from the East China Sea Laboratory of China University of Geosciences (Wuhan) is "A high-precision multi-time series water depth inversion method, computing equipment and storage medium". The patent number is CN202311489467.2. By constructing a water depth inversion model, multiple single-time series water depth inversion result data sets are generated, and multiple time-series water depth inversion data sets are composed for traversal. The maximum outlier method is used to eliminate invalid error data sets, and the water depth distribution histogram generated by the valid data set is Gaussian fitted to obtain the water depth distribution curve. The least squares method is used to fit the water depth aquatic distribution curve parameters to further generate the water depth value of the target area.

[0010] However, the above results do not fully utilize the inherent spectral characteristics of water bodies, which will further increase the error of inverted water depth. Summary of the Invention

[0011] In response to the above problems, the present invention discloses a multispectral remote sensing random optimal ratio logarithmic water depth inversion method.

[0012] The present invention is achieved by adopting the following technical solutions:

[0013] Multispectral remote sensing random optimal ratio logarithmic water depth inversion method, including:

[0014] S1: Since the water body signal is generally less than 5%, the effects of solar spectrum, total upward scattering rate, solar zenith angle, global total transmittance, vertical irradiance and total spherical albedo are comprehensively considered to construct an atmospheric coupled radiation calibration correction model to correct the atmospheric and water surface reflectance, see formula (1):

[0015] (1)

[0016] Where, For the Band; For the Grayscale value of the band; For the The absolute calibration gain value of the band; For the The offset value of the band; The true surface reflectance after correction for atmospheric coupled radiation calibration; is the atmospheric correction factor calculated by this model.

[0017] in, The coefficient can be obtained through formula (2):

[0018] (2)

[0019] Where, (funct filter) function filter, (sol spect) is the solar spectrum, (upward total scattering) is the total upward scattering rate, is the solar zenith angle, (total global transmittance) is the global total transmittance, (downward total scattering) is the total downward scattering rate, (direct irradiance) is the vertical irradiance, (total spherical albedo) is the total spherical albedo.

[0020] S2: According to the transit time of satellite remote sensing images, a difference quotient interpolation tide height correction model is constructed to obtain the instantaneous tide height, and the tide height correction is performed using formula (3):

[0021] (3)

[0022] Where, Indicates instantaneous water depth data, Indicates the measured water depth data, It represents the difference between mean sea level and tidal datum. Indicates the tidal height value. Through the difference quotient interpolation transformation, we can get:

[0023]

[0024] Where, The transit time of satellite remote sensing images , It is earlier than the 0th hour when the satellite remote sensing image passes through. It is the first hour when the satellite remote sensing image passes by. It is the second hour when the satellite remote sensing image passes by. Later than the satellite remote sensing image passing time, Later than the satellite remote sensing image passing At this moment, the tide height data point is known , To get Calculate The calculated difference quotients of each order at each known point are as follows:

[0025]

[0026] S3: When extracting water body information, near-infrared band pixels may confuse ground object information. Combined with the low reflectivity of the short-wave infrared band on the water surface, the FNDWI normalized difference water index is constructed to accurately extract the water body range:

[0027] (4)

[0028] Where, and The FNDWI is the reflectivity of the green band and short infrared band in satellite remote sensing images, respectively. The threshold of 0 is used to determine whether it is a water body, that is, if FNDWI>0, it is a water body, and then a water and land mask is performed.

[0029] S4: Construct a dataset based on satellite transit images and water depth data:

[0030]

[0031] Where, Indicates the Secondary input data characteristics; Represents the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tide height correction;

[0032] S5: Construct the optimal ratio logarithmic model for water depth inversion. Optimize the model parameters through randomized algorithms.

[0033] S5-1: Logarithmic Odds Model:

[0034] (5)

[0035] Where, is the scaling constant of the ratio to water depth, It's the green band. It's the blue band, Is the water depth The water reflectivity in the green band is Is the water depth The water reflectivity in the blue band is It is the offset of water depth relative to the horizontal plane of 0m.

[0036] S5-2: Construct a random optimal value model, find the optimal model parameter data set, and then find the A and B values ​​with the smallest error:

[0037] (6)

[0038] further:

[0039] (7)

[0040] Where, is the water depth value of the training set, >0 is the probability of obtaining the minimum difference between A and B, is the time required to find the minimum difference, and It is the time required to find the maximum or minimum difference between A and B.

[0041] S5-3: Construct extreme gradient boosting method to optimize model parameters. Based on S5-2 dataset , where Indicates the Secondary input data characteristics; It represents the corrected water depth value of the target area corresponding to the remote sensing data characteristics, and optimizes the model parameters. The specific formula is shown in (8):

[0042] (8)

[0043] Where, is the water depth prediction value of the model; is the loss function, usually the mean square error; is a regularization term that controls the complexity of the model.

[0044] For each tree , then:

[0045] (9)

[0046] Where, is the number of leaf nodes in the tree, is the weight of the leaf node, and is a hyperparameter that controls the complexity of the model. Gradient boosting involves adding a new tree each time to fit the residual of the current model; the output of the new tree is a fit to the residual, thereby updating the model:

[0047] (10)

[0048] Compared with the prior art, the present invention has significant improvements in that:

[0049] 1. As the water body signal is generally less than 5%, the present invention comprehensively considers the influence of solar spectrum, total upward scattering rate, solar zenith angle, global total transmittance, vertical irradiance and total spherical albedo, and constructs an atmospheric coupled radiation calibration correction model to correct the atmospheric and water surface reflectivity;

[0050] 2. The difference quotient interpolation tidal height correction model constructed by the present invention can realize the calculation of satellite time tidal height by using more than two data;

[0051] 3. Considering that near-infrared band pixels may confuse ground object information, combined with the low reflectivity of the water surface in the short-wave infrared band, the present invention constructs the FNDWI normalized difference water index to extract the water body range;

[0052] 4. The logarithmic ratio model provided by the present invention utilizes the attenuation ratio of the blue-green band in underwater propagation, which can better invert the water depth and is insensitive to the bottom sediment.

[0053] 5. To address the problem of large data volumes, the random algorithm model provided by the present invention extracts sub-datasets, reducing the time complexity of data extraction;

[0054] 6. To address the difficulty in inverting model parameters, the extreme gradient boosting method provided by the present invention updates and iteratively optimizes model parameters. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0056] Figure 2This is the experimental area location of the present invention

[0057] Figure 4 The water depth data of the present invention

[0058] FIG3 is the Landsat-9 satellite image data of the present invention

[0059] Figure 5 This is the tide height correction process of the present invention

[0060] Figure 6 Scatter plot of water depth inverted by classic model and true water depth

[0061] Figure 7 Scatter plot of the inverted water depth and the true water depth of the present invention Specific implementation plan

[0062] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, preferred embodiments are given below and the specific implementation of the present invention is further described in detail with reference to the accompanying drawings.

[0063] Example:

[0064] Combine Figure 1 , describing the specific process of this embodiment.

[0065] S1: Combine Figure 2 , indicating that the experimental area is located in the western coastal area of ​​Weizhou Island, Beihai City, Guangxi Zhuang Autonomous Region, with a latitude range of 20°03'N-20°01'N and a longitude range of 109°03'E-109°05'E.

[0066] S2: Data and its preprocessing

[0067] S2-1: Image data and its preprocessing

[0068] Combine Figure 3 , indicating that the image used in the experimental area is Landsat-9 remote sensing image, the product type is Collection 2 Level 1 (parameters see Table 1), the image shooting time is April 9, 2022 3:10 (UTC), and the bands 1-7 and band 10 used are shown in Table 2 (parameters).

[0069] Table 1. Landsat 9 satellite orbit parameters

[0070]

[0071] Table 2. Landsat 9 satellite payload technical specifications

[0072]

[0073] Based on formulas (1)-(2), atmospheric coupled radiation calibration is performed in the IDL environment and the results are returned to ENVI. The calibration coefficients of each band are read in, the DN value (Digital Number) is converted to radiance value, and the calibration results are saved.

[0074] Based on the atmospheric coupled radiation calibration correction model, 、 、 After the three coefficients are obtained, the input text file is created by formatting and editing the values ​​and then output. Save the input value file as a txt file (in.txt) and put it in the same directory as the model execution file (main.exe). Enter the command line: main<in.txt> out.txt, call the model calculation to get the output result file out.txt.

[0075] The output file contains four parts: the first part lists the input parameters; the second part gives the direct, scattered and ambient radiation components in the atmospheric radiation transfer; the third part gives the influence of various atmospheric components on the uplink and downlink radiation transfer processes;

[0076] The fourth part of the output file includes the three coefficient values ​​required for atmospheric correction of remote sensing images ( 、 、 ) .

[0077] Apply atmospheric correction coefficients for all Landsat-9 bands to the radiance image using the IDL language. Compile and execute the IDL program to obtain the true surface radiance file.

[0078] S2-2: Obtaining water depth ENCs electronic nautical charts

[0079] Combine Figure 4 , indicating that the water depth data of the experimental area of ​​this embodiment was obtained from the website https: / / www.gebco.net / (parameters are shown in Table 3), resampled to 30 meters using ArcGis, and the water depth data was cropped in ENVI to match the size of the experimental area;

[0080] Table 3. Overview of DBMs dataset

[0081]

[0082] S2-3: Obtain tide height data and correct tide height

[0083] Combine Figure 5, indicating that the tide table of the experimental area of ​​this embodiment is obtained through the website. Based on formula (3), according to the transit time of the Landsat-9 satellite remote sensing image, the corresponding instantaneous tide level data of the water depth points in the measured water depth dataset are calculated through the difference quotient interpolation tide height correction model (part of the parameters are shown in Table 4), and the tide height is corrected to match the dataset.

[0084] Table 4. Water depth data points of some training sets and their satellite transit time

[0085]

[0086] S3: Based on formula (4), the satellite remote sensing data is masked and the threshold 0 is used to determine whether it is a water body, that is, FNDWI>0 is a water body, and then the land and water mask separation is performed. After masking, the water body pixel value in the study area is 1, and the non-water body pixel value is 0.

[0087] S4: Based on formula (5), the blue band and green band of the Landsat-9 satellite remote sensing image corresponding to the image data point are selected as the input characteristic spectral parameters, and the data is initially input in the IDL environment; the measured water depth data after tidal correction of the target area is used as the secondary input data feature to obtain the data sample set , Indicates the Secondary input data characteristics; Indicates the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tide height correction;

[0088] S4-1: Based on formulas (6)-(7), a random optimal value model is constructed to find the optimal model parameter data set, and then the A and B values ​​with the minimum error are obtained.

[0089] S4-2: Based on the random optimal value model, we further optimize the model parameters using the extreme gradient boosting method. We input the dataset and train the optimized model. After traversing the dataset, we find the optimal parameter combination for the model. We also use regularization to control the model complexity to prevent overfitting.

[0090] S5: Combine Figure 6 and 7 , verifying the accuracy of the present invention. The present invention carried out the inversion depth experiment of the Stumpf model in the same area and image data. Figure 6 and 7 It can be seen that compared with the Stumpf model, the R 2 The accuracy is improved by 0.592 and the RMSE is improved by 0.841m.

[0091] The above embodiments are intended to illustrate the present invention only and are not intended to limit the present invention. Persons skilled in the art may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions fall within the scope of the present invention, and the scope of patent protection of the present invention shall be defined by the claims.

[0092] The technical contents not described in detail in this invention are all well-known technologies.

Claims

1. Multispectral remote sensing random optimal ratio logarithmic water depth inversion method, characterized by The following steps are involved: S1: Since the water body signal is generally less than 5%, the effects of solar spectrum, total upward scattering rate, solar zenith angle, global total transmittance, vertical irradiance and total spherical albedo are comprehensively considered to construct an atmospheric coupled radiation calibration correction model to correct the atmospheric and water surface reflectance, see formula (1): (1)； Where, For the Band; For the Grayscale value of the band; For the The absolute calibration gain value of the band; For the The offset value of the band; The true surface reflectance after correction for atmospheric coupled radiation calibration; The atmospheric correction factor calculated for this model; in, The coefficient can be obtained through formula (2): (2)； Where, (funct filter) function filter, (sol spect) is the solar spectrum, (upwardtotal scattering) is the total upward scattering rate, is the solar zenith angle, (total global transmittance) is the global total transmittance, (downward total scattering) is the total downward scattering rate, (direct irradiance) is the vertical irradiance, (total spherical albedo) is the total spherical albedo; S2: According to the transit time of satellite remote sensing images, a difference quotient interpolation tide height correction model is constructed to obtain the instantaneous tide height, and the tide height correction is performed using formula (3): (3)； Where, Indicates instantaneous water depth data, Indicates the measured water depth data, It represents the difference between mean sea level and tidal datum. represents the tidal height value; Through the difference quotient interpolation transformation, we can get: ； Where, The transit time of satellite remote sensing images , It is earlier than the 0th hour when the satellite remote sensing image passes through. It is the first hour when the satellite remote sensing image passes by. It is the second hour when the satellite remote sensing image passes by. Later than the satellite remote sensing image passing time, Later than the satellite remote sensing image passing At this moment, the tide height data point is known , To get Calculate The calculated difference quotients of each order at each known point are as follows: ； S3: When extracting water body information, near-infrared band pixels may confuse ground object information. Combined with the low reflectivity of the short-wave infrared band on the water surface, the FNDWI normalized difference water index is constructed to accurately extract the water body range: (4)； Where, and The reflectances of the green band and short infrared band in satellite remote sensing images are respectively; the threshold 0 is used to determine whether it is a water body, that is, FNDWI>0 is a water body, and then the land and water mask is performed; S4: Construct a dataset based on satellite transit images and water depth data: ； Where, Indicates the Secondary input data characteristics; Represents the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tide height correction; S5: Construct the optimal ratio logarithmic model for water depth inversion; optimize the model parameters through random algorithm; S5-1: Logarithmic Odds Model: (5)； Where, is the scaling constant of the ratio to water depth, It's the green band. It's the blue band, Is the water depth The water reflectivity in the green band is Is the water depth The water reflectivity in the blue band is is the offset relative to the water depth at 0m horizontal plane; S5-2: Construct a random optimal value model, find the optimal model parameter data set, and then find the A and B values ​​with the smallest error: (6)； further: (7)； Where, is the water depth value of the training set, >0 is the probability of obtaining the minimum difference between A and B, is the time required to find the minimum difference, and It is the time required to find the maximum or minimum difference between A and B; S5-3: Constructing extreme gradient boosting method to optimize model parameters; based on S5-2 dataset , where Indicates the Secondary input data characteristics; It represents the corrected water depth value of the target area corresponding to the remote sensing data characteristics, and optimizes the model parameters. The specific formula is shown in (8): (8)； Where, is the water depth prediction value of the model; is the loss function, usually the mean square error; is a regularization term that controls the complexity of the model; For each tree , then: (9)； Where, is the number of leaf nodes in the tree, is the weight of the leaf node, and It is a hyperparameter that controls the complexity of the model. Gradient boosting involves adding a new tree each time to fit the residual of the current model. The output of the new tree is a fit to the residual, thereby updating the model. (10)。

Images