A sea surface wind speed retrieval method based on marine radar images

Abstract

The application discloses a sea surface wind speed inversion method based on a marine radar image and belongs to the technical field of ocean remote sensing. The application solves the problem of low wind speed inversion precision of the existing method. The technical scheme of the application is as follows: step 1, original radar image sets are obtained by reading radar files offline, each original radar image is processed to obtain the full distance integral area of the counter-wind direction corresponding to each original radar image, and then a function relationship between the wind speed and the full distance integral area of the counter-wind direction is established according to the wind speed when each original radar image is collected; step 2, the full distance integral area of the counter-wind direction of the radar image to be measured is calculated; and step 3, the full distance integral area of the counter-wind direction calculated in step 2 is substituted into the function of the wind speed and the full distance integral area of the counter-wind direction established in step 1 to obtain the sea surface wind speed when the radar image read in step 2 is collected. The method can be applied to sea surface wind speed inversion.

Description

Technical Field

[0001] This invention belongs to the field of marine remote sensing technology, specifically relating to a method for inverting sea surface wind speed based on marine radar images. Background Technology

[0002] Sea surface wind field is an important parameter for studying the interaction between the sea and the atmosphere. By regulating heat, water vapor, air-sea flux, and particulate matter, the sea surface wind field modulates the coupling between the atmosphere and the ocean, thereby maintaining global and regional climate. Currently, sea surface wind speed is mainly obtained through field sensors (such as anemometers) located on the mast of ships. However, anemometers are easily affected by atmospheric turbulence and airflow disturbances caused by the superstructure of the ship. Even if the anemometer is installed in an unobstructed location, the error in wind parameter measurement may be very high (Thornhill E, Wall A, McTavish S, et al. Ship anemometer bias management[J]. Ocean Engineering, 2020, 216: 107843.). In recent decades, X-band marine radar has become one of the effective methods for obtaining sea surface wind field information due to its advantages of high resolution and timely feedback (Wang Hui, Lu Zhizhong. Sea surface wind direction inversion algorithm based on wavenumber energy spectrum[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2014, 42(12): 96-100.). Furthermore, since most ships are already equipped with X-band marine radar, which has a relatively small additional cost compared to other remote sensing instruments, these studies have laid the foundation for marine radar to extract sea surface wind field information.

[0003] Currently, there are two types of methods for retrieving wind speed from X-band marine radar image sequences. The first type of method retrieves wind speed by establishing an empirical model function that correlates a reference wind speed with backscattering intensity information in the radar image. The total variance of the background noise of the radar image was compared with the reference wind speed, and a linear function between the two was established (Hatten H, Seemann J, Horstmann J, et al. Azimuthal dependence of the radar cross section and the spectral background noise of a nautical radar at grazing incidence[C] / / IGARSS'98. Sensing and Managing the Environment. 1998 IEEE International Geoscience and Remote Sensing. Symposium Proceedings. (Cat. No. 98 CH36174). IEEE, 1998, 5: 2490-2492. and Izquierdo P, Soares C G. Analysis of sea waves and wind from X-band radar[J]. Ocean Engineering, 2005, 32(11-12): 1404-1419.). In 2012, Lund discovered a cubic polynomial relationship between echo intensity and reference wind speed through the study of marine radar images. He then established a cubic polynomial function between the two using least squares fitting to invert wind speed (Lund B, Graber HC, Romeiser R. Wind Retrieval From Shipborne Nautical X-Band Radar Data[J].IEEE Transactions on Geoscience & Remote Sensing, 2012, 50(10): 3800-3811.).Subsequently, in 2013, Vicen-Bueno inverted wind speed from a third-order polynomial continuous function, which depends on the maximum range distance of a pre-selected intensity value (Vicen-Bueno R, Horstmann J, Terril E, et al. Real-time ocean wind vector retrieval from marine radar image sequences acquired at grazing angle[J]. Journal of Atmospheric and Oceanic Technology, 2013, 30(1): 127-139.). For Decca and Furono radars, in 2014, Liu proposed using a hyperbolic fitting method for quadratic fitting to extract sea surface wind speed information from measured marine radar data (Liu Y, Huang W, Gill EW, et al. Comparison of algorithms for wind parameters extraction from shipborne X-band marine radar images[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 8(2): 896-906.). In 2015, Chen et al. proposed an empirical wind speed function for shore-based radar data. They first established the relationship between radar cross section (RCS) and wind speed using the probability distribution function method, and then added the significant wave height parameter to improve the inversion accuracy. However, this model is only applicable to wind and wave conditions (Chen Z, He Y, Zhang B, et al. Determination of nearshore sea surface wind vector from marine X-band radar images[J]. Ocean Engineering, 2015, 96(mar.1):79-85.). Furthermore, different radar models require modeling of sea surface wind and RCS. Therefore, the model function method is not universally applicable.Another method for deriving wind speed from X-band ocean radar image sequences uses neural networks (Dankert H, Horstmann J, Rosenthal W. Ocean wind fields retrieved from radar-image sequences[J]. Journal of Geophysical Research: Oceans, 2003, 108(C11). and Dankert H, Horstmann J. Amarine radar wind sensor[J]. Journal of Atmospheric and Oceanic Technology, 2007, 24(9):1629-1642.). This method parameterizes the relationship between wind speed and normalized radar cross section, and incorporates sea state information and atmospheric parameters as inputs to obtain a model. However, neural networks are physically undefined, therefore they cannot provide a definitive model for radar imaging mechanisms. Furthermore, the learning process of neural networks converges slowly, and the number of hidden layers and nodes is difficult to determine, making the network unsuitable for certain radar data and difficult to meet engineering requirements.

[0004] Existing methods for wind speed retrieval using marine radar imagery generally suffer from drawbacks, requiring highly detailed radar data. A 360° unobstructed radar image is necessary to obtain radar backscatter information correlated with wind speed. Furthermore, radar backscatter information is affected by stationary targets, which can influence the backscatter information of the selected area, thus reducing the accuracy of wind speed retrieval.

[0005] Furthermore, radar echo intensity is modulated by distance. The distance dependence of echo intensity on the ocean surface can be shown in radar images. Radar echo intensity decreases with increasing distance. Not only does the echo signal from a distant sea surface decrease with increasing distance, but there are also more pixels with zero intensity values. This makes long-distance backscattering information unusable. However, long-distance backscattering information is more sensitive under high wind speeds, while short-distance backscattering information is more sensitive under low wind speeds. Therefore, backscattering information at all distances needs to be considered. Studies have shown that any radial distance dependence in radar images is the result of the grazing angle dependence of ocean surface echo intensity. The echo intensity and radial distance have a cubic attenuation relationship (Lund B, Graber HC, Romeiser R. Wind Retrieval From Shipborne Nautical X-Band Radar Data[J].IEEE Transactions on Geoscience & Remote Sensing,2012,50(10):3800-3811.). Therefore, when retrieving wind information from marine radar images, the effect of range attenuation cannot be ignored.

[0006] In summary, due to the presence of obstructions in actual radar data, the influence of fixed targets on radar backscattering information, and the modulation of radar echo intensity by distance, the existing methods for wind speed inversion still suffer from low accuracy. Summary of the Invention

[0007] The purpose of this invention is to solve the problem of low wind speed inversion accuracy in existing methods, and to propose a sea surface wind speed inversion method based on marine radar images.

[0008] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0009] A method for inverting sea surface wind speed based on marine radar images, the method specifically includes the following steps:

[0010] Step 1: Read the radar files offline to obtain the original radar image set, and process each original radar image to obtain the full range integral area in the upwind direction corresponding to each original radar image;

[0011] Then, based on the wind speed when each original radar image was acquired, a functional relationship between wind speed and the area of ​​the entire range of the upwind direction was established.

[0012] Step 2: Read the radar image of the wind speed to be measured online and calculate the area of ​​the radar image of the wind speed to be measured over the entire distance in the upwind direction.

[0013] Step 3: Substitute the total area of ​​the upwind direction calculated in Step 2 into the function of wind speed and total area of ​​the upwind direction established in Step 1 to obtain the sea surface wind speed when the radar image was acquired in Step 2.

[0014] The beneficial effects of this invention are:

[0015] This invention utilizes the range attenuation of echo intensity from marine radar images to invert wind speed, overcoming the impact of range attenuation on wind speed retrieval. It overcomes interference from fixed targets and outliers on the sea surface by employing histogram statistics, making it applicable to sea areas affected by fixed objects. The attenuation function value in the upwind direction is calculated based on the dependence of wind direction and azimuth, and then the area of ​​the upwind direction is calculated using the attenuation function value to invert wind speed. This overcomes the impact of radar image field-of-view obstruction on wind speed retrieval, allowing for the inversion of sea surface wind speed from obstructed radar images.

[0016] Experiments show that the method of this invention can significantly improve the accuracy of wind speed inversion, with a root mean square error of only 2.4 m / s. Attached Figure Description

[0017] Figure 1 This is a flowchart of a sea surface wind speed inversion method based on marine radar images according to the present invention;

[0018] Figure 2 This is the complete original radar image in polar coordinates.

[0019] Figure 3 This is the normalized original radar image in polar coordinates.

[0020] Figure 4 The histogram of echo intensity at a distance of 7.5m is shown.

[0021] Figure 5 This represents one-dimensional ideal attenuation data from the original radar image.

[0022] Figure 6 This is the ideal attenuation model fitted based on the ideal attenuation data;

[0023] Figure 7 This is the cosine function fitted based on the attenuation level component;

[0024] Figure 8 This is a curve showing the area integral over the entire distance in the headwind direction.

[0025] Figure 9 The graph shows the wind speed as a linear function of the area integraled over the entire distance in the headwind direction.

[0026] Figure 10This is a comparison chart of wind speeds measured by two different methods and a reference wind speed. Detailed Implementation

[0027] The present application will now be described in further detail with reference to specific embodiments and accompanying drawings. Obviously, the described embodiments are merely a part of the embodiments of the present invention, and not all of them. Other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are all within the scope of protection of the present invention.

[0028] Specific Implementation Method 1: Combination Figure 1 This embodiment describes a method for retrieving sea surface wind speed based on marine radar images. The method specifically includes the following steps:

[0029] Step 1: Read the radar files offline to obtain the original radar image set, and process each original radar image to obtain the full range integral area in the upwind direction corresponding to each original radar image;

[0030] Then, based on the wind speed when each original radar image was acquired, a functional relationship between wind speed and the area of ​​the entire range of the upwind direction was established.

[0031] Step 2: Read the radar image of the wind speed to be measured online, and calculate the full-range integral area of ​​the radar image of the wind speed to be measured in the upwind direction using the method in Step 1.

[0032] Step 3: Substitute the total area of ​​the upwind direction calculated in Step 2 into the function of wind speed and total area of ​​the upwind direction established in Step 1 to obtain the sea surface wind speed when the radar image was acquired in Step 2.

[0033] This invention utilizes the dependence of echo intensity in radar images on the radial distance from the sea surface to derive an ideal attenuation model for radar images. After excluding obstructed areas, the attenuation function for each azimuth angle is obtained through the ideal attenuation model, and the horizontal attenuation component, i.e., the relative average echo intensity, is calculated. The azimuth dependence of the grayscale level is used to invert the wind direction. The attenuation model for the upwind direction is obtained through the horizontal attenuation component, and a cubic polynomial relationship between the integral area of ​​the attenuation model and the reference wind speed is established to invert the wind speed.

[0034] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that the specific process of step 1 is as follows:

[0035] Step 1.1: Read the radar file offline to obtain the original radar image, and record the azimuth, radial distance and echo intensity information of the original radar image. After performing co-channel interference suppression processing on the recorded echo intensity information, normalize the processed echo intensity information.

[0036] Step 1.2: Based on the normalized echo intensity information, calculate the area of ​​the whole distance integral in the headwind direction;

[0037] Step 1.3: Repeat steps 1.1 to 1.2 to obtain the full-range integrated area in the upwind direction corresponding to H original radar images;

[0038] The least squares fitting method is used to establish a functional relationship between the wind speed and the area of ​​the whole distance integral in the upwind direction when the original radar image is acquired.

[0039] The other steps and parameters are the same as in Specific Implementation Method 1.

[0040] Specific Implementation Method 3: This implementation method differs from Specific Implementation Method 1 or 2 in that the co-frequency interference suppression processing of the recorded echo intensity information is performed using a filtering algorithm.

[0041] Other steps and parameters are the same as in specific implementation method one or two.

[0042] Specific Implementation Method Four: This implementation method differs from Specific Implementation Methods One to Three in that the specific process of step 1.2 is as follows:

[0043] Step 1.2.1: Calculate the threshold T in histogram statistics:

[0044] T = k × M

[0045] Where: M is the number of azimuth lines in the original radar image, that is, the number of pixels at the same radial distance in the polar coordinate system of the original radar image, and k is an empirical parameter;

[0046] Step 1.2.2: In the polar coordinate system of the original radar image, select a group of pixels that are equidistant from the center of the original radar image (i.e., the position of the radar antenna in the image). Based on the normalized echo intensity of the selected pixels, store the selected pixels in m histogram panes with echo intensity ranging from 0 to 1.

[0047] If the number of pixels in a histogram pane is less than the threshold T, the normalized echo intensity of the pixels in that histogram pane is considered an outlier and is excluded; if the number of pixels in a histogram pane is not less than the threshold T, the normalized echo intensity of the pixels in that histogram pane is considered a normal value and is retained.

[0048] The largest echo intensity value is selected from the remaining normalized echo intensity values, and the selected echo intensity value is used as the echo intensity value at the current radial distance (i.e., the radial distance from the selected pixel to the center of the original radar image under the current condition).

[0049] Step 1.2.3: Repeat step 1.2.2 until the radial distance values ​​from each pixel to the center of the original radar image have been traversed, and the echo intensity values ​​at each radial distance are obtained respectively.

[0050] Step 1.2.4: Apply least squares fitting to fit the echo intensity values ​​at each radial distance to the ideal distance attenuation model, and determine the regression parameters of the ideal distance attenuation model;

[0051] Step 1.2.5: Calculate the attenuation horizontal component C in the upwind direction. wind ;

[0052] Step 1.2.6: Attenuation of the horizontal component C based on the headwind direction wind The attenuation function value in the upwind direction is calculated using the regression parameters of the ideal distance attenuation model;

[0053] Then, calculate the total area S along the headwind direction based on the attenuation function value in the headwind direction. wind .

[0054] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0055] Specific Implementation Method Five: This implementation method differs from Specific Implementation Methods One to Four in that the selection of the largest echo intensity value from the retained echo intensity values ​​specifically involves:

[0056]

[0057] Where r is the radial distance from the radar antenna to the sea surface, i.e., the radial distance from the selected pixel to the center of the original radar image, {X(r, θ)} is the set of normalized echo intensity values ​​retained at the current radial distance r, and θ is the azimuth angle of the original radar image. ra (r) is the maximum echo intensity value selected at the current radial distance r.

[0058] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0059] Specific Implementation Method Six: This implementation method differs from Specific Implementation Methods One to Five in that the specific process of step 1.2.4 is as follows:

[0060] The formula for the ideal distance decay model D(r) is as follows:

[0061]

[0062] Where: b0 and b1 are regression parameters, and D(r) is the ideal distance attenuation function model value corresponding to radial distance r, that is, the echo intensity fitting value at radial distance r.

[0063] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0064] Specific Implementation Method Seven: This implementation method differs from Specific Implementation Methods One to Six in that the specific process of step 1.2.5 is as follows:

[0065] Step 1) Initialize the weights and threshold δ:

[0066] The initial weights are:

[0067]

[0068] Where, r n ω(r) is the radial distance from the radar antenna to the nth discrete point on the sea surface, where n = 1, 2, 3, ..., p, and p is the number of discrete points corresponding to each azimuth angle. n ) is the initial weight corresponding to the nth discrete point;

[0069] Initialize the value of the threshold δ to 1;

[0070] Step 2) Initialize the azimuth angle θ to 0°;

[0071] Step 3) In the original radar image, select the normalized echo intensity of all discrete points (each pixel corresponds to a discrete point) at the azimuth angle θ, and reset the weight of the discrete points with a normalized echo intensity of 0 to zero.

[0072] The attenuation level component C at the current azimuth angle θ is calculated by minimizing the following formula using a nonlinear least squares method. θ :

[0073]

[0074] Where X(r) n ,θ) is the radial distance r at the azimuth angle θ. n Normalized echo intensity at point D(r) n () is the radial distance r at the azimuth angle θ. n The echo intensity fitting value at the location; r n = n×l, where l is the radar's range resolution;

[0075] Step 4), based on the attenuation level component C calculated in Step 3). θ After updating the weights, update the threshold δ, and then use the updated weights and threshold δ to return to step 3) to calculate the updated attenuation level component C. θ ;

[0076] Based on the updated attenuation level component C θCalculate the final weights and update the threshold δ again to obtain the final threshold δ. Calculate the final attenuation level component C based on the final weights and the final threshold δ. θ ;

[0077] Step 5) Define the clockwise direction as the direction of azimuth increase (that is, starting from 0°, the value of the azimuth in the clockwise direction gradually increases), let the azimuth θ = Δθ + θ, where Δθ is the azimuth increment, and return to step 3);

[0078] The process continues until the azimuth angle θ reaches 2π, at which point the final attenuation horizontal component C at each azimuth angle is obtained. θ ;

[0079] Step 6) Apply the least squares method to the final attenuation level component C θ As a function of azimuth, it is fitted to the cosine function:

[0080] C θ =a0+a1 cos 2 (0.5(θ-w))

[0081] Where a0 and a1 are regression parameters, and w is the angle of the peak point of the fitted function;

[0082] Step 7) Calculate the attenuation horizontal component C in the headwind direction based on regression parameters a0 and a1. wind :

[0083] C wind =a0+a1

[0084] The other steps and parameters are the same as those in one of the specific implementation methods one to six.

[0085] Specific Implementation Method Eight: This implementation method differs from Specific Implementation Methods One to Seven in that the specific process of step 1.2.6 is as follows:

[0086]

[0087] Where, D(r) wind It is the attenuation function value in the direction of the headwind;

[0088]

[0089] Among them, S wind It is the area of ​​the whole distance integral in the headwind direction, r max It is the maximum detection range of the radar.

[0090] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0091] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One through Eight in that the specific process of establishing the functional relationship between wind speed and the area of ​​the entire range of the headwind direction when acquiring the original radar image is as follows:

[0092] V=β1S wind +β0

[0093] Where V is the wind speed, and β0 and β1 are regression parameters.

[0094] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0095] Specific Implementation Method Ten: This implementation method differs from Specific Implementation Methods One to Nine in that the attenuation level component C calculated in step 3) is... θ The specific process for updating the weights is as follows:

[0096]

[0097] Where, ω′(r n ) is the updated weight.

[0098] Based on the updated attenuation level component C θ The method for calculating the final weights is the same as that in this embodiment.

[0099] The threshold δ is updated as follows: the value of the threshold δ after each update is half of the threshold δ used in the previous execution of step 3). That is, the initial threshold δ = 1, the value of the threshold δ after the first update is 0.5, and the final value of the threshold δ is 0.25.

[0100] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0101] Example

[0102] The method for wind speed inversion of marine radar images based on radar image range attenuation, proposed in this invention, will be further described in detail below with reference to the accompanying drawings.

[0103] Implementation process as follows Figure 1 As shown, the process can be divided into the following steps: First, establish the linear function relationship between wind speed and the area of ​​the entire range in the upwind direction offline; Second, read the radar file online to obtain the original radar image, remove co-channel interference and normalize it, and calculate the area of ​​the entire range in the upwind direction; Third, calculate the sea surface wind speed.

[0104] The embodiment of this invention uses an X-band marine radar, which is installed on the forward side of the main mast of the marine surveillance vessel. It acquires data using short pulses, monitoring a radial distance range of 4.2 km, with an angular resolution of approximately 1 degree. Each image acquisition time is approximately 2.7 seconds, and 32 images are stored as a time series. Each radar image has approximately 3000 lines, with 560 points on each line, and a range resolution of 7.5 m. The wind reference data used in this invention comes from the anemometer data of the ship during sea trials. Similar to the marine radar, this anemometer is also installed on the main mast, and wind parameters are recorded in minutes to verify the accuracy of the radar-derived wind speed.

[0105] Combination Figures 1 to 10 The specific implementation steps of this invention are as follows:

[0106] The first step is to establish an offline linear function relationship between wind speed (wind speed reference data comes from the anemometer data of the ship during the sea trial) and the integral area over the entire distance in the upwind direction.

[0107] Step 1.1: Offline reading of radar files to obtain the original radar image, removal of co-channel interference, and normalization processing. Step 1.1 includes the following steps:

[0108] Step 1.1.1: Load the radar file using a radar image processing program and record information such as the radar image acquisition time, azimuth, radial distance, and echo intensity.

[0109] Step 1.1.2: Perform co-channel interference suppression processing on the original radar image using the selected filtering algorithm. In this embodiment, median filtering is selected for co-channel interference suppression. Specifically, the echo intensity value of each pixel is replaced by the median of the echo intensities of the other 8 pixels within a 3×3 neighborhood window of that pixel.

[0110] Step 1.1.3: Normalize the radar image data. In this embodiment, the radar image data is generated by digitizing the radar backscattering information and storing it as a 14-bit grayscale depth image sequence. That is, the digitized backscattering intensity range is 0 to 8192. Therefore, the intensity of each pixel is divided by 8192 to normalize it to an intensity range of 0 to 1. Figure 2 The image is an unnormalized raw radar image. Figure 3 The original radar image after removing co-channel interference and normalization.

[0111] Step 1.2: Calculate the area of ​​the entire range in the upwind direction. In the polar coordinate system of the radar image obtained in Step 1.1, ideal attenuation data is obtained through histogram statistics, and least squares fitting is applied to fit the attenuation data to the ideal range attenuation model to determine the regression parameters. Then, the attenuation horizontal component of each azimuth is determined by fitting and comparing the ideal range attenuation model with the intensity of radar image pixels at each azimuth. The dependence of wind direction and azimuth on the maritime radar image under horizontal polarization is used to obtain the fitting curve and determine the attenuation horizontal component in the upwind direction. Finally, the attenuation function in the upwind direction is determined and the area of ​​the entire range in the upwind direction is calculated.

[0112] Step 1.2 specifically includes the following steps:

[0113] Step 1.2.1: Calculate the threshold T in histogram statistics using formula (1). In this invention, 360 angle data points were uniformly selected from 3000 angle data points, M is 360, and the empirical parameter k is selected as 0.01. Therefore, the threshold T = 3.6.

[0114] The formula for calculating the threshold T is:

[0115] T=k×M (1)

[0116] In the formula:

[0117] M represents the number of azimuth lines in the radar image, that is, the number of pixels at the same distance;

[0118] k represents an empirical parameter;

[0119] Step 1.2.2: Starting from the center of the radar image, select one-dimensional data points with the same distance and store them in m histogram panes with an intensity range of 0 to 1. In this embodiment, m is selected as 256, where panes with a number less than the threshold T are considered outliers and excluded. Figure 4 As shown, this is a histogram intensity value statistics at a distance of 7.5m. Pixels containing panes with a count less than the threshold of 3.6 are considered outliers and excluded.

[0120] Step 1.2.3: For pixels with the same distance r determined in Step 1.2.2, select the pixel with the highest intensity value as the value at this distance. Figure 4 We selected an intensity of 0.1521 as the intensity value at 7.5m. Then, using the image center as the origin, we repeated steps 1.2.2 to 1.2.3 at larger distances until the farthest distance was reached, thus determining the ideal distance attenuation data, as shown below. Figure 5 The image shows the final extracted ideal distance attenuation data.

[0121] The highest intensity value X at a distance r raThe formula for calculating (r) is as follows:

[0122] X ra (r)=max θ {X(r, θ)}, θ∈[0, 2π) (2)

[0123] In the formula:

[0124] r represents the radial distance from the radar antenna to the sea surface, with a range of 0-4.2 km.

[0125] θ represents the azimuth angle of the radar image;

[0126] X represents the normalized radar image intensity;

[0127] After r has traversed all distance values, the ideal distance decay data is obtained.

[0128] Step 1.2.4: Select an existing ideal distance decay function model D(r), and apply least squares fitting to fit the ideal distance decay data to the ideal distance decay model, thereby determining the regression parameters b0 = 0.1714 and b1 = 0.8288 for D(r). The curve of the fitted ideal distance decay model is shown below. Figure 6 As shown.

[0129] The formula for the ideal distance decay model D(r) is as follows:

[0130]

[0131] In the formula:

[0132] r represents the radial distance from the radar antenna to the sea surface;

[0133] b0 represents the regression parameter;

[0134] b1 represents the regression parameters;

[0135] Step 1.2.5: Initialize the weighting function ω and the threshold δ = 1.

[0136] ω(r n The initialization calculation formula for ) is as follows:

[0137]

[0138] In the formula:

[0139] r n r is the radial distance from the radar antenna to the nth discrete point on the sea surface. n = n×l, where l is the radar's range resolution, which is 7.5m in this embodiment;

[0140] n = 1, 2, 3, ..., p, where p is the number of discrete points corresponding to each azimuth angle, which is 560 in this embodiment;

[0141] Step 1.2.6: Starting from 0° and moving clockwise to 2π in the radar image, sequentially select one-dimensional data X(r, θ) at each azimuth angle. Then, assign weights ω(r, θ) to the echo intensity values ​​of 0 in X(r, θ). n Set the value to zero, and calculate the attenuation horizontal component C for each azimuth angle by minimizing the following formula using a nonlinear least squares method. θ .

[0142] C θ The calculation formula is as follows:

[0143]

[0144] Where n represents the discrete point position of the radial distance;

[0145] Step 1.2.7: Preliminarily estimate C based on Step 1.2.6. θ Then use iterative constraints on C θ The values ​​are optimized by updating the weighting function ω and the threshold δ during the iteration process, where ω is one of 560 weight values ​​ranging from 0km to 4.2km. n ) is one of the weight values. The weight ω is updated using formula (6), and the threshold δ is updated to half of its original value after each iteration. The updated parameters ω and δ are input into formula (5) for the next iteration. The entire process requires 3 iterations to obtain the final attenuation level component C. θ The attenuation horizontal component C at each azimuth angle θ It exhibits a cosine function dependence on the radar illumination direction, with a peak value in the windward direction. For example... Figure 7 As shown, 170° to 330° is the fixed obstruction area of ​​the ship's mast.

[0146] The update formula for the weighted function ω is as follows:

[0147]

[0148] Step 1.2.8: Apply the least squares method to reduce the attenuation level component C. θ The azimuth angle is fitted to a cosine function. First, a fixed obstruction region is removed. The data used in this invention exhibits an approximately 160° obstruction region between 170° and 330°, therefore data from this region is excluded. Then, the attenuation level component C from the removed obstruction region is... θ Fit to formula (7), as follows Figure 7As shown, the parameters obtained are a0 = 0.3772, a1 = -0.0405, and w = -24.5185.

[0149] The cosine function is as follows:

[0150] C θ =a0+a1 cos 2 (0.5(θ-w)) (7)

[0151] θ represents the azimuth angle of the radar image;

[0152] a0 represents the regression parameter;

[0153] a1 represents the regression parameter;

[0154] w represents the angle of the peak point of the fitted function;

[0155] Step 1.2.9: Calculate the attenuation horizontal component C in the headwind direction based on the fitting function obtained from formula (7) in step 1.2.8. wind =0.3772-0.0405=0.3367.

[0156] The calculation formula is as follows:

[0157] C wind =a0+a1 (8)

[0158] In the formula:

[0159] C wind The horizontal component representing the attenuation in the headwind direction;

[0160] a0 represents the regression parameter;

[0161] a1 represents the regression parameter;

[0162] Step 1.2.10: Attenuate the horizontal component C in the headwind direction. wind Substituting into formula (9), we obtain the attenuation function D(r) in the headwind direction. wind C wind =0.3367, b0=0.1714, b1=0.8288, as Figure 8 As shown.

[0163] The calculation formula is as follows:

[0164]

[0165] In the formula:

[0166] D(r) wind The attenuation function representing the direction of headwind;

[0167] C windThe horizontal component representing the attenuation in the headwind direction;

[0168] b0 represents the regression parameter;

[0169] b1 represents the regression parameters;

[0170] Step 1.2.11: Calculate the total area S along the upwind direction. wind In this embodiment, the distance is 0-4.2km, then r max = 4.2km, S is calculated wind =0.284.

[0171] The calculation formula is as follows:

[0172]

[0173] D(r) wind The attenuation function representing the direction of headwind;

[0174] S wind The area represented by the total distance integral in the direction of headwind;

[0175] r max This represents the maximum detection range of the radar;

[0176] r represents the radial distance from the radar antenna to the sea surface;

[0177] Step 1.3: Establish a linear function relationship between wind speed and the area of ​​the whole distance integral in the headwind direction.

[0178] Step 1.3.1: Repeat steps 1.1 and 1.2 to obtain the integral area of ​​H radar image data. In this invention, H = 2000.

[0179] Step 1.3.2, as follows Figure 9 As shown, the reference wind speed V is established using least squares fitting as the integration area S. wind The result is a linear function, yielding β1 = 76.35 and β0 = -10.23.

[0180] The calculation formula is as follows:

[0181] V=β1S wind +β0 (11)

[0182] V represents wind speed;

[0183] S wind The area represented by the total distance integral in the direction of headwind;

[0184] β1 represents the regression parameter;

[0185] β0 represents the regression parameter;

[0186] The second step involves online reading of the radar file to obtain the original radar image, removing co-channel interference, and normalizing it, using the same method as in step 1.2 for offline reading of the radar file to obtain the original radar image. The area integral over the entire range in the upwind direction is then calculated, using the same method as in step 1.2 for offline calculation of the area integral over the entire range in the upwind direction.

[0187] Step 3: Calculate the sea surface wind speed. This includes the following steps:

[0188] The S of the desired radar image wind Substitute the fitted formula (11) from step 1.3.2 into the formula to calculate the sea surface wind speed.

[0189] The sea surface wind speed obtained in this embodiment is: V = 76.35 × 0.284 - 10.23 = 11.45 m / s, which differs from the reference wind speed of 12.1 m / s by 0.65 m / s.

[0190] The performance verification of the algorithm of this invention was conducted based on measurement data from sea trials. During the experimental test, the radar's monitoring range was within 4.2 km, the acquisition time for each image was approximately 2.7 seconds, and 32 images were stored as a time series. The radar image bus had approximately 3000 lines, with 560 points on each line, and a range resolution of 7.5 m. The wind reference data used in this invention came from the anemometer data from the sea trials. Similar to the marine radar, this anemometer was also installed on the ship's main mast, and wind parameters were recorded in minutes to verify the accuracy of the radar's wind field inversion. This invention randomly selected 2000 sets of data from March 7th to March 8th, 2012 to establish the relationship between wind speed and radar images, and selected 15914 sets of data from March 7th to March 8th, 2012 to verify the accuracy of the invention.

[0191] To facilitate better analysis, the method of this invention was compared with the single-curve fitting method that directly superimposes and averages the echo intensity at each azimuth angle, such as... Figure 10 As shown in Table 1, the overall results are as follows:

[0192] Table 1 Overall Calculation Results of Wind Speed

[0193]

[0194] As shown in Table 1, the wind speed inversion method based on the distance attenuation method of this invention has a better effect than the single curve fitting method of direct averaging, and can effectively invert the wind speed of radar images.

[0195] The method for wind speed inversion from marine radar images based on the horizontal component of radar image range attenuation proposed in this invention improves the accuracy of wind speed inversion. It not only overcomes the influence of sea state on radar echoes, but also makes full use of the range attenuation characteristics of radar images, resulting in a significant improvement in the accuracy of wind speed inversion.

[0196] The above examples of this invention are merely illustrative of the computational model and process of this invention, and are not intended to limit the implementation of this invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of this invention are still within the scope of protection of this invention.

Claims

1. A method for inverting sea surface wind speed based on marine radar images, characterized in that, The method specifically includes the following steps: Step 1: Read the radar files offline to obtain the original radar image set, and process each original radar image to obtain the full range integral area in the upwind direction corresponding to each original radar image; Then, based on the wind speed when each original radar image was acquired, a functional relationship between wind speed and the area of ​​the entire range of the upwind direction was established. The specific process of step 1 is as follows: Step 1.1: Read the radar file offline to obtain the original radar image, and record the azimuth, radial distance and echo intensity information of the original radar image. After performing co-channel interference suppression processing on the recorded echo intensity information, normalize the processed echo intensity information. Step 1.2: Based on the normalized echo intensity information, calculate the area of ​​the whole distance integration in the headwind direction; The specific process of step 1.2 is as follows: Step 1.2.1: Calculate the threshold T in histogram statistics: Where: M is the number of azimuth lines in the original radar image, that is, the number of pixels at the same radial distance in the polar coordinate system of the original radar image, and k is an empirical parameter; Step 1.2.2: In the polar coordinate system of the original radar image, select a group of pixels that are equidistant from the center of the original radar image. Based on the normalized echo intensity of the selected pixels, store the selected pixels in m histogram panes with echo intensity ranging from 0 to 1. If the number of pixels in a histogram pane is less than the threshold T, the normalized echo intensity of the pixels in that histogram pane is considered an outlier and is excluded; if the number of pixels in a histogram pane is not less than the threshold T, the normalized echo intensity of the pixels in that histogram pane is considered a normal value and is retained. The largest echo intensity value is selected from the remaining normalized echo intensity values, and the selected echo intensity value is used as the echo intensity value at the current radial distance. The largest echo intensity value is selected from the retained normalized echo intensity values, specifically: in, It is the radial distance from the radar antenna to the sea surface, that is, the radial distance from the selected pixel to the center of the original radar image. At the current radial distance The set of normalized echo intensity values ​​that are retained below It is the azimuth angle of the original radar image. It is the selected current radial distance The maximum echo intensity value; Step 1.2.3: Repeat step 1.2.2 until the radial distance values ​​from each pixel to the center of the original radar image have been traversed, and the echo intensity values ​​at each radial distance are obtained respectively. Step 1.2.4: Fit the echo intensity values ​​at each radial distance to the ideal distance attenuation model and determine the regression parameters of the ideal distance attenuation model; Step 1.2.5: Calculate the attenuation horizontal component in the upwind direction. ; The specific process of step 1.2.5 is as follows: Step 1) Initialize weights and thresholds : The initial weights are: in, It is the radial distance from the radar antenna to the nth discrete point on the sea surface. , It is the number of discrete points corresponding to each azimuth angle. It is the initial weight corresponding to the nth discrete point; threshold The value is initialized to 1; Step 2) Initialize azimuth angle It is 0°; Step 3) Select the azimuth angle from the original radar image. The normalized echo intensity of all discrete points on the graph is calculated, and the weights of discrete points with a normalized echo intensity of 0 are reset to zero. Calculate the current azimuth angle The attenuation level component below : in, It is the azimuth angle. Up, radial distance Normalized echo intensity at the location, It is the azimuth angle. Above, radial distance The fitted value of the echo intensity at that location; Step 4) Based on the attenuation level component calculated in Step 3) After updating the weights, then adjust the threshold. Perform an update, utilizing the updated weights and thresholds. Return to step 3) Calculate the updated attenuation level component. ; Based on the updated attenuation level components Calculate the final weights and update the thresholds again. Obtain the final threshold Based on the final weight and the final threshold Calculate the final attenuation level component ; Step 5) Define the clockwise direction as the azimuth increase direction, and set the azimuth angle... , (Increment the azimuth angle and return to step 3). Until the azimuth Reaching 2 At the end of the time, the final attenuation horizontal component at each azimuth angle is obtained. ; Step 6) Convert the final attenuation level component As a function of azimuth, it is fitted to the cosine function: in, and These are regression parameters. It is the angle of the peak point of the fitted function; Step 7) Based on regression parameters and Calculate the attenuation horizontal component in the upwind direction : Step 1.2.6: Attenuation of the horizontal component based on the headwind direction The attenuation function value in the upwind direction is calculated using the regression parameters of the ideal distance attenuation model; Then, calculate the total area of ​​the upwind integral based on the attenuation function value in the upwind direction. ; The specific process of step 1.2.6 is as follows: in, It is the attenuation function value in the direction of the headwind. and These are regression parameters; in, It is the area integral over the entire distance in the headwind direction. This is the maximum detection range of the radar; Step 1.3: Repeat steps 1.1 to 1.2 to obtain the full-range integrated area in the upwind direction corresponding to H original radar images; The least squares fitting method was used to establish a functional relationship between the wind speed and the area of ​​the whole distance integral in the upwind direction when the original radar image was acquired. Step 2: Read the radar image of the wind speed to be measured online and calculate the area of ​​the radar image of the wind speed to be measured over the entire distance in the upwind direction. Step 3: Substitute the total area of ​​the upwind direction calculated in Step 2 into the function of wind speed and total area of ​​the upwind direction established in Step 1 to obtain the sea surface wind speed when the radar image was acquired in Step 2.

2. The method for inverting sea surface wind speed based on marine radar images according to claim 1, characterized in that, The co-frequency interference suppression processing of the recorded echo intensity information uses a filtering algorithm.

3. The method for inverting sea surface wind speed based on marine radar images according to claim 2, characterized in that, The process of selecting the largest echo intensity value from the retained normalized echo intensity values ​​specifically involves: in, It is the radial distance from the radar antenna to the sea surface, that is, the radial distance from the selected pixel to the center of the original radar image. At the current radial distance The set of normalized echo intensity values ​​that are retained below It is the azimuth angle of the original radar image. It is the selected current radial distance The maximum echo intensity value.

4. The method for inverting sea surface wind speed based on marine radar images according to claim 3, characterized in that, The specific process of step 1.2.4 is as follows: in: and These are regression parameters. Radial distance The corresponding ideal distance decay function model value, i.e., radial distance The fitted value of the echo intensity at that location.

5. The method for inverting sea surface wind speed based on marine radar images according to claim 4, characterized in that, The specific process of step 1.2.6 is as follows: in, It is the attenuation function value in the direction of the headwind; in, It is the area integral over the entire distance in the headwind direction. It is the maximum detection range of the radar.

6. The method for inverting sea surface wind speed based on marine radar images according to claim 5, characterized in that, The specific process for establishing the functional relationship between wind speed and the area of ​​the entire range in the upwind direction when acquiring the original radar image is as follows: + in, It's wind speed. and These are regression parameters.

7. The method for inverting sea surface wind speed based on marine radar images according to claim 6, characterized in that, The attenuation level component calculated based on step 3) The specific process for updating the weights is as follows: in, This is the updated weight.

Images