CN112013822A - Multispectral remote sensing water depth inversion method based on improved GWR model - Google Patents

Multispectral remote sensing water depth inversion method based on improved GWR model Download PDF

Info

Publication number
CN112013822A
CN112013822A CN202010711999.6A CN202010711999A CN112013822A CN 112013822 A CN112013822 A CN 112013822A CN 202010711999 A CN202010711999 A CN 202010711999A CN 112013822 A CN112013822 A CN 112013822A
Authority
CN
China
Prior art keywords
water depth
model
inversion
point
data
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.)
Pending
Application number
CN202010711999.6A
Other languages
Chinese (zh)
Inventor
胡娜
余利剑
李叶繁
刘修国
花卫华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan Zhitu Yunqi Technology Co ltd
Original Assignee
Wuhan Zhitu Yunqi Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Wuhan Zhitu Yunqi Technology Co ltd filed Critical Wuhan Zhitu Yunqi Technology Co ltd
Priority to CN202010711999.6A priority Critical patent/CN112013822A/en
Publication of CN112013822A publication Critical patent/CN112013822A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C13/00—Surveying specially adapted to open water, e.g. sea, lake, river or canal
    • G01C13/008—Surveying specially adapted to open water, e.g. sea, lake, river or canal measuring depth of open water
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/55—Specular reflectivity
    • 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
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/02—Systems using the reflection of electromagnetic waves other than radio waves
    • G01S17/06—Systems determining position data of a target
    • 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
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/1793—Remote sensing

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Electromagnetism (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Analytical Chemistry (AREA)
  • Pathology (AREA)
  • Immunology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biochemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Hydrology & Water Resources (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

The invention provides a multispectral remote sensing water depth inversion method based on an improved GWR model, which comprises the following steps: s1, preprocessing the image to obtain reflectivity data; s2, corresponding the actually measured water depth data with the reflectivity data according to the position information; s3, constructing a first water depth inversion model to perform water depth inversion on the whole research area, and obtaining water depth data corresponding to each sample point in the research area; s4, constructing a second water depth inversion model based on the improved GWR model to perform pixel-by-pixel coefficient calibration to realize water depth inversion; and S5, comparing the inversion results of the two models by using the verification points to verify the precision, and outputting a water depth value obtained by inverting the second water depth inversion model after the verification is finished. The invention has the beneficial effects that: vertical heterogeneity is considered, normalized vertical water depth information and horizontal coordinates are added into calculation of an optimal bandwidth and a space weight matrix, and the method is suitable for sea areas with large underwater terrain variation.

Description

Multispectral remote sensing water depth inversion method based on improved GWR model
Technical Field
The invention relates to the technical field of space remote sensing, in particular to a multispectral remote sensing water depth inversion method based on an improved GWR model.
Background
The seawater depth is an important hydrological factor and is essential basic geographic space information for activities such as marine navigation, marine environment treatment, marine engineering construction, marine resource development and the like. Although the traditional water depth measurement method is high in precision, the problems of low measurement efficiency and easiness in missing small obstacle targets generally exist, and the traditional water depth measurement method is easily influenced by factors such as measurement time, weather and waves.
The physical basis of the water depth remote sensing is that light rays have certain penetrating capacity on a water body, and a water depth inversion model algorithm can be obtained by a simplified radiation transmission model, wherein the simplified radiation transmission model is as follows:
Li＝Li∞+CiRBi(λ)e-kifZ
in the formula, LiRepresenting the radiation value, L, of the i-th band received by the sensori∞Represents the radiation value of the wave band in a deep water region, CiDenotes a constant, k, related to solar irradiance, atmospheric and water surface transmittance, water surface refraction, and the likeiRepresenting the diffuse attenuation coefficient, RBiThe reflectivity of the wave band at the bottom of the water body is shown, lambda represents the wavelength, f represents the path length of the water body (usually 2), and Z represents the water depth value.
At present, a plurality of water depth measurement methods based on multispectral remote sensing images are proposed, and in a commonly used semi-empirical semi-theoretical model for multispectral remote sensing water depth inversion, a single-band model, a dual-band model and a multiband model are all converted by applying a semi-theoretical semi-empirical concept based on the radiation transmission model, namely:
Figure BDA0002596902090000011
wherein Z represents a depth value of water, a0、aiAnd (i is 1, …, n) represents the coefficient to be determined of the model, n represents the number of wave bands participating in inversion, and i represents the ith wave band and represents an error term. However, most of the methods need the assumed conditions that the atmosphere is uniform, the seabed sediment is single, the water quality is uniform, and the water attenuation coefficients are equal and difficult to meet, but the water attenuation coefficients caused by the change of the sediment type, the water quality and the water depth cause the problem of spatial heterogeneity, so that the fitting model coefficients at different positions are different; the semi-empirical semi-theoretical model only reflects the relation between the spatial average of the spectral radiation and the water depth, the water depth obtained by the inversion of the global model still generates a large error, and the single water depth inversion model cannot ensure the water depth inversion accuracy under the influence of the spatial heterogeneity problem.
A common method for solving the problem of spatial heterogeneity is to adopt a Geographical Weighted Regression (GWR), and compared with a global model, the GWR model considers the problem of spatial heterogeneity caused by water quality and substrate variation in the water body in the horizontal direction, and embeds the spatial structure of data into a common linear regression water depth inversion model, and establishes different spatial relationships at different geographical positions, and the model formula is as follows:
Figure BDA0002596902090000021
wherein (x, y) represents the horizontal coordinate of the sample point, Z (x, y) represents the water depth at the sample point (x, y), a0(x,y)、ai(x, y) represents the undetermined coefficient of the sample point (x, y) in the GWR model,ithe error term for the ith band is indicated.
The GWR model estimates model parameters of all sample points in the data set, thereby obtaining a better water depth inversion result. However, in a sea area with obvious underwater terrain variation, vertical spatial heterogeneity becomes large, and the variation of the vertical water depth causes different water attenuation coefficients to influence model coefficients, so that a GWR model considering only horizontal heterogeneity cannot obtain a good inversion result in the sea area with large underwater terrain variation.
Disclosure of Invention
In view of the above, the invention provides a multispectral remote sensing water depth inversion method based on an improved GWR model, and aims to solve the problem that the existing GWR model lacks consideration of vertical heterogeneity problem, so that the inversion accuracy of water depth in sea areas with large vertical heterogeneity is low.
The invention provides a multispectral remote sensing water depth inversion method based on an improved GWR model, which comprises the following steps:
s1, performing image preprocessing on the multispectral remote sensing image of the research area to obtain water surface reflectivity data;
s2, acquiring actually measured water depth data, and corresponding the actually measured water depth data to the reflectivity data of the multispectral remote sensing image according to the position information;
s3, constructing a first water depth inversion model, and obtaining water depth data of the whole research area by using the first water depth inversion model;
s4, constructing a second water depth inversion model based on the improved GWR model:
Figure BDA0002596902090000031
in the formula (x)j,yj,zj) Represents the spatial coordinates of the real measuring point j, where zjRepresenting water depth data obtained at the actual measuring point j according to the first water depth inversion model; zj(x, y, z) represents the water depth of the sample point (x, y, z) obtained by inverting the second water depth inversion model according to the real-time point j; a is0(xj,yj,zj)、ai(xj,yj,zj) Representing undetermined coefficients of actual measurement points j in the second water depth inversion model; xi＝ln[ρi-ρi∞]，ρiDenotes the reflectivity, p, of the i-th bandi∞The reflectivity of the wave band in a deep water area is shown, i is 1, …, n is n, and n is the number of the wave bands participating in inversion;j2an error term representing the second water depth inversion model;
taking the logarithm value of the reflectivity on the multispectral image corresponding to the actual measurement point, the actual measurement point water depth and the actual measurement point position information as modeling data of the second water depth inversion model, taking the logarithm value of the reflectivity of each wave band corresponding to each sample point in the research area, the water depth data obtained based on the first water depth inversion model and the corresponding position information as input data of the second water depth inversion model, and performing pixel-by-pixel coefficient calibration to realize water depth inversion;
and S5, comparing inversion results of the first water depth inversion model and the second water depth inversion model by using the verification points, so as to perform precision verification, and outputting the water depth value obtained by the second water depth inversion model after the verification is completed.
Further, in step S3, a first water depth inversion model is constructed based on the GWR model:
Figure BDA0002596902090000041
in the formula (x)j,yj) Horizontal coordinate, Z, representing real measuring point jj(x, y) represents the water depth at the sample point (x, y) obtained by inversion from the actual point j, a0(xj,yj)、ai(xj,yj) Representing undetermined coefficients of actual measurement points j in the first water depth inversion model;j1an error term representing the first water depth inversion model;
and taking the logarithmic value of the reflectivity of each wave band of the actually measured water depth data and the position information thereof as modeling data of the first water depth inversion model, and performing water depth inversion on the whole research area according to a weighted least square method to obtain water depth data corresponding to each sample point in the research area.
Further, in step S1, the image preprocessing process includes:
carrying out radiometric calibration on the multispectral remote sensing image, and converting a pixel value DN into a spectral radiance value L at the top of the atmosphere, wherein the pixel value L is DN, absCalfactor/delta lambda, wherein absCalfactor represents an absolute calibration factor, and delta lambda represents the effective width of a wave band;
atmospheric correction is carried out on the radiance image of each wave band obtained by radiometric calibration, and reflectivity data of each wave band is obtained;
and extracting a water body part image from the reflectivity data by adopting an improved normalized water body index, and finally removing solar flare to finish image preprocessing.
Further, the specific process of performing water depth inversion on the whole research area according to the weighted least square method in step S3 is as follows: for any sample point (x, y) in the research area, coefficient calibration is carried out by using a weighted least square method to obtain a model coefficient:
aj(x,y)＝(XTWj(x,y)X)-1XTWj(x,y)Z
in the formula, aj(X, y) represents a regression coefficient of the actual measurement point j at the sample point (X, y), where X ═ X1,…,Xn]TZ represents a vector formed by actual-point water depth data, Wj(x, y) represents the weight matrix of the real measurement point j at sample point (x, y).
Further, the weight matrix WjThe diagonal elements of (x, y) are a function of the coordinates (x, y) of the sample points to ensure that real points closer to the sample points have greater weight.
Further, in step S4, the specific process of performing water depth inversion includes:
s41, carrying out normalization processing on the plane coordinate data and the vertical water depth data, and estimating the optimal bandwidth by adopting a cross validation method:
Figure BDA0002596902090000051
wherein CV represents the sum of squares of the cross-validation residuals, (x, y, z) represents the three-dimensional space coordinates of the regression point k, scale (. cndot.) represents the normalization function, N represents the total number of data calculated by the regression parameters, and b represents the bandwidth，ZkRepresents the measured water depth value at the regression point k,
Figure BDA0002596902090000052
representing a fitted water depth value which does not include regression points to carry out regression parameter calculation per se; the bandwidth corresponding to the minimum CV value is the optimal bandwidth;
and S42, calculating the space weight by adopting a bi-square function according to the optimal bandwidth b:
Figure BDA0002596902090000053
in the formula (d)kjRepresents the distance, d, between the real point j and the regression point kbRepresenting the distance of the inversion point to the b-th nearest neighbor;
s43, taking any sample point (x, y, z) in the multispectral remote sensing image as an inversion point, calculating the optimal bandwidth and weight of each corresponding wave band by utilizing the steps S41 and S42 to obtain a weight matrix, and further estimating the model coefficient by adopting a weighted least square method:
aj(x,y,z)＝(XTWj(x,y,z)X)-1XTWj(x,y,z)Z
in the formula, aj(X, y, z) represents a regression coefficient of the measured point j at the sample point (X, y, z), where X is [ X ═ X1,…,Xn]TZ represents a vector formed by actual-point water depth data, Wj(x, y, z) represents a weight matrix of the measured point j at the sample point (x, y, z), whose diagonal elements are obtained according to step S42.
Further, the distance dkjComprises the following steps:
Figure BDA0002596902090000054
wherein (x)k,yk,zk)、(xj,yj,zj) And respectively representing the three-dimensional space coordinates of the regression point k and the actual measurement point j.
Further, in step S5, the precision verification includes: and counting the judgment coefficients of the first water depth inversion model and the second water depth inversion model, and the root mean square error, the average absolute error and the average relative error of the inversion result, and performing model precision comparison analysis.
The technical scheme provided by the invention has the beneficial effects that: not only the heterogeneity in the horizontal direction is considered, but also the heterogeneity in the vertical direction of the sea area is considered, and the optimal bandwidth and the space weight matrix are calculated by using the space three-dimensional distance to replace the two-dimensional plane distance; meanwhile, in the sea area with large underwater topography change, the method of water depth zoning is used, so that when the water depth inversion model of the shallow water area or the deep water area is established, the corresponding actual measurement water depth data of the deep water area or the shallow water area does not participate in the establishment of the corresponding model, and the inversion accuracy of the model is improved.
Drawings
Fig. 1 is a flowchart of a multispectral remote sensing water depth inversion method based on an improved GWR model according to an embodiment of the present invention;
FIG. 2 is a WorldView-2 remote sensing image map provided by an embodiment of the present invention;
fig. 3 is a comparison graph of inversion results of the first water depth inversion model and the second water depth inversion model at the verification point according to the embodiment of the present invention;
fig. 4 is a comparison graph of inversion results of the first water depth inversion model and the second water depth inversion model at the verification point in the shallow water region according to the embodiment of the present invention;
fig. 5 is a comparison graph of inversion results of the first water depth inversion model and the second water depth inversion model at the verification point in the deep water region according to the embodiment of the invention
Fig. 6 is a final water depth inversion result diagram provided by the embodiment of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention will be further described with reference to the accompanying drawings.
Referring to fig. 1, an embodiment of the present invention provides a multispectral remote sensing water depth inversion method based on an improved GWR model, and the method is used for researching a ganquan island located in south China sea and an area near the ganquan island, and includes the following steps:
s1, multispectral remote sensing data preprocessing: the multispectral remote sensing data participating in water depth inversion is subjected to image preprocessing, a WorldView-2 remote sensing image is taken as an example (data shown in figure 2), and the image preprocessing comprises radiometric calibration, atmospheric correction, water and land separation, solar flare correction and the like.
Radiometric calibration is to convert the pixel value (Digital Number, DN) of the remote sensing image into the spectral radiance value of the top of the atmosphere, specifically, firstly, convert the pixel DN value into band-integrated radiance (L)0DN absCalFactor (unit W m)-2·sr-1) Then calculating the spectral radiance L ═ L0[ Delta ] lambda (unit W.m.)-2·sr-1·μm-1) Wherein absCalfactor represents an absolute scaling factor and Δ λ represents the effective width of the band (μm), both of which are available in the metadata file (. IMD) of the image.
After obtaining the radiance image of each wave band through radiometric calibration, obtaining remote sensing reflectivity data of each wave band by adopting a FLAASH atmospheric correction, 6S atmospheric correction or MODTRAN atmospheric correction method; next, in order to extract the water depth information more effectively and prevent the ground objects at the boundary between water and land from affecting the water depth inversion result, the non-sea area needs to be removed, in this embodiment, an improved normalized water body index (MNDWI) is used to extract the water body part image, where the improved normalized water body index is:
Figure BDA0002596902090000071
wherein Green represents the reflectivity of a Green waveband, and MIR represents the reflectivity of a mid-infrared waveband.
Since solar flare caused by sea surface waves and the like affects the real emissivity of the water body, solar flare removal processing needs to be performed on the reflectivity data of the sea surface, and the flare removal is performed by using a flare correction method based on a near-infrared band in the embodiment.
S2, acquiring actually measured water depth data: in this embodiment, a SHOAL-3000 airborne laser radar depth measurement system of Optech, canada, is used to measure the longitude and latitude coordinates, and the measured water depth data after tidal correction is corresponded to the multispectral remote sensing image through the position information to obtain a reflectivity data set on the multispectral remote sensing image corresponding to the measured point.
S3, constructing a first water depth inversion model: when the problem of spatial heterogeneity caused by substrate and water quality difference at different geographic positions is solved, compared with a global model, a GWR model is more feasible and effective, and can detect the spatial change of model coefficients, wherein the GWR model is the extension of a water depth inversion linear model, and the regression coefficients are obtained by performing local weighted regression on modeling points in a certain range around the inversion points and are different along with the change of the geographic positions; therefore, in this embodiment, a first water depth inversion model is constructed based on a GWR model, and the model formula is as follows:
Figure BDA0002596902090000081
in the formula (x)j,yj) Horizontal coordinate, Z, representing real measuring point jj(x, y) represents the water depth at the sample point (x, y) obtained by inversion from the actual point j, a0(xj,yj)、ai(xj,yj) Representing undetermined coefficients of actual measurement points j in the first water depth inversion model; xi＝ln[ρi-ρi∞]，ρiDenotes the reflectivity, p, of the i-th bandi∞The reflectivity of the wave band in a deep water area is shown, i is 1, …, n is n, and n is the number of the wave bands participating in inversion;j1an error term representing the first water depth inversion model;
and taking the logarithm value of the reflectivity on the multispectral image corresponding to the actual measurement point, the actual measurement point water depth and the actual measurement point position information as modeling data of the first water depth inversion model, taking the water depth of the position where each pixel (sample point) is located in the whole research area as a predicted value, solving the first water depth inversion model according to a weighted least square method, and obtaining the water depth inversion value of each pixel in the whole research area as an input value for improving the optimal bandwidth and weight matrix calculation in the GWR model. The solving process is as follows:
obtaining model coefficients according to the weighted least squares:
aj(x,y)＝(XTWj(x,y)X)-1XTWj(x,y)Z，
in the formula, aj(X, y) represents a regression coefficient of the actual measurement point j at the sample point (X, y), Z represents a vector composed of the actual measurement point water depth data, and X ═ X1,…,Xn]T，Wj(x, y) represents a weight matrix of the real measurement point j at the sample point (x, y) for ensuring that the real measurement point closer to the sample point has a larger weight, and the expression of the weight matrix is as follows:
Figure BDA0002596902090000082
in the formula, wj1,wj2,…,wjnRespectively representing the weight of the corresponding band observation point j.
S4, constructing a second water depth inversion model based on the improved GWR model: in order to consider the influence of vertical water depth information in parameter estimation, the vertical water depth is added into the calculation of a space weight matrix and an optimal bandwidth, so that a GWR model is improved, the vertical water depth information is used as a third dimension of space, namely, the distance between geographic space positions influencing parameter estimation values is improved from a two-dimensional plane distance to a three-dimensional space distance, and the problem of spatial heterogeneity in a spatial relationship is more comprehensively considered; in addition, in order to enable the vertical water depth information to play a better role in the calculation of the bandwidth and the weight matrix, the horizontal coordinate and the vertical water depth information need to be subjected to normalization processing, so that the problems caused by different spatial resolutions of remote sensing data, different distances of actually measured water depth sample points and regional differences are solved, and the applicability of the method is improved.
Specifically, the formula for obtaining the second water depth inversion model based on the improved GWR model is as follows:
Figure BDA0002596902090000091
wherein (x)j,yj,zj) Representing the spatial coordinates, z, of the actual measurement point jjRepresenting water depth data obtained at the actual measuring point j according to the first water depth inversion model; zj(x, y, z) represents the water depth of the sample point (x, y, z) obtained by inverting the second water depth inversion model according to the real-time point j; a is0(xj,yj,zj)、ai(xj,yj,zj) Representing undetermined coefficients of actual measurement points j in the second water depth inversion model;j2an error term representing the second water depth inversion model;
and taking the logarithm value of the reflectivity on the multispectral image corresponding to the actual measurement point, the actual measurement point water depth and the actual measurement point position information as modeling data of the second water depth inversion model, taking the logarithm value of the reflectivity of each wave band corresponding to each pixel of the whole research area, the water depth inversion value obtained based on the first water depth inversion model and the position information thereof as input data of the second water depth inversion model, meanwhile, carrying out normalization processing on three-dimensional coordinate data when calculating a bandwidth and a weight matrix, and finally carrying out coefficient calibration on pixel by pixel to realize water depth inversion.
The specific process of step S4 includes:
and S41, calculating the optimal bandwidth. This embodiment adopts Cross Validation (CV) method to estimate the optimal bandwidth:
Figure BDA0002596902090000092
where CV represents the sum of squares of the cross-validation residuals, N represents the total number of data for regression parameter calculations, b represents the bandwidth, and Z represents the sum of the squares of the cross-validation residualskRepresents the regression point (x)k,yk) The measured water depth value of the position (A),
Figure BDA0002596902090000105
numbers representing the periphery of regression points alone, excluding the regression points themselves, in the estimation of regression parametersCalculating a fitted water depth value according to the regression parameters; and drawing the different bandwidths and the CV values thereof into a trend line, wherein the bandwidth corresponding to the minimum CV value is the optimal bandwidth.
In the improved GWR model, the calculation of the optimal bandwidth also needs to perform normalization processing on the plane coordinate data and the vertical water depth data, so the CV calculation formula is as follows:
Figure BDA0002596902090000101
wherein scale (. cndot.) represents a normalization function,
Figure BDA0002596902090000102
representing fitted water depth values excluding the regression points themselves after normalization of the coordinate data, CV (scale (x), scale (y), scale (z)) representing the sum of squared cross-validation residuals after normalization of the coordinate data; and finding out the self-adaptive optimal bandwidth as the input of the spatial weight matrix by a cross validation method.
And S42, calculating the space weight. This example uses the bi-square function to compute the spatial weights:
Figure BDA0002596902090000103
in the formula (d)kjRepresents the distance, d, between the real point j and the regression point kbRepresenting the distance from the inversion point to the b-th nearest neighbor point, namely the adaptive bandwidth; within the bandwidth range, the closer the distance, the higher the weight, and the weight outside the bandwidth range is 0; in the second water depth inversion model, the distance between sample points is taken as a three-dimensional spatial distance, and thus:
Figure BDA0002596902090000104
in the formula (x)k,yk,zk)、(xj,yj,zj) Respectively representing a regression point k and an actual measurement pointSpatial coordinates of j, wherein zk、zjThe water depths at the regression point k and the actual measurement point j are respectively represented.
It should be noted that, because the magnitude of the difference between the horizontal coordinate and the vertical water depth is large, the vertical water depth is far smaller than the horizontal distance, and if the three-dimensional coordinate information is directly used for calculation, the influence of the vertical water depth information on the model coefficient is small. In order to provide the same important position for the horizontal coordinate and the vertical water depth in the weight matrix calculation, the two-dimensional coordinate in the horizontal direction and the vertical water depth value are normalized, such as the most value normalization, the standard score normalization, the median normalization, and the like, preferably, the present embodiment adopts the standard score normalization method, and therefore, the distance formula is:
Figure BDA0002596902090000111
wherein scale (·) represents a data normalization function.
And S43, estimating the model coefficient by using a weighted least square method to obtain a second water depth inversion model. Specifically, for any sample point (x, y, z) in the multispectral remote sensing image, the sample point is used as an inversion point, the optimal bandwidth and weight of each corresponding wave band are calculated by steps S41 and S42, a weight matrix is obtained, and the model coefficient is further determined as follows:
aj(x,y,z)＝(XTWj(x,y,z)X)-1XTWj(x,y,z)Z
in the formula, aj(X, y, z) represents a regression coefficient of the measured point j at the sample point (X, y, z), where X is [ X ═ X1,…,Xn]TZ represents a vector formed by actual-point water depth data, Wj(x, y, z) represents a weight matrix of the measured point j at the sample point (x, y, z), whose diagonal elements are obtained according to step S42.
S5, verifying the inversion accuracy of water depth: comparing and verifying the water depth inversion results of the first water depth inversion model and the second water depth inversion model by using the verification points, and counting the decision coefficient of the models and the root mean square error, the average absolute error and the average relative error of the inversion results; meanwhile, model precision comparison analysis is carried out on the statistical average absolute error and the average relative error of different water depth sections, and therefore precision verification is carried out. And outputting an inversion result of the second water depth inversion model after the verification is finished.
Specifically, referring to table 1, the determination coefficients (R) of the water depth inversion results of the OLS model (multiband model), the GWR model (first water depth inversion model), and the improved GWR model (second water depth inversion model) are counted2) The inversion accuracy of the three models is quantitatively evaluated according to Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Relative Error (MRE). From the overall modeling situation, the determination coefficient of each model is more than 0.80, wherein R of the improved GWR model2The maximum value is 0.991, which is 6.445% higher than that of a GWR model, and the result accuracy of the water depth inversion based on the improved GWR model is higher; from the estimation deviation situation, the RMSE of the three models is below 2m, wherein the RMSE of the improved GWR model is 0.627m, the RMSE is reduced by 56.848% compared with the GWR model, the MAE is 0.589m, the RMSE is reduced by 47.032% compared with the GWR model, the MRE is 29.299%, and the RMSE is reduced by 70.617% compared with the GWR model, which shows that the overall estimation error of the improved GWR model is smaller, and the result is better than that of the GWR model.
TABLE 1 comparison of inversion results for WorldView-2 data
Figure BDA0002596902090000121
Referring to fig. 3, a and b in fig. 3 respectively show the correlation between the inversion values and the measured values of the GWR model and the improved GWR model, and it can be seen that the scatter diagram of the improved GWR model is reduced by some points with larger deviation compared with the scatter diagram of the GWR model, and the correlation coefficient between the inversion water depth and the measured water depth obtained by improving the GWR model is 0.982, which is improved by 4.025% compared with the GWR model.
Whether analyzed from model fitting scenarios or estimation errors, the improved GWR model outperforms the GWR model.
In order to better evaluate the water depth inversion effect of the improved GWR model, please refer to table 2, and perform precision analysis on the model inversion results of different water depth regions. In the research area of the Ganquan island, the precision of the improved GWR model is improved in two different water depth ranges compared with the GWR model. For the shallow water region of the Ganquan island, the determination coefficients of the two models are both more than 0.8 from the overall modeling condition, and the determination coefficient of the improved GWR model is improved by 8.558% compared with the GWR model, which indicates that the inversion result of the improved GWR model is more excellent; from the estimation bias situation, the RMSE of the two models is below 0.4m, and the RMSE of the improved GWR model is reduced by 32.544 percent compared with the RMSE of the GWR model; the MAE of the improved GWR model was 0.199m, 24.335% lower compared to the GWR model, the MRE was 68.582% lower than the GWR model, 39.406% lower compared to the GWR model. In addition, for the deep water area of the Ganquan island, the inversion result shows that the situation is similar to that of the shallow water area, and the comparison result of the two models shows that the RMSE, MAE and MRE of the improved GWR model are respectively 0.540m, 0.596m and 8.310 percent and the performance is better than that of the GWR model no matter from the aspects of model establishment and estimated deviation.
TABLE 2 inversion results of water depth of two models at different water depth sections
Figure BDA0002596902090000122
Referring to fig. 4 and 5, portions a and b in fig. 4 represent the correlation between the inversion values and the actual values of the shallow water GWR model and the improved GWR model, respectively, and portions a and b in fig. 5 represent the correlation between the inversion values and the actual values of the deep water GWR model and the improved GWR model, respectively. As can be seen from the scatter diagram, the improved GWR model has smaller deviation of inversion results and higher correlation coefficient.
After the verification is completed, a final water depth inversion result of the second water depth inversion model is output, as shown in fig. 6.
According to another embodiment of the invention, the first water depth inversion model can also be a multiband model.
While the present invention has been described with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, which are illustrative and not restrictive, and it will be apparent to those skilled in the art that various changes and modifications can be made therein without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (8)

1. The multispectral remote sensing water depth inversion method based on the improved GWR model is characterized by comprising the following steps:
s1, performing image preprocessing on the multispectral remote sensing image of the research area to obtain water surface reflectivity data;
s2, acquiring actually measured water depth data, and corresponding the actually measured water depth data to the reflectivity data of the multispectral remote sensing image according to the position information;
s3, constructing a first water depth inversion model, and obtaining water depth data of the whole research area by using the first water depth inversion model;
s4, constructing a second water depth inversion model based on the improved GWR model:
Figure FDA0002596902080000011
in the formula (x)j,yj,zj) Represents the spatial coordinates of the real measuring point j, where zjRepresenting water depth data obtained at the actual measuring point j according to the first water depth inversion model; zj(x, y, z) represents the water depth of the sample point (x, y, z) obtained by inverting the second water depth inversion model according to the real-time point j; a is0(xj,yj,zj)、ai(xj,yj,zj) Representing undetermined coefficients of actual measurement points j in the second water depth inversion model; xi＝ln[ρi-ρi∞]，ρiDenotes the reflectivity, p, of the i-th bandi∞The reflectivity of the wave band in a deep water area is shown, i is 1, …, n is n, and n is the number of the wave bands participating in inversion;j2an error term representing the second water depth inversion model;
taking the logarithm value of the reflectivity on the multispectral image corresponding to the actual measurement point, the actual measurement point water depth and the actual measurement point position information as modeling data of the second water depth inversion model, taking the logarithm value of the reflectivity of each wave band corresponding to each sample point in the research area, the water depth data obtained based on the first water depth inversion model and the corresponding position information as input data of the second water depth inversion model, and performing pixel-by-pixel coefficient calibration to realize water depth inversion;
and S5, comparing inversion results of the first water depth inversion model and the second water depth inversion model by using the verification points, so as to perform precision verification, and outputting the water depth value obtained by the second water depth inversion model after the verification is completed.
2. The method for multispectral remote sensing water depth inversion based on an improved GWR model as claimed in claim 1, wherein in step S3, a first water depth inversion model is constructed based on the GWR model:
Figure FDA0002596902080000021
in the formula (x)j,yj) Horizontal coordinate, Z, representing real measuring point jj(x, y) represents the water depth at the sample point (x, y) obtained by inversion from the actual point j, a0(xj,yj)、ai(xj,yj) Representing undetermined coefficients of actual measurement points j in the first water depth inversion model;j1an error term representing the first water depth inversion model;
and taking the logarithmic value of the reflectivity of each wave band of the actually measured water depth data and the position information thereof as modeling data of the first water depth inversion model, and performing water depth inversion on the whole research area according to a weighted least square method to obtain water depth data corresponding to each sample point in the research area.
3. The method for multispectral remote water depth inversion based on an improved GWR model as claimed in claim 1, wherein in step S1, the image preprocessing comprises:
carrying out radiometric calibration on the multispectral remote sensing image, and converting a pixel value DN into a spectral radiance value L at the top of the atmosphere, wherein the pixel value L is DN, absCalfactor/delta lambda, wherein absCalfactor represents an absolute calibration factor, and delta lambda represents the effective width of a wave band;
atmospheric correction is carried out on the radiance image of each wave band obtained by radiometric calibration, and reflectivity data of each wave band is obtained;
and extracting a water body part image from the reflectivity data by adopting an improved normalized water body index, and finally removing solar flare to finish image preprocessing.
4. The method for multispectral remote sensing water depth inversion based on an improved GWR model as claimed in claim 1, wherein the specific process of performing water depth inversion on the whole research area according to the weighted least squares method in step S3 is as follows: for any sample point (x, y) in the research area, coefficient calibration is carried out by using a weighted least square method to obtain a model coefficient:
aj(x,y)＝(XTWj(x,y)X)-1XTWj(x,y)Z
in the formula, aj(X, y) represents a regression coefficient of the actual measurement point j at the sample point (X, y), where X ═ X1,…,Xn]TZ represents a vector formed by actual-point water depth data, Wj(x, y) represents the weight matrix of the real measurement point j at sample point (x, y).
5. The improved GWR model-based multispectral remote water depth inversion method as claimed in claim 4, wherein the weight matrix W isjThe diagonal elements of (x, y) are a function of the coordinates (x, y) of the sample points to ensure that real points closer to the sample points have greater weight.
6. The method for multispectral remote sensing water depth inversion based on an improved GWR model as claimed in claim 1, wherein in step S4, the specific process for performing water depth inversion is as follows:
s41, carrying out normalization processing on the plane coordinate data and the vertical water depth data, and estimating the optimal bandwidth by adopting a cross validation method:
Figure FDA0002596902080000031
wherein CV represents the sum of squares of the cross-validation residuals, (x, y, Z) represents the three-dimensional space coordinates of the regression point k, scale (. cndot.) represents the normalization function, N represents the total number of data calculated by the regression parameters, b represents the bandwidth, and Z represents the bandwidthkRepresents the measured water depth value at the regression point k,
Figure FDA0002596902080000033
representing a fitted water depth value which does not include regression points to carry out regression parameter calculation per se; the bandwidth corresponding to the minimum CV value is the optimal bandwidth;
and S42, calculating the space weight by adopting a bi-square function according to the optimal bandwidth b:
Figure FDA0002596902080000032
in the formula (d)kjRepresents the distance, d, between the real point j and the regression point kbRepresenting the distance of the inversion point to the b-th nearest neighbor;
s43, taking any sample point (x, y, z) in the multispectral remote sensing image as an inversion point, calculating the optimal bandwidth and weight of each corresponding wave band by utilizing the steps S41 and S42 to obtain a weight matrix, and further estimating the model coefficient by adopting a weighted least square method:
aj(x,y,z)＝(XTWj(x,y,z)X)-1XTWj(x,y,z)Z
in the formula, aj(X, y, z) represents a regression coefficient of the measured point j at the sample point (X, y, z), where X is [ X ═ X1,…,Xn]TZ represents the direction of water depth data composition of actual measurement pointsAmount, Wj(x, y, z) represents a weight matrix of the measured point j at the sample point (x, y, z), whose diagonal elements are obtained according to step S42.
7. The improved GWR model-based multi-spectral remote water depth inversion method as claimed in claim 6, wherein said distance dkjComprises the following steps:
Figure FDA0002596902080000041
wherein (x)k,yk,zk)、(xj,yj,zj) And respectively representing the three-dimensional space coordinates of the regression point k and the actual measurement point j.
8. The method for multispectral remote water depth inversion based on an improved GWR model as claimed in claim 1, wherein in step S5, the accuracy verification comprises:
and counting the judgment coefficients of the first water depth inversion model and the second water depth inversion model, and the root mean square error, the average absolute error and the average relative error of the inversion result, and performing model precision comparison analysis.
CN202010711999.6A 2020-07-22 2020-07-22 Multispectral remote sensing water depth inversion method based on improved GWR model Pending CN112013822A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010711999.6A CN112013822A (en) 2020-07-22 2020-07-22 Multispectral remote sensing water depth inversion method based on improved GWR model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010711999.6A CN112013822A (en) 2020-07-22 2020-07-22 Multispectral remote sensing water depth inversion method based on improved GWR model

Publications (1)

Publication Number Publication Date
CN112013822A true CN112013822A (en) 2020-12-01

Family

ID=73498852

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010711999.6A Pending CN112013822A (en) 2020-07-22 2020-07-22 Multispectral remote sensing water depth inversion method based on improved GWR model

Country Status (1)

Country Link
CN (1) CN112013822A (en)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113326470A (en) * 2021-04-11 2021-08-31 桂林理工大学 Remote sensing water depth inversion tidal height correction method
CN113639716A (en) * 2021-07-29 2021-11-12 北京航空航天大学 Depth residual shrinkage network-based water depth remote sensing inversion method
CN114459438A (en) * 2022-01-10 2022-05-10 山东科技大学 Method for judging validity of high-resolution multispectral water depth inversion data based on spectral roughness information
CN114758254A (en) * 2022-06-15 2022-07-15 中国地质大学（武汉） Dual-band unsupervised water depth inversion method and system
CN115060656A (en) * 2022-05-11 2022-09-16 自然资源部第一海洋研究所 Satellite remote sensing water depth inversion method based on sparse prior actual measurement points
CN115371642A (en) * 2022-08-22 2022-11-22 河海大学 Geological statistics inversion method for hydrogeological parameter dynamic change characteristics
CN117789056A (en) * 2024-02-27 2024-03-29 杭州蚁联传感科技有限公司 Remote sensing data processing method and device with solar flare and storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105445751A (en) * 2015-11-18 2016-03-30 国家***第一海洋研究所 Shallow water area water depth ratio remote sensing inversion method
CN105469393A (en) * 2015-11-18 2016-04-06 国家***第一海洋研究所 Shallow water depth multi-temporal remote sensing image inversion method based on decision fusion
CN105627997A (en) * 2015-12-23 2016-06-01 国家***第一海洋研究所 Multi-angle remote sensing water depth decision fusion inversion method
CN105651263A (en) * 2015-12-23 2016-06-08 国家***第海洋研究所 Shallow sea water depth multi-source remote sensing fusion inversion method
CN105865424A (en) * 2016-04-13 2016-08-17 中测新图(北京)遥感技术有限责任公司 Nonlinear model-based multispectral remote sensing water depth inversion method and apparatus thereof
CN109543356A (en) * 2019-01-07 2019-03-29 福州大学 Consider the ocean interior temperature-salinity structure remote sensing inversion method of Space atmosphere
CN110046771A (en) * 2019-04-25 2019-07-23 河南工业大学 A kind of PM2.5 concentration prediction method and apparatus

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105445751A (en) * 2015-11-18 2016-03-30 国家***第一海洋研究所 Shallow water area water depth ratio remote sensing inversion method
CN105469393A (en) * 2015-11-18 2016-04-06 国家***第一海洋研究所 Shallow water depth multi-temporal remote sensing image inversion method based on decision fusion
CN105627997A (en) * 2015-12-23 2016-06-01 国家***第一海洋研究所 Multi-angle remote sensing water depth decision fusion inversion method
CN105651263A (en) * 2015-12-23 2016-06-08 国家***第海洋研究所 Shallow sea water depth multi-source remote sensing fusion inversion method
CN105865424A (en) * 2016-04-13 2016-08-17 中测新图(北京)遥感技术有限责任公司 Nonlinear model-based multispectral remote sensing water depth inversion method and apparatus thereof
CN109543356A (en) * 2019-01-07 2019-03-29 福州大学 Consider the ocean interior temperature-salinity structure remote sensing inversion method of Space atmosphere
CN110046771A (en) * 2019-04-25 2019-07-23 河南工业大学 A kind of PM2.5 concentration prediction method and apparatus

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
ANDRZEJ CHYBICKI: "Three-Dimensional Geographically Weighted Inverse Regression (3GWR) Model for Satellite Derived Bathymetry Using Sentinel-2 Observations", 《MARINE GEODESY》 *
刘源 等: "基于地理加权回归模型的水深反演算法研究", 《海洋学研究》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113326470A (en) * 2021-04-11 2021-08-31 桂林理工大学 Remote sensing water depth inversion tidal height correction method
CN113326470B (en) * 2021-04-11 2022-08-16 桂林理工大学 Remote sensing water depth inversion tidal height correction method
CN113639716A (en) * 2021-07-29 2021-11-12 北京航空航天大学 Depth residual shrinkage network-based water depth remote sensing inversion method
CN114459438A (en) * 2022-01-10 2022-05-10 山东科技大学 Method for judging validity of high-resolution multispectral water depth inversion data based on spectral roughness information
CN114459438B (en) * 2022-01-10 2024-02-02 山东科技大学 Method for judging validity of high-resolution multispectral water depth inversion data
CN115060656A (en) * 2022-05-11 2022-09-16 自然资源部第一海洋研究所 Satellite remote sensing water depth inversion method based on sparse prior actual measurement points
CN115060656B (en) * 2022-05-11 2024-06-07 自然资源部第一海洋研究所 Satellite remote sensing water depth inversion method based on sparse priori real measurement points
CN114758254A (en) * 2022-06-15 2022-07-15 中国地质大学（武汉） Dual-band unsupervised water depth inversion method and system
CN115371642A (en) * 2022-08-22 2022-11-22 河海大学 Geological statistics inversion method for hydrogeological parameter dynamic change characteristics
CN117789056A (en) * 2024-02-27 2024-03-29 杭州蚁联传感科技有限公司 Remote sensing data processing method and device with solar flare and storage medium
CN117789056B (en) * 2024-02-27 2024-05-07 杭州蚁联传感科技有限公司 Remote sensing data processing method and device with solar flare and storage medium

Similar Documents

Publication Publication Date Title
CN112013822A (en) Multispectral remote sensing water depth inversion method based on improved GWR model
CN111583214B (en) Sea surface wind speed inversion method based on RBF neural network and based on marine radar image
Huang et al. Discharge estimation in high-mountain regions with improved methods using multisource remote sensing: A case study of the Upper Brahmaputra River
CN109059796B (en) Shallow sea water depth multispectral satellite remote sensing inversion method for water depth control point-free area
Ai et al. Convolutional neural network to retrieve water depth in marine shallow water area from remote sensing images
Su et al. Prediction of water depth from multispectral satellite imagery—the regression Kriging alternative
CN108303044B (en) Leaf area index obtaining method and system
Xu et al. Deriving highly accurate shallow water bathymetry from Sentinel-2 and ICESat-2 datasets by a multitemporal stacking method
Kanno et al. Shallow water bathymetry from multispectral satellite images: Extensions of Lyzenga's method for improving accuracy
Jawak et al. Spectral information analysis for the semiautomatic derivation of shallow lake bathymetry using high-resolution multispectral imagery: A case study of Antarctic coastal oasis
CN114201732A (en) Sentinel-2A image-based shallow sea water depth inversion method
CN113639716A (en) Depth residual shrinkage network-based water depth remote sensing inversion method
CN114297938A (en) Optical shallow water bottom depth inversion method based on neural network
Coleman et al. Holes in the ocean: Filling voids in bathymetric lidar data
CN111561916B (en) Shallow sea water depth uncontrolled extraction method based on four-waveband multispectral remote sensing image
CN112381033A (en) Satellite remote sensing quantitative evaluation method for intertidal zone erosion and deposition conditions caused by offshore wind turbine
Lumban-Gaol et al. Satellite-derived bathymetry using convolutional neural networks and multispectral sentinel-2 images
Berezowski et al. Flooding extent mapping for synthetic aperture radar time series using river gauge observations
CN117274831B (en) Offshore turbid water body depth inversion method based on machine learning and hyperspectral satellite remote sensing image
Albanai Mapping Kuwait bathymetry using passive multispectral remote sensing
CN116817869B (en) Submarine photon signal determination method using laser radar data
Ashphaq et al. Analysis of univariate linear, robust-linear, and non-linear machine learning algorithms for satellite-derived bathymetry in complex coastal terrain
Kao et al. Determination of shallow water depth using optical satellite images
CN115060656B (en) Satellite remote sensing water depth inversion method based on sparse priori real measurement points
CN111650128A (en) High-resolution atmospheric aerosol inversion method based on surface reflectivity library

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