CN103969643B - One carries out X-band pathfinder inverting ocean wave parameter method based on novel wave dispersion relation band filter - Google Patents

Abstract

The invention belongs to Remote sensing of ocean wave technical field, be specifically related to that a kind of sea clutter image utilizing pathfinder to obtain carries out ocean wave parameter inverting carries out X-band pathfinder inverting ocean wave parameter method based on novel wave dispersion relation band filter.The present invention includes: radar image data collection；Radar image pretreatment；To the image sequence application Fourier transformation under cartesian coordinate system, obtain the three-dimensional wave number frequency image spectrum of radar image；Ocean wave spectrum information retrieval；Wave Information inverting.The wave filter of the present invention remains the impact of ship's speed dispersion relations, and the bandwidth efficiently solving Conventional wide band bandpass filter can increase this problem with the increase of movement velocity so that can be filtered in the case of radar platform is with Ship Motion；The band of the Novel Filter of the present invention leads in deriving in border and does not measure approximation to any, reduces calculating error, will not produce impact in border logical on band, and bandwidth calculation is more accurate, improves wave inversion accuracy.

Description

One carries out X-band pathfinder inverting ocean wave parameter method based on novel wave dispersion relation band filter

Technical field

The invention belongs to Remote sensing of ocean wave technical field, be specifically related to a kind of sea clutter image utilizing pathfinder to obtain and carry out wave Parametric inversion carry out X-band pathfinder inverting ocean wave parameter method based on novel wave dispersion relation band filter.

Background technology

The echo of X-band radar is analyzed, it is possible to obtain the parameter of wave.But except containing in this echo Ocean-wave Signal is outer possibly together with substantial amounts of noise, such as co-channel interference, background noise, sleet interference and object echo etc..These are made an uproar Sound will have a strong impact on the parametric inversion of wave, make result distortion.So needing to remove these noises letter from the echo of radar Number and extract useful Ocean-wave Signal.Owing to this signal meets the dispersion relation of gravitational wave, it is possible to utilize dispersion relation Structure 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,ω)

Wave meets this theory of dispersion relation of gravitational wave and was proposed in 1978 by LeBlond and Mysak.See reference Document [1] (Paul H.LeBlond, Lawrence A.Mysak.Waves in the ocean [M] .1978, Elsevier.)

Assuming that in selected analyzed area sea wave field and flow field be space uniform and the time stable, then the color of single order gravitational wave Dissipate relationship description as follows

${\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 being one Flow relative to the surface of antenna platform of radar speedWhereinSo introduce a Doppler frequency shift Dispersion relation becomes:

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

Then

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

Above-mentioned for radar with sea surface there is not relative velocity with there is relative velocity under dispersion relation curve.Curve a is flow velocity Being dispersion relation curve when 0, it is symmetrical about the longitudinal axis.Curve b is the dispersion relation curve that there is Doppler frequency shift.Contour Line represents the distribution of energy.

Show according to document, the most existing band filter master based on dispersion relation extracted for sea clutter image Ocean-wave Signal It is divided into two big classes, i.e. arrowband and broadband.In the development history of band filter, occur that the wave filter of narrowband versions the earliest, Subsequent to solve the deficiency of narrow band filter and occur in that broadband filter.

The principle of narrow band filter is to find suitable method to calculate the velocity vector in dispersion relation formula (3)So color Scattered relation is known.Owing to the Wave energy overwhelming majority is distributed near dispersion relation curve, therefore by bent in dispersion relation Narrow bandpass is set near line, just can obtain real Wave Information, filter noise jamming.

Narrow bandpass dispersion relation wave filter the earliest is that Reichert, Nieto et al. in 1999 propose.Reichert et al. is respectively Two kinds of similar band filter definition of thought are given, as shown in formula (4) (5) in document [2-4].

Owing to sea clutter image spectrum is obtained by 3 d-dem Fourier transformation, therefore wave numberAnd frequencies omega is all the point after discretization. In formula (4), I (k_x,k_y, ω) and it is sea clutter image spectrum；F(k_x,k_y, ω) and it is that the signal pattern after bandpass filtering is composed, this wave filter To all of wave numberCarry out traveling through the frequencies omega obtaining correspondence, and only retain the energy at this as Wave energy, other parts As noise filtering [2] [3].In formula (5), d is the depth of water,Obtained by estimating stream, ω_pRepresent the frequency of a certain frequency chip.This filter Ripple device is equally to all of wave numberCarry out traversal and obtain corresponding frequencies omega, if ω is ∈ [ω_p-1,ω_p+1], then in reserve frequency sheet Energy [4] under this wave number on p, p-1 or p, p+1 sheet.Therefore, this wave filter band is logical more logical than band described in (4) formula the widest. See reference document [2-4] (Jose Carols Nieto Borge, Konstanze Reichert, Jurgen Dittmer.Use of nautical radar as a wave monitoring instrument[J].Coastal Engineering,1999,37:331-342P.Konstanze Reichert, Jose Carols Nieto Borge, Jurgen Dittmer.WaMoS2:An operational wave monitoring system.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, above-mentioned wave filter is improved by the red grade of Shen Ji, it is contemplated that think after high-order dispersion relation that Wave energy is not only It is distributed on 1st order chromatic dispersion relation curve.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 m is the order of dispersion relation, ω_mFrequency for m order ripple.Visible as m=1 formula (6) be single order color Dissipate relation equation (3), and the ω that high-order dispersion relation is asked for_mCan increase.High-order dispersion relation border is utilized to generate new wave filter, This wave filter is defined as follows

Wherein,For the wave filter upper bound, high-order dispersion relation curve finding limit obtain, than formula (4) (5) In result increase；For wave filter lower bound, identical with the result in formula (4) (5). This wave filter still falls within narrow band filter, needs to estimate flow rate and direction before bandpass filtering, and equally to each discretization Wave numberTravel through.See reference document [5] (Shen Jihong, Li Ying, Shujuan WANG, base in Dai Yuntao, Wang Xiao Oriolus chinensis diffusus .X band radar Wave filter in high-order dispersion relation designs [J]. Chinese navigation, and 2012,4 (35): 4-7.)

Being more than to have the current situation of narrow band filter in document at present, its filter effect places one's entire reliance upon and estimates stream precision.By The factor impact such as limited in image resolution ratio, spectral power distribution is uneven, noise characteristic is complicated, the most still lacks accurately Estimating the technological means of stream, estimate stream effect the most not ideal enough, flow rate and direction estimation difference has had a strong impact on the performance of narrow band filter, Real not dispersion relation curve based on narrow bandpass can cause wave filter to remain in a large number being not belonging to the noise of wave and filter True Wave Information.

For this problem, first Rune devised based on Peak Flow Rate U in 2002_maxBroadband-pass filter [6].The party Method need not accurately estimate stream, solves the deficiency that narrow band filter is maximum.The method has become as the main flow being widely used at present Wave filter.Considering in the case of deep water, this filter method is defined as foloows:

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_maxPlunge into the commercial sea the Peak Flow Rate that surface is likely to be breached for normal condition. This filter edges formula B_n、B_PIt is derived by by the basic dispersion relation formula (3) of band flow velocity.To portion in derivation Subitem has taken approximation, first when the depth of water is more than or equal to 1/2nd wavelength, and Ke YisheSecondly because's Value is less than 1,Value be more thanValue, and wave wavelength is the biggest, the two gap is the biggest.Therefore for most while simplifying calculating Reduce error possibly, can be byBring formula (3) into, then obtain

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

Therefore when cos θ=1,

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

When cos θ=-1,

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

They distinguish lower boundary and the coboundary of correspondence wave filter.Simultaneously take account of the impact of frequency resolution and wavenumber resolution Just B can be obtained_n、B_pThe border of wave filter.See reference document [6] (Rune Gangeskar.Ocean current estimated from X-band radar sea surface images[J].IEEE transactions on geosciences and remote sensing, 2002,40:783-791.)

2012, above-mentioned broadband based on dispersion relation filtering approach was carried out improving [7] again by Jia Ruicai et al..The method with Maximum wave number k assumed_mIt is different from above-mentioned filtering mode.This wave filter is defined as

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_mComputational methods be first calculate spectrum I (k_x,k_y, ω) in maximum C, then obtain all more than 0.9C of spectrum The wave number modulus value of point, then by wave number modulus value arithmetic average, i.e. obtain | km |.The same B of this wave filter_n、B_pWave filter is similar, with master Guided wave number k_mInstead of ω²/ g so that filter bandwidht amplifies, but and not as ω²/ g is accurate.See reference document [7] (Yuan Jiangxi South, 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 both broadband filters all by frequencies omega is traveled through, ask for wave number modulus value corresponding for this ω the upper bound and After next time, the wave number under each ω is screened.This wave filter that ω is traveled through and the filter that wave number k is traveled through Ripple utensil has identical time complexity, is 0 (n²).But, bring into from analysis above, both wave filter Approximation will cause this wave filter to there are two defects: one is to be only applicable to the depth of water more than or equal to 1/2nd wavelength case；Two are Only under big wavelength, this filter filtering effect is the most accurate.

In addition, two kinds of wideband filtered methods all cannot directly use under platform kinestate, works as seen from from algorithm defines When radar platform is kept in motion, the bandwidth of wave filter will significantly be amplified, and does not has filter effect so that inverting is tied Contain substantial amounts of noise jamming in Guo, affect precision.With B_n、B_pAs a example by wave filter, when radar platform is static, U_max=3m/s, along with the increase of ship's speed, U_maxAlso should be the biggest, under band is logical to broaden therewith, and the pace of change of coboundary is more than Border.It addition, the U in this wave filter_maxAlso lead to bandwidth under kinestate increase, it is impossible to play filter effect.

As can be seen here, although broadband filter need not estimate stream, greatest hope flow velocity U is used_maxSolve and estimate stream under narrow bandpass not The problem accurately brought, but the band of wave filter is logical in a dynamic state can infinitely be amplified, it is impossible to play filter effect.

For this problem, the present invention devises a kind of novel dispersion relation band filter.Wave number is carried out time by this wave filter Go through, and belong to the category in broadband, but the bandwidth efficiently solving again Conventional wide band bandpass filter can be with the increasing of movement velocity Increase this problem greatly.By the movement velocity of radar platform in view of carrying in logical accounting equation, make the bandwidth under kinestate Identical with static state.Under kinestate, the increase of ship's speed can cause Wave energy to follow dispersion relation producing skew, and newly carry logical in add Add this side-play amount so that new band is logical has the advantages that to follow Energy distribution.And this wave filter inherit broadband filter without The feature of stream need to be estimated in advance.Prove that this novel dispersion relation band filter is applicable to sea clutter under dynamic environment by experimental data The filtering of image.

Summary of the invention

It is an object of the invention to provide and a kind of promote carrying out based on novel wave dispersion relation band filter of wave inversion accuracy X-band pathfinder inverting ocean wave parameter method.

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 the boat of synchronous recording boats and ships at that time To with speed of a ship or plane information；

(2) radar image pretreatment: selected distance stem to angle is, pixel is that the rectangular area of N1*N1 is as analysis Subimage sequence, rectangular area is carried out medium filtering, uses the method for closest approach interpolation to obtain image cartesian coordinate；

(3) to the image sequence application Fourier transformation under cartesian coordinate system, the three-dimensional wave number frequency image of radar image is obtained Spectrum；

(4) ocean wave spectrum information retrieval: use dispersion relation band filter to be filtered processing to three-dimensional wave number frequency image spectrum, Obtain three-dimensional wave image spectrum；Frequency domain in composing this image is integrated, and obtains ocean waves wave number image spectrum；Use and adjust Modulation trnasfer function MTF obtains ocean waves spectrum.

(5) Wave Information inverting: utilize ocean wave spectrum Inversion Calculation to obtain ocean wave parameter, including significant wave height, peak period, ripple Peak-to-peak to.

Described step (2) including:

(2.1) selected distance stem to angle isPixel is that the rectangular area of N1*N1 is as the subimage sequence analyzed；

(2.2) analyzed area of N width image is rotated as follows:

Wherein,It it is the frame heart angle of piece image；The course corresponding for N width image and the angle in the first width course；It it is the frame heart angle after the i-th width image rotation；

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

(2.4) for the every bit in radar original coordinate system, closest approach interpolation method is utilized to obtain the cartesian coordinate of image.

Described step (4) including:

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

Wherein,For the frame heart angle chosen, u_ShipFor ship's speed；

(4.2) to the three-dimensional wave number frequency energy spectrum I (k asked in step (3)_x,k_y, ω) and it is filtered 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、ω_pUpper and lower border for filter band；U_maxIt is maximum flow of water speed；It is The amount of knowing；

ω_n、ω_pDrawn by process calculated below:

Therefore dispersion relation

Sea surface velocity maximum can be with value as U_max, the band of wave filter leads to border and is

Due toThen

${\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 beneficial effects of the present invention is: owing to the most still lacking the technological means accurately estimating stream, estimate stream effect the most inadequate Ideal, flow rate and direction estimation difference has had a strong impact on the performance of narrow band filter, and engineering cannot use.The wave filter of the present invention Inherit broadband filter without the advantage estimating stream so that filter effect is not affected by estimating stream precision, has the highest engineering and uses It is worth and application prospect；When radar platform is kept in motion, Conventional wide band filtering method all cannot be directly at platform motion shape Using under state, the bandwidth of its wave filter is amplified rapidly along with ship's speed increases, algorithm complete failure under the high speed of a ship or plane so that inverting Result contains the interference of substantial amounts of random background noise.The wave filter of the present invention remains the impact of ship's speed dispersion relations, has The bandwidth solving to effect Conventional wide band bandpass filter can increase this problem with the increase of movement velocity so that can be at radar Platform is filtered with in the case of Ship Motion；The band of the Novel Filter of the present invention lead in deriving in border not to any measure near Like value, reducing calculating error, will not produce impact in border logical on band, bandwidth calculation is more accurate, improves wave inverting essence Degree.

Accompanying drawing explanation

Fig. 1 flow velocity be zero with flow velocity non-zero in the case of dispersion relation curve chart and its Energy distribution；

B under Fig. 2 kinestate_n、B_pFilter effect is inconspicuous；

Fig. 3 detailed description of the invention flow chart；

Fig. 4 analyzed area position view；

Fig. 5 closest approach interpolation schematic diagram；

Fig. 6 ship's speed coordinate system transition diagram；

Under Fig. 7 operational configuration, surface, inverting sea of the present invention significant wave height contrasts with WAMOS significant wave height；

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

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

Figure 10 frequency is sea clutter image spectrum during 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).

Detailed description of the invention

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

The invention discloses one utilizes novel wave dispersion relation band filter to carry out X-band pathfinder inverting ocean wave parameter Method, the ocean wave parameter of Inversion Calculation include significant wave height, peak period, crest peak to etc., described method comprises radar map The collection of picture, radar image pretreatment, sea clutter image spectrum obtain, Wave Information extracts and five parts of sea information inverting. Obtain sea clutter image by Fourier transformation to compose, and use novel band filter by Wave energy from sea clutter image is composed Extract, lay a good foundation for follow-up inverting.Compared with traditional method, filtering method disclosed in this invention can be accurate Extraction Ocean-wave Signal.Utilize inclined to value and WAVEX system measurement of crest peak that sail measurement data inverting obtains Difference is only 7.02o, and period of wave, the relative error of value was only 2%, illustrates that the method not only customer service traditional method estimates stream inaccurate Defect and compensate for conventional filter and cannot be applicable to the deficiency of 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 ship at that time to and ship's speed.

Step 2, Image semantic classification.Selected distance stem to angle isPixel is that the rectangular area of N1*N1 is as analyzing Subimage sequence.This analyzed area is carried out medium filtering, and uses the method for closest approach interpolation to obtain image cartesian coordinate.

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

Step 4, Wave Information extracts.The Novel belt bandpass filter of the use present invention design three-dimensional wave number to obtaining in the 3rd step Frequency image spectral filter, obtains three-dimensional wave image spectrum.

Step 5, Wave Information inverting.Frequency domain in the three-dimensional wave image spectrum of the 4th step is integrated, obtains Two-dimensional Sea Wave image spectrum.Modulation transfer function (MTF) (MTF) is used to obtain ocean wave spectrum.Utilize ocean wave spectrum inverting obtain wave height, crest peak to, The ocean wave parameters such as period of wave.

Described step 2 comprises the following steps:

Step 2.1, selected distance stem to angle isPixel is that the rectangular area of Nx*Ny is as the subimage sequence analyzed；

Step 2.2, rotates as follows to the analyzed area of N width image:

Wherein,It it is the frame heart angle of piece image；The course corresponding for N width image and the angle in the first width course；It it is the frame heart angle after the i-th 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 coordinate of image.

Described step 4 comprises the following steps:

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

Wherein,For the frame heart angle chosen, u_ShipFor ship's speed.

Step 4.2, to the three-dimensional wave number frequency energy spectrum I (k asked in step 3_x,k_y, ω) and it is filtered 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、ω_pUpper and lower border for filter band；U_maxThe maximum flow of water speed assumed that, can be set it；It it is known quantity.

ω in described step 4.2_n、ω_pDrawn by process calculated below:

Shown in basic dispersion relation equation such as formula (4), whereinTherefore dispersion relation equation becomes

Assuming that normal condition plunges into the commercial sea surface velocity maximum can be with value as U_max, the band of wave filter the most of the present invention leads to border and is

Due toThen

${\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, correspond to N/2 frequencies omega_i(i=1...N). Each frequency chip ω_iIn have N_x*N_yIndividualEachCan be corresponding be decomposed into unique k_x、k_y, i.e. (k_x,k_y,ω_i)。 It is to say, it is eachThere is N number of frequencies omega_iCorrespond.The resolution d ω of frequency is defined as follows:

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

Wherein, time scale T is the total time that N width image is corresponding.N_x*N_yPixel for selected analyzed area.

The present invention in selected analyzed area with wavenumber resolution for step-length to k_x、k_yTravel through, resolution dk of wave number_x、dk_y It is 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 selected.According to formula (18), (19) to eachObtain corresponding ω_n、ω_pIf, ω_n≤ω_i≤ω_pThen retaining this energy, the energy of other parts is set to 0.

The novel dispersion relation band filter being applicable to proposition of the present invention under dynamic environment below in conjunction with accompanying drawing is made further Detailed description.Detailed description of the invention flow chart is shown in Fig. 3.Wherein, the first step is image information collecting；Second step is to radar Image carries out pretreatment；3rd step obtains wave number frequency spectrum by three-dimensional Fourier transform；4th step is to utilize the band invented herein Wave number frequency spectrum is filtered processing by bandpass filter；5th step to the 7th step is for being calculated wave composing from filtered image Wave height, crest peak to the parameters such as period of wave.

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

Second step, carries out Image semantic classification to selected analyzed area.

(1) choose away from 75 ° of directions of stem, the actual range of distance radar platform be 600 meters, real space a size of Rectangular area (pixel resolution is 7.5m, i.e. 128*128 the picture element) pie graph of 960m*960m is as analytical sequence.Point Analysis regional location is as shown in Figure 4.Owing to during navigation, stem has certain swing, cause each image selects frame position slightly Variant.In order to the analyzed area making selection is constant, need on the basis of the course of the first width therein, revolve 32 width images Turn.This makes the central angle selecting frame can swing near 75 °.Formula is as follows

Wherein,It is the frame heart angle of piece image, is set as 75 ° herein；It it is the course that 32 width images are corresponding Angle with the first width course；It it is the frame heart angle after the i-th 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, (s is t) that radar image pixel point is in polar coordinate position (s, image echo strength value t) to g；F'(r, θ) for scheming after filtering As the gray value at polar coordinate position (r, θ)；Putting the pixel point at (r, θ) place centered by N (r, θ), (s t) takes centered by (r, θ) 8 points.

The template center of 3 × 3 template median filters is overlapped with certain pixel position of polar coordinate image, by itself and surrounding 8 The echo strength value of adjacent picture elements point arranges, and the echo strength value taking centre is assigned to the pixel of center, and template traversal is complete Width radar image obtains the image sequence after medium filtering.

(3) closest approach interpolation method is utilized to obtain image cartesian coordinate.By the every bit in radar original coordinate system with polar coordinate Form be expressed as (r, θ, z), be expressed as after being converted to Cartesian coordinate (x, y, z).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 set to (x_i,y_i).Utilizing closest approach interpolation is exactly to allow the every bit in rectangle frame all in sector region Find from its nearest polar coordinate (r a little obtaining its correspondence_i,θ_i), and the echo strength of this point in sector is assigned to rectangle Point (x in frame_i,y_i).Appoint and take any and be designated as (x₀,y₀), the polar coordinate mode obtaining 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 coordinate₀,θ₀) corresponding the returning of point Point (the x that intensity of wave value is assigned in rectangle frame₀,y₀) just complete the Coordinate Conversion of radar image analyzed area.Fig. 5 is closest approach Interpolation schematic diagram.

3rd step, to image sequence use Fourier transformation, becomes wave number frequency domain by the image of spatial temporal domain, obtains three-dimensional Wave number frequency image is composed.

To sea clutter sequence η (x, y, t) carry out 3 d-dem Fourier transformation, 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 selected herein is that the length of side is L_x*L_yRectangle frame, wherein L_x=L_y=960m；Time scale T is the total time that 32 width images are corresponding；k_xAnd k_yRespectively For the wave number in 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)$

4th step, utilizes the band filter based on dispersion relation invented herein to be filtered.

(1) read the ship's speed data obtained in first step information gathering, be transformed into the cartesian coordinate system at analyzed area place Under.Converting schematic diagram according to flow velocity, Fig. 6, transform mode is as follows

Wherein,Frame heart angle for choosing brings 75 ° into, u_ShipFor ship's speed.Owing to stem is to towards x-axis positive direction, be then equivalent to Water is to x-axis negative direction stream.Therefore in Fig. 6 ship's speed towards x-axis negative direction.

(2) to three the wave number frequency energy spectrum I (k asked in the 3rd step_x,k_y, ω) and it is filtered operation.In order to will be relative to thunder The surface velocity of sea in real time reaching platform is incorporated in the border of dispersion relation band filter, and the border making wave filter can be with flow velocity Change and change, design band filter according to the dispersion relation of wave

ω_nIt is the frequency lower boundary of wave filter, ω_pIt is the frequency coboundary of wave filter,WithIt is known terms, U_maxIt is normal The maximal rate that under situation, current can reach is set as 3m/s herein.AndThenTherefore ω_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, correspond to 32 frequencies omega_i(i=1...32).ω herein_iWith frequency discrimination Rate d ω is that step-length is incremented by, and the resolution d ω of frequency is defined as follows:

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

Wherein, time scale T is the total time that 32 width images are corresponding.

Due to the real space a size of 960m*960m of analyzed area rectangle frame, range resolution ratio is set to 7.5m, the most herein Individual frequencies omega_iIn have 128*128

128*128 is the pixel of selected analyzed area.

The present invention in selected analyzed area with wavenumber resolution for step-length to k_x,k_yTravel through, k_x,k_ySpan be (-0.4189,0.4123).Resolution dk of wave number_x、dk_yIt is 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 selected.According to formula (18), (19) To eachObtain corresponding ω_n、ω_pIf, ω_n≤ω_i≤ω_pThen retaining this energy, the energy of other parts is set to 0.

5th step, three-dimensional wave number frequency energy spectrum I (k_x,k_y, ω) filtered after obtain wave 3-D view spectrum F (k_x,k_y, ω), by it It is integrated frequencies omega obtaining two dimensional image spectrum F (k_x,k_y).Owing to Fourier transformation is about the symmetry of initial point, if Whole ω volume integration will lose phase velocity directional information, therefore only ω ＞ 0 (or ω ＜ 0) region is integrated, is defined as follows:

$I(k_x,k_y)={\∫}_{\mathrm{\ω}>0}E(k_x,k_y,\mathrm{\ω})\mathrm{d\ω}---\left(28\right)$

6th step, uses 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, can exist certain between two dimensional image spectrum and the ocean wave spectrum of reality Difference, present non-linear relation, show to become apparent from the region that wave number is bigger.Accordingly, it would be desirable to by modulation transmission letter Image spectrum is converted and obtains ocean wave spectrum by number MTF accordingly.Computing formula is as follows:

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

E (k in formula_x,k_y) it is ocean waves, | M (k_x,k_y)|²For modulation transfer function (MTF).

The theoretical derivation form that modulation transfer function (MTF) is the most corresponding at present, but utilize great number tested data matching to obtain.Under Face formula gives the modulation transfer function (MTF) definition of experience:

|M(k_x,k_y)|²≈k^-β(30)

Wherein β is empirical coefficient, typically in about 1.2 values, is widely used.

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

(1) inversion method of wave height

The X-band radar that the present invention uses calculate the method for significant wave height be square root based on significant wave height and signal to noise ratio linearly The principle of relation realizes, i.e.

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

Wherein, A, B are the coefficients of fitting a straight line, and SNR is signal to noise ratio.Undetermined coefficient A and B are to obtain wave significant wave height Crucial empirical parameter, needs to be demarcated by long-term substantial amounts of test data.

The definition of signal to noise ratio and calculation 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 background noise energy.Can be by following in deep water situation SIG and BGN Formula defines:

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

$\mathrm{BGN}=\underset{1}{\overset{N_x}{\mathrm{\Σ}}}\underset{1}{\overset{N_y}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}I(k_x,k_y,{\mathrm{\ω}}_i)\·\mathrm{d\ω}{\mathrm{dk}}_y{\mathrm{dk}}_x-\underset{1}{\overset{N_x}{\mathrm{\Σ}}}\underset{1}{\overset{N_y}{\mathrm{\Σ}}}I(k_x,k_y){\mathrm{dk}}_y{\mathrm{dk}}_x---\left(34\right)$

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

Energy spectrum of wave number under rectangular coordinate is as follows with the transformational relation of the energy spectrum of wave number under polar coordinate

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 obtained is between [-π, π].

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

$E(f,\mathrm{\θ})=E(k,\mathrm{\θ})\frac{\mathrm{dk}}{\mathrm{df}}=\frac{k}{2\mathrm{\π}}\frac{\∂k}{\∂\mathrm{\ω}}E(k,\mathrm{\θ})---\left(36\right)$

The direction θ integration being further directed in directional spectrum E (f, θ), it is possible to obtain one-dimensional frequency spectrum, i.e.

$S\left(f\right)={\∫}_0^{2\mathrm{\π}}E(f,\mathrm{\θ})\mathrm{d\θ}---\left(37\right)$

It is further directed to the frequency f integration in directional spectrum E (f, θ), it is possible to obtain one-dimensional directional spectrum, i.e.

$S\left(\mathrm{\θ}\right)={\∫}_0^{f_{\max }}E(f,\mathrm{\θ})\mathrm{df}---\left(38\right)$

After obtaining S (θ) and S (f), the cycle of wave can be calculated, the parameter such as wave direction.For one-dimensional directional spectrum S (θ), The coordinate assuming spectral peaks position is θ, then have the wave direction to be

Pdir=θ

Obtain obtaining peak period T from one-dimensional spectrum S (f)_pWith main ripple wavelength X_P。

T_p=1/f_p

Wherein, f_pFor the frequency that the maximum of S (f) is corresponding, can be calculated by formula below

$f_p=\frac{{\∫}_{f_1}^{f_2}S\left(f\right)\mathrm{fdf}}{{\∫}_{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

Substitute into dispersion relation, can try to achieve

${\mathrm{\λ}}_P=\frac{g}{2\mathrm{\π}\·{f_p}^2}---\left(40\right)$

The band filter based on dispersion relation being applicable under dynamic environment that the application present invention proposes is in November, 2011 22-29 Day is filtered processing at the underway dynamic X-band radar echo of Bohai Offshore.The antenna height of X-band radar is 40 Rice, average swing circle is 2.39 wonderful, is operated in burst mode, and radius of investigation is about 2 km marine sites.Within every 3 minutes, gather One group of data, often group comprises 32 width images.The significant wave height reference value as wave height of on-the-spot WAMOS system offer is provided, WAMOS output in every 2 minutes is with this result；The period of wave provided using WAVEX radar and crest peak are to as cycle and crest Peak to reference value, every 4 minutes of WAVEX output is with this result.Experimental image regional center is positioned at 75 degree of directions of stem, face Long-pending 960*960 square metre, this marine site mean depth 20 meters.For reducing error as far as possible, this experiment obtains in choosing 20 minutes Inverting data are to wave height, and period of wave and crest peak are to carrying out arithmetic average.And and WAMOS, after WAVEX20 minute is average Numerical value compare.Fig. 7 is that the application present invention significant wave height that inverting obtains under dynamic environment records significant wave with WAMOS High comparison diagram, it is seen that its tendency is essentially identical, standard root-mean-square error shown in table two is 0.18 standard 0.5m being less than setting； The crest wind direction that Fig. 8 obtains by applying inverting of the present invention and WAVEX are surveyed contrast, and deviation shown in table two is 7.02 ° and is less than setting Standard 20 °；The Periods that Fig. 9 obtains by applying inverting of the present invention and WAVEX are surveyed contrast, relative standard shown in table two Error is 2% standard 10% being less than setting.Therefore this test is effective, and illustrates to have reached preferable filter effect.

Choosing a piece of observation filter effect of certain in 32 frequency chips, this sheet frequency is 0.74rad/s.Figure 10 is under this frequency Sea clutter energy spectrum, the Ocean-wave Signal image obtained after the band filter of present invention design is composed such as Figure 11.Comparison Figure 12 Visible current main flow wave filter B_n、B_POwing to being provided with bigger U under kinestate_maxCannot play filter effect, And filter filtering designed by the present invention is respond well, remain the energy of the real number of major part and wave.

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

Table one: experiment parameter is arranged

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 deviation	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 (2)

1. carrying out an X-band pathfinder inverting ocean wave parameter method based on novel wave dispersion relation band filter, it is special Levy and be:

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

(2) radar image pretreatment: selected distance stem to angle isPixel is that the rectangular area of N1*N1 is as analysis Subimage sequence, rectangular area is carried out medium filtering, uses the method for closest approach interpolation to obtain image cartesian coordinate；

(3) to the image sequence application Fourier transformation under cartesian coordinate system, the three-dimensional wave number frequency image of radar image is obtained Spectrum；

(4) ocean wave spectrum information retrieval: use dispersion relation band filter to be filtered processing to three-dimensional wave number frequency image spectrum, Obtain three-dimensional wave image spectrum；Frequency domain in composing this image is integrated, and obtains ocean waves wave number image spectrum；Use and adjust Modulation trnasfer function MTF obtains ocean waves spectrum；

(5) Wave Information inverting: utilize ocean wave spectrum Inversion Calculation to obtain ocean wave parameter, including significant wave height, peak period, ripple Peak-to-peak to；Described step (2) including:

(2.1) selected distance stem to angle isPixel is that the rectangular area of N1*N1 is as the subimage sequence analyzed；

(2.2) analyzed area of N width image is rotated as follows:

Wherein,It it is the frame heart angle of piece image；The course corresponding for N width image and the angle in the first width course；It it is the frame heart angle after the i-th width image rotation；

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

(2.4) for the every bit in radar original coordinate system, closest approach interpolation method is utilized to obtain the cartesian coordinate of image.

It is anti-that one the most according to claim 1 carries out X-band pathfinder based on novel wave dispersion relation band filter Drill ocean wave parameter method, it is characterised in that described step (4) including:

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

Wherein,For the frame heart angle chosen, u_ShipFor ship's speed；

(4_.2) to the three-dimensional wave number frequency energy spectrum I (k asked in step (3)_x,k_y, ω) and it is filtered operation

Wherein,

ω_n、ω_pUpper and lower border for filter band；U_maxIt is maximum flow of water speed；It is The amount of knowing；k_xAnd k_yIt is respectively the wave number in x and y direction,D is the depth of water；

ω_n、ω_pDrawn by process calculated below:

Therefore dispersion relation

Sea surface velocity maximum can be with value as U_max, the band of wave filter leads to border and is

Due toThen