CN103969643A - Method for conducting X wave band navigation radar wave parameter inversion through band-pass filter based on novel wave dispersion relation - Google Patents

Abstract

The invention belongs to the field of wave remote sensing technology, and particularly relates to a method for conducting X wave band navigation radar wave parameter inversion through a band-pass filter based on a novel wave dispersion relation, wherein sea clutter images are obtained through navigation radar so that wave parameter inversion can be conducted. The method comprises the steps that radar image data are collected; radar images are preprocessed; Fourier transform is conducted on an image sequence under a Cartesian coordinate system, so that a three-dimensional wave number frequency image spectrum of the radar images is obtained; wave spectrum information is extracted; wave information inversion is achieved. Influence on the dispersion relation by the ship speed is kept by the filter, the problem that the band pass of a traditional broadband-pass filter is increased along with increase of the movement speed is effectively solved, and therefore filtering can be conducted under the condition that a radar platform is moved along with movement of a ship. In band pass boundary derivation of the novel filter, an approximate value of any quantity is avoided, calculation errors are reduced, influence on the band pass boundary is avoided, band width calculation is more accurate, and wave inversion accuracy is improved.

Description

One is carried out X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter

Technical field

The invention belongs to Remote sensing of ocean wave technical field, be specifically related to a kind of utilize that sea clutter image that navar obtains carries out ocean wave parameter inverting carry out X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter.

Background technology

Echo to X-band radar is analyzed, and can obtain the parameter of wave.But in this echo, except containing Ocean-wave Signal, also contain a large amount of noises, as co-channel interference, ground unrest, sleet interference and object echo etc.These noises will have a strong impact on the parametric inversion of wave, make result distortion.So need to remove these noise signals from the echo of radar and extract useful Ocean-wave Signal.Because this signal meets the dispersion relation of gravity wave, thus can utilize dispersion relation to construct wave filter, by signal extraction out.Concrete filtering flow process is as follows:

Sea clutter image spectrum I (k _x, k _y, ω) and → dispersion relation filtering → signal pattern spectrum F (k _x, k _y, ω)

This theory of dispersion relation that wave meets gravity wave is proposed in 1978 by LeBlond and Mysak.Document [1] (Paul H.LeBlond, Lawrence A.Mysak.Waves in the ocean[M] .1978, Elsevier.) sees reference

Suppose that in selected analyzed area sea wave and flow field are that space uniform and time are stable, the dispersion relation of single order gravity wave is described below

${\mathrm{\ω}}_0^2=g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)---\left(1\right)$

Wherein, d is the depth of water, and g is local gravitational acceleration, for wave number, λ is wave wavelength.If there is the stream of the surface with respect to antenna platform of radar speed wherein the dispersion relation of introducing so a Doppler shift becomes:

$\mathrm{\ω}={\mathrm{\ω}}_0+\stackrel{\→}{k}\·\stackrel{\→}{u}---\left(2\right)$

?

$\mathrm{\ω}=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+\stackrel{\→}{k}\·\stackrel{\→}{u}---\left(3\right)$

Above-mentioned is that radar does not exist relative velocity and has the dispersion relation curve under relative velocity with surface, sea.Curve a is that flow velocity is the dispersion relation curve of 0 o'clock, and it is about longitudinal axis symmetry.Curve b is the dispersion relation curve that has Doppler shift.Level line represents the distribution of energy.

Show according to document, the existing bandpass filter based on dispersion relation of extracting for sea clutter image Ocean-wave Signal is mainly divided into two large classes, i.e. arrowband and broadband at present.In the development history of bandpass filter, what occur the earliest is the wave filter of arrowband form, has occurred broadband filter subsequently in order to solve the deficiency of narrow band filter.

The principle of narrow band filter is to find suitable method to calculate the velocity vector in dispersion relation formula (3) dispersion relation is known like this.Because the Wave energy overwhelming majority is distributed near dispersion relation curve, therefore by narrow bandpass being set near dispersion relation curve, just can obtain real Wave Information, filtering noise disturbs.

Narrow bandpass dispersion relation wave filter is the earliest that the people such as Reichert, Nieto in 1999 propose.The people such as Reichert have provided respectively two kinds of bandpass filter definition that thought is similar in document [2-4], as shown in formula (4) (5).

Because sea clutter image spectrum is obtained by 3 d-dem Fourier transform, therefore wave number and frequencies omega is all the points after discretize.In formula (4), I (k _x, k _y, ω) and be sea clutter image spectrum; F (k _x, k _y, ω) and be the signal pattern spectrum after bandpass filtering, this wave filter is to all wave numbers travel through and obtain corresponding frequencies omega, and the energy that only retains this place is as Wave energy, other parts are as noise filtering [2] [3].In formula (5), d is the depth of water, obtain ω by estimating stream _prepresent the frequency of a certain frequency chip.This wave filter is equally to all wave numbers travel through and obtain corresponding frequencies omega, if ω ∈ is [ω _p-1, ω _p+1], the energy [4] under this wave number on p, p-1 or p, p+1 sheet in reserve frequency sheet.Therefore, this wave filter band is logical logical slightly wider than band described in (4) formula.Document [2-4] (Jose Carols Nieto Borge sees reference, Konstanze Reichert, Jurgen Dittmer.Use of nauticalradar as a wave monitoring instrument[J] .Coastal Engineering, 1999, 37:331-342P.KonstanzeReichert, Jose Carols Nieto Borge, Jurgen Dittmer.WaMoS2:An operational wave monitoringsystem.Proceedings-Oceanology International ' 98, Brighton, 1998.K.Reichert, K.Hessner, J.C.Nieto Borge, ea tl.WaMoS2:A radar based wave and current monitoring system[J] .ISOPE ' 99, Brest, 1999, 3:1-5P.)

Within 2012, the red grade of Shen Ji is improved above-mentioned wave filter, has considered to think that Wave energy is not only distributed on 1st order chromatic dispersion relation curve after high-order dispersion relation.High-order dispersion relation formula is

${\mathrm{\ω}}_m=m\sqrt{\frac{g\left|\stackrel{\→}{k}\right|}{m}\·\tanh \left(\frac{d\left|\stackrel{\→}{k}\right|}{m}\right)}+\stackrel{\→}{k}\·\stackrel{\→}{u}m=\mathrm{1,2},...\left(6\right)$

Here the order that m is dispersion relation, ω _mfor m order wave frequency.Visible in the time of m=1 formula (6) be 1st order chromatic dispersion relation equation (3), and the ω that asks for of high-order dispersion relation _mcan increase.Utilize high-order dispersion to be related to that border generates new wave filter, this filter definition is as follows

Wherein, for the wave filter upper bound, ask the limit to obtain by high-order dispersion relation curve, increase than the result in formula (4) (5); for wave filter lower bound, with coming to the same thing in formula (4) (5).This wave filter still belongs to narrow band filter, need to before bandpass filtering, estimate flow rate and direction, and is the wave number to each discretize equally travel through.See reference document [5] (Shen Jihong, Li Ying, Shujuan WANG, Dai Yuntao, the design of filter [J] based on high-order dispersion relation in Wang Xiao warbler .X band radar. Chinese navigation, 2012,4 (35): 4-7.)

Be more than the current situation of narrow band filter at present existing document, its filter effect places one's entire reliance upon and estimates stream precision.Due to the factors such as image resolution ratio is limited, spectral power distribution is inhomogeneous, noisiness is complicated impact, still lack at present the technological means of accurately estimating stream both at home and abroad, estimate stream effect all not ideal enough, flow rate and direction evaluated error has had a strong impact on the performance of narrow band filter, and the real not dispersion relation curve based on narrow bandpass can cause wave filter to retain not belonging to the noise of wave in a large number and filter true Wave Information.

For this problem, first Rune designed based on Peak Flow Rate U in 2002 _maxbroadband-pass filter [6].The method does not need accurately to estimate stream, has solved the deficiency of narrow band filter maximum.The method has become the main flow wave filter being widely used at present.Consider that this filter method is defined as foloows the in the situation that of deep water:

Wherein

$\left|B_n\right|\frac{{(\mathrm{\ω}-\frac{\mathrm{\Δ\ω}}{2}-U_{\max }\frac{{\mathrm{\ω}}^2}{g})}^2}{g}-\frac{\mathrm{\Δk}}{\sqrt{2}}$ $\left|B_P\right|\frac{{(\mathrm{\ω}-\frac{\mathrm{\Δ\ω}}{2}-U_{\max }\frac{{\mathrm{\ω}}^2}{g})}^2}{g}+\frac{\mathrm{\Δk}}{\sqrt{2}}$

Δ ω is frequency resolution, and Δ k is wavenumber resolution, U _maxfor the normal condition Peak Flow Rate that surface may reach of plunging into the commercial sea.This wave filter border formula B _n, B _pderive and obtain by the basic dispersion relation formula (3) with flow velocity.In derivation, subitem is got to approximate value, first in the time that the depth of water is more than or equal to 1/2nd wavelength, Ke Yishe secondly due to value be less than 1, value be greater than value, and wave wavelength is larger, and the two gap is larger.Therefore reduce as much as possible error for simplifying when calculating, can be by bring formula (3) into, obtain so

$\mathrm{\ω}=\sqrt{g\left|\stackrel{\→}{k}\right|}+\frac{{\mathrm{\ω}}^2}{g}\·\left|\stackrel{\→}{u}\right|\·\cos \mathrm{\θ}---\left(9\right)$

Therefore in the time of cos θ=1,

${\left|\stackrel{\→}{k}\right|}_{\min }=\frac{{(\mathrm{\ω}-\frac{{\mathrm{\ω}}^2}{g}\·|\stackrel{\→}{u}\left|\right)}^2}{g}---\left(10\right)$

In the time of cos θ=-1,

${\left|\stackrel{\→}{k}\right|}_{\max }=\frac{{(\mathrm{\ω}+\frac{{\mathrm{\ω}}^2}{g}\·|\stackrel{\→}{u}\left|\right)}^2}{g}---\left(11\right)$

Their are respectively corresponding lower boundary and the coboundary of wave filter.The impact of simultaneously considering frequency resolution and wavenumber resolution just can obtain B _n, B _pthe border of wave filter.See reference document [6] (Rune Gangeskar.Ocean current estimated fromX-band radar sea surface images[J] .IEEE transactions on geosciences and remote sensing, 2002,40:783-791.)

2012, the people such as Jia Ruicai carried out improving [7] to the above-mentioned wide band-pass filtering method based on dispersion relation again.The method is with the maximum wave number k of hypothesis _mbe different from above-mentioned filtering mode.This filter definition is

Wherein,

$\left|k_P\right|\frac{{(\mathrm{\ω}-\frac{\mathrm{\Δ\ω}}{2}-U_{\max }|k_m\left|\right)}^2}{g}+\frac{\mathrm{\Δk}}{\sqrt{2}}$ $\left|k_n\right|\frac{{(\mathrm{\ω}-\frac{\mathrm{\Δ\ω}}{2}-U_{\max }|k_m\left|\right)}^2}{g}-\frac{\mathrm{\Δk}}{\sqrt{2}}$

K _mcomputing method be first to calculate spectrum I (k _x, k _y, ω) in maximal value C, then obtain spectrum value be greater than 0.9C wave number mould value a little, then by wave number mould value arithmetic mean, obtain | km|.The same B of this wave filter _n, B _pwave filter is similar, with leading wave number k _mreplace ω ²/ g, amplifies filter bandwidht, but and not as ω ²/ g is accurate.See reference document [7] (Yuan Gannan, Jia Ruicai, Shujuan WANG, Dai Yuntao. dispersion relation Design of Bandpass [J]. Central China University of Science and Technology's journal (natural science edition), 2012,3 (40): 36-30.)

Existing these two kinds of broadband filters, all by frequencies omega is traveled through, are asked for the upper bound of the wave number mould value that this ω is corresponding and after the next time, the wave number under each ω are screened.This wave filter that ω is traveled through has identical time complexity with the wave filter that wave number k is traveled through, and is 0 (n ²).But from analysis above, the approximate value of bringing in these two kinds of wave filters will cause this wave filter to have two defects: the one, be only applicable to the depth of water and be more than or equal to 1/2nd wavelength situations; The 2nd, only under large wavelength, this filter filtering effect is comparatively accurate.

In addition, two kinds of wideband filtered methods all cannot directly be used under platform motion state, visible in the time that Texas tower is kept in motion from algorithm definition, the bandwidth of wave filter will significantly be amplified, do not have filter effect, make to have comprised a large amount of noise in inversion result, affect precision.With B _n, B _pwave filter is example, in the time that Texas tower is static, and U _max=3m/s, along with the increase of ship's speed, U _maxalso should be larger, band is logical to broaden thereupon, and the pace of change of coboundary is greater than lower boundary.In addition, the U in this wave filter _maxalso will cause bandwidth under motion state to increase, can not play filter effect.

As can be seen here, although broadband filter does not need to estimate stream, use greatest hope flow velocity U _maxsolve and under narrow bandpass, estimated the inaccurate problem of bringing of stream, but led to and can infinitely be amplified at the band of dynamically lower wave filter, cannot play filter effect.

For this problem, the present invention has designed a kind of novel dispersion relation bandpass filter.This wave filter travels through wave number, and belongs to the category in broadband, and the bandwidth that has still effectively solved again traditional broadband-pass filter can increase this problem with the increase of movement velocity.The movement velocity of Texas tower is considered with in logical accounting equation, made the bandwidth under motion state identical with static state.Under motion state, the increase of ship's speed can cause Wave energy to follow dispersion relation generation skew, and new band has added this side-play amount in leading to, and makes the newly logical feature with the energy distribution of following of band.And this wave filter has been inherited broadband filter without the feature of estimating in advance stream.This novel dispersion relation bandpass filter of digital proof is applicable to the filtering of sea clutter image under dynamic environment by experiment.

Summary of the invention

The object of the present invention is to provide a kind of promote wave inversion accuracy carry out X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter.

The object of the present invention is achieved like this:

(1) radar image data collection: gather marine site, N width space clutter consecutive image, and synchronous recording course and the speed of a ship or plane information of boats and ships at that time;

(2) radar image pre-service: selected distance stem to angle is , medium filtering, as the subimage sequence of analyzing, is carried out to rectangular area in the rectangular area that pixel is N1*N1, uses the method for closest approach interpolation to obtain image Cartesian coordinates;

(3), to the image sequence application Fourier transform under cartesian coordinate system, obtain the three-dimensional wave number frequency image spectrum of radar image;

(4) ocean wave spectrum information extraction: use dispersion relation bandpass filter to carry out filtering processing to three-dimensional wave number frequency image spectrum, obtain three-dimensional wave image spectrum; Frequency field in this image spectrum is carried out to integration, obtain ocean waves wave number image spectrum; Use modulation transfer function to obtain ocean waves spectrum.

(5) Wave Information inverting: utilize ocean wave spectrum Inversion Calculation to obtain ocean wave parameter, comprise significant wave height, crest cycle, crest peak to.

Described step (2) comprising:

(2.1) selected distance stem to angle is pixel is that the rectangular area of N1*N1 is as the subimage sequence of analysis;

(2.2) analyzed area of N width image is rotated in the following manner:

Wherein, it is the frame heart angle of piece image; for the angle in course corresponding to N width image and the first width course; it is the frame heart angle after i width image rotation;

(2.3) the 2-D nonlinear smoothing medium filtering of radar image analytical sequence application 3 × 3 templates of step (2.1) being chosen;

(2.4), for the every bit in radar original coordinate system, utilize closest approach interpolation method to obtain the Cartesian coordinates of image.

Described step (4) comprising:

(4.1) speed of a ship or plane data of record in step (1) are transformed under the cartesian coordinate system at analyzed area place:

Wherein, for the frame heart angle of choosing, u _shipfor ship's speed;

(4.2) to the three-dimensional wave number frequency energy spectrum I (k asking in step (3) _x, k _y, ω) and carry out filtering operation

Wherein,

${\mathrm{\ω}}_n=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

${\mathrm{\ω}}_p=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

ω _n, ω _pfor the upper and lower border of filter band; U _maxmaximum flow of water speed; it is known quantity;

ω _n, ω _pdrawn by following computation process:

therefore dispersion relation

Sea surface velocity maximum can value be U _max, the logical border of band of wave filter is

Due to ?

${\mathrm{\ω}}_n=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

${\mathrm{\ω}}_p=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max }.$

Beneficial effect of the present invention is: owing to still lacking at present the technological means of accurately estimating stream both at home and abroad, estimate stream effect all not ideal enough, flow rate and direction evaluated error has had a strong impact on the performance of narrow band filter, in engineering, cannot use.Wave filter of the present invention has been inherited broadband filter without the advantage of estimating stream, makes filter effect do not estimated stream Accuracy, has very high engineering use value and application prospect; In the time that Texas tower is kept in motion, tradition wideband filtered method all cannot directly be used under platform motion state, the bandwidth of its wave filter is amplified rapidly along with ship's speed increases, and algorithm complete failure under the high speed of a ship or plane makes to have comprised a large amount of random ground unrests in inversion result and disturbs.Wave filter of the present invention has retained the impact of ship's speed dispersion relations, and the bandwidth that has effectively solved traditional broadband-pass filter can increase this problem with the increase of movement velocity, makes to carry out filtering in Ship Motion situation at Texas tower; In the band logical border derivation of Novel Filter of the present invention, not to any approximate value that measures, reduce the error of calculation, can be to not being with logical border to exert an influence, bandwidth calculation is more accurate, has promoted wave inversion accuracy.

Brief description of the drawings

Fig. 1 flow velocity be zero with flow velocity non-zero situation under dispersion relation curve map with and energy distribution;

B under Fig. 2 motion state _n, B _pfilter effect is not obvious;

Fig. 3 embodiment process flow diagram;

Fig. 4 analyzed area position view;

Fig. 5 closest approach interpolation schematic diagram;

Fig. 6 ship's speed coordinate system conversion schematic diagram;

The inverting surperficial significant wave height in sea of the present invention and the contrast of WAMOS significant wave height under Fig. 7 operational configuration;

Under Fig. 8 operational configuration inverting of the present invention sea surface wave peak-to-peak to WAVEX crest peak to contrast;

Under Fig. 9 operational configuration, inverting sea surface wave cycle of the present invention and WAVEX contrast period of wave;

Sea clutter image spectrum when Figure 10 frequency is 0.74rad/s;

Figure 11 is through B _n, B _pfiltered wave image spectrum (U _max=10m/s);

Figure 12 is through ω _n, ω _pfiltered wave image spectrum (U _max=3m/s).

Embodiment

Below in conjunction with accompanying drawing, the present invention is described further:

The invention discloses a kind of method of utilizing novel wave dispersion relation bandpass filter to carry out X-band navar inverting ocean wave parameter, the ocean wave parameter of Inversion Calculation comprise significant wave height, crest cycle, crest peak to etc., the collection that described method comprises radar image, radar image pre-service, sea clutter image spectrum are obtained, Wave Information extracts and five parts of sea information inverting.Obtain sea clutter image spectrum by Fourier transform, and use novel bandpass filter that Wave energy is extracted from sea clutter image spectrum, for follow-up inverting is laid a good foundation.Compare with classic method, filtering method disclosed in this invention can extract Ocean-wave Signal accurately.Utilizing the crest peak that the inverting of sail measurement data obtains is only 7.02o to the deviation of value and WAVEX systematic survey value, the relative error of value is only 2% period of wave, illustrate the method not only customer service classic method estimate the inaccurate defect of stream but also made up conventional filter the deficiency that cannot be applicable to wave parametric inversion under dynamic environment.

The present invention includes:

Step 1, data acquisition.Gather marine site, N width space clutter consecutive image, and synchronous recording at that time ship to and ship's speed.

Step 2, image pre-service.Selected distance stem to angle is pixel is that the rectangular area of N1*N1 is as the subimage sequence of analysis.This analyzed area is carried out to medium filtering, and use the method for closest approach interpolation to obtain image Cartesian coordinates.

Step 3, to the image sequence application Fourier transform under cartesian coordinate system, obtains three-dimensional wave number frequency image spectrum.

Step 4, Wave Information extracts.Use the Novel belt bandpass filter of the present invention's design to the three-dimensional wave number frequency image spectral filter obtaining in the 3rd step, obtain three-dimensional wave image spectrum.

Step 5, Wave Information inverting.Frequency field in the three-dimensional wave image spectrum of the 4th step is carried out to integration, obtain ocean waves image spectrum.Use modulation transfer function (MTF) to obtain ocean wave spectrum.Utilize ocean wave spectrum inverting to obtain wave height, crest peak to ocean wave parameters such as, periods of wave.

Described step 2 comprises the following steps:

Step 2.1, selected distance stem to angle is pixel is that the rectangular area of Nx*Ny is as the subimage sequence of analysis;

Step 2.2, is rotated in the following manner to the analyzed area of N width image:

Wherein, it is the frame heart angle of piece image; for the angle in course corresponding to N width image and the first width course; it is the frame heart angle after i width image rotation.

Step 2.3, the 2-D nonlinear smoothing medium filtering of the radar image analytical sequence that step 2.1 is chosen application 3 × 3 templates;

Step 2.4, for the every bit in radar original coordinate system, utilizes closest approach interpolation method to obtain the Cartesian coordinates of image.

Described step 4 comprises the following steps:

Step 4.1, is transformed into the speed of a ship or plane data of record in step 1 under the cartesian coordinate system at analyzed area place.Transform mode is as follows

Wherein, for the frame heart angle of choosing, u _shipfor ship's speed.

Step 4.2, to the three-dimensional wave number frequency energy spectrum I (k asking in step 3 _x, k _y, ω) and carry out filtering operation

Wherein,

${\mathrm{\ω}}_n=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

${\mathrm{\ω}}_p=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

ω _n, ω _pfor the upper and lower border of filter band; U _maxbe the maximum flow of water speed of hypothesis, can set it; it is known quantity.

ω in described step 4.2 _n, ω _pdrawn by following computation process:

Basic dispersion relation equation is as shown in formula (4), wherein therefore dispersion relation equation becomes

Suppose that the normal condition surface velocity maximum of plunging into the commercial sea can value be U _max, the logical border of the band of wave filter of the present invention is

Due to ?

${\mathrm{\ω}}_n=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max }---\left(18\right)$

${\mathrm{\ω}}_p=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max }---\left(19\right)$

From described step 3, the present invention selects N width image to carry out inverting, corresponding N/2 frequencies omega _i(i=1...N).Each frequency chip ω _iin have N _x* N _yindividual each can be corresponding be decomposed into unique k _x, k _y, i.e. (k _x, k _y, ω _i).That is to say, each there is N frequencies omega _iin contrast should.The resolution d ω of frequency is defined as follows:

$\mathrm{d\ω}=\frac{2\mathrm{\π}}{T}$

Wherein, time scale T is T.T. corresponding to N width image.N _x* N _yfor the pixel of selected analyzed area.

The present invention in selected analyzed area taking wavenumber resolution as step-length is to k _x, k _ytravel through the resolution dk of wave number _x, dk _ybe defined as follows:

${\mathrm{dk}}_x=\frac{2\mathrm{\π}}{L_x},{\mathrm{dk}}_y=\frac{2\mathrm{\π}}{L_y}$

Wherein, L _x, L _ythe length of side of the analyzed area rectangle frame of selecting.According to formula (18), (19) to each obtain corresponding ω _n, ω _pif, ω _n≤ ω _i≤ ω _pretain this energy, the energy of other parts is set to 0.

The novel dispersion relation bandpass filter being applicable under dynamic environment the present invention being proposed below in conjunction with accompanying drawing is described in further detail.Embodiment process flow diagram is shown in Fig. 3.Wherein, the first step is image information collecting; Second step is that radar image is carried out to pre-service; The 3rd step is obtained wave number frequency spectrum by three-dimensional Fourier transform; The 4th step is to utilize the bandpass filter of invention herein to carry out filtering processing to wave number frequency spectrum; The 5th step to the seven steps be the wave height that calculates wave from filtered image spectrum, crest peak to the parameter such as period of wave.

The first step, gathers 32 width continuous sea clutter images, and recording its total length of time is T (approximately 1.5 minutes), synchronous recording course (stem to), the speed of a ship or plane (relative velocity) value.

Second step, carries out image pre-service to selected analyzed area.

(1) choosing apart from 75 ° of directions of stem, apart from the actual range of Texas tower is rectangular area (pixel resolution is 7.5m, i.e. 128*128 picture element) the composing images analytical sequence that 600 meters, real space are of a size of 960m*960m.Analyzed area position as shown in Figure 4.Because stem in navigation process has certain swing, what cause every width image selects frame position slightly variant.For the analyzed area that makes to select constant, need to be by 32 width images taking the course of the first width wherein as benchmark is rotated.This makes to select the central angle of frame near 75 °, to swing.Formula is as follows

Wherein, be the frame heart angle of piece image, be set as 75 ° herein; be course that 32 width images are corresponding and the angle in the first width course; it is the frame heart angle after i width image rotation.

(2) the 2-D nonlinear smoothing medium filtering to selected radar image analytical sequence application 3 × 3 templates

$f^{\′}(r,\mathrm{\θ})=\underset{(s,t)\∈(r,\mathrm{\θ})}{\mathrm{median}}\left\{g(s,t)\right\}---\left(21\right)$

In formula, g (s, t) is for radar image pixel point is in the image echo strength value of polar coordinates position (s, t); F'(r, θ) for image after filtering is at the gray-scale value of polar coordinates position (r, θ); The pixel point that centered by N (r, θ), point is located at (r, θ), (s, t) gets 8 points centered by (r, θ).

The template center of 3 × 3 template median filters is overlapped with certain pixel position of polar coordinate image, itself and the echo strength value of 8 adjacent picture elements points are around arranged, echo strength value in the middle of getting is assigned to the pixel of center, and template traversal full width radar image obtains the image sequence after medium filtering.

(3) utilize closest approach interpolation method to obtain image Cartesian coordinates.Every bit in radar original coordinate system is expressed as to (r, θ, z) with polar form, is expressed as (x, y, z) after being converted to Cartesian coordinate.Have

$\left\{\begin{array}{c}x=r^{*}\cos \mathrm{\θ}\\ y=r*\sin \mathrm{\θ}\end{array}\right.---\left(22\right)$

The coordinate of rectangle frame point is made as (x _i, y _i).Utilizing closest approach interpolation is exactly to allow every bit in rectangle frame all in sector region, find from its nearest its corresponding polar coordinates (r that a bit obtains _i, θ _i), and by the echo strength assignment of this point on fan-shaped to the point (x in rectangle frame _i, y _i).Appoint and get any and be designated as (x ₀, y ₀), the polar coordinate mode that obtains its closest approach is:

$\left\{\begin{array}{c}{r}_0=\mathrm{round}\left(\mathrm{sqrt}({x_0}^2+{y_0}^2)\right)\\ {\mathrm{\θ}}_0=\mathrm{round}\left(\mathrm{rem}(\arctan (y_0,x_0)+2\mathrm{\π},2\mathrm{\π})\right)\end{array},-,-,-,\left(23\right)\right.$

Wherein, rem () is MOD function, and round () is to closest approach bracket function.By (r under polar coordinates ₀, θ ₀) put corresponding echo strength value assignment to the point (x in rectangle frame ₀, y ₀) just complete the coordinate conversion of radar image analyzed area.Fig. 5 is closest approach interpolation schematic diagram.

The 3rd step, to image sequence use Fourier transform, becomes wave number frequency field by the image in space time territory, obtains three-dimensional wave number frequency image spectrum.

Extra large clutter sequence η (x, y, t) is carried out to 3 d-dem Fourier transform, as follows:

$F(k_x,k_y,\mathrm{\ω})=\underset{0}{\overset{L_y}{\mathrm{\Σ}}}\underset{0}{\overset{L_x}{\mathrm{\Σ}}}\underset{0}{\overset{T}{\mathrm{\Σ}}}\mathrm{\η}(x,y,t)\exp [-2\mathrm{\πi}(k_{x}x/L_x+k_{y}y/L_y+\mathrm{\ωt}/T)---\left(24\right)$

Wherein, F (k _x, k _y, ω) and for obtaining three-dimensional wave number-frequency image spectrum after conversion.The analyzed area of selecting is herein that the length of side is L _x* L _yrectangle frame, wherein L _x=L _y=960m; Time scale T was 32 T.T.s corresponding to width image; k _xand k _ybe respectively the wave number of x and y direction, ω is frequency.

Three-dimensional wave number frequency energy spectrum is defined by formula below:

$I(k_x,k_y,\mathrm{\ω})=\frac{1}{L_x{L}_{y}T}{\left|F(k_x,k_y,\mathrm{\ω})\right|}^2---\left(25\right)$

The 4th step, utilizes the bandpass filter based on dispersion relation of invention herein to carry out filtering.

(1) read the ship's speed data of obtaining in first step information acquisition, be transformed under the cartesian coordinate system at analyzed area place.Transform schematic diagram according to flow velocity, Fig. 6, transform mode is as follows

Wherein, for bringing 75 ° into, u in the frame heart angle of choosing _shipfor ship's speed.Because stem is to towards x axle positive dirction, is equivalent to water and flows to x axle negative direction.Therefore in Fig. 6 ship's speed towards x axle negative direction.

(2) to three wave number frequency energy spectrum I (k that ask in the 3rd step _x, k _y, ω) and carry out filtering operation.For the real-time extra large surface velocity with respect to Texas tower is incorporated in the border of dispersion relation bandpass filter, the border of wave filter can be changed with the change of flow velocity, according to the dispersion relation design bandpass filter of wave

ω _nthe frequency lower boundary of wave filter, ω _pthe frequency coboundary of wave filter, with known terms, U _maxthat the maximal rate that under normal condition, current can reach is set as 3m/s herein.And ? therefore ω _n, ω _pbecome

${\mathrm{\ω}}_n=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

${\mathrm{\ω}}_p=\sqrt{g\left|\stackrel{\→}{k}\right|\tanh \left(\right|\stackrel{\→}{k}\left|d\right)}+k_x\·u_x+k_y\·u_y-\left|\stackrel{\→}{k}\right|\·U_{\max };$

The present invention selects 32 width images to carry out inverting, corresponding 32 frequencies omega _i(i=1...32).ω herein _iincrease progressively taking frequency resolution d ω as step-length, the resolution d ω of frequency is defined as follows:

$\mathrm{d\ω}=\frac{2\mathrm{\π}}{T}$

Wherein, time scale T was 32 T.T.s corresponding to width image.

Because the real space of analyzed area rectangle frame is of a size of 960m*960m, range resolution is decided to be 7.5m herein, therefore each frequencies omega _iin have 128*128

128*128 is the pixel of selected analyzed area.

The present invention in selected analyzed area taking wavenumber resolution as step-length is to k _x, k _ytravel through k _x, k _yspan be (0.4189,0.4123).The resolution dk of wave number _x, dk _ybe defined as follows:

${\mathrm{dk}}_x=\frac{2\mathrm{\π}}{L_x},{\mathrm{dk}}_y=\frac{2\mathrm{\π}}{L_y}$

Wherein, L _x=960m, L _y=960m is the length of side of the analyzed area rectangle frame of selection.According to formula (18), (19) to each obtain corresponding ω _n, ω _pif, ω _n≤ ω _i≤ ω _pretain this energy, the energy of other parts is set to 0.

The 5th step, three-dimensional wave number frequency energy spectrum I (k _x, k _y, ω) and after filtering, obtain wave 3-D view spectrum F (k _x, k _y, ω), it is carried out to integration to frequencies omega and can obtain two dimensional image spectrum F (k _x, k _y).Because Fourier transform is about the symmetry of initial point, if will lose phase velocity directional information at whole ω territory integration, therefore only integration is carried out in ω > 0 (or ω < 0) region, be defined as follows:

$I(k_x,k_y)={\&Integral;}_{\mathrm{\&omega;}>0}E(k_x,k_y,\mathrm{\&omega;})\mathrm{d\&omega;}---\left(28\right)$

The 6th step, is used modulation transfer function that image spectrum is converted to ocean wave spectrum

Due to the non-linearity of X-band marine radar image-forming mechanism, between two dimensional image spectrum and actual ocean wave spectrum, can there is certain difference, present nonlinear relationship, show more obviously in the larger region of wave number.Therefore, need to carry out corresponding conversion to image spectrum by modulation transfer function and obtain ocean wave spectrum.Computing formula is as follows:

$E(k_x,k_y)={\left|M(k_x,k_y)\right|}^2\&CenterDot;F(k_x,k_y)---\left(29\right)$

E (k in formula _x, k _y) be ocean waves, | M (k _x, k _y) | ²for modulation transfer function.

Modulation transfer function does not have corresponding theoretical resolvant at present, but utilizes great number tested data matching to obtain.Below formula the modulation transfer function definition of experience:

|M(k _x,k _y)| ²≈k ^-β(30)

Wherein β is experience factor, generally, in 1.2 left and right values, is widely used.

The 7th step, ask for wave height, crest peak to and period of wave.

(1) inversion method of wave height

The method of the X-band radar calculating significant wave height that the present invention uses is that the linear principle of square root based on significant wave height and signal to noise ratio (S/N ratio) realizes,

$H_s=A+B\sqrt{\mathrm{SNR}}---\left(31\right)$

Wherein, A, B is the coefficient of fitting a straight line, SNR is signal to noise ratio (S/N ratio).Undetermined coefficient A and B are the crucial empirical parameters that obtains wave significant wave height, need to demarcate by long-term a large amount of test figure.

The definition of signal to noise ratio (S/N ratio) and account form have a variety of, in this article, are defined as follows:

$\mathrm{SNR}=\frac{\mathrm{SIG}}{\mathrm{BGN}}---\left(32\right)$

Wherein, SIG is ocean wave spectrum energy, and BGN is ground unrest energy.Can be defined by formula below at deep water situation SIG and BGN:

$\mathrm{SIG}=\underset{1}{\overset{N_x}{\mathrm{\&Sigma;}}}\underset{1}{\overset{N_y}{\mathrm{\&Sigma;}}}E(k_x,k_y){\mathrm{dk}}_x{\mathrm{dk}}_y---\left(33\right)$

$\mathrm{BGN}=\underset{1}{\overset{N_x}{\mathrm{\&Sigma;}}}\underset{1}{\overset{N_y}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}I(k_x,k_y,{\mathrm{\&omega;}}_i)\&CenterDot;\mathrm{d\&omega;}{\mathrm{dk}}_y{\mathrm{dk}}_x-\underset{1}{\overset{N_x}{\mathrm{\&Sigma;}}}\underset{1}{\overset{N_y}{\mathrm{\&Sigma;}}}I(k_x,k_y){\mathrm{dk}}_y{\mathrm{dk}}_x---\left(34\right)$

(2) crest peak to the inversion method of period of wave

The transformational relation of the energy spectrum of wave number under energy spectrum of wave number and polar coordinates under rectangular coordinate is as follows

E(k,θ)＝E(k _x,k _y)k (35)

Wherein

$k=\sqrt{k_x^2+k_y^2},$ θ＝arctan2(k _y,k _x)+2π。

The angular range obtaining is between [π, π].

Ocean wave spectrum E (k, θ) is converted to directional wave spectra E (f, θ):

$E(f,\mathrm{\&theta;})=E(k,\mathrm{\&theta;})\frac{\mathrm{dk}}{\mathrm{df}}=\frac{k}{2\mathrm{\&pi;}}\frac{\&PartialD;k}{\&PartialD;\mathrm{\&omega;}}E(k,\mathrm{\&theta;})---\left(36\right)$

Further, for the direction θ integration in directional spectrum E (f, θ), just can obtain one dimension frequency spectrum,

$S\left(f\right)={\&Integral;}_0^{2\mathrm{\&pi;}}E(f,\mathrm{\&theta;})\mathrm{d\&theta;}---\left(37\right)$

Further, for the frequency f integration in directional spectrum E (f, θ), just can obtain one dimension directional spectrum,

$S\left(\mathrm{\&theta;}\right)={\&Integral;}_0^{f_{\max }}E(f,\mathrm{\&theta;})\mathrm{df}---\left(38\right)$

Obtain S (θ) and S (f) after, can calculate cycle of wave, the parameters such as wave direction.For one dimension directional spectrum S (θ), suppose that the coordinate of wave spectrum peak value position is θ, there is wave direction to be

Pdir＝θ

From one dimension spectrum S (f), obtain obtaining crest cycle T _pwith main ripple wavelength X _p.

T _p＝1/f _p

Wherein, f _pfor frequency corresponding to maximal value of S (f), can be calculated by formula below

$f_p=\frac{{\&Integral;}_{f_1}^{f_2}S\left(f\right)\mathrm{fdf}}{{\&Integral;}_{f_1}^{f_2}S\left(f\right)\mathrm{df}}---\left(39\right)$

Here f ₁and f ₂determined by 80% of the ceiling capacity in energy spectrum.Again by

ω _p＝2π/T _P，k _p＝2π/λ _P

Substitution dispersion relation, can try to achieve

${\mathrm{\&lambda;}}_P=\frac{g}{2\mathrm{\&pi;}\&CenterDot;{f_p}^2}---\left(40\right)$

The bandpass filter based on dispersion relation being applicable under dynamic environment that application the present invention proposes is carried out filtering processing to 22-29 day in November, 2011 at the underway dynamic X-band radar echo of Bohai Offshore.The antenna height of X-band radar is 40 meters, and average swing circle is 2.39 wonderful, is operated in short pulse pattern, and radius of investigation is about 2 km marine sites.Within every 3 minutes, gather one group of data, every group comprises 32 width images.Use significant wave height that on-the-spot WAMOS system the provides reference value as wave height, WAMOS output in every 2 minutes is with this result; The period of wave providing using WAVEX radar and crest peak to as cycle and crest peak to reference value, WAVEX output in every 4 minutes is with this result.Experimental image regional center is positioned at stem 75 and spends direction, area 960*960 square metre, 20 meters of this marine site mean depths.For reduce error as far as possible, this experiment is chosen the inverting data that obtain in 20 minutes to wave height, and period of wave and crest peak are to carrying out arithmetic mean.And and WAMOS, within WAVEX20 minute, the numerical value after average compares.Fig. 7 is the comparison diagram that application the present invention's significant wave height that inverting obtains under dynamic environment and WAMOS record significant wave height, and its tendency is basic identical as seen, and the root-mean-square error of standard shown in table two is the 0.18 standard 0.5m lower than setting; Fig. 8 is surveyed contrast by application inverting of the present invention the crest wind direction and the WAVEX that obtain, deviation shown in table two be 7.02 ° lower than 20 ° of the standards of setting; Fig. 9 is surveyed contrast by application the inverting of the present invention Periods and the WAVEX that obtain, and relative standard deviation shown in table two is 2% standard 10% lower than setting.Therefore this test is effective, and explanation has reached desirable filter effect.

Choose certain a slice observation filter effect in 32 frequency chips, this sheet frequency is 0.74rad/s.Figure 10 is the extra large clutter energy spectrum under this frequency, and the Ocean-wave Signal image obtaining after the bandpass filter of the present invention's design is composed as Figure 11.The visible current main flow wave filter B of comparison Figure 12 _n, B _powing to being provided with larger U under motion state _maxcannot play filter effect, and the designed filter filtering of the present invention is respond well, has retained the energy of most of real number and wave.

By above-mentioned analysis and example calculation, the bandpass filter of the present invention's design is simple and effective.

Table one: experiment parameter setting

Parameter	Numerical value
		Frequency resolution Δ ω	0.082s -1
Wavenumber resolution Δ k	0.0065m -1
		Space x is to resolution ax x	7.5m
Space y is to resolution ax y	7.5m
		Nonlinear energy correction index β	1.2

Table two: navigation experimental error statistics

Error statistics	Wave height	Crest peak to	Period of wave
				Hrb average	1.34	66.59	6.18
WAMOS/WAVEX average	1.34	69.2	6.23
				Mean difference	0	2.05	0.05
Mean square deviation	0.03	49.28	0.03
				Standard root-mean-square error-deviation	0.18	7.02	0.18
Relative standard deviation			2.00％

Claims (3)

1. carry out an X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter, it is characterized in that:

(1) radar image data collection: gather marine site, N width space clutter consecutive image, and synchronous recording course and the speed of a ship or plane information of boats and ships at that time;

(2) radar image pre-service: selected distance stem to angle is pixel be the rectangular area of N1*N1 as the subimage sequence of analyzing, medium filtering is carried out in rectangular area, use the method for closest approach interpolation to obtain image Cartesian coordinates;

(3), to the image sequence application Fourier transform under cartesian coordinate system, obtain the three-dimensional wave number frequency image spectrum of radar image;

(4) ocean wave spectrum information extraction: use dispersion relation bandpass filter to carry out filtering processing to three-dimensional wave number frequency image spectrum, obtain three-dimensional wave image spectrum; Frequency field in this image spectrum is carried out to integration, obtain ocean waves wave number image spectrum; Use modulation transfer function to obtain ocean waves spectrum;

(5) Wave Information inverting: utilize ocean wave spectrum Inversion Calculation to obtain ocean wave parameter, comprise significant wave height, crest cycle, crest peak to.

2. one according to claim 1 is carried out X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter, it is characterized in that, described step (2) comprising:

(2.1) selected distance stem to angle is pixel is that the rectangular area of N1*N1 is as the subimage sequence of analysis;

(2.2) analyzed area of N width image is rotated in the following manner:

Wherein, it is the frame heart angle of piece image; for the angle in course corresponding to N width image and the first width course; it is the frame heart angle after i width image rotation;

(2.3) the 2-D nonlinear smoothing medium filtering of radar image analytical sequence application 3 × 3 templates of step (2.1) being chosen;

(2.4), for the every bit in radar original coordinate system, utilize closest approach interpolation method to obtain the Cartesian coordinates of image.

3. one according to claim 1 and 2 is carried out X-band navar inverting ocean wave parameter method based on novel wave dispersion relation bandpass filter, it is characterized in that, described step (4) comprising:

(4.1) speed of a ship or plane data of record in step (1) are transformed under the cartesian coordinate system at analyzed area place:

Wherein, for the frame heart angle of choosing, u _shipfor ship's speed;

(4.2) to the three-dimensional wave number frequency energy spectrum I (k asking in step (3) _x, k _y, ω) and carry out filtering operation

Wherein,

${\mathrm{\&omega;}}_n=\sqrt{g\left|\stackrel{\&RightArrow;}{k}\right|\tanh \left(\right|\stackrel{\&RightArrow;}{k}\left|d\right)}+k_x\&CenterDot;u_x+k_y\&CenterDot;u_y-\left|\stackrel{\&RightArrow;}{k}\right|\&CenterDot;U_{\max };$

${\mathrm{\&omega;}}_p=\sqrt{g\left|\stackrel{\&RightArrow;}{k}\right|\tanh \left(\right|\stackrel{\&RightArrow;}{k}\left|d\right)}+k_x\&CenterDot;u_x+k_y\&CenterDot;u_y-\left|\stackrel{\&RightArrow;}{k}\right|\&CenterDot;U_{\max };$

ω _n, ω _pfor the upper and lower border of filter band; U _maxmaximum flow of water speed; it is known quantity;

ω _n, ω _pdrawn by following computation process:

therefore dispersion relation

Sea surface velocity maximum can value be U _max, the logical border of band of wave filter is

Due to ?

${\mathrm{\&omega;}}_n=\sqrt{g\left|\stackrel{\&RightArrow;}{k}\right|\tanh \left(\right|\stackrel{\&RightArrow;}{k}\left|d\right)}+k_x\&CenterDot;u_x+k_y\&CenterDot;u_y-\left|\stackrel{\&RightArrow;}{k}\right|\&CenterDot;U_{\max };$

${\mathrm{\&omega;}}_p=\sqrt{g\left|\stackrel{\&RightArrow;}{k}\right|\tanh \left(\right|\stackrel{\&RightArrow;}{k}\left|d\right)}+k_x\&CenterDot;u_x+k_y\&CenterDot;u_y-\left|\stackrel{\&RightArrow;}{k}\right|\&CenterDot;U_{\max }.$