A method for retrieving rainfall intensity levels

Abstract

The present invention discloses a method for inverting rainfall intensity levels. First, radar images under different rainfall intensities are collected offline, a suitable rainfall detection area is selected, and the corresponding wavenumber energy spectrum image is obtained by using two-dimensional discrete Fourier transform (2D-DFT). The characteristic parameters of the energy spectrum images under different rainfall intensities are extracted and fed into the constructed GA-BPNN model for training to obtain a model that meets the experimental requirements. Then, the radar image to be detected is selected, and its wavenumber energy spectrum image and characteristic parameters are calculated. Finally, the characteristic parameters are fed into the GA-BPNN model obtained offline to obtain the rainfall intensity level of the image to be detected. The present invention realizes the requirement of real-time acquisition of marine rainfall information during ship navigation and offshore operations.

Description

Technical Field

[0001] The present invention belongs to the technical field of marine remote sensing under marine environmental conditions, and relates to a method for retrieving rainfall intensity levels, in particular to a method for retrieving rainfall intensity levels based on a combination of wave number energy spectrum (WES) and a genetic algorithm optimized backpropagation neural network (GA-BPNN). The type of marine radar applicable to this method is an X-band marine radar. Background Art

[0002] In a complex marine environment, marine radars are widely used in ship navigation to guide ships to avoid marine obstacles and ensure the navigation safety of ships. Rainfall intensity is an important meteorological parameter of the marine environment and is one of the necessary meteorological element information for the continuous real-time takeoff and landing safety guarantee of carrier-based aircraft under the conditions of marine power carriers. Therefore, while performing the mission of safety guarantee, the marine radar on the ship takes into account the detection of rainfall information, makes full use of the rich rain modulation information in the marine radar image, and obtains rainfall information in real time, which has important practical value for saving ship resources, controlling equipment costs, improving equipment convenience and effectiveness, etc.

[0003] Currently, in the research on obtaining maritime rainfall information using marine radar, the representative achievements are as follows: In 2012, Lund et al. found that there are certain differences in the zero-pixel percentage between rainfall and rainless radar images, and for the first time proposed using the zero-pixel percentage parameter to detect whether there is rainfall in the radar image (Lund B, Graber H C, Romeiser R. Wind Retrieval From Shipborne Nautical X-Band Radar Data[J]. IEEE Transactions on Geoscience & Remote Sensing, 2012, 50(10): 3800-3811). Huang et al. proposed the high clutter pixel (HCD) method and the support vector machine (SVM) method for rainfall radar image detection (Huang W, Liu X, Gill E W. An Empirical Mode Decomposition Method for Sea Surface Wind Measurements From X-Band Nautical Radar Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(11): 6218-6227. Chen X, Huang W, Zhao C, et al. Rain Detection From X-Band Marine Radar Images: A Support Vector Machine-Based Approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58). Lu et al. proposed the correlation feature method and the wave texture difference (WTD) method for rainfall radar image detection (Zheng Y, Shi Z, Lu Z, et al. A Method for Detecting Rainfall From X-Band Marine Radar Images[J]. IEEE Access, 2020, 8: 19046-19057. Sun L, Lu Z, Wang H, et al. A Wave Texture Difference Method for Rainfall Detection Using X-Band Marine Radar[J]. Journal of Sensors, 2022).These studies mainly focus on using marine radar images for rainfall image detection, but there is little research on real-time inversion of rainfall intensity. Lu Z, Sun L, Zhou Y. A Method for Rainfall Detection and Rainfall Intensity Level Retrieval from X-Band Marine Radar Images[J]. Applied Sciences, 2021, 11(4): 1565 proposed an RZE method based on the occluded area of radar images for rainfall image detection and rainfall intensity level inversion. However, this method uses the least squares method for high-order fitting to invert rainfall intensity, and there is a certain degree of randomness in the selection and processing of data. Therefore, the accuracy of this method is limited. The present invention discloses an algorithm for inverting rainfall intensity level by combining WES and GA-BPNN, which transforms the size of spatial rainfall into the distribution problem of energy in the frequency domain, making it easier to extract parameter characteristics under different rainfall intensities. At the same time, an optimized neural network model is used to predict the rainfall intensity level, greatly improving the accuracy and robustness of the algorithm, and also providing a new research idea for real-time acquisition of rainfall information from marine radars. Summary of the Invention

[0004] 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 rainfall intensity level by combining wave number energy spectrum (WES) and genetically algorithm optimized backpropagation neural network (GA-BPNN), so as to realize real-time and accurate acquisition of rainfall intensity information during ship navigation.

[0005] To solve the above technical problem, a method for inverting rainfall intensity level of the present invention includes the following steps:

[0006] Step 1: Offline obtain a GA-BPNN model applicable to rainfall intensity level inversion

[0007] Offline collect multiple radar images under different rainfall intensities, screen the images, select the rainfall detection area, obtain the Cartesian image using the nearest neighbor interpolation method, and then obtain the wave number energy spectrum corresponding to the Cartesian image through two-dimensional discrete Fourier transform; divide the energy spectrum into different wave number segments, and count the proportion of wave number energy in each wave number segment to the total energy under different rainfall intensities as characteristic parameters, and send the characteristic parameters into the constructed GA-BPNN model for training and verification to obtain a model that meets the experimental requirements;

[0008] Step 2: Read the radar image to be detected, select the rainfall detection area, and calculate the Cartesian image using the same nearest neighbor interpolation method as in Step 1;

[0009] Step 3: Calculate the wavenumber energy spectrum using the same 2D-DFT as in Step 1, and calculate the energy proportion in different wavenumber segments as characteristic parameters;

[0010] Step 4: Input the characteristic parameters into the GA-BPNN model trained in Step 1 to obtain the rainfall intensity level of the image to be detected.

[0011] Further, the specific method for obtaining the Cartesian image using the nearest neighbor interpolation method in Step 1 is as follows:

[0012]

[0013] In the formula, r is the distance between a certain point in the radar image and the center point, and θ is the angle formed with the bow direction.

[0014] Further, the specific method for obtaining the wavenumber energy spectrum corresponding to the Cartesian image through two-dimensional discrete Fourier transform in Step 1 is as follows:

[0015] Step 1.1.1: The number of full waves appearing in a 2π length is defined as the wavenumber k, satisfying:

[0016]

[0017] In the formula, (k x ,k y ) are the coordinates in the energy spectrum domain, that is, the wavenumber components of k on the x and y axes, with the unit of rad / m; λ is the wavelength of the sea surface stripes presented in the radar image;

[0018] Step 1.1.2: Perform 2D-DFT on the Cartesian image z = f(x, y), and its mathematical model can be expressed as:

[0019]

[0020] In the formula, N is the number of sampling points of the sub-image z, x, y = 1, 2,..., N are frequency variables, F(k x ,k y ) are the Fourier coefficients of the sub-image z, (k x ,k y ) are the wavenumber domain coordinates k x and k y are the components of the wavenumber vector in the x and y directions respectively. Taking the vector absolute value of F(k x ,k y ) gives the wavenumber energy spectrum Specifically:

[0021]

[0022] In the formula, Re(F(kx , k y ), and Im(F(k x , k y )) are the real and imaginary parts of F(k x , k y ).

[0023] Furthermore, in step 1, the energy spectrum is divided into different wavenumber segments, and the proportion of the wavenumber energy in each wavenumber segment under different rainfall intensities to the total energy is counted as the characteristic parameter, specifically as follows:

[0024] Step 1.2.1: Let λ d represent the minimum value of the ocean wave wavelength scale, and λ t represent the maximum value of the ocean wave wavelength scale. Then the lower limit of the wavenumber of the ocean wave energy spectrum is:

[0025]

[0026] The upper limit of the wavenumber of the ocean wave energy spectrum is:

[0027]

[0028] The wavenumber energy spectrum is divided into three wavenumber segments: low, medium, and high. Among them, the energy spectrum of the low wavenumber segment is The energy spectrum of the ocean wave wavenumber segment (medium wavenumber segment) is The energy spectrum of the high wavenumber segment is

[0029] Step 1.2.2: The energy of each wavenumber segment is expressed as:

[0030]

[0031] In the formula, En L , En M , En H are the energies of the low, medium, and high wavenumber segments respectively. The total energy of the wavenumber energy spectrum is:

[0032]

[0033] Therefore, the proportions of the energies of the low, medium, and high wavenumber segments to the total energy are:

[0034]

[0035] In the formula, P L , P M , P H are the energy proportions of the low, medium, and high wavenumber segments respectively, and they are used as characteristic parameters.

[0036] Furthermore, in step 1, the characteristic parameters are sent into the constructed GA - BPNN model for training and verification, specifically as follows:

[0037] Step 1.3.1: Initialize the BPNN model, send the feature parameter data into the network model, the data is divided into a training set and a validation set, and the model generates initial weights and thresholds;

[0038] Step 1.3.2: Use GA to perform binary encoding on the generated initial values, and randomly generate N individuals as the initial population;

[0039] Step 1.3.3: Establish a fitness function Fitness:

[0040]

[0041] In the formula, and respectively represent the i-th predicted output and the actual output obtained by the k-th training sample through the neural network, Nt and Ns respectively represent the number of training set samples and the number of output nodes;

[0042] Step 1.3.4: Perform selection, crossover, and mutation operations on the initial population;

[0043] Step 1.3.5: After several generations of selection, crossover, and mutation operations, if the calculation error reaches the allowable error, stop the calculation, and replace the initial weights and thresholds in the BPNN with the calculated optimal weights and thresholds;

[0044] Step 1.3.6: Adjust the weights and thresholds through the training of the BPNN. If the set number of training times is completed or the set accuracy is reached, the training terminates, and the accuracy of the model prediction data is output.

[0045] Advantages of the present invention: The present invention discloses a rainfall intensity level inversion method combining WES and GA - BPNN, which solves the problem of real - time and accurate acquisition of rainfall intensity information during ship navigation in actual ship production and life, improves the inversion efficiency, and further realizes the conservation of ship resources, control of equipment costs, and improvement of equipment convenience and timeliness. Description of the Drawings

[0046] Figure 1 is the determination of the rainfall detection area in the radar image;

[0047] Figure 2 is the wave number energy spectrum (moderate rain);

[0048] Figure 3 is the divided low, medium, and high wave number segments;

[0049] Figure 4 is the construction of GA - BPNN;

[0050] Figure 5It is an iterative process of fitness;

[0051] Figure 6 It is the training performance of GA-BPNN;

[0052] Figure 7(a) is the inversion result of GA-BPNN, specifically the comparison between the predicted rainfall intensity level and the actual rainfall intensity level;

[0053] Figure 7(b) is the error graph of Figure 7(a);

[0054] Figure 8 It is the determination of the accuracy of GA-BPNN;

[0055] Figure 9 It is the wavenumber energy spectrum (light rain);

[0056] Figure 10 It is the simultaneously recorded rainfall rate;

[0057] Figure 11 It is Figure 10 the classification result of rainfall intensity levels;

[0058] Figure 12(a) is the inversion result of the RZE method, specifically the comparison between the predicted rainfall intensity level and the actual rainfall intensity level;

[0059] Figure 12(b) is the error graph of Figure 12(a);

[0060] Figure 13 It is the determination of the accuracy of the RZE method;

[0061] Figure 14 It is the method flow chart. Specific implementation mode

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

[0063] In conjunction with Figure 14 , the present invention includes the following steps:

[0064] Step 1: Offline obtain a GA-BPNN model applicable to rainfall intensity level inversion

[0065] Offline collect multiple radar images under different rainfall intensities, screen the images. Select a suitable rainfall detection area, obtain a Cartesian image using the nearest neighbor interpolation method, and then obtain the wavenumber energy spectrum corresponding to the Cartesian image through two-dimensional discrete Fourier transform (2D-DFT). Divide the energy spectrum into different wavenumber segments, and count the proportion of the wavenumber energy in each wavenumber segment under different rainfall intensities to the total energy as a characteristic parameter. Feed the characteristic parameter into the constructed GA-BPNN model for training and verification to obtain a model that meets the experimental requirements.

[0066] Step 2: Read the radar image to be detected, select a suitable rainfall detection area, and calculate the Cartesian image using the nearest neighbor interpolation method.

[0067] Step 3: Calculate the wave number energy spectrum using 2D-DFT, and calculate the energy proportion in different wave number segments as characteristic parameters.

[0068] Step 4: Send the characteristic parameters into the GA-BPNN model trained in Step 1 to obtain the rainfall intensity level of the image to be detected.

[0069] The GA-BPNN model meeting the experimental requirements in Step 1 specifically includes:

[0070] Step 1.1, Obtaining the Cartesian image. Use the nearest neighbor interpolation method for the selected rainfall detection area to obtain the Cartesian image, and the calculation method is

[0071]

[0072] where r is the distance between a certain point in the radar image and the center point, and θ is the angle formed with the bow direction.

[0073] Step 1.2, Calculate the wave number energy spectrum

[0074] Step 1.2.1, The number of full waves appearing in a 2π length is defined as the wave number k, and its calculation formula is

[0075]

[0076] where (k x ,k y ) are the coordinates in the energy spectrum domain, that is, the wave number components of k on the x and y axes, with the unit of rad / m; λ is the wavelength of the sea surface stripes presented in the radar image.

[0077] Step 1.2.2, Perform 2D-DFT on the Cartesian image z = f(x, y), and its mathematical model can be expressed as:

[0078]

[0079] where N is the number of sampling points of the sub-image z, x, y = 1, 2,..., N are frequency variables, F(k x ,k y ) are the Fourier coefficients of the sub-image z, (k x ,k y ) are the wave number domain coordinates k x and k y are the components of the wave number vector in the x and y directions respectively. For F(k x ,k y)Take the absolute value of the vector to obtain the wavenumber energy spectrum

[0080]

[0081] where Re(F(k x ,k y )) and Im(F(k x ,k y )) are the real and imaginary parts of F(k x ,k y ).

[0082] Step 1.3, Obtaining characteristic parameters in the energy spectrum

[0083] Step 1.3.1, Since ocean waves have a certain wavelength range, let λ d represent the minimum value of the ocean wave wavelength scale, and λ t represent the maximum value of the ocean wave wavelength scale. Then the lower limit of the wavenumber of the ocean wave energy spectrum can be expressed as

[0084]

[0085] The upper limit of the wavenumber of the ocean wave energy spectrum can be expressed as

[0086]

[0087] In this way, the wavenumber energy spectrum is divided into three wavenumber segments: low, medium, and high. Among them, the energy spectrum of the low wavenumber segment is The energy spectrum of the ocean wave wavenumber segment (medium wavenumber segment) is The energy spectrum of the high wavenumber segment is

[0088] Step 1.3.2, The energy of each wavenumber segment can be expressed as

[0089]

[0090] where En L , En M , En H are the energies of the low, medium, and high wavenumber segments respectively. The total energy of the wavenumber energy spectrum is

[0091]

[0092] Therefore, the proportions of the energies of the low, medium, and high wavenumber segments in the total energy can be expressed as

[0093]

[0094] where P L , P M , P HThey are the energy proportions in the low, medium, and high wavenumber bands respectively, which are used as characteristic parameters.

[0095] Step 1.4: Train and validate the constructed GA-BPNN model

[0096] Step 1.4.1: Initialize the BPNN model. Send the characteristic parameter data into the network model (three-fourths of the data is used as the training set, and the rest is used as the validation set), and the model generates the initial weights and thresholds.

[0097] Step 1.4.2: Use GA to perform binary encoding on the generated initial values. Randomly generate N individuals as the initial population.

[0098] Step 1.4.3: Establish the fitness function Fitness

[0099]

[0100] In the formula and respectively represent the i-th predicted output and the actual output obtained by the k-th training sample through the neural network. Nt and Ns represent the number of training set samples and the number of output nodes respectively.

[0101] Step 1.4.4: Perform selection, crossover, and mutation operations. Among them, selection refers to the operation of selecting excellent individuals from the population and eliminating inferior individuals, and the roulette wheel selection method is often used. Crossover represents replacing or recombining parts of the structures of two parents to generate new individuals, which can be expressed by the formula

[0102]

[0103] Here, the real number crossover method is used to perform crossover on two chromosomes a k and b k to obtain a' k and b' k , r is a random number between [0, 1], k = 1, 2... l, and l is the chromosome length.

[0104] The basic content of the mutation operator is to change the gene values at certain sites of the individual strings in the population. When the mutation operator mutates the k-th gene of individual a with a certain probability, the mutation operation is

[0105]

[0106] Here, a max and a min are respectively the upper and lower bounds of the gene value, G is the current iteration number, G maxis the maximum number of iterations, and r1 and r2 are random numbers between [0, 1]. To enhance the global search ability in the initial stage of iteration and the local search ability in the later stage of iteration, we also adopt an adaptive adjustment of the number of mutation bits, and the method can be expressed as: Here represents rounding up, and L max represents the preset maximum number of mutation bits.

[0107] Step 1.4.5, one generation refers to a process of completing selection, crossover, and mutation. After several generations, if the calculation error reaches the allowable error, the calculation is stopped, and the calculated optimal weights and thresholds are used to replace the initial weights and thresholds in the BPNN.

[0108] Step 1.4.6, fine-tune the weights and thresholds through the training of the BPNN. If the network reaches the training goal (completes the set number of training times or reaches the set accuracy), the training is terminated, and the accuracy of the model prediction data is output.

[0109] The specific method for calculating the nearest neighbor interpolation of the Cartesian image in Step 2 includes:

[0110]

[0111] In the formula, r is the distance between a certain point in the radar image and the center point, and θ is the angle formed with the bow direction.

[0112] The specific calculation of the characteristic parameters in Step 3 includes:

[0113] Step 3.1, calculate the wavenumber energy spectrum of the Cartesian image. Perform 2D-DFT on the sub-image z = f(x, y), and its mathematical model can be expressed as:

[0114]

[0115] In the formula, N is the number of sampling points of the sub-image z, x, y = 1, 2,..., N are frequency variables, F(k x , k y ) is the Fourier coefficient of the sub-image z, (k x , k y ) is the wavenumber domain coordinate k x and k y are the components of the wavenumber vector in the x and y directions respectively. Take the vector absolute value of F(k x , k y ) to obtain the wavenumber energy spectrum

[0116]

[0117] In the formula, Re(F(k x , ky )) and Im(F(k x , k y )) are the real and imaginary parts of F(k x , k y ).

[0118] Step 3.2: Divide the energy spectrum into three wave number segments according to the wave number segment where the ocean wave is located in Step 1.3. The proportion of each wave number segment to the total wave number segment is used as a parameter variable. That is, the energy proportions of the low, medium, and high wave number segments to the total energy can be expressed as

[0119]

[0120] In the formula, P L , P M , P H are the energy proportions of the low, medium, and high wave number segments respectively, and they are used as characteristic parameters. En L , En M , En H are the energies of the low, medium, and high wave number segments respectively, and En is the total energy of the energy spectrum.

[0121] The GA - BPNN model in Step 4 is the network model trained in Step 1.

[0122] The following gives an embodiment in combination with specific parameters:

[0123] The marine radar used in the embodiment of the present invention is a standard RM - 1290 X - band navigation radar, operating in the short - pulse mode. The operating frequency of the radar is 9320 - 9500 MHz. The rotation speed of the radar antenna is about 26 r / min, and the pulse repetition frequency is about 1300 Hz. The radial resolution of the collected radar image is 7.5 m, the azimuth resolution is 0.1°, there are 600 sampling points in the range direction, and the actual observation range is 4.5 km.

[0124] Table 1 Technical parameters of X - band marine radar

[0125]

[0126] Combined with the accompanying drawings, the specific experimental steps of the present invention are as follows:

[0127] Step 1: Obtain the GA - BPNN model applicable to rainfall intensity level inversion offline.

[0128] Step 1.1, Offline collect radar images under different rainfall intensities. According to the rainfall recorded by the rain gauge synchronously, screen 1000 radar images for each of the four rainfall intensity levels from no rain to light rain, light rain, moderate rain, heavy rain to rainstorm. Select a suitable rainfall intensity inversion area for these images. After analyzing the images, the red fan-shaped area is the sea wave area without fixed objects such as mountains and poles. Therefore, the radial distance from 1000m to 1960m and the angle starting from 135° in the radar image are selected as the rainfall intensity inversion area, as Figure 1 shown. Use the nearest neighbor interpolation method to map the pixel intensity values in the selected area to the Cartesian coordinate system to obtain the corresponding Cartesian image. The calculation method can be described as follows:

[0129]

[0130] where r is the distance between a point in the radar image and the center point, and θ is the angle formed with the bow direction. The formula indicates that the point (r,θ) in polar coordinates is mapped to (x0,y0) in Cartesian coordinates, and the obtained Cartesian image is 960m*960m, that is, 128*128 pixel points.

[0131] Step 1.2, Calculate the wave number energy spectrum. Perform a two-dimensional discrete Fourier transform (2D-DFT) on the obtained Cartesian image z = f(x,y), and its formula can be expressed as

[0132]

[0133] where N is the number of sampling points of the sub-image z, x,y = 1,2,…,N are frequency variables, F(k x ,k y ) is the Fourier coefficient of the sub-image z, (k x ,k y ) is the wave number domain coordinate, k x and k y are the components of the wave number vector in the x and y directions respectively. Take the vector absolute value of F(k x ,k y ) to obtain the wave number energy spectrum as Figure 2 shown, and its formula can be expressed as

[0134]

[0135] where Re(F(k x ,k y )) and Im(F(k x ,k y )) are the real and imaginary parts of F(k x ,k y ).

[0136] Step 1.3: Divide the wave number energy spectrum into several wave number segments, and take the energy proportion within each wave number segment as the characteristic parameter variable. Since the common wave length range of ocean waves is 50 - 200 meters, and the lower limit of the wave number of the energy spectrum of ocean waves can be expressed as

[0137]

[0138] The upper limit of the wave number of the energy spectrum of ocean waves can be expressed as

[0139]

[0140] Therefore, in this example, the low, medium, and high wave number segments are 0 - 0.0314 rad / m, 0.0314 rad / m - 0.1257 rad / m, and greater than 0.1257 rad / m respectively, as Figure 3 shown. Calculate the energy proportion of wave numbers in each wave number segment as

[0141]

[0142] In the formula, En L , En M and En H are the energies of the low, medium, and high wave number segments respectively, P L , P M and P H are the corresponding energy proportions of the wave number segments, and En is the total energy. In this calculation, P L , P M and P H are 0.469, 0.341, and 0.19 respectively, and the corresponding rainfall intensity is 0.6 mm / 10 mins. It can be seen from Table 2 that the rainfall intensity level is moderate rain.

[0143] Table 2 Classification standard of rainfall intensity level

[0144]

[0145] Step 1.4: The constructed GA - BPNN model includes two parts: the GA and BPNN algorithms, as Figure 4 shown. The GA part includes the encoding of the initial values of the neural network, the establishment of the fitness function, the selection, crossover, and mutation operations, and the update of the population. The BPNN algorithm part includes the calculation of errors and the fine - tuning of weights and thresholds to ensure the best prediction output.

[0146] Since the energy spectrum has P L , P M and P HThree characteristic parameters are used to output one of the four rainfall intensity levels (no rain to light rain, light rain, moderate rain, heavy rain to rainstorm). Therefore, the input layer of the BPNN part has 3 nodes and the output layer has 1 node. The number of hidden layer nodes can be calculated according to the empirical formula.

[0147]

[0148] Where l, n, and m represent the nodes of the hidden layer, input layer, and output layer respectively, and β is a constant between 0 and 10. After calculation, the number of hidden layer nodes is 8. Therefore, the structure of the BPNN model is 3-8-1. Here, the training function of the network is Trainlm, and the number of iterations, learning rate, and minimum training target error are 100, 0.01, and 0.001 respectively. The parameters of the genetic algorithm are: the number of iterations is 50, the population size is 10, the crossover probability is 0.3, and the mutation probability is 0.1.

[0149] The constructed GA-BPNN model is trained and verified. To avoid the impact of training, the selection of data categories should be approximately equal, and samples of different categories should be input crosswise. The 4000 groups of sample data screened as the training set and test set (3000 groups are randomly selected for the training set, and the rest are for the test set) are sent into the constructed GA-BPNN model. The main training process is as follows: BPNN generates initial random weights and thresholds, which are sent to the GA part for binary encoding and the establishment of a fitness function. New populations are generated through selection, crossover, and mutation operations. When the average fitness of the population iteration remains unchanged, the optimal weights and thresholds are obtained. They are sent back to the BPNN part to replace the initial weights and thresholds. The error is calculated by the real output and the predicted output, so as to fine-tune the weights and thresholds. When the training target meets the set number of training times or reaches the set algorithm accuracy, the network model training ends, and at this time, a neural network model that meets the algorithm requirements is obtained.

[0150] Figure 5 and Figure 6 are the statistical results of the rainfall intensity level inversion of the GA-BPNN model for this time. Among them, Figure 5 The abscissa is the generation number of evolution, and the ordinate is the fitness. The average fitness reaches 344 after 30 generations from about 413, and the GA-BPNN model tends to be stable. At this time, the GA searches out the optimal weights and thresholds. Figure 6 is the corresponding network training performance. The abscissa represents the number of training iterations, and the ordinate represents the root mean square error (MSE). As can be seen from the figure, after 9 iterations of the network with the optimal weights and thresholds replaced, the root mean square error remains basically unchanged, and the network reaches convergence. At this time, the MSE is about 0.038954.

[0151] Figures 7(a) and 7(b) show the detection results of the test set (1000 groups) under this training. Among them, Figure 7(a) is the comparison between the predicted rainfall intensity and the true rainfall intensity. The abscissa is the number of data, and the ordinate is the rainfall intensity, which is divided into no rain or light rain, light rain, moderate rain, heavy rain to rainstorm, and marked as 1, 2, 3, and 4 respectively. The red hollow dots are the true rainfall intensity levels, and the blue solid dots are the predicted rainfall intensity levels. It can be found that most of the red and blue dots coincide, which can be more intuitively reflected by Figure 7(b). After statistical analysis, the inversion accuracy of the rainfall intensity level is 97.2%. To avoid contingency, 100 model trainings are carried out, and the accuracy rates are statistically analyzed, as Figure 8 shown. It can be seen from the figure that the accuracy rates are mainly distributed between 96.5% and 99%. Therefore, the average value of these values is taken as the final accuracy rate of the algorithm, and its value is 97.4%. At this time, the GA-BPNN model meets the accuracy requirements of the experiment and can be used for the inversion of the rainfall intensity level.

[0152] Step 2: Select a radar image to be detected, select the rainfall detection area according to Step 1.1, and calculate the Cartesian image using the nearest neighbor interpolation method. The calculation method is

[0153]

[0154] where r is the distance between a certain point in the radar image and the center point, and θ is the angle formed with the bow direction. The formula indicates that the point (r, θ) in polar coordinates is mapped to (x0, y0) in Cartesian coordinates, and the obtained Cartesian image is 960m * 960m, that is, 128 * 128 pixel points.

[0155] Step 3: Calculate the wavenumber energy spectrum using 2D-DFT. Let the Cartesian image be z = f(x, y), and the calculation method is

[0156]

[0157] where N is the number of sampling points of the sub-image z, x, y = 1, 2,..., N are frequency variables, F(k x , k y ) is the Fourier coefficient of the sub-image z, (k x , k y ) is the wavenumber domain coordinate k x and k y are the components of the wavenumber vector in the x and y directions respectively. Take the vector absolute value of F(k x , k y ) to obtain the wavenumber energy spectrum As Figure 9 shown, its formula can be expressed as

[0158]

[0159] Where Re(F(k x ,k y )) and Im(F(k x ,k y )) are the real and imaginary parts of F(k x ,k y ). The energy spectrum is divided into three wavenumber bands of low, medium and high according to step 1.3, with ranges of 0 - 0.0314 rad / m, 0.0314 rad / m - 0.1257 rad / m and greater than 0.1257 rad / m respectively. The characteristic parameters are calculated using the formula for the proportion of energy in the wavenumber band

[0160]

[0161] Where, En L , En M and En H are the energies of the low, medium and high wavenumber bands respectively, P L , P M and P H are the proportions of the corresponding wavenumber band energies, and En is the total energy. The values of the three parameters P L , P M and P H are 0.378, 0.222 and 0.4 respectively.

[0162] Step 4, input the values of the three parameters into the GA - BPNN model trained in step 1.4, and the output value of the model is 2 (light rain), and the prediction result is consistent with the real rainfall intensity recorded synchronously by the rain gauge, and the prediction result is correct.

[0163] To prove the superiority of the algorithm disclosed in this patent, the RZE method based on the occluded area of the radar image is compared with the rainfall intensity level inversion algorithm in this paper. Using the same 4000 groups of data, 1000 groups of data are randomly selected from them, and the real rainfall intensities recorded synchronously are as Figure 10 shown.

[0164] To compare more clearly with the proposed method, the real rainfall intensities are marked as 1, 2, 3 and 4 in the same way as in Table 2, as Figure 11As shown. According to the method in the literature "Lu Z, Sun L, Zhou Y. A Method for Rainfall Detection and Rainfall Intensity Level Retrieval from X-Band Marine Radar Images[J]. Applied Sciences, 2021, 11(4): 1565.", the rainfall intensity is retrieved using the offline fitting formula. The retrieved rainfall intensity level is compared with the true rainfall intensity level, and the error map is calculated, as shown in Figures 12(a) and 12(b). Repeat the same experiment by randomly selecting 1000 groups from 4000 groups of data, a total of 100 times, and the results are as Figure 13 shown. As shown in the figure, the accuracy rate is mainly distributed between 82% and 87%, and the final accuracy rate is determined to be 84.5%. Therefore, the method proposed in this application improves the accuracy rate of rainfall intensity level retrieval by about 12.9% compared with the RZE method. A rainfall intensity level retrieval algorithm based on the combination of WES and GA-BPNN proposed by the present invention converts the spatial rainfall into the distribution of energy in the wavenumber domain, and uses an optimized neural network model to more accurately retrieve the rainfall intensity level. Compared with the rainfall intensity recorded synchronously by the rain gauge, this algorithm has a high retrieval accuracy.

Claims

1. A method for inverting rainfall intensity levels, characterized in that: Step 1: Obtain an offline GA-BPNN model applicable to inverting rainfall intensity levels Collect multiple radar images offline under different rainfall intensities, screen the images, select the rainfall detection area, obtain the Cartesian image using the nearest neighbor interpolation method, and then obtain the wave number energy spectrum corresponding to the Cartesian image through two-dimensional discrete Fourier transform; divide the energy spectrum into different wave number segments, and count the proportion of the wave number energy in each wave number segment under different rainfall intensities in the total energy as a characteristic parameter, and send the characteristic parameter into the constructed GA-BPNN model for training and verification to obtain a model that meets the experimental requirements; Step 2: Read the radar image to be detected, select the rainfall detection area, and calculate the Cartesian image using the same nearest neighbor interpolation method as in Step 1; Step 3: Calculate the wave number energy spectrum using the same 2D-DFT as in Step 1, and calculate the proportion of the energy in different wave number segments as a characteristic parameter; Step 4: Send the characteristic parameter into the GA-BPNN model trained in Step 1 to obtain the rainfall intensity level of the image to be detected.

2. The rainfall intensity level inversion method according to claim 1, characterized in that: The specific method for obtaining the Cartesian image using the nearest neighbor interpolation method in Step 1 is: In the formula, r is the distance between a certain point in the radar image and the center point, and θ is the angle formed with the bow direction.

3. A rainfall intensity level inversion method according to claim 1, characterized in that: The specific method for obtaining the wave number energy spectrum corresponding to the Cartesian image through two-dimensional discrete Fourier transform in Step 1 is: Step 1.1.1: Define the number of full waves appearing in a 2π length as the wave number k, satisfying: In the formula, (k x , k y ) are the coordinates in the energy spectrum domain, that is, the wave number components of k on the x and y axes, with the unit of rad / m; λ is the wavelength of the sea surface stripes presented in the radar image; Step 1.1.2: Perform 2D-DFT on the Cartesian image z = f(x, y), and its mathematical model can be expressed as: where N is the number of sampling points of the sub-image z, x, y = 1, 2, …, N are frequency variables, and F(k x , k y ) is the Fourier coefficient of the sub-image z, (k x , k y ) is the coordinate in the wavenumber domain, k x and k y are the components of the wavenumber vector in the x and y directions respectively. Taking the vector absolute value of F(k x , k y ) gives the wavenumber energy spectrum Specifically: where Re(F(k x ,k y )) and Im(F(k x ,k y )) are the real and imaginary parts of F(k x ,k y ).

4. A rainfall intensity level inversion method according to claim 3, characterized in that: The specific method for dividing the energy spectrum into different wave number segments in Step 1 and counting the proportion of the wave number energy in each wave number segment under different rainfall intensities in the total energy as a characteristic parameter is: Step 1.2.1: Use λ d to represent the minimum value of the ocean wave wavelength scale, and λ t to represent the maximum value of the ocean wave wavelength scale. Then the lower limit of the wavenumber of the ocean wave energy spectrum is: The upper limit of the wave number of the energy spectrum of ocean waves is: The wave number energy spectrum is divided into three wave number segments: low, medium, and high. Among them, the energy spectrum of the low wave number segment is: L(k x ,k y ), The energy spectrum in the medium wave number band, i.e., the wave number band of ocean waves, is as follows: M(k x ,k y ), The high wavenumber band energy spectrum is H(k x ,k y ), Step 1.2.2: The energy of each wave number segment is expressed as: where En L , En M , and En H are the energies in the low, medium, and high wavenumber bands respectively, and the total energy of the wavenumber energy spectrum is: Therefore, the proportion of the energy of the low, medium, and high wave number segments in the total energy is: Wherein, P L , P M , P H are the energy ratios of the low, medium, and high wavenumber bands respectively, and are used as characteristic parameters.

5. A rainfall intensity level inversion method according to claim 1, characterized in that: The specific method for sending the characteristic parameter into the constructed GA-BPNN model for training and verification in Step 1 is: Step 1.3.1: Initialize the BPNN model, send the characteristic parameter data into the network model, and the data is divided into a training set and a verification set, and the model generates initial weights and thresholds; Step 1.3.2: Use GA to perform binary encoding on the generated initial values, and randomly generate N individuals as the initial population; Step 1.3.3: Establish a fitness function Fitness: wherein, and respectively represent the i-th predicted output and the actual output obtained by the k-th training sample through the neural network, Nt and Ns respectively represent the number of training set samples and the number of output nodes; Step 1.3.4: Perform selection, crossover, and mutation operations on the initial population; Step 1.3.5: After several generations of selection, crossover, and mutation operations, if the calculation error reaches the allowable error, stop the calculation, and replace the initial weights and thresholds in BPNN with the calculated optimal weights and thresholds; Step 1.3.6: Adjust the weights and thresholds through the training of BPNN. If the set number of training times is completed or the set accuracy is reached, the training terminates, and the accuracy of the model prediction data is output.

Images