CN102662164B - Sea surface current information extraction method based on X-band radar image and particle swarm optimization - Google Patents

Abstract

The invention provides a sea surface current information extraction method based on an X-band radar image and particle swarm optimization, belonging to the technical field of sea wave parameter inversion. The method firstly initializes the parameters such as, positions and speed of the particles and the like, and provides a ship speed initialization method and an flow estimating initialization method to determine initial positions of the particles; and then, takes the initial positions as optimal positions where the particles pass by, based on calculating the value of the fitness function of each particle, selects the position of a particle with a minimum value as the current optimal position of a swarm; and then, updates the speed and the positions of the particles in the next generation, and calculates the values of the fitness function of the updated positions, and updates the optimal position of each particle and the optimal position of the swarm according to the values of the fitness function of the updated positions, performs constant iteration and stops the constant iteration until reaching a stop condition. By adopting the method, the accuracy of sea wave information estimation can be improved, and the possibility of being trapped in a local optimum in sea wave inversion process can be greatly reduced.

Description

Sea surface flow information extraction method based on X-band radar image and particle swarm optimization

Technical Field

The invention belongs to the technical field of sea wave parameter inversion, and particularly relates to a sea surface flow information extraction method based on X-band radar images and particle swarm optimization.

Background

Sea surface currents, also known as ocean currents, are large-scale, relatively stable flows of seawater formed by internal changes and the influence of various external forces. With the development of scientific technology, the prior art can utilize ocean currents to select routes, generate electricity, catch fish and the like. Ocean currents have profound influences and restriction on formation and changes of various geological processes, physical processes, chemical processes, climate and weather over the ocean, so that understanding and mastering of ocean current information has important significance on fishery, marine transportation, pollution discharge, military affairs and the like.

In 1985, Yong et al proposed the use of Least Squares (LSM) iterationsThe method for estimating ocean current parameters by using a generation approximation method has the main idea that theoretical values are used

And from the energy spectrum

Calculate to get the truth value

The square root of the difference between reaches a minimum. The method makes it possible to obtain sea surface flow information through an X-band radar image. In 2001, Senetil considers the higher harmonic caused by nonlinearity of sea wave imaging and the frequency spectrum confusion phenomenon caused by low sampling, provides an iterative least square sea surface flow estimation algorithm, and provides a sea surface flow error estimation model, thereby greatly improving the estimation accuracy. In 2002, Gangeskar improves a calculation method of a sea surface flow field, spectral energy is used as weight to participate in iterative operation of flow velocity, a derivation method is also based on a least square method, an obtained equation can directly calculate x and y components of the flow velocity, and the size and direction of the sea surface flow can be obtained through direction conversion.

Although the sea surface flow estimation method based on the least square strategy has a certain accuracy, the optimal estimation strategy of the least square method is fixed (i.e. the square of the measured value and the theoretical value is the minimum), so the estimation accuracy is influenced to a certain extent.

Disclosure of Invention

The invention provides a sea surface flow information extraction method based on X-band radar images and particle swarm optimization, aiming at the problem that the estimation precision is influenced due to the fact that an estimation strategy is fixed by a sea surface flow estimation method based on a least square strategy in the prior art. According to the method, the particle swarm optimization method is used for replacing a least square strategy to estimate the sea surface flow, a flexible fitness function form is adopted, the problem that the estimation precision is influenced due to the fact that an estimation strategy is fixed in the traditional method is effectively solved, accuracy of sea flow information estimation is improved, and in addition, due to the global optimization capability of the particle swarm optimization, the possibility that the sea flow falls into local optimization in the sea flow inversion process is greatly reduced.

The invention provides a sea surface flow information extraction method based on X-band radar images and particle swarm optimization, which specifically comprises the following steps:

initializing parameters, wherein the parameters comprise a particle position, a particle speed, a maximum value and a minimum value of an inertia weight, a self-learning factor, a global learning factor, a maximum iteration number and an iteration stop threshold, and specifically comprise the following steps:

location of particles

Wherein,

the method comprises the steps that the size of the x-y direction component of the sea surface flow of the ith particle is represented by i, 1,2, a, n, wherein n represents the ith particle, n is more than or equal to 1 and is the population number of the particle, l represents the population algebra, and the initial l is 1; velocity of particles

r_x，r_yRandom numbers obeying a standard normal distribution; maximum value w of inertial weight_maxAnd a minimum value w_min：0＜w_max≤1，0≤w_min≤w_max(ii) a Self learning factor c₁Is greater than 0; global learning factor c₂Is greater than 0; maximum number of iterations l_max(ii) a The iteration stop threshold value xi is more than or equal to 2.

Step two, firstly, regarding each particle, the particle position of the particle when the l-th generation is 1 as the current optimal position of the particle

$P_i^b=(U_{x,i}^l,U_{y,i}^l)=U_i^l,$ i＝1，2，...，n；

Then, a fitness function of each particle of the current generation is determined

The value of (c):

$F\left(U_i^l\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-w^{\′}\left(\stackrel{\→}{k^{j^{\′},k^{\′}}}\right)|}^{\mathrm{\μ}}\·A)$

$=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B|}^{\mathrm{\μ}}\·A)$

wherein, i is 1,2_i′An actual value representing the frequency of the wave of the radar image after Fourier transformation,

expressed in ocean currents of

The theoretical value of the corresponding wave frequency,

is a wave number vector, w_i′、

Obtaining 32 or 64 continuous X-band radar images through three-dimensional Fourier transform; n is₃Number of X-band radar images, n₃32 or n₃64; at the same position on 32 or 64 consecutive radar images, a rectangular frame is used to select a part of the radar image, n₁，n₂Respectively equal to the number of horizontal and vertical pixel points of the radar image in the rectangular frame; d is the water depth; g is the acceleration of gravity; mu > 0 is the order of the error function; a > 0 is an energy weighting quantity, and an X-band radar image energy spectrum is set asWhen in use $A=\mathrm{fft}{(w_{i^{\′}},k_x^{j^{\′},k^{\′}},k_y^{j^{\′},k^{\′}})}^f$ When, represents the value of the difference

Carrying out f power energy weighting; p is the order of the dispersion relation, is a function of i ', j ', k ', determines the difference e (p, i ', j ', k ') under three values of p being 0,1 and 2, and selects the p value corresponding to the maximum difference, wherein the formula of the difference e (p, i ', j ', k ') is as follows:

$e(p,i^{\′},j^{\′},k^{\′})=w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B$

wherein B represents the influence of aliasing effect, and B is m.w_NParameter offloor (x) denotes taking x down integers, w_NRepresenting a cutoff frequency of the radar image;

finally, selecting the particle with the minimum fitness function value, and recording the position of the particle as the best position P found by the population when the generation l is 1^g。

Step three, updating the displacement of each particle in the (l + 1) th generation

$\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}}=\mathrm{\ω}\·\stackrel{\→}{{\mathrm{\υ}}_i^l}+c_1{r}_1^l\·\stackrel{\→}{U_i^l{P}_i^b}+c_2{r}_2^l\·\stackrel{\→}{U_i^l{P}^g},$ i＝1，2，...，n

Wherein the inertial weight

u is more than or equal to 0 and is a monotonous control quantity, and u is taken as 1;

is [0,1 ] of the first generation]The intervals follow evenly distributed random numbers.

Step four, updating the position of each particle in the (l + 1) th generation

$U_i^{l+1}=U_i^l+\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}},$ i＝1，2，...，n

Wherein the position of each particle of the (l + 1) th generation

Step five, firstly, determining the first generation and the second generation of each particleFitness function

i＝1，2，...，n：

$F\left(U_i^{l+1}\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^{l+1}-k_y^{j^{\′},k^{\′}}{U}_{y,i}^{l+1}+B|}^{\mathrm{\μ}}\·A)$

Then, the optimal position of each particle is updated, specifically: and for each particle, comparing the value of the fitness function of the current optimal position of the particle with the value of the fitness function of the particle in the (l + 1) th generation, and updating the position corresponding to the value of the smaller fitness function as the optimal position of the particle.

Step six, updating the optimal position P of the population^gSpecifically, the method comprises the following steps: is provided withFor the current optimum position of each particle

The position with the minimum fitness function value is compared with the fitness function value F (P) of the optimal position of the current population^g) Andfitness function value of

If it is

Then order

Otherwise, P is maintained^gThe value of (a) is not changed.

Step seven, updating l to l +1, if the iteration number l after updating is larger than the maximum iteration number l_maxOr the best position P of the population^gIf the xi iterations are not changed, the method is ended, otherwise, the step III is executed continuously.

The method has the advantages and positive effects that:

1. according to the sea surface flow information extraction method based on the X-band radar image and the particle swarm optimization, the global optimization capability of the particle swarm algorithm is utilized, the probability that two-dimensional variables of iterative solution fall into local optimal solution is effectively reduced, and the possibility that the two-dimensional variables fall into local optimal solution in the sea flow inversion process is greatly reduced;

2. according to the sea surface flow information extraction method based on the X-band radar image and the particle swarm optimization, due to the fact that the adaptability function form of the particle swarm algorithm is high in controllability, the problem of inaccurate solution caused by unnecessary difference weighting can be reduced as much as possible, and therefore the estimation accuracy of the sea surface flow information is improved;

3. according to the sea surface flow information extraction method based on the X-band radar image and particle swarm optimization, due to the fact that the adaptability function form of the particle swarm algorithm is high in flexibility, when the influence of aliasing effect and high-order dispersion relation is considered, the method is more convenient, the evaluation strategy of the aliasing effect and the high-order dispersion relation influence can be conveniently and visually introduced into the adaptability function, and the estimation precision of the sea surface flow information is improved;

4. compared with the traditional method, the sea surface flow information extraction method based on the X-band radar image and particle swarm optimization has stronger robustness, cannot be accurately solved due to the influence of low-frequency noise like the traditional method, and is more suitable for practical application in fishery, shipping, pollution discharge, military and the like.

Drawings

FIG. 1 is a flowchart illustrating the overall steps of the sea surface flow information extraction method of the present invention;

FIG. 2 is a diagram illustrating a method for extracting sea surface flow information according to the present invention

A schematic diagram of the time flow velocity inversion result;

FIG. 3 is a diagram illustrating a method for extracting sea surface flow information according to the present invention

A schematic diagram of the time flow direction inversion result;

fig. 4 is a schematic diagram of a flow velocity inversion result when μ ═ 1 in the sea surface flow information extraction method of the present invention;

fig. 5 is a schematic diagram of a flow direction inversion result when μ ═ 1 in the sea surface flow information extraction method of the present invention;

fig. 6 is a schematic diagram illustrating a flow velocity inversion result when μ ═ 2 in the sea surface flow information extraction method of the present invention;

fig. 7 is a schematic diagram of a flow direction inversion result when μ ═ 2 in the sea surface flow information extraction method of the present invention;

FIG. 8 is a schematic diagram of a flow velocity inversion result of a sea surface flow estimation method based on a least square strategy;

fig. 9 is a schematic diagram of a flow direction inversion result of a sea surface flow estimation method based on a least square strategy.

Detailed Description

The technical solution of the present invention will be described in detail with reference to the accompanying drawings and specific embodiments.

The initialization of sea surface flow is an important environment, the accuracy of the final iteration result is directly influenced by the quality of the initialization, and the method provides two initialization methods.

The first initialization method comprises the following steps: and a ship speed initialization method, wherein the ship speed is used as an initialization average value, and an empirical value is used as an initialization range to initialize the particles. Set the ship speed as

V′_x，V′_yRespectively representing the ship speed values in the x direction and the y direction, initializing the particle position

Wherein:

$\left\{\begin{array}{c}{U}_{x,i}^l=V_x^{\′}+(r^{\′}-0.5)\·b_1\\ U_{y,i}^l=V_y^{\′}+(r^{\′\′}-0.5)\·b_1\end{array}\right.---\left(1\right)$

wherein,

the component size of the ith particle sea surface flow in the x and y directions is represented by i, which is 1, 2. l represents a population algebra, and initial l is 1; r ', r' are [0,1]Random numbers with evenly distributed intervals; b₁The initial range of the ship speed is more than 0, and the empirical value is b₁3, this experienceThe value is obtained according to the variation range of the flow velocity true value when the ship speed is zero.

The second initialization method comprises the following steps: and (3) a flow estimation initialization method, which adopts least square to calculate the flow value as an initialization average value under the condition of not considering the high-order dispersion relation, and initializes the particles by taking the empirical value as an initialization range. Let the least square estimate be

V″_x，V″_yRespectively representing the component magnitudes in the x-direction and the y-direction, the particle position is initialized $U_i^l=(U_{x,i}^l,U_{y,i}^l),$ Wherein

$\left\{\begin{array}{c}{U}_{x,i}^l=V_x^{\text{\′\′}}+(r^{\′}-0.5)\·b_2\\ \stackrel{}{U_{y,i}^l}=V_y^{\′\′}+(r^{\′\′}-0.5)\·b_2\end{array}\right.---\left(2\right)$

Wherein, i is 1,2, n, which represents the ith particle, and n is more than or equal to 1 and is the population number of the particles; l represents a population algebra, and initial l is 1; r ', r' are [0,1]Random numbers with evenly distributed intervals; b₂> 0, initializing a range for estimating the flow, with an empirical value of b₂The empirical value is obtained from the difference between the estimated flow value and the true flow rate value, 1.

The fitness function of the particle swarm algorithm is used in the method, and the fitness function in the method is determined as follows: inversion of sea by using X-band radar imageThe surface flow requires three-dimensional Fourier transformation of the radar image, and the wavenumber energy spectrum of each frequency obtained from the transformation result is fft (w, k)_x，k_y) Wherein w represents frequency, k_xDenotes the wave number, k, in the x direction_yRepresenting the wavenumber in the y-direction, the dispersion relation function in the absence of noise is:

$w^{\′}\left(\stackrel{\→}{k}\right)=(P+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k}\left|\right|d}{p+1}\right)}+\stackrel{\→}{k}\stackrel{\→}{U}---\left(3\right)$

wherein

A theoretical value representing the corresponding wave frequency when the ocean current is U,

is a vector of wave numbers and is,

for sea surface flow velocity, d is water depth, g is gravitational acceleration, and p is 0,1, 2.

And fitting the energy spectrum by using the dispersion relation function to obtain the size of the sea surface flow. Since the corresponding spectral energy is already very small when the order p ≧ 3 of the spectrum, the influence of higher-order waves of order p ≧ 3 can be ignored, and only the dispersion relation when order p is 0,1,2 is used for fitting. Fitness function of particle swarm optimizationCan be set as follows:

$F\left(U_i^l\right)=\underset{i\′=1}{\overset{n_3/2+1}{\∑}}\underset{j\′=1}{\overset{n_2}{\∑}}\underset{k\′=1}{\overset{n_1}{\∑}}({\text{|}{w}_{i^{\′}}-w^{\′}\left(\stackrel{\→}{k^{j^{\′},k^{\′}}}\right)|}^{\mathrm{\μ}}\·A)---\left(4\right)$

wherein w_i′The actual value representing the frequency of the wave after the fourier transform,

is a wave number vector, w_i′，

Obtaining 32 (or 64) continuous X-band radar images through three-dimensional Fourier transform; n is₃Is the number of X-band radar images, n₃32 or n₃As can be seen from the nature of the fourier transform, 64,the spectrum obtained after fourier transformation is symmetric about the origin, so only half of the transformed results need to be fitted, i.e. i' is 1,2₃/2+1. Because the data volume of the whole radar image is very large, the calculation amount is overlarge, and because all data are not necessary to be calculated, the frame selection processing is carried out at the same position on 32 or 64 continuous radar images, a rectangular frame is adopted to select a part of the whole radar image so as to reduce the calculation amount and improve the calculation speed, and n is n₁，n₂The number of the horizontal and vertical pixel points of the selected rectangular radar image is equal to the number of the horizontal and vertical pixel points of the selected rectangular radar image; mu > 0 is the order of the error function, when mu is 2, the error function is the target function of the least square algorithm, and the difference value is obtainedIs squared to the difference

The method increases the influence of the energy points which are not on the dispersion relation curved surface on the fitness function, and reduces the influence of the energy points which are on the dispersion relation curved surface on the fitness function. When the energy of the noise point is large, the sea wave energy point is easily submerged in the noise and the sea current information cannot be extracted, so that the sea surface current inversion is inevitably inaccurate, but the problem is inevitable in the least square algorithm, the particle swarm algorithm can be adopted to randomly adjust the value of the mu value, and the unnecessary difference value is effectively prevented

Inaccurate calculation caused by weighting; a > 0 is an energy weighting amount, which indicates that the difference is energy weighted when a is 1, and when a is fft (w, k)_x，k_y)^fTime means that the difference is weighted with energy to the power of f. In addition, to take into account the influence of energy points on the higher order dispersion relation, it is necessary to operate on the current operation point (k)_x，k_y) The order of the p value is determined, and then the sum operation described in the formula (4) is participated in.

Due to the limitation of sampling frequency, the bandwidth of the radar image is limited, so that aliasing effect can occur when the wave frequency exceeds the bandwidth of the radar image, and the cutoff frequency of the radar image is set as w_NThen the frequency of the radar image is distributed in [0, w ]_N]Interval, when the frequency out of range [0, w ] is calculated according to the dispersion relation function_N]When it is needed to

And correcting by the following formula:

$w^{\′}\left(\stackrel{\→}{k}\right)=(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k}\left|\right|}{\text{p+1}}\tanh \left(\frac{\left|\right|\stackrel{\→}{k}\left|\right|d}{p+1}\right)}+\stackrel{\→}{k}\stackrel{\→}{U}-B---\left(5\right)$

wherein B is m.w_N，

The function floor (x) represents the fitness function obtained by taking the whole of x and substituting the formula (5) into the formula (4)Comprises the following steps:

$F\left(U_i^l\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B|}^{\mathrm{\μ}}\·A)---\left(6\right)$

where p is the order of the dispersion relation, a function of i ', j', k ', making the difference e (p, i', j ', k') equal

$e(p,i^{\′},j^{\′},k^{\′})=w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′}{,k}^{\′}}{U}_{y,i}^l+B---\left(7\right)$

Let e (p ', i', j ', k') -max { e (0, i ', j', k '), e (1, i', j ', k'), e (2, i ', j', k ') } then p ═ p'.

The invention provides a particle swarm optimization-based sea surface flow information extraction method, which specifically comprises the following steps as shown in figure 1:

the method comprises the following steps: and initializing parameters.

Initializing the particle position based on the formula (1) of the ship speed initialization method or the formula (2) of the estimated flow initialization method

Wherein

Is the magnitude of the x, y direction component of the sea surface flow, i ═1,2, a.n, which represents the ith particle, wherein n is more than or equal to 1 and is the population quantity; l ═ 1 denotes the population algebra; initialization shift

r_x、r_yRandom numbers obeying a standard normal distribution; maximum value of inertial weight: omega is more than 0_max1, minimum value of inertial weight: omega is more than or equal to 0_min≤ω_max(ii) a The self-learning factor and the global learning factor are respectively c₁＞0，c₂Is greater than 0; maximum number of iterations is l_max(ii) a The iteration stop threshold value xi is more than or equal to 2.

Step two: and determining the particles with the minimum fitness value under the current algebra l by utilizing a fitness function, recording the particles as the best positions found by the population, and recording the best positions of the particles.

Order toi＝1，2，...，n，

The optimum position passed by the particle in the 1 st generation is shown. Application of equation (6) to determine fitness function

The values are:

$F\left(U_i^l\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B|}^{\mathrm{\μ}}\·A)---\left(8\right)$

wherein, w_i′Is the frequency of the waves and is,

is a wave number vector, w_i′，

Is obtained by three-dimensional Fourier transform of 32 (or 64) continuous X-band radar images, n₃32 (or n)₃64) the spectrum obtained after the fourier transform is symmetric about the origin, and therefore only half of the transformed result needs to be calculated (i.e. i' 1, 2.. n)₃/2+1)，n₁，n₂The number of the horizontal and vertical pixel points of the selected rectangular radar image is equal to the number of the horizontal and vertical pixel points of the selected rectangular radar image; d is the water depth, g is the acceleration of gravity, p is the order of the dispersion relation, a function on i ', j', k ', the difference e (p, i, j, k') is calculated according to equation (7) as

$e(p,i^{\′},j^{\′},k^{\′})=w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′}{,k}^{\′}}{U}_{y,i}^l+B---\left(9\right)$

Where p is 0,1,2, and then point-aligning

And judging the order of the corresponding dispersion relation: let e (p ', i', j ', k') -max { e (0, i ', j', k '), e (1, i', j ', k'), e (2, i ', j', k ') }, then p ═ p'; mu > 0 is the order of the error function; a > 0 is an energy weighting quantity, and an X-band radar image energy spectrum is set as

Then when a is 1, it means that the difference is energy weighted, when

Time means that the energy weighting of the difference is performed to the power of f; b is the term of aliasing effect B ═ m.w_N，

floor (x) indicates that the integer is removed for x.

Taking the particle with the minimum fitness value, and setting the position of the particle as U_i′Let us order

Record the best position found by the population when the l-1 generation.

Step three: particle renewal Using equation (10)

Shift in the l +1 th generation

i＝1，2，...，n。

$\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}}=\mathrm{\ω}\·\stackrel{\→}{{\mathrm{\υ}}_i^l}+c_1{r}_1^l\·\stackrel{\→}{U_i^l{P}_i^b}+c_2{r}_2^l\·\stackrel{\→}{U_i^l{P}^g}---\left(10\right)$

Wherein the inertial weight

u is more than or equal to 0 and is a monotonous control quantity, and u is generally 1;

is [0,1 ] in the first generation]The intervals follow evenly distributed random numbers.

Step four: particle renewal Using equation (11)

In the position of the next generation

i＝1，2，...，n：

$U_i^{l+1}=U_i^l+\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}}---\left(11\right)$

Wherein the position of each particle of the (l + 1) th generation

Step five: and determining the fitness function value of each particle in the next generation, and updating the optimal position of each particle.

Calculation of each particle by applying equation (6)

Fitness function in the next generation

The (i ═ 1, 2.., n) values are:

$F\left(U_i^l\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B|}^{\mathrm{\μ}}\·A)---\left(12\right)$

for each particle, comparing the value of the fitness function obtained in the (l + 1) th generation with the value of the fitness function of the stored optimal position of the particle, and setting the position with smaller fitness function value as the optimal position of the particle, namely if

Then orderOtherwise, the current optimal position of the particle is kept unchanged.

Step six: and updating the optimal position of the whole population. Is provided withFor the current optimum position of each particle

The position with the minimum fitness value is compared with the fitness function value of the current best position of the population

The fitness function value of if

Then order

If it is

The current optimal position of the population is kept unchanged.

Step seven: if the iteration number l is larger than the maximum iteration number l, taking l as l +1_maxOr P^gAnd if no change occurs in the continuous xi iterations, stopping the operation, otherwise, turning to the step three.

Example (b):

the radar image data used in the experiment is a sea surface echo radar image acquired on a sea island in Tan county, Tan, Fujian province in 2010 at 10 months and 23 days, the average water depth of a radar imaging area is about 28 meters, the maximum flow velocity value of sea waves is about 1.5m/s, and the wind direction of the imaging area is mainly southwest wind and northeast wind. 32 radar images obtained in succession are selected as a time series of radar images required for data processing. The height of the radar is 40m, and the image covers halfThe diameter is about 2km, the average rotation speed of the antenna is 27 revolutions per minute, namely the rotation period is about 2.25s, and the total time length is 72 seconds. For each image in the radar image sequence, because the wave distribution is affected by the terrain and is not uniform when approaching the shore, through the effective echo area analysis of the radar image, a more ideal sea area selection rectangular frame is selected in the 75-degree direction. The minimum distance between the rectangular area and the radar is 600 meters, 128 multiplied by 128 pixel points are taken, and if the distance resolution of the radar is 7.5 meters, namely the actual sea area corresponding to each pixel point is 7.5 meters, the size of the selected rectangular area is 960 multiplied by 960 m. It is emphasized that the selection of the rectangular frame in each image of the time series of images of the radar echo signal must be at the same position to ensure the consistency of the data. From the Fourier transform theory, the spectral resolution takes values as follows: Δ k_x＝Δk_y＝0.0065m^-1，Δw＝0.0872s^-1，Δk_x、Δk_yFor the wave number after fourier transform, Δ w is the frequency interval after fourier transform, and the Nyquist (Nyquist) cutoff wave number and cutoff frequency in the x, y directions are: k is a radical of_Nx＝k_Ny＝0.4187m^-1，w_N＝1.3956s^-1。

When a real radar image sequence is inverted, due to the fact that a large-area real-time sea surface flow truth value is difficult to obtain, an inversion result cannot be compared with the truth value, and therefore evaluation indexes need to be defined by users according to the properties of the sea surface flow to verify the quality of an experimental result. Physical analysis shows that the change of the sea surface flow field in a certain area in time is relatively small, and the uniformity in space is high. Therefore, the evaluation index of the problem can be defined according to the characteristic that the sea surface flow velocity and the flow direction change rate in the adjacent time are low.

Definition of evaluation index: the average value of sea surface flow difference values obtained by inverting two groups of radar images adjacent in time (each group of radar images consists of 32 continuous radar images) is expressed by the formula:

$\frac{1}{N-1}\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}|u_i-u_{i+1}|---\left(13\right)$

where N is the number of groups of radar images, u_iAnd (4) obtaining the flow velocity or flow direction value for the inversion of the ith group of radar images. The smaller the evaluation index is, the better the calculation result is.

Order learning factor c₁＝1.4962，c₂1.4962; upper and lower limits omega of inertial weight_max＝0.9，ω_min0.4; maximum number of iterations l _max100; the population number n is 15; the iteration stop threshold value xi is 10; initializing flow speed by adopting flow estimation initialization method, and estimating flow initialization range b₂＝1。

In order to verify the flexibility of the fitness function of the particle swarm algorithm and the superiority of the sea surface flow estimation method based on the least square strategy, different error weighted values are adopted, namely the error weighted value mu in the formula (6) is adopted, different values of mu are compared, the influence on the flow velocity inversion result is exerted, and meanwhile, the result is compared with the sea surface flow estimation method RLSM result based on the least square strategy.

Respectively take

Mu-1, mu-2, the results obtained with the method of the invention are shown in fig. 2, 3, 4, 5, 6, 7, sea surface using least squares based strategyThe flow estimation method results are shown in fig. 8 and 9. The abscissa in fig. 2 to 9 represents the number N of groups of radar images participating in the simulation. The results obtained for evaluation index e are shown in Table 1, and it can be seen that when

The inversion result of the method is ideal when mu is 1, but the fluctuation of the inversion result is larger when mu is 2, which directly leads to the reduction of the precision, the fitness function adopted by the sea surface flow estimation method based on the least square strategy is the case when mu is 2, comparing the corresponding evaluation index e when mu is 2 in fig. 6 and 7 with those in fig. 8 and 9 and table 1, it can be seen that the dispersion of the method of the present invention with respect to the result obtained by the sea surface current estimation method based on the least squares strategy is substantially the same when μ is 2, the method not only can reproduce the operation mode of the sea surface flow estimation method of the traditional least square strategy, but also can flexibly improve the fitness function on the basis, therefore, the calculation accuracy is improved, and the experiment shows that the method has obvious superiority compared with the sea surface flow estimation method based on the least square strategy.

TABLE 1 comparison of the Performance of the method of the present invention with RLSM algorithm under different weights μ

Compared with the traditional sea surface flow estimation method based on the least square strategy, the sea surface flow information extraction method provided by the invention has stronger flexibility, and due to the flexibility of the selection of the fitness function, the error weight mu can be properly selected in the sea surface flow inversion process, so that the problem of inaccurate calculation caused by the meaningless error weighting of the traditional least square algorithm is avoided. The accurate fitness function enables the sea surface flow to be calculated more accurately, the noise interference resistance is higher, and the sea surface flow is more visual and reliable when the influence of the dispersion relation and the aliasing effect is considered.

Claims (4)

1. A sea surface flow information extraction method based on X-band radar images and particle swarm optimization is characterized by comprising the following steps:

initializing parameters, wherein the parameters comprise a particle position, a particle speed, a maximum value and a minimum value of an inertia weight, a self-learning factor, a global learning factor, a maximum iteration number and an iteration stop threshold, and specifically comprise the following steps:

location of particles

Wherein,

the method comprises the steps that the size of the x-y direction component of the sea surface flow of the ith particle is represented by i, 1,2, a, n, wherein n represents the ith particle, n is more than or equal to 1 and is the population number of the particle, l represents the population algebra, and the initial l is 1; velocity of particles

r_x，r_yRandom numbers obeying a standard normal distribution; maximum value w of inertial weight_maxAnd a minimum value w_min：0＜w_max≤1，0≤w_min≤w_max(ii) a Self learning factor c₁Is greater than 0; global learning factor c₂Is greater than 0; maximum number of iterations l_max(ii) a The iteration stop threshold value xi is more than or equal to 2;

step two, firstly, regarding each particle, the particle position of the particle when the l-th generation is 1 as the current optimal position of the particle $P_i^b:P_i^b=(U_{x,i}^l,U_{y,i}^l)=U_i^l,i=\mathrm{1,2},...,n;$

Then, a fitness function of each particle of the current generation is determined

The value of (c):

$F\left(U_i^l\right)=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-w^{\′}\left(\stackrel{\→}{k^{j^{\′},k^{\′}}}\right)|}^{\mathrm{\μ}}\·A)$

$=\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B|}^{\mathrm{\μ}}\·A)$

wherein, i is 1,2_i′An actual value representing the frequency of the wave of the radar image after Fourier transformation,expressed in ocean currents ofThe theoretical value of the corresponding wave frequency,

is a wave number vector, w_i′、

Obtaining 32 or 64 continuous X-band radar images through three-dimensional Fourier transform; n is₃Number of X-band radar images, n₃32 or n₃64; at the same position on 32 or 64 consecutive radar images, a rectangular frame is used to select a part of the radar image, n₁，n₂Respectively equal to the number of horizontal and vertical pixel points of the radar image in the rectangular frame; d is the water depth; g is the acceleration of gravity; mu > 0 is the order of the error function; a > 0 is an energy weighting quantity, and an X-band radar image energy spectrum is set as

When in use $A=\mathrm{fft}{(w_{i^{\′}},k_x^{j^{\′},k^{\′}},k_y^{j^{\′},k^{\′}})}^f$ When, represents the value of the difference

Carrying out f power energy weighting; p is the order of the dispersion relation, is a function of i ', j ', k ', and determines the difference e (p, i ', j ', k ') under the three values of p =0,1 and 2, and selects the p value corresponding to the maximum difference and the difference e (p, i ', j ', k ')The formula is as follows:

$e(p,i^{\′},j^{\′},k^{\′})=w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^l-k_y^{j^{\′},k^{\′}}{U}_{y,i}^l+B$

wherein B represents the influence of aliasing effect, and B is m.w_NParameter of

floor (x) denotes taking x down integers, w_NRepresenting a cutoff frequency of the radar image;

finally, selecting the particle with the minimum fitness function value, and recording the position of the particle as the best position P found by the population when the generation l is 1^g；

Step three, updating the displacement of each particle in the (l + 1) th generation

$\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}}=\mathrm{\ω}\·\stackrel{\→}{{\mathrm{\υ}}_i^l}+c_1{r}_1^l\·\stackrel{\→}{U_i^l{P}_i^b}+c_2{r}_2^l\·\stackrel{\→}{U_i^l{P}^g},i=\mathrm{1,2},...,n$

Wherein the inertial weight

u is more than or equal to 0 and is a monotonous control quantity, and u is taken as 1;is [0,1 ] of the first generation]Random numbers with evenly distributed intervals;

step four, updating each particle in the first stepPosition of +1 generation

$U_i^{l+1}=U_i^l+\stackrel{\→}{{\mathrm{\υ}}_i^{l+1}},i=\mathrm{1,2},...,n$

Wherein the position of each particle of the (l + 1) th generation

Step five, firstly, determining the fitness function of each particle in the (l + 1) th generation

$F\left(U_i^{l+1}\right)\text{=}\underset{i^{\′}=1}{\overset{n_3/2+1}{\mathrm{\Σ}}}\underset{j^{\′}=1}{\overset{n_2}{\mathrm{\Σ}}}\underset{k^{\′}=1}{\overset{n_1}{\mathrm{\Σ}}}({|w_{i^{\′}}-(p+1)\sqrt{\frac{g\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|}{p+1}\tanh \left(\frac{\left|\right|\stackrel{\→}{k^{j^{\′},k^{\′}}}\left|\right|d}{p+1}\right)}-k_x^{j^{\′},k^{\′}}{U}_{x,i}^{l+1}-k_y^{j^{\′},k^{\′}}{U}_{y,i}^{l+1}+B|}^{\mathrm{\μ}}\·A)$

Then, the optimal position of each particle is updated, specifically: for each particle, comparing the value of the fitness function of the current optimal position of the particle with the value of the fitness function of the particle in the (l + 1) th generation, and updating the position corresponding to the value of the smaller fitness function as the optimal position of the particle;

step six, updating the optimal position P of the population^gSpecifically, the method comprises the following steps: is provided with

For the current optimum position of each particleThe position with the minimum fitness function value is compared with the fitness function value F (P) of the optimal position of the current population^g) Andfitness function value of

If it is

Then orderOtherwise, P is maintained^gThe value of (d) is unchanged;

step seven, updating l to l +1, if the iteration number l after updating is larger than the maximum iteration number l_maxOr the best position P of the population^gIf the xi iterations are not changed, the method is ended, otherwise, the step III is executed continuously.

2. The method for extracting sea surface flow information based on X-band radar image and particle swarm optimization according to claim 1, wherein the particle position in the first step has two initialization methods, which are respectively:

the first method comprises the following steps: the ship speed initialization method is to initialize the particles by taking the ship speed as an initialization average value, and the formula (1) is as follows:

$\left\{\begin{array}{c}{U}_{x,i}^l=V_x^{\′}+(r^{\′}-0.5)\·b_1\\ U_{y,i}^l=V_y^{\′}+(r^{\′\′}-0.5)\·b_1\end{array}\right.---\left(1\right)$

wherein the particle positions are initialized

Initial 1, boat speed

V_x′,V_y'represents the ship speed values in the x direction and the y direction respectively, and r', r '' are [0,1]Random numbers whose intervals follow an even distribution, b₁For initializing the range of the ship's speed, b₁＞0；

And the second method comprises the following steps: the stream estimation initialization method is characterized in that a least square stream estimation value obtained by applying a least square under the condition of not considering a high-order dispersion relation is used as an initialization average value to initialize particles, and the formula (2) is as follows:

$\left\{\begin{array}{c}{U}_{x,i}^l=V_x^{\′\′}+(r^{\′}-0.5)\·b_2\\ U_{y,i}^l=V_y^{\′\′}+(r^{\′\′}-0.5)\·b_2\end{array}\right.---\left(2\right)$

wherein the least squares estimate

V_x′′,V_y'' indicates the x-direction and the y-direction respectivelyThe magnitude of the component in the direction, r ', r ' ' is [0,1 ]]Random numbers whose intervals follow an even distribution, b₂Initializing a range for estimating a flow, b₂＞0。

3. The method for extracting sea surface flow information based on X-band radar image and particle swarm optimization according to claim 2, wherein the ship speed initialization range b is₁Has an empirical value of 3.

4. The method for extracting sea surface flow information based on X-band radar image and particle swarm optimization according to claim 2, wherein the flow estimation initialization range b is₂Has an empirical value of 1.

Images