Multispectral remote sensing random optimal ratio logarithmic water depth inversion method

Abstract

Aiming at the situation that a water body signal is generally lower than 5%, the invention provides a multispectral remote sensing random optimal ratio logarithmic water depth inversion method, which comprehensively considers the influence of a solar spectrum, an upward total scattering rate, a solar zenith angle, a global total transmittance, a vertical irradiance and a total spherical reflection illuminance, and constructs an atmospheric coupled radiation calibration correction model to correct atmospheric and water surface reflectance. The quotient interpolation tidal height correction model corrects the tidal height, the FNDWI normalized difference water body index accurately extracts the water body, the random optimal model effectively extracts the sub-data set, the extreme gradient lifting method optimizes the model parameters, the optimal ratio logarithmic model is obtained, and water depth inversion is carried out. Compared with the inversion water depth of a Stumpf model, the R2 precision of the inversion water depth of the model is improved by 0.592, and the RMSE is improved by 0.841 m.

Description

Technical Field

[0001] The 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] Depth data provides important data support for marine activities such as navigation ships, underwater topography and landform monitoring, and coastal ecosystem monitoring, and plays an important role in the field of marine research. Traditional depth measurement methods are mainly based on shipborne single-beam and multi-beam sonar and airborne lidar. Although traditional sonar depth measurement and airborne lidar methods are highly accurate, they have become a problem that modern researchers in the field of ocean depth have troubled due to their high costs. Therefore, the depth data obtained by traditional methods is very limited, which limits the acquisition of depths in large areas and multiple islands and reefs. In recent years, inversion of depth based on multispectral remote sensing images has become an effective method for quickly and efficiently obtaining large-scale, high-resolution depth data. Multispectral remote sensing inversion depth uses the measured depth as a known number, combined with sea surface spectral information, to construct and train the inversion model. At present, the acquisition of measured depth data mainly relies on traditional shipborne sonar and airborne lidar. However, the range of measured depth data obtained by these traditional methods is limited. 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 images. By preprocessing images, determining inversion factors, building and verifying multiple models, the blue, green and red bands with high correlation are selected in combination with Sentinel-2A images to select the optimal model for water depth inversion. It has the characteristics of wide coverage, fast update and low cost, and effectively makes 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, which 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, which realizes the refined extraction of underwater terrain photons and improves 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 sea water depth inversion based on physical constraint neural network", has 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 moment of imaging, and solves the problem of large errors in depth inversion and terrain inversion in high turbidity waters.

[0007] The invention patent "A method for water depth inversion using satellite remote sensing images" applied for by Yang Zhixiang and others from China Railway Water Resources and Hydropower Planning and Design Group Co., Ltd. has patent number CN202311666078.2. It is based on the water area representativeness index and selects 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, "A method for measuring water depth inversion based on satellite-borne active and passive remote sensing information during high and low tides", has patent number CN202311666078.2. It substitutes multi-spectral data at high tide and low tide into the constructed water depth inversion model to obtain water depth inversion results in deep water areas, providing 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 titled "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 view of the above problems, the present invention discloses a multi-spectral remote sensing random optimal ratio logarithmic water depth inversion method.

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

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

[0014] S1: Since the water body signal is generally less than 5%, the influence 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] In the formula, For the Band; is the gray value of the i-th band; For the The absolute calibration gain value of the band; For the The offset value of the band; The true reflectivity of the ground surface after calibration correction for atmospheric coupled radiation; is the atmospheric correction factor calculated by this model.

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

[0018] (2)

[0019] In the formula, (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 tidal height correction model is constructed to obtain the instantaneous tidal height, and the tidal height is corrected using formula (3):

[0021] (3)

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

[0023]

[0024] In the formula, The transit time of satellite remote sensing images , It is earlier than the 0th hour when the satellite remote sensing image passes by. 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 various orders 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] In the formula, and They are the reflectance of the green band and short infrared band in satellite remote sensing images respectively. The threshold 0 is used to determine whether it is a water body, that is, FNDWI>0 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] In the formula, Indicates Secondary input data characteristics; It represents the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tidal height correction;

[0032] S5: Construct the optimal ratio logarithmic model for water depth inversion. Parameter optimization model through random algorithm;

[0033] S5-1: Logarithmic Odds Model:

[0034] (5)

[0035] In the formula, is the scaling constant of the ratio to water depth, It is the green band. It is the blue band, It is water depth Right The water reflectivity of the band, It is the offset of water depth relative to the horizontal plane of 0m.

[0036] S5-2: Construct a random optimal value model to 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] In the formula, 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 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] In the formula, 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] In the formula, 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 is to add 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 satellite time tidal height calculation by using more than two data;

[0051] 3. In view of the fact that near-infrared band pixels may confuse ground object information, combined with the low reflectivity of the water surface 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 characteristics of the attenuation ratio of the blue-green band in underwater propagation, which can better invert the water depth and is insensitive to the bottom quality;

[0053] 5. To address the problem of large data volume, the random algorithm model provided by the present invention extracts sub-datasets, thereby 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 iterates the optimization of model parameters. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

[0058] FIG. 3 is 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 real water depth

[0061] Figure 7 Scatter plot of the inverted water depth and the real water depth of the present invention DETAILED DESCRIPTION

[0062] In order to make the purpose, 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 in conjunction with the accompanying drawings.

[0063] Example:

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

[0065] S1: Combination Figure 2 , indicating that the experimental area is located in the coastal area on the west side 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] Combination 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 are shown in Table 1), the image was taken at 3:10 am on April 9, 2022 (UTC), and the bands 1-7 and band 10 were used (parameters are shown in Table 2).

[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 correction 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 a radiation brightness value, and the calibration results are saved.

[0074] Based on the atmospheric coupled radiation calibration correction model, , , The three coefficients are obtained by formatting and editing the values ​​to create an input text file and output it. Save the input value file as a txt file (in.txt), put it in the same directory as the model execution file (main.exe), and enter the following command in the command line: main<in.txt> out.txt, call the model to calculate and 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 upward and downward radiation transfer process;

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

[0077] The atmospheric correction coefficients of all bands of Landsat-9 are applied to the radiance image using IDL language to achieve atmospheric correction. Compile and execute the IDL program to obtain the real radiance file of the surface.

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

[0079] Combination Figure 4 , indicating that the water depth data of the experimental area of ​​this embodiment is obtained from the website https: / / www.gebco.net / (parameters are shown in Table 3), ArcGis is used to resample to 30 meters, and the water depth data consistent with the size of the experimental area is cut out in ENVI;

[0080] Table 3. Overview of DBMs dataset

[0081]

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

[0083] Combination 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 point in the measured water depth data set are calculated through the difference quotient interpolation tide height correction model (part of the parameters are shown in Table 4), and the tide height correction is performed to match the data set.

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

[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 initial data input is performed under the IDL environment; the measured water depth data after tidal correction in the target area is used as the secondary input data feature to obtain the data sample set , Indicates Secondary input data characteristics; It represents the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tidal 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 smallest error are obtained.

[0089] S4-2: Based on the model of randomly taking the optimal value, further optimize the extreme gradient boosting method to optimize the model parameters. Input the data set, train the optimized model, traverse to obtain the best parameter combination of the model, and use regularization to control the complexity of the model to prevent overfitting.

[0090] S5: Combination 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 only used to illustrate the present invention, but not to limit the present invention. A person 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 also belong to the scope of the present invention, and the patent protection scope of the present invention shall be defined by the claims.

[0092] The technical contents not described in detail in the present invention are all 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 influence 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)； In the formula, For the Band; is the gray value of the i-th band; For the The absolute calibration gain value of the band; For the The offset value of the band; The true reflectivity of the ground surface after calibration correction for atmospheric coupled radiation; The atmospheric correction factor calculated for this model; in, The coefficient can be obtained through formula (2): (2)； In the formula, (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 tidal height correction model is constructed to obtain the instantaneous tidal height, and the tidal height is corrected using formula (3): (3)； In the formula, represents instantaneous water depth data, Indicates the measured water depth data, It is the difference between mean sea level and tidal height datum. represents the tidal height value; Through the difference quotient interpolation transformation, we get: ； In the formula, The transit time of satellite remote sensing images , It is earlier than the 0th hour when the satellite remote sensing image passes by. 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 various orders 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)； In the formula, and They are the reflectance of the green band and the short infrared band in the satellite remote sensing image 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 water and land mask is performed; S4: Construct a dataset based on satellite transit images and water depth data: ； In the formula, Indicates Secondary input data characteristics; It represents the target area value corresponding to the remote sensing data feature, that is, the measured water depth data after tidal 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)； In the formula, is the scaling constant of the ratio to water depth, It is the green band. It is the blue band, It is water depth Right The water reflectivity of the band, is the offset relative to the water depth of 0m horizontal plane; S5-2: Construct a random optimal value model to find the optimal model parameter data set, and then find the A and B values ​​with the smallest error: (6)； further: (7)； In the formula, 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: Construct extreme gradient boosting method to optimize model parameters; based on S5-2 dataset , where Indicates 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)； In the formula, 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)； In the formula, 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 is to add 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