A method for inverting significant wave height from marine radar images

Abstract

The present invention discloses a method for inverting effective wave height from a marine radar image. First, offline buoy data and radar are used to obtain the relationship between the wave crest period and the average zero-crossing period of the wave by least square fitting, and then a two-dimensional wave spectrum is obtained by 3D-FFT spectrum analysis. The wave crest period of the wave is obtained from the two-dimensional wave spectrum, and then the average zero-crossing period is calculated. The root mean square steepness of the wave is calculated by the shadow method, and finally the effective wave height of the wave is calculated by the average zero-crossing period and the root mean square steepness parameter. The present invention can independently invert the effective wave height parameters of the wave through radar data, so that the process of measuring the effective wave height parameters of the wave by radar is completely separated from the support of the buoy, and does not need to rely on other instruments, thereby increasing the applicability and independence of the algorithm.

Description

Technical Field

[0001] The invention belongs to the field of inverting sea wave parameters by remote sensing means, and relates to a method for inverting effective wave height from a marine radar image, in particular to a method for inverting effective wave height from a marine radar image based on a combination strategy. Background Art

[0002] In recent years, the use of X-band marine radar to perform real-time remote sensing of ocean waves has become a research hotspot in the field of ocean monitoring. Currently, the two main methods for inverting significant wave height using X-band marine radar images are spectral analysis and shadow method.

[0003] The spectral analysis method is an algorithm that comprehensively utilizes the time and space information of sea clutter images. This method performs FFT transformation, filtering, integration and other operations on continuous sea clutter images to ultimately obtain the frequency directional spectrum and signal-to-noise ratio of the waves. The wave direction and frequency information can be further obtained from the frequency directional spectrum of the waves. The effective wave height information can be calculated through the linear relationship between the signal-to-noise ratio and the effective wave height.

[0004] The advantage of the spectrum analysis method is that it can obtain multiple parameters of the waves at the same time - wave direction, frequency, and effective wave height. The influence of occlusion on the spectrum analysis method can be ignored, which is an advantage that many methods cannot match. However, the coefficient between the signal-to-noise ratio and the effective wave height in the spectrum analysis method needs to be determined through buoy experiments, which increases the engineering application cost of the spectrum analysis method. In addition, the spectrum method requires multiple consecutive images, and the algorithm calculation time is long, which also limits the application of this method.

[0005] In 2014, Norwegian oceanographer Rune

[13] summarized previous research and clearly pointed out that sea clutter is mainly generated by shadow modulation, and proposed a shadow method that does not require calibration parameters to calculate effective wave height. Compared with the traditional spectral analysis method, the biggest advantage of the shadow method is that the calculation of the shadow method is completely based on theoretical derivation and does not require buoy experiments to calibrate the relevant coefficients. In addition, the shadow method only requires one radar image for calculation, and the algorithm has low time cost, which enhances the applicability of the algorithm. However, in the current research process, the average zero-crossing period of the waves in the shadow method is provided by the buoy, which limits the application and development of the shadow method. Summary of the invention

[0006] In view of the above-mentioned prior art, the technical problem to be solved by the present invention is to provide a method for inverting effective wave height from marine radar images based on a combination strategy of 3D-FFT spectral analysis and shadow method, so as to solve the problems of long algorithm calculation time and high engineering application cost of 3D-FFT spectral analysis method and dependence of wave height parameter calculation on buoy data in shadow statistics method.

[0007] In order to solve the above technical problems, a method for inverting significant wave height from a marine radar image of the present invention comprises the following steps:

[0008] Step 1: Using offline buoy data and radar, the relationship between the wave crest period and the average zero-crossing period is obtained by least squares fitting;

[0009] Step 2: Select a research area, obtain the two-dimensional wave spectrum of the selected research area through 3D-FFT spectrum analysis method, and obtain the wave crest period of the wave from the two-dimensional wave spectrum;

[0010] Step 3: Substitute the wave crest period obtained in step 2 into the relationship between the wave crest period and the average zero-crossing period obtained in step 1 to calculate the average zero-crossing period T m02 ;

[0011] Step 4: Select a research area in the original radar image, divide the original radar image into m sub-areas in the radial dimension in the selected research area; calculate the threshold τ to distinguish between shadow areas and non-shadow areas S ; According to τ S Calculate the shadow ratio function of each sub-area; then calculate the root mean square wave steepness σ according to the shadow ratio function RMS ;

[0012] Step 5: By averaging the zero-crossing period T m02 and the root mean square steepness σ RMS Calculate the effective wave height H of the ocean wave S , H S satisfy:

[0013]

[0014] Where g is the acceleration due to gravity.

[0015] Furthermore, the relationship between the wave crest period and the average zero-crossing period of the ocean wave is specifically:

[0016] T m02 =aT M +b

[0017] Where, T M is the peak period, T m02 is the average zero-crossing period;

[0018] The average zero-crossing period sequence of the buoy in the offline experiment is T m02i , i = 1, 2, 3...n, and the corresponding peak period sequence obtained by 3D-FFT spectrum analysis using radar data is recorded as T Mi , i=1,2,3...n, according to the least squares method, we get:

[0019]

[0020]

[0021] In the formula, and are the average values ​​of the average zero-crossing period sequence and the peak period sequence respectively.

[0022] Furthermore, step 2 includes:

[0023] Step 2.1: Through the nearest point interpolation, the echo intensity value in the polar coordinates of the original radar image is converted to the corresponding point in the Cartesian coordinate system. The coordinates of a point in the original radar image are recorded as (r0, θ0), and the corresponding Cartesian coordinates after conversion are recorded as (x0, y0). The radar image sequence obtained after interpolation is recorded as η(x, y, t). The conversion formula is:

[0024]

[0025] Where r0 is the radial distance of a point in the original radar image, θ0 is the azimuth of a point in the original radar image, x0 and y0 are the horizontal and vertical coordinates of the Cartesian coordinate system, sqrt(), round() and rem() are the square root function, integer function and remainder function respectively;

[0026] Step 2.2: The size of the selected research area is L x *L y The Cartesian frame, L x ,L y is the range of the Cartesian frame along the x and y directions, the total time of the radar image sequence is T, and the radar image sequence η(x, y, t) is transformed into a three-dimensional image spectrum by a three-dimensional Fourier transform. Then, the image sequence is transformed from the time-space domain to the wavenumber frequency domain using the discrete form of the transformation formula:

[0027]

[0028] In the formula, k x and k y are the wave numbers in the x and y directions respectively, and ω is the angular frequency. The resolution is calculated as follows:

[0029]

[0030] Step 2.3: Assume that the theoretical frequency corresponding to the wave number is ω0, and the real frequency ω is expressed by the dispersion relation, then the sum of the squares of the interpolation between the two is:

[0031]

[0032] Taking the energy value in the three-dimensional image spectrum as the weight, the objective function is:

[0033]

[0034] When the partial derivatives of SSE with respect to are all zero, SSE has a minimum solution, and the flow velocity is calculated as:

[0035]

[0036] Where h is the water depth, g is the gravitational acceleration, and F is the three-dimensional image spectrum F(k x ,k y ,ω), sea surface laminar flow Wave number Among them, k x ,k y are the wave numbers in the x and y directions respectively;

[0037] Step 2.4: Filter the three-dimensional image spectrum using a bandpass filter based on the dispersion relation. The filter used is a bandpass filter based on the maximum flow velocity:

[0038]

[0039] In the formula, E(k x ,k y ,ω) is the three-dimensional image spectrum obtained after filtering, where

[0040]

[0041]

[0042] Where k is the wave number, dω is the frequency resolution, dk = (dk x =dk y ) is the wave number resolution, U max is the maximum flow rate, when B n When the value is negative, let B n =0;

[0043] Then the three-dimensional image spectrum E(k x ,k y ,ω) and perform frequency ω-dimensional integration to obtain the two-dimensional image spectrum I(k x ,k y ):

[0044] I(k x ,k y )=∫ ω＞0 E(k x ,k y ,ω)dω

[0045] Step 2.5: Use the empirical modulation transfer function to transform the two-dimensional image spectrum I(k x ,k y ) is converted into a two-dimensional wave spectrum E(k x ,k y ):

[0046] E(k x ,k y )=|M(k x ,k y )| 2 ·I(k x ,k y )

[0047] The empirical modulation transfer function is:

[0048] |M(k x ,k y )| 2 ≈k -β

[0049] Among them, β is the empirical coefficient;

[0050] Step 2.6: Through the two-dimensional wave spectrum E(k x ,k y ) Using the dispersion relation to calculate the peak period T of the ocean wave M , the calculation formula is:

[0051]

[0052] Where g is the acceleration due to gravity, k is m is the wave spectrum E(k x ,k y ) is the wave number at the point where the energy value is the largest, and the calculation formula is:

[0053]

[0054] Among them, k xm , k ym are the horizontal and vertical coordinates of the point with the maximum wave number.

[0055] Furthermore, in step 4, the threshold τ for distinguishing the shadow area from the non-shadow area is calculated. S include:

[0056] Step 4.2.1: Image edge detection: by performing convolution operations on the original radar image in t adjacent different directions based on the differential operator, the edge gradient image I is obtained. Ei (r,θ):

[0057]

[0058] Where I(r,θ) is the original radar image, r is the distance, and θ is the azimuth. represents the convolution operation, I Ei (r,θ) is the edge gradient image obtained after convolution operation with different difference operators, H i In (r,θ), i=1,2,3...t represents the pixel difference operator in t different directions;

[0059] Step 4.2.2: Subtract the difference image I Ei The pixel points with the intensity values ​​in the first N% of (r,θ) are taken as the boundary pixels of the shadow area and the non-shadow area to obtain the sub-edge image I Ti (r,θ), where i=1,2,3...t, the specific formula is:

[0060]

[0061] Step 4.2.3: I Ti (r,θ) are superimposed to obtain the original edge image I T (r,θ), where i = 1, 2, 3...t:

[0062]

[0063] Then use the edge comparison method to transform the original edge image I T Each point in (r,θ) is related to the threshold τ F Compare to remove isolated noise points and get edge image I F (r,θ):

[0064]

[0065] Step 4.2.4: For each isolated 1 in the edge image, find the pixel η of the original radar image at that position, and establish an image histogram function F through the correspondence between the edge image and the original radar image. H (η), for F H (η) Taking the mode, we can get the threshold τ to distinguish between shadow area and non-shadow area S :

[0066] τ S =mode(F H (η)).

[0067] Furthermore, in step 4, according to τ S The shadow ratio function for each sub-area is calculated as follows:

[0068] Step 4.3.1: Using threshold τ S The original radar image I T(r,θ) is divided into shadow area and non-shadow area, and the shadow image I s (r,θ) is:

[0069]

[0070] Among them, I(r,θ) is the original radar image, I s (r,θ) is the shadow image after segmentation;

[0071] Step 4.3.2: Calculate the grazing angle γ of each sub-region i , the grazing angle corresponding to the center point of the sub-region is the grazing angle of the sub-region. Assume that the radar antenna height is H, the radar image distance resolution is R, the number of reference points at the starting position is h0, and the number of relative points at the start of the i-th sub-region is p i , the relative number of points at the end is q i , where i = 1, 2, 3...m, then the calculation formula for the grazing angle of the i-th sub-area is:

[0072]

[0073] Step 4.3.3: Calculate the shadow ratio function S(γ i ), where i = 1, 2, 3...m, then the calculation formula of the shadow ratio function of the i-th sub-area is:

[0074]

[0075] Where P i Represents the shadow image I corresponding to the i-th sub-region s The number of pixels with a pixel value of 1 in (r,θ), Q i Represents the shadow image I corresponding to the i-th sub-region s The number of all pixels in (r,θ).

[0076] Furthermore, the root mean square steepness calculated according to the shadow ratio function in step 4 includes:

[0077] Step 4.4.1: Shadow ratio function S(γ i ) and lighting ratio L(γ i )satisfy:

[0078] L(γ i )=1-S(γ i )

[0079] The least squares method is used to fit the illumination ratio L(γ i ) and the grazing angle γ i The relationship between , where i = 1, 2, 3 ... m, is used to estimate the root mean square wave steepness σRMS , the Smith function is:

[0080]

[0081] In the formula,

[0082]

[0083] μ i =tan(γ i )

[0084] Among them, erfc() is the Gaussian error function. The erfc() function is derived from the Gaussian surface probability density function. The expression of the Gaussian surface probability density function is:

[0085]

[0086] Where ξ is the height of the Gaussian sea level from the horizontal plane, υ is the root mean square of the height deviation, and Δξ is the height increment of the Gaussian surface when Δτ is added in the horizontal axis direction.

[0087] Furthermore, the marine radar is an X-band marine radar.

[0088] The beneficial effects of the present invention are as follows: the present invention designs an effective wave height combined inversion algorithm based on marine radar images, especially X-band marine radar images, which combines 3D-FFT spectrum analysis method and shadow statistics method. This algorithm makes up for the shortcomings of long calculation time and high engineering application cost of the 3D-FFT spectrum analysis method and the dependence of wave height parameter calculation in the shadow statistics method on buoy data, i.e., the need for external reference equipment for calibration, thereby increasing the applicability and independence of the algorithm.

[0089] Compared with the existing technology, its advantages are:

[0090] (1) The present invention uses the shadow method to calculate the average steepness of the waves, and calculates the effective wave height of the waves by the average zero-crossing period and the average steepness parameter. It does not require buoy experiments to calibrate related coefficients. The algorithm has low time cost, which can reduce R&D costs and shorten the R&D cycle.

[0091] (2) The present invention combines the spectral analysis method and the shadow statistics method to independently invert the effective wave height parameters of the waves through radar data, so that the process of measuring the effective wave height parameters of the waves by radar is completely independent of the support of the buoy and no longer requires the help of other instruments, thereby improving the engineering applicability of the algorithm and reducing the additional engineering cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] Figure 1 It is the fitting diagram of radar peak period and buoy average zero-crossing period;

[0093] Figure 2 This is a schematic diagram of the study area selection of 3D-FFT spectrum analysis method;

[0094] Figure 3 It is a schematic diagram of Cartesian box selection of 3D-FFT spectrum analysis method;

[0095] Figure 4 It is a flow chart of obtaining ocean wave parameters using 3D-FFT spectrum analysis method;

[0096] Figure 5 is the two-dimensional wave spectrum corresponding to a significant wave height of 2.0 m;

[0097] Figure 6 It is a line chart of average zero-crossing period;

[0098] Figure 7 This is a schematic diagram of region segmentation using the shadow method;

[0099] Figure 8 It is a line chart of effective wave height;

[0100] Figures 9(a) and 9(b) are comparisons of the radar-inverted effective wave height and the buoy effective wave height;

[0101] Fig.10 It is a flow chart of the present invention. DETAILED DESCRIPTION

[0102] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0103] Combination Fig.10 , the present invention comprises the following steps:

[0104] Step 1, offline data processing, obtaining the relationship between the wave crest period and the average zero-crossing period of the wave by least squares fitting, including the following steps:

[0105] Step 1.1, using offline buoy data and radar, the relationship between the wave crest period and the average zero-crossing period is obtained by least squares fitting. Assume that there is a linear relationship between the two, and record the corresponding linear relationship as:

[0106] T m02 =aT M +b

[0107] Where, T M is the peak period, T m02 is the average zero-crossing period.

[0108] The average zero-crossing period sequence of the buoy in the offline experiment is T m02i, i = 1, 2, 3...n, and the corresponding peak period sequence obtained by 3D-FFT spectrum analysis using radar data is recorded as T Mi , i=1,2,3...n, the relevant knowledge of least squares method is:

[0109]

[0110]

[0111] In the formula, and are the average values ​​of the average zero-crossing period sequence and the peak period sequence respectively.

[0112] The most commonly used parameter for evaluating linear relationships is the correlation coefficient r. The value range of r is [0,1], which reflects the closeness of the linear relationship between the two sequences. The closer r is to 1, the better the fitting effect is, and vice versa. The definition of r is:

[0113]

[0114] Step 2, obtaining a two-dimensional wave spectrum by 3D-FFT spectrum analysis method, and obtaining the wave crest period of the wave from the two-dimensional wave spectrum, including the following steps:

[0115] Step 2.1, select the study area size as L x *L y , the spatial resolution of the radar is R, so the pixel size of the study area is M x *M y .

[0116] Step 2.2, obtain a two-dimensional wave spectrum by 3D-FFT spectrum analysis method, and obtain the wave crest period of the wave from the two-dimensional wave spectrum. The 3D-FFT spectrum analysis method is divided into image preprocessing, 3D-FFT, dispersion relation estimation, bandpass filtering, integration, empirical modulation transfer function, calculation of wave crest period and other steps, and the step 2.2 includes the following steps:

[0117] Step 2.2.1, image preprocessing: Through the nearest point interpolation, the echo intensity value in the polar coordinates of the original radar image is converted to the corresponding point in the Cartesian coordinate system. The coordinates of a point in the original radar image are recorded as (r0, θ0), and the corresponding Cartesian coordinates after conversion are recorded as (x0, y0). The radar image sequence obtained after interpolation is recorded as η(x, y, t), and the conversion formula is:

[0118]

[0119] Where r0 is the radial distance of a point in the original radar image, θ0 is the azimuth of a point in the original radar image, x0 and y0 are the abscissa and ordinate values ​​of the Cartesian coordinate system, sqrt(), round() and rem() are the square root function, integer function and remainder function, respectively.

[0120] Step 2.2.2, perform a three-dimensional Fourier transform on the radar image sequence η(x, y, t) interpolated in step 2.2.1 to obtain a three-dimensional image spectrum F(k x ,k y ,ω), the transformation formula is:

[0121]

[0122] The size of the selected research area is L x *L y The Cartesian frame, L x ,L y is the range of the Cartesian frame along the x and y directions, and the total time of the radar image sequence is T, then the above transformation formula can be written as:

[0123]

[0124] The discrete form of the transformation formula is used to transform the image sequence from the time-space domain to the wavenumber frequency domain. Its discrete form is:

[0125]

[0126] In the formula, k x and k y are the wave numbers in the x and y directions respectively, and ω is the angular frequency. The resolution is calculated as follows:

[0127]

[0128] Where, L x ,L y is the extent of the Cartesian box selected for study along the x and y directions, and T is the total time of the radar image sequence.

[0129] Step 2.2.3, dispersion relation estimation: Assuming that the theoretical frequency corresponding to the wave number is ω0, and the real frequency ω is expressed by the dispersion relation, the sum of the squares of the interpolation between the two can be expressed as:

[0130]

[0131] The energy value in the three-dimensional image spectrum is used as the weight, so the objective function is:

[0132]

[0133] When the partial derivatives of SSE with respect to are all zero, SSE has a minimum solution:

[0134]

[0135]

[0136] The flow velocity can then be calculated as:

[0137]

[0138] Where h is the water depth, g is the gravitational acceleration, and F is the three-dimensional image spectrum F(k x ,k y ,ω), sea surface laminar flow Wave number Among them, k x ,k y are the wave numbers in the x and y directions respectively.

[0139] Step 2.2.4, bandpass filtering and integration: The three-dimensional image spectrum obtained in step 2.2.2 is filtered using a bandpass filter based on the dispersion relation. The filter used is a bandpass filter based on the maximum flow velocity:

[0140]

[0141] In the formula, E(k x ,k y ,ω) is the three-dimensional image spectrum obtained after filtering, F(k x ,k y ,ω) is the three-dimensional image spectrum obtained by step 2.2.2, where

[0142]

[0143]

[0144] Where k is the wave number, dω is the frequency resolution, dk = (dk x =dk y ) is the wave number resolution, U max is the maximum flow rate, when B n When the value is negative, let B n =0.

[0145] in,

[0146]

[0147] Where, L x ,L yis the extent of the Cartesian box selected for study along the x and y directions, and T is the total time of the radar image sequence.

[0148] Then the three-dimensional image spectrum E(k x ,k y ,ω) and perform frequency ω-dimensional integration to obtain the two-dimensional image spectrum I(k x ,k y ), the integral formula is:

[0149] I(k x ,k y )=∫ ω＞0 E(k x ,k y ,ω)dω

[0150] Step 2.2.5, empirical modulation transfer function: The two-dimensional image spectrum I(k x ,k y ) is converted into a two-dimensional wave spectrum E(k x ,k y ), the conversion formula is:

[0151] E(k x ,k y )=|M(k x ,k y )| 2 ·I(k x ,k y )

[0152] The empirical modulation transfer function is obtained by fitting experimental data. The commonly used empirical modulation transfer function is in the form of:

[0153] |M(k x ,k y )| 2 ≈k -β

[0154] Among them, β is the empirical coefficient.

[0155] Step 2.2.6, calculate the wave crest period: Use the two-dimensional wave spectrum E(k x ,k y ) Using the dispersion relation to calculate the peak period T of the ocean wave M , the calculation formula is:

[0156]

[0157] Where g is the acceleration due to gravity, k is m is the wave spectrum E(k x ,ky ) is the wave number at the point where the energy value is the largest, and the calculation formula is:

[0158]

[0159] Among them, k xm , k ym are the horizontal and vertical coordinates of the point with the maximum wave number.

[0160] Step 3, calculating the average zero-crossing period, includes the following steps:

[0161] Step 3.1, from step 1, the relationship between the wave crest period and the average zero-crossing period can be obtained. From step 2, the wave crest period can be calculated by 3D-FFT spectrum analysis method, and the wave crest period is substituted into the relationship to calculate the average zero-crossing period.

[0162] Step 4, calculating the root mean square steepness of the waves by the shadow method, comprises the following steps:

[0163] Step 4.1, select the study area. Step 4.1 includes the following steps:

[0164] Step 4.1.1, select a sector area with an azimuth angle range of θ1-θ2 and a radial distance range of r1-r2 in the original radar image as the study area, and divide the original radar image into m sub-areas in the radial dimension in the selected study area.

[0165] Step 4.2, calculating the threshold for distinguishing the shadow area from the non-shadow area, step 4.2 includes the following steps:

[0166] Step 4.2.1, image edge detection, by performing convolution operations on the original radar image in t adjacent different directions based on the differential operator, the edge gradient image I is obtained Ei (r,θ), the calculation formula is:

[0167]

[0168] Where I(r,θ) is the original radar image, r is the distance, and θ is the azimuth. represents the convolution operation, I Ei (r,θ) is the edge gradient image obtained after convolution operation with different difference operators, H i In (r,θ), i=1,2,3...t, which represents the pixel difference operator in t different directions.

[0169] Step 4.2.2: The differential image I obtained in step 4.2.1. Ei The pixel points with the intensity values ​​in the first N% of (r,θ) are taken as the boundary pixels of the shadow area and the non-shadow area to obtain the sub-edge image ITi (r,θ), where i=1,2,3...t, the specific formula is:

[0170]

[0171] Step 4.2.3, the sub-edge image I obtained in step 4.2.2 Ti (r,θ) are superimposed to obtain the original edge image I T (r,θ), where i = 1, 2, 3...t:

[0172]

[0173] Then use the edge comparison method to transform the original edge image I T Each point in (r,θ) is related to the threshold τ F Compare to remove isolated noise points and get edge image I F (r,θ), the specific formula is:

[0174]

[0175] Step 4.2.4, for each isolated 1 in the edge image, find the pixel η of the original radar image at that position, and establish an image histogram function F through the correspondence between the edge image and the original radar image H (η), for F H (η) Take the majority to get the shadow threshold τ S .

[0176] τ S =mode(F H (η))

[0177] Step 4.3, calculating the shadow ratio function of each sub-area, step 4.3 includes the following steps:

[0178] Step 4.3.1, use the shadow threshold τ obtained in step 4.2 S The original radar image I T (r,θ) is divided into shadow area and non-shadow area, and the shadow image I s (r,θ) is:

[0179]

[0180] Among them, I(r,θ) is the original radar image, I s (r,θ) is the segmented shadow image.

[0181] Step 4.3.2, calculate the grazing angle γ of each sub-area i, the grazing angle corresponding to the center point of the sub-region is the grazing angle of the sub-region. Assume that the radar antenna height is H, the radar image distance resolution is R, the number of reference points at the starting position is h0, and the number of relative points at the start of the i-th sub-region is p i , the relative number of points at the end is q i , where i = 1, 2, 3...m, then the calculation formula for the grazing angle of the i-th sub-area is:

[0182]

[0183] Step 4.3.3, calculate the shadow ratio function S(γ i ), that is, calculate the ratio of the number of pixels with a pixel value of 1 in the shadow image obtained in step 4.2.2 corresponding to each sub-region to the number of all pixels in the sub-region where it is located, where i = 1, 2, 3...m, then the calculation formula of the shadow ratio function of the i-th sub-region is:

[0184]

[0185] Where P i Represents the shadow image I corresponding to the i-th sub-region s The number of pixels with a pixel value of 1 in (r,θ), Q i Represents the shadow image I corresponding to the i-th sub-region s The number of all pixels in (r,θ).

[0186] Step 4.4, calculating the root mean square wave steepness, step 4.4 includes the following steps:

[0187] Step 4.4.1, the shadow ratio function S(γ i ) and lighting ratio L(γ i ) has the following relationship:

[0188] L(γ i )=1-S(γ i )

[0189] The least squares method is used to fit the illumination ratio L(γ i ) and the grazing angle γ i The relationship between , where i = 1, 2, 3 ... m, is used to estimate the root mean square wave steepness σ RMS , the Smith function is:

[0190]

[0191] In the formula,

[0192]

[0193] μi =tan(γ i )

[0194] Among them, erfc() is the Gaussian error function. The erfc() function is derived from the Gaussian surface probability density function. The expression of the Gaussian surface probability density function is:

[0195]

[0196] Where ξ is the height of the Gaussian sea level from the horizontal plane, υ is the root mean square of the height deviation, and Δξ is the height increment of the Gaussian surface when Δτ is added in the horizontal axis direction.

[0197] Step 5, calculating the effective wave height of the wave by averaging the zero-crossing period and the root mean square wave steepness parameter, step 5 comprises the following steps:

[0198] Step 5.1, calculate the effective wave height using the relationship between the effective wave height and the root mean square wave steepness and the average zero-crossing period. The specific formula is:

[0199]

[0200] Where g is the acceleration due to gravity, H S is the effective wave height, σ RMS is the root mean square wave steepness, T m02 is the average zero-crossing period.

[0201] An embodiment is given below in combination with specific parameters.

[0202] The present invention can be specifically divided into the following steps: the first step is offline data processing, and the relationship between the wave crest period and the average zero-crossing period is obtained by least square fitting; the second step is to obtain a two-dimensional wave spectrum by 3D-FFT spectrum analysis, and obtain the wave crest period of the wave from the two-dimensional wave spectrum; the third step is to calculate the average zero-crossing period; the fourth step is to calculate the root mean square steepness of the wave by the shadow method; and the fifth step is to calculate the effective wave height of the wave by the average zero-crossing period and the root mean square steepness parameter.

[0203] The marine radar used in the implementation of the present invention is an X-band navigation radar, which works in short pulse mode. The echo data is digitized and stored in polar coordinates by line. The time interval between two adjacent stored lines is less than 1ms. The radar antenna scans one circle in about 2.5s. The total number of lines in a radar image is 2048, and there are 2048 pixel points on each line. The azimuth resolution is about 0.1° and the radial resolution is about 2.5m.

[0204] The main technical parameters of the marine radar used in the present invention are shown in Table 1:

[0205] Table 1 Technical parameters of marine radar

[0206]

[0207] Combination Fig.10 , the specific experimental steps of the present invention are:

[0208] The first step is offline data processing, which uses the least squares method to fit the relationship between the wave crest period and the average zero-crossing period, including the following steps:

[0209] Step 1.1, using offline buoy data and radar, the relationship between the wave crest period and the average zero-crossing period is fitted by the least squares method. Figure 1 This is the fitting diagram of the radar wave crest period and the average zero-crossing period of the buoy in this example. The horizontal axis in the figure is the wave crest period calculated by the radar, and the vertical axis is the average zero-crossing period of the wave measured by the buoy. The fitting relationship is:

[0210] T m02 =0.6136T M +2.538

[0211] Where, T M is the peak period, T m02 is the average zero-crossing period. The correlation coefficient of the fitting result in this example is 0.75, which shows that the linear model can better reflect the relationship between the two.

[0212] The second step is to obtain a two-dimensional wave spectrum through 3D-FFT spectrum analysis method, and obtain the wave crest period of the wave from the two-dimensional wave spectrum, including the following steps:

[0213] Step 2.1, select the study area size as L x *L y , the spatial resolution of the radar is R, so the pixel size of the study area is M x *M y In this example, L x =L y =960m, R=7.5m, M x =M y =128, Figure 2 Select a schematic diagram for your study area. Figure 3 Select the diagram for the Cartesian box, Figure 3 The black area is the selected Cartesian box.

[0214] Step 2.2, obtain the two-dimensional wave spectrum through 3D-FFT spectrum analysis method, and obtain the wave crest period from the two-dimensional wave spectrum. The flow chart of obtaining wave parameters by 3D-FFT spectrum analysis method is shown in Figure 4 , step 2.2 comprises the following steps:

[0215] Step 2.2.1, through the nearest point interpolation, convert the echo intensity value in the polar coordinates of the original radar image to the corresponding point in the Cartesian coordinate system. The coordinates of a point in the original radar image are marked as (r0, θ0), and the corresponding Cartesian coordinates after conversion are marked as (x0, y0). The radar image sequence obtained after interpolation is recorded as η(x, y, t), and the conversion formula is:

[0216]

[0217] Where r0 is the radial distance of a point in the original radar image, θ0 is the azimuth of a point in the original radar image, x0 and y0 are the abscissa and ordinate values ​​of the Cartesian coordinate system, sqrt(), round() and rem() are the square root function, integer function and remainder function, respectively.

[0218] Step 2.2.2, perform a three-dimensional Fourier transform on the radar image sequence η(x, y, t) obtained by the nearest point interpolation in step 2.2.1 to obtain a three-dimensional image spectrum F(k x ,k y ,ω), the transformation formula is:

[0219]

[0220] The size of the selected research area is L x *L y The Cartesian frame, L x ,L y is the range of the selected Cartesian box along the x and y directions, and the total time of the radar image sequence is T, then the above transformation formula can be written as:

[0221]

[0222] The discrete form of the transformation formula is used to transform the image sequence from the time-space domain to the wavenumber frequency domain. Its discrete form is:

[0223]

[0224] In the formula, k x and k y are the wave numbers in the x and y directions respectively, and ω is the angular frequency. The resolution is calculated as follows:

[0225]

[0226] Among them, L x ,L y is the extent of the selected Cartesian box in the x and y directions, and T is the total time of the radar image sequence. x =L y=960m, T=90s.

[0227] Step 2.2.3, estimate the ocean current information based on the dispersion principle and the least squares principle. The ocean current calculation formula is:

[0228]

[0229] Where ω0 is the theoretical frequency corresponding to the wave number, h is the water depth, g is the gravitational acceleration, and F is the three-dimensional image spectrum F(k x ,k y ,ω), sea surface laminar flow Wave number Among them, k x ,k y are the wave numbers in the x and y directions respectively.

[0230] Step 2.2.4, filtering the three-dimensional image spectrum obtained in step 2.2.2 using a bandpass filter based on dispersion relations, the filter used is a bandpass filter based on maximum flow velocity:

[0231]

[0232] In the formula, E(k x ,k y ,ω) is the three-dimensional image spectrum obtained after filtering, F(k x ,k y ,ω) is the three-dimensional image spectrum obtained by step 2.2.2, where

[0233]

[0234]

[0235] Where k is the wave number, dω is the frequency resolution, dk = (dk x =dk y ) is the wave number resolution, U max is the maximum flow rate, when B n When the value is negative, let B n =0.

[0236] in,

[0237]

[0238] Where, L x ,L y is the extent of the selected Cartesian box in the x and y directions, and T is the total time of the radar image sequence. x =L y =960m, T=90s.

[0239] Then the three-dimensional spectrum E(k x ,k y ,ω) and perform frequency ω-dimensional integration to obtain the two-dimensional image spectrum I(k x ,k y ), the integral formula is:

[0240] I(k x ,k y )=∫ ω＞0 E(k x ,k y ,ω)dω

[0241] Step 2.2.5, the two-dimensional image spectrum I(k x ,k y ) is converted into a two-dimensional wave spectrum E(k x ,k y ), the conversion formula is:

[0242] E(k x ,k y )=|M(k x ,k y )| 2 ·I(k x ,k y )

[0243] The empirical modulation transfer function is obtained by fitting experimental data. The commonly used empirical modulation transfer function is in the form of:

[0244] |M(k x ,k y )| 2 ≈k -β

[0245] Here, β is an empirical coefficient, and in this example, β is set to 1.2.

[0246] Step 2.2.6, calculate the wave crest period, using the two-dimensional wave spectrum E(k x ,k y ) Using the dispersion relation to calculate the peak period T of the ocean wave M , the calculation formula is:

[0247]

[0248] Where g is the acceleration due to gravity, k is m is the wave spectrum E(k x ,k y ) is the wave number at the point where the energy value is the largest, and the calculation formula is:

[0249]

[0250] Among them, k xm , k ym are the horizontal and vertical coordinates of the point with the maximum wave number.

[0251] This example selects a set of data with an effective wave height of 2.0m, and then uses the 3D-FFT spectrum analysis method to obtain the corresponding two-dimensional wave spectrum. The results are shown in Figure 5 , the corresponding peak period is calculated to be 8.21s. It can be seen from the figure that the two-dimensional wave spectrum can better characterize the peak period of the waves.

[0252] The third step is to calculate the average zero-crossing period, including the following steps:

[0253] Step 3.1: The relationship between the wave crest period and the average zero-crossing period can be obtained from the first step. The wave crest period can be obtained by 3D-FFT spectrum analysis in the second step. The wave crest period is substituted into the relationship to calculate the average zero-crossing period. The calculation results are shown in Figure 6 In the figure, the solid line is the average zero-crossing period of the buoy output, and the dotted line is the average zero-crossing period of the waves obtained by using the radar data through the 3D-FFT spectrum analysis algorithm. This example selects a set of data with an effective wave height of 2.0m. The corresponding wave crest period is calculated to be 8.21s in the second step. The corresponding average zero-crossing period is 7.57s when the wave crest period is substituted into the relationship.

[0254] The fourth step is to calculate the root mean square steepness of the waves by the shadow method, which includes the following steps:

[0255] Step 4.1, select the study area. Step 4.1 includes the following steps:

[0256] Step 4.1.1, select a sector area with an azimuth angle range of θ1-θ2 and a radial distance range of r1-r2 in the original radar image as the study area, and divide the original radar image into m sub-areas in the selected study area. In this example, m = 50, θ1 = 120°, θ2 = 230°, r1 = 200m, r2 = 2500m, and the shadow method area segmentation diagram is shown in Figure 7 .

[0257] Step 4.2, calculating the threshold for distinguishing the shadow area from the non-shadow area, step 4.2 includes the following steps:

[0258] Step 4.2.1, image edge detection, by performing convolution operations on the original radar image in t adjacent different directions based on the differential operator, the edge gradient image I is obtained Ei (r,θ), the calculation formula is:

[0259]

[0260] Where I(r,θ) is the original radar image, r is the distance, = is the azimuth, represents the convolution operation, I Ei (r,θ) is the edge gradient image obtained after convolution operation with different difference operators, H i (r,θ) represents the pixel point difference operator in t different directions, where i=1,2,3...t, and t=8 in this example, specifically:

[0261]

[0262] Step 4.2.2: The differential image I obtained in step 4.2.1. Ei The pixel points with the intensity values ​​in the first N% of (r,θ) are taken as the boundary pixels of the shadow area and the non-shadow area to obtain the sub-edge image I Ti (r,θ), where i=1,2,3...t, the specific formula is:

[0263]

[0264] In this example, t=8, and N is 10.

[0265] Step 4.2.3, the sub-edge image I obtained in step 4.2.2 Ti (r,θ) are superimposed, where i = 1, 2, 3...t, to obtain the original edge image I T (r,θ):

[0266]

[0267] In this example, t=8.

[0268] Then use the edge comparison method to transform the original edge image I T Each point in (r,θ) is related to the threshold τ F Compare to remove isolated noise points and get edge image I F (r,θ), the specific formula is:

[0269]

[0270] In this example, τ F The value is 5.

[0271] Step 4.2.4, for each isolated 1 in the edge image, find the pixel η of the original radar image at that position, and establish an image histogram function F through the correspondence between the edge image and the original radar image H (η), for F H (η) Take the majority to get the shadow threshold τ S .

[0272] τ S =mode(F H (η))

[0273] Step 4.3, calculating the shadow ratio function of each sub-area, step 4.2 includes the following steps:

[0274] Step 4.3.2: Use the shadow threshold τ obtained in step 4.1 S The original radar image I T (r,θ) is divided into shadow area and non-shadow area, and the shadow image I s (r,θ) is:

[0275]

[0276] Among them, I(r,θ) is the original radar image, I s (r,θ) is the segmented shadow image.

[0277] Step 4.3.3, calculate the grazing angle γ of the sub-region i , the grazing angle corresponding to the center point of the sub-region is the grazing angle of the sub-region. Assume that the radar antenna height is H, the radar image distance resolution is R, the number of reference points at the starting position is h0, and the number of relative points at the start of the i-th sub-region is p i , the relative number of points at the end is q i , i = 1, 2, 3...m, then the calculation formula for the grazing angle of the i-th sub-area is:

[0278]

[0279] In this example, the radar antenna height H is 45 m, the radar image distance resolution R is 7.5 m, the number of starting position reference points h0 is 25, and m=50.

[0280] Step 4.3.4, calculate the shadow ratio function S(γ i ), that is, calculate the ratio of the number of pixels with a pixel value of 1 in the shadow image obtained in step 4.2.2 corresponding to each sub-region to the number of all pixels in the sub-region where it is located, i = 1, 2, 3...m, then the calculation formula of the shadow ratio function of the i-th sub-region is:

[0281]

[0282] Where P i Represents the shadow image I corresponding to the i-th sub-region s The number of pixels with a pixel value of 1 in (r,θ), Q i Represents the shadow image I corresponding to the i-th sub-region s The number of all pixels in (r,θ), in this example m=50.

[0283] Step 4.4, calculating the root mean square wave steepness, step 4.4 includes the following steps:

[0284] Step 4.4.1, the shadow ratio function S(γ i ) and lighting ratio L(γ i ) has the following relationship:

[0285] L(γ i )=1-S(γ i )

[0286] The least squares method is used to fit the illumination ratio L(γ i ) and the grazing angle γ i The relationship between i = 1, 2, 3...m, in this example m = 50, the estimated root mean square wave steepness σ RMS , the Smith function is:

[0287]

[0288] In the formula,

[0289]

[0290] μ i =tan(γ i )

[0291] Among them, erfc() is the Gaussian error function. The erfc() function is derived from the Gaussian surface probability density function. The expression of the Gaussian surface probability density function is:

[0292]

[0293] Where ξ is the height of the Gaussian sea level from the horizontal plane, υ is the root mean square of the height deviation, and Δξ is the height increment of the Gaussian surface when Δτ is added in the horizontal axis direction.

[0294] The fifth step is to calculate the effective wave height of the waves by averaging the zero-crossing period and the root mean square wave steepness parameter, including the following steps:

[0295] Step 5.1, calculate the effective wave height using the relationship between the effective wave height and the root mean square wave steepness and the average zero-crossing period. The specific formula is:

[0296]

[0297] Where g is the acceleration due to gravity, H S is the effective wave height, σ RMS is the root mean square wave steepness, T m02 is the average zero-crossing period.

[0298] In this example, a set of data with an effective wave height of 2.0m is selected. The corresponding peak period is calculated to be 8.21s in step 2. The corresponding average zero-crossing period is calculated to be 7.57s by substituting the peak period into the relationship obtained by fitting in step 1 in step 3. The root mean square steepness is calculated to be 2.38 in step 4. The effective wave height is calculated to be 2.80m in step 5.

[0299] The effective wave height inversion method based on X-band marine radar images proposed in this invention is applied to conduct experimental analysis on a large amount of radar data acquired in the East China Sea in 2014 and the sea condition information of the relevant time period. This experiment selects the shore-based data collected in Pingtan, Fujian from 18:00 on November 16, 2014 to 0:00 on November 19, 2014 for the experiment. The shadow method, spectral analysis method and shadow statistics method combined algorithm are used for comparative experiments. The calculation results of the effective wave height series by the shadow method, spectral analysis method and shadow statistics method combined algorithm are shown in Figure 8 The figure shows the effective wave heights obtained by the buoy, shadow method and combined algorithm. Figure 8 It can be seen that, taking the effective wave height data output by the buoy as the true value, the effective wave height data calculated by the shadow method and the combined algorithm are relatively high overall. Figure 9(a) and Figure 9(b) are scatter plots comparing the effective wave height obtained by the shadow method and the combined algorithm with the effective wave height of the buoy, where Figure 9(a) is a scatter plot of the effective wave height of the shadow method and the effective wave height of the buoy, and Figure 9(b) is a scatter plot of the effective wave height of the combined algorithm and the effective wave height of the buoy.

[0300] In the present invention, firstly, the average zero-crossing period T outputted by the radar 3D-FFT analysis method and the buoy is respectively m02 For comparative analysis, the wave period of the buoy is taken as the reference truth value, and the average deviation BIAS and the correlation coefficient CORR are used as the parameters of the evaluation algorithm. i is the effective wave height sequence calculated by the combined algorithm of the invention, Y i is the effective wave height true value output by the buoy, where i=1,2,...N, then:

[0301]

[0302]

[0303] Table 2 can be obtained by calculating the above formula:

[0304] Table 2 Algorithm performance parameters

[0305]

[0306] The experimental results show that although the performance of the combined algorithm is lower than that of the shadow method, the combined algorithm can independently use radar data to obtain the effective wave height parameters of waves, making the process of using radar to measure the effective wave height parameters of waves completely independent of the support of buoys and no longer requiring the help of other instruments, which increases the applicability and independence of the algorithm.

Claims

1. A method for inverting significant wave height from marine radar images, characterized in that: The following steps are involved: Step 1: Using offline buoy data and radar, the relationship between the wave crest period and the average zero-crossing period is obtained by least squares fitting; Step 2: Select a research area, obtain the two-dimensional wave spectrum of the selected research area through 3D-FFT spectrum analysis method, and obtain the wave crest period of the wave from the two-dimensional wave spectrum; Step 3: Substitute the wave crest period obtained in step 2 into the relationship between the wave crest period and the average zero-crossing period obtained in step 1 to calculate the average zero-crossing period T m02 ; Step 4: Select a research area in the original radar image, divide the original radar image into m sub-areas in the radial dimension in the selected research area; calculate the threshold τ to distinguish between shadow areas and non-shadow areas S ; According to τ S Calculate the shadow ratio function of each sub-area; then calculate the root mean square wave steepness σ according to the shadow ratio function RMS ; Step 5: By averaging the zero-crossing period T m02 and the root mean square steepness σ RMS Calculate the effective wave height H of the ocean wave S , H S satisfy: Where g is the acceleration due to gravity.

2. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: The relationship between the wave crest period and the average zero-crossing period of the ocean wave is specifically: T m02 =aT M +b Where, T M is the peak period, T m02 is the average zero-crossing period; The average zero-crossing period sequence of the buoy in the offline experiment is T m02i , i = 1, 2, 3...n, and the corresponding peak period sequence obtained by 3D-FFT spectrum analysis using radar data is recorded as T Mi , i=1,2,3...n, according to the least squares method, we get: In the formula, and are the average values ​​of the average zero-crossing period sequence and the peak period sequence respectively.

3. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: The step 2 comprises: Step 2.1: Through the nearest point interpolation, the echo intensity value in the polar coordinates of the original radar image is converted to the corresponding point in the Cartesian coordinate system. The coordinates of a point in the original radar image are recorded as (r0, θ0), and the corresponding Cartesian coordinates after conversion are recorded as (x0, y0). The radar image sequence obtained after interpolation is recorded as η(x, y, t). The conversion formula is: Where r0 is the radial distance of a point in the original radar image, θ0 is the azimuth of a point in the original radar image, x0 and y0 are the horizontal and vertical coordinates of the Cartesian coordinate system, sqrt(), round() and rem() are the square root function, integer function and remainder function respectively; Step 2.2: Select the study area size as L x *L y The Cartesian frame, L x ,L y is the range of the Cartesian frame along the x and y directions, the total time of the radar image sequence is T, and the radar image sequence η(x, y, t) is transformed into a three-dimensional image spectrum by a three-dimensional Fourier transform. Then, the image sequence is transformed from the time-space domain to the wavenumber frequency domain using the discrete form of the transformation formula: In the formula, k x and k y are the wave numbers in the x and y directions respectively, and ω is the angular frequency. The resolution is calculated as follows: Step 2.3: Assuming that the theoretical frequency corresponding to the wave number is ω0, and the real frequency ω is expressed by the dispersion relation, the sum of the squares of the interpolation between the two is: Taking the energy value in the three-dimensional image spectrum as the weight, the objective function is: When the partial derivatives of SSE with respect to are all zero, SSE has a minimum solution, and the flow velocity is calculated as: Where h is the water depth, g is the gravitational acceleration, and F is the three-dimensional image spectrum F(k x ,k y ,ω), sea surface laminar flow Wave number Among them, k x ,k y are the wave numbers in the x and y directions respectively; Step 2.4: Filter the three-dimensional image spectrum using a bandpass filter based on the dispersion relation. The filter used is a bandpass filter based on the maximum flow velocity: In the formula, E(k x ,k y ,ω) is the three-dimensional image spectrum obtained after filtering, where Where k is the wave number, dω is the frequency resolution, dk = (dk x =dk y ) is the wave number resolution, U max is the maximum flow rate, when B n When the value is negative, let B n =0; Then the three-dimensional image spectrum E(k x ,k y ,ω) and perform frequency ω-dimensional integration to obtain the two-dimensional image spectrum I(k x ,k y ): I(k x ,to y )=∫ ω＞0 E(k x ,to y ,ω)dω Step 2.5: Use the empirical modulation transfer function to transform the two-dimensional image spectrum I(k x ,k y ) is converted into a two-dimensional wave spectrum E(k x ,k y ): E(k x ,k y )=|M(k x ,k y )| 2 ·I(k x ,k y ) The empirical modulation transfer function is: |M(k x ,k y )| 2 ≈k -β Among them, β is the empirical coefficient; Step 2.6: Through the two-dimensional wave spectrum E(k x ,k y ) Using the dispersion relation to calculate the peak period T of the ocean wave M , the calculation formula is: Where g is the acceleration due to gravity, k is m is the wave spectrum E(k x ,k y ) is the wave number at the point where the energy value is the largest, and the calculation formula is: Among them, k xm , k ym are the horizontal and vertical coordinates of the point with the maximum wave number.

4. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: Step 4 calculates the threshold τ for distinguishing between shadow area and non-shadow area S include Step 4.2.1: Image edge detection: by performing convolution operations on the original radar image in t adjacent different directions based on the differential operator, the edge gradient image I is obtained. Ei (r, θ): Where I(r, θ) is the original radar image, r is the distance, and θ is the azimuth. represents the convolution operation, I Ei (r, θ) is the edge gradient image obtained after convolution operation with different difference operators, H i In (r, θ), i=1, 2, 3...t, representing the pixel difference operator in t different directions; Step 4.2.2: Subtract the difference image I Ei The pixel points with the intensity values ​​in the first N% of (r, θ) are taken as the boundary pixels of the shadow area and the non-shadow area to obtain the sub-edge image I Ti (r, θ), where i = 1, 2, 3...t, the specific formula is: Step 4.2.3: I Ti (r, θ) are superimposed to obtain the original edge image I T (r, θ), where i = 1, 2, 3, ... t: Then use the edge comparison method to transform the original edge image I T Each point in (r, θ) is equal to the threshold τ F Compare to remove isolated noise points and get edge image I F (r, θ): Step 4.2.4: For each isolated 1 in the edge image, find the pixel η of the original radar image at that position, and establish an image histogram function F through the correspondence between the edge image and the original radar image. H (η), for F H (η) Taking the mode, we can get the threshold τ to distinguish between shadow area and non-shadow area S : t S =mode(F H (or)).

5. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: Step 4: According to τ S The shadow ratio function for each sub-area is calculated as follows: Step 4.3.1: Using threshold τ S The original radar image I T (r, θ) is divided into shadow area and non-shadow area, and the shadow image I s (r, θ) is: Among them, I(r, θ) is the original radar image, I s (r, θ) is the shadow image after segmentation; Step 4.3.2: Calculate the grazing angle γ of each sub-region i , the grazing angle corresponding to the center point of the sub-region is the grazing angle of the sub-region. Assume that the radar antenna height is H, the radar image distance resolution is R, the number of reference points at the starting position is h0, and the number of relative points at the start of the i-th sub-region is p i , the relative number of points at the end is q i , where i = 1, 2, 3...m, then the calculation formula for the grazing angle of the i-th sub-area is: Step 4.3.3: Calculate the shadow ratio function S(γ i ), where i = 1, 2, 3...m, then the calculation formula of the shadow ratio function of the i-th sub-area is: Where P i Represents the shadow image I corresponding to the i-th sub-region s The number of pixels with a pixel value of 1 in (r, θ), Q i Represents the shadow image I corresponding to the i-th sub-region s The number of all pixels in (r, θ).

6. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: The step 4 of calculating the root mean square wave steepness according to the shadow ratio function includes: Step 4.4.1: Shadow ratio function S(γ i ) and lighting ratio L(γ i )satisfy: L(γ i )=1-S(γ i ) The least squares method is used to fit the illumination ratio L(γ i ) and the grazing angle γ i The relationship between , where i = 1, 2, 3...m, is used to estimate the root mean square wave steepness σ RMS , the Smith function is: In the formula, m i =tan(γ i ) Among them, erfc() is the Gaussian error function. The erfc() function is derived from the Gaussian surface probability density function. The expression of the Gaussian surface probability density function is: Where ξ is the height of the Gaussian sea level from the horizontal plane, υ is the root mean square of the height deviation, and Δξ is the height increment of the Gaussian surface when Δτ is added in the horizontal axis direction.

7. The method for inverting significant wave height from marine radar images according to claim 1, characterized in that: The navigation radar is an X-band navigation radar.

Images