CN102073037B - Iterative current inversion method based on adaptive threshold selection technique - Google Patents

Abstract

The present invention provides an iterative current inversion method based on an adaptive threshold selection technique, which comprises the following steps: (1) collecting radar images and carrying out A/D (analogue/digital) conversion to obtain 32 original radar images; (2) obtaining images within Cartesian coordinates and then obtaining an image spectrum through Fourier transformation; (3), determining initial estimates, only taking influence of 0-order secondary waves into consideration, and obtaining a rough flow through inversion; and (4) iterating fitting flow information by adopting a average weighted least squares method based on the adaptive threshold selection technique, taking influence of the 0-order secondary waves and 1-order secondary waves into consideration. By adopting the method, wave signals can be automatically identified, wave energy can be accurately separated from the image spectrum, the insufficient empirical value can be made up, the current trend can be better reflected, and the purposes of fewer errors and better stability can be achieved.

Description

Iteration ocean current inversion method based on the adaptive threshold selecting technology

Technical field

The present invention relates to a kind of remote sensing technology, particularly technological based on the ocean remote sensing of X-band navar.

Background technology

Flow measurement appearance (current meter, acoustics ocean current profile instrument) can be measured the stream information on surface, sea; But it receives the destruction of artificial or disaster easily; So only use at present in the Yu Haiyang scientific investigation, China does not also form real-time and effective businessization monitoring capability to the flow field, sea at present.

X-band radar not only can be kept watch on mobile target, and it is operated in short pulse pattern following time, and it can measure ocean current and wave information, and the measurement of flow velocity (size, direction) is the gordian technique of X-band radar wave monitoring.Compare with other sea sensor (wave buoy, flow measurement appearance), the X-band radar investigative range is wider, and is safe, safeguards advantages such as simple, can be used as the substitute of wave buoy and flow measurement appearance.Domestic existing how tame unit carries out Primary Study with regard to X-band radar flow measurement problem.At present, ocean current inversion technique commonly used realizes as follows:

(1) (x, y t) arrive wave number frequency field F through 3 dimension discrete Fourier transformations for the time of radar image and spatial sequence g ⁽³⁾(k _x, k _y, ω), wherein, x, y represent two components in radar image space respectively, t representes the time shaft of radar image, k _x, k _yBe respectively the x and the y component of wave number, ω is the frequency of wave.

(2) because the wave number and the frequency of wave itself meet dispersion relation, as follows:

$S\left(k\right)=\sqrt{g\left|k\right|\tanh \left(\right|k\left|h\right)}+k\·u---\left(1\right)$

Wherein g is an acceleration of gravity, k=(k _x, k _y) be wave number, h is the depth of water, u is a flow velocity.

In the image spectrum, select an energy threshold C; Require energy value to be generally dozens of greater than the number of the discrete point of C; Use these conventional point,, use the least square fitting of energy weighting to obtain rough flow velocity according to the characteristic of wave number of wave own and frequency.

(3) at first according to the flow velocity that obtains, confirm 0,1 rank dispersion relation.Rule of thumb choose energy threshold C then ₁, require C ₁Much smaller than C, choose greater than C ₁The spread in energy point, judge respectively that then it still is 1 order ripple that these points belong to 0 rank, simulate flow velocity according to the energy weighted least-squares method, new flow velocity is confirmed 0,1 new rank dispersion relation, repeat this step, will obtain constantly flow velocity accurately.

$w_p=(p+1)\sqrt{\frac{\mathrm{gk}}{p+1}\tanh \left(\frac{\mathrm{kh}}{p+1}\right)}+k\·u---\left(2\right)$

Wherein, p=0,1 o'clock respectively corresponding 0,1 rank dispersion relation.

Mainly there is following deficiency in this method: 1. the extra large subsurface runoff that measures of X-band radar should be an average flow information in the zone of analyzing, and the stream that the least square method match of energy weighting obtains, and tends to the stream of the big ripple of energy; 2. rely on the method applicability of experience selected threshold bad, the empirical value of different waters, sea situation is different, needs to obtain through a large amount of experiments, and versatility is bad.

Summary of the invention

The object of the present invention is to provide a kind of really to the flow velocity MEASURING QUASI, versatility is good, can improve the iteration ocean current inversion method based on the adaptive threshold selecting technology that X-band radar is measured the precision of extra large subsurface runoff.

The objective of the invention is to realize like this:

(1) collection of radar image with D be converted to 32 original radar images,

The X-band navar is operated under the short pulse, and the capillary wave generation Bragg diffraction on the electromagnetic wave of emission and sea, back scattering are received machine and receive; Vision signal converts digital signal into through A/D, and the vision signal between two stem signals is formed a width of cloth radar image, selects the continuous radar image g (x of 32 width of cloth; Y, t), wherein; X, y are respectively two components in radar image space, and t is the time component;

(2) obtain the image under the Cartesian coordinates and it is obtained image spectrum as Fourier transform;

Rectangle frame need is analyzed in selection in radar image, use the closest approach interpolation realize in the rectangle frame radar image by polar coordinates to the Cartesian coordinates conversion table be shown g1 (x, y, t); (x, y t) become the wave number frequency field to radar image sequence g1 in the rectangle frame after discrete 3 dimension Fourier transforms;

$\mathrm{\γ}(k,w)={\∫}_0^{L_X}{\∫}_0^{L_Y}{\∫}_0^{T}g(x,y,t)\exp \left[i(k_{x}x+k_{y}y-\mathrm{wt})\right]\mathrm{dxdydt}$

K=(k wherein _x, k _y), L _x, L _y, T is respectively the 180 degree fuzzy problems of considering actual conditions behind length and width and the seasonal effect in time series total length of time three-dimensional Fourier transform of rectangular area and eliminating spectrum, only keeps the part of w＞0 so the energy spectrum of 3 d image and is:

$F^{\left(3\right)}(k,w)=\frac{1}{L_X{L}_{Y}T}{\mathrm{\γ}}^2(k,w);$

(3) influence of 0 order ripple is only considered in initial valuation, and inverting obtains rough stream;

In the wave number frequency field, choose the bigger point of dozens of energy value, according to the dispersion relation between Ocean-wave Signal wave number and frequency, use the average weighted least square method, simulate stream information, basic mathematical model is following:

In order to draw subsurface runoff u, with formula

${\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&omega;}}_i-S\left(k_i\right))}^2$

To u _x, u _yAsk the single order partial derivative, and be zero:

$\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2/\&PartialD;u_x=0$

$\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2/\&PartialD;u_y=0$

Following formula is write as the form of matrix:

$\left(\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{yx}}& D_{\mathrm{yy}}\end{array}\right)\left(\begin{array}{c}{u}_x\\ u_y\end{array}\right)=\left(\begin{array}{c}{b}_x\\ b_y\end{array}\right)$

Be abbreviated as:

Du＝b

The element of matrix D does

$D_{\mathrm{xx}}=\underset{i=1}{\overset{\mathrm{<figure-callout\; id="2"\; label="n\; k\; ix"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">n</figure-callout>}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}^{}{}_2^{},D_{\mathrm{xy}}=D_{\mathrm{yx}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}{k}_{\mathrm{iy}}$

$D_{\mathrm{yy}}=\underset{i=1}{\overset{\mathrm{<figure-callout\; id="2"\; label="n\; k\; iy"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">n</figure-callout>}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{iy}}^{}{}_2^{}$

The element of vector b does

$b_x=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}{w}_{\mathrm{iD}},$ $b_y=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{iy}}{w}_{\mathrm{iD}}$

Wherein, w _ID=w _i-ζ (k _i) in order to calculate the u that effluents, it is D that D is necessary for nonsingular matrix (det (D) ≠ 0) ^-1Exist, at this moment

u＝D ^-1·b；

(4) based on the adaptive threshold selecting technology, consider the influence of 0,1 order ripple, use average weighted least square method iterative fitting stream information;

1. use the adaptive threshold selecting technology, obtain threshold value C _It, greater than C _ItEnergy just comprised 0 order and 1 order Wave energy, number is N ₁, the stream that utilizes initial valuation to obtain is by formula

Calculate 0 order and 1 order wave frequency;

2. judge this N that surveys ₁Individual energy point meets 0 order or 1 order wave dispersion relation, if | w _i-w (k _i) |＜| w _i-w ₁(k _i) |, then this frequency meets 0 order ripple; If | w _i-w (k _i) |＞| w _i-w ₁(k _i) |, then this frequency meets 1 order ripple, with judge good data based different minimal value function suc as formula

Use average least square method, match obtains new subsurface runoff;

3. use new subsurface runoff and make up BPF., obtain new wave signal to noise ratio (S/N ratio), use the adaptive threshold selecting technology, obtain new threshold value;

The subsurface runoff that 4. will newly obtain; Substitution formula

obtains 0 new order and 1 order wave frequency; Repeat above-mentioned steps, will obtain constantly accurate subsurface runoff.

Said adaptive threshold selecting technology is:

1. in the short period of time, the variation of wave on time and space is less, and the signal to noise ratio (S/N ratio) of a preceding sequence chart is the signal to noise ratio snr of this sequence chart in approximate thinking;

2. the image spectrum F that obtains by 3 dimension Fourier transforms ⁽³⁾(k w) obtains the Wave energy P that contains in the image spectrum with signal to noise ratio snr

$P=\frac{\mathrm{SNR}+1}{\mathrm{SNR}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_w}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},w_k)\&CenterDot;{\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y\mathrm{\&Delta;w}$

3. energy points all in the image spectrum is sorted according to from big to small order, the energy of ordering back i sampled point is P _i,, therefore can make because Wave energy point is generally big than the ground unrest energy point

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i$

Obtain the energy point number m that is used for LSM by following formula, the energy value of m sampled point is P _m, P then _mBe iteration energy threshold C _It

Main technical points of the present invention is:

(1) initial valuation obtains rough flow velocity (size, direction).Suppose that the wave space is even in the rectangle frame, the time is stable, and then the wave model meets Gaussian distribution.When not having under the flow rate conditions, the dispersion relation of extra large capillary wave and gravity wave does

$\mathrm{\&zeta;}\left(k\right)\&ap;\sqrt{g\left|k\right|\tanh \left(\right|k\left|h\right)}---\left(3\right)$

Wherein, ζ (k) is a wave frequency, and k is a wave-number vector, and h is the depth of water, and g is an acceleration of gravity.If a subsurface runoff u with respect to radar is arranged, just introduced the Doppler shift item in the frequency, dispersion relation becomes:

S(k)＝ζ+k·u＝ζ+|k||u|cosθ(4)

Wherein, S (k) is theoretical wave frequency, and ku is a Doppler shift, u=(u _x, u _y) be the relative velocity between wave field and radar antenna platform; Comprise Texas tower movement velocity (such as ship's speed) and subsurface runoff speed vector, θ is an angle between wave number k and the u, when θ=90; The Doppler shift item is 0; Have only stream u on the direction of wave number k, could influence Doppler shift when important, when Texas tower was static, u referred to extra large subsurface runoff.

Can obtain rough stream information through the average weighted least square method.The minimal value function definition of average weighted least square method is suc as formula shown in (5):

${\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&omega;}}_i-S\left(k_i\right))}^2---\left(5\right)$

Wherein, ω _iBe the wave frequency of radar observation, n is the energy point number that is used for least square method.

(2) the automatic selecting technology of threshold value obtains threshold value C _It

Big when wind speed, when wave height is higher, the echo of radar image is stronger, the signal to noise ratio (S/N ratio) of Ocean-wave Signal is bigger, the quality of the data of this moment is better.Otherwise the Ocean-wave Signal signal to noise ratio (S/N ratio) is less, and the quality of data is poor.The quality of the signal to noise ratio (S/N ratio) judgment data quality through Ocean-wave Signal, selected threshold C according to this _It

At first provide the definition of wave image energy signal to noise ratio (S/N ratio):

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

$\mathrm{SIG}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}}){\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y---\left(7\right)$

$\mathrm{BGN}=\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_{\mathrm{\&omega;}}}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_k){\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y\mathrm{\&Delta;\&omega;}---\left(8\right)$

$-\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}{F}^{\left(2\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}}){\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y$

$F^{\left(2\right)}\left(k\right)=2{\&Integral;}_{\mathrm{\&omega;}>0}{F}^{\left(3\right)}(k_x,k_y,\mathrm{\&omega;})\mathrm{\&delta;}(\mathrm{\&omega;}-{\mathrm{\&omega;}}_0)\mathrm{d\&omega;}---\left(9\right)$

Wherein, SIG is the wave spectrum energy, and BGN is a ground unrest, F ⁽³⁾Be the 3 d image spectrum, F ⁽²⁾(k) be 2 dimension image spectrums, N _Kx, N _Ky, N _ωBe the scope of spectrum, Δ k _x, Δ k _yBe wavenumber resolution, Δ ω is a frequency resolution, δ (ω-ω ₀) be BPF..

The basic realization of threshold value automatically selecting method is in conjunction with accompanying drawing 7 explanations.

1. in the short period of time, it promptly is stable in time that the signal to noise ratio (S/N ratio) of wave changes less, and the signal to noise ratio (S/N ratio) that therefore obtains before the hypothesis is the signal to noise ratio snr of this sequence chart.

2. the image spectrum F that obtains by 3 dimension Fourier transforms ⁽³⁾(k w) can obtain the Wave energy P that contains in the image spectrum with signal to noise ratio snr.

$P=\frac{\mathrm{SNR}+1}{\mathrm{SNR}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_w}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},w_k)\&CenterDot;{\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y\mathrm{\&Delta;w}---\left(10\right)$

3. energy points all in the image spectrum is sorted according to from big to small order, the energy of ordering back i sampled point is P _i,, therefore can make because Wave energy point is generally big than the ground unrest energy point

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i---\left(11\right)$

Can obtain the energy point number m that is used for LSM by following formula, the energy value of m sampled point is P _m, P then _mBe iteration energy threshold C _It

(3) based on the iteration valuation of the automatic selecting technology of threshold value.

At this, we only consider the situation of 0,1 order harmonic wave.

1. the application self-adapting selection of threshold is technological, selected threshold C _iT can think greater than C _ItEnergy partly comprised 0 order and 1 order Wave energy, number is N ₁The stream that utilizes initial valuation to obtain calculates 0 order and 1 order wave frequency by formula (12) respectively.

$w_p=(p+1)\sqrt{\frac{\mathrm{gk}}{p+1}\tanh \left(\frac{\mathrm{kh}}{p+1}\right)}+k\&CenterDot;u---\left(12\right)$

The dispersion relation of p=0,1 respectively corresponding 0,1 order ripple.

2. judge this N that actual measurement obtains ₁Individual energy point belongs to 0 order or 1 order ripple, if | w _i-w (k _i) |＜| w _i-w ₁(k _i) |, then this frequency meets 0 order ripple; If | w _i-w (k _i) |＞| w _i-w ₁(k _i) |, then this frequency meets 1 order ripple.With judging good data based different minimal value function, use least square method, can obtain new subsurface runoff.

3. use new subsurface runoff and make up BPF., obtain new Ocean-wave Signal signal to noise ratio (S/N ratio), application self-adapting selection of threshold technology obtains new threshold value.

The subsurface runoff that 4. will newly obtain in the substitution formula (12), obtains 0 new order and 1 order wave dispersion relation, repeats above-mentioned steps, will obtain constantly accurate subsurface runoff.

The stream inverting new method advantage of choosing based on automatic threshold compared with prior art of invention:

(1) different wind speed and zone down in the radar sequence chart Ocean-wave Signal be different, mainly be to embody through signal to noise ratio (S/N ratio).The energy of wave in the image sequence can be calculated by the gross energy of signal to noise ratio (S/N ratio) and image, threshold value can be drawn thus.The present invention proposes the automatic selecting technology of threshold value based on above thinking, it can discern the Ocean-wave Signal quality automatically, and through iterative technique, from the image spectrum, isolates Wave energy accurately, has remedied the deficiency of empirical value.

(2) can know from Fig. 1,2, table 1, when flow velocity is big (＞0.3m/s), the more former methods and results of stream inversion result of the inventive method has two aspect advantages: the trend that 1. can better reflect stream; 2. has better stability of littler sum of errors.

(3) flow velocity hour (＜0.2m/s), signal a little less than, from weak signal, extracting stream information is a difficult problem.The more former method of the inversion result of method of the present invention has very big improvement shown in Fig. 3,4, table 1.

Description of drawings

When Fig. 1 flow velocity was big, method, former method and spot sensor flow rate result were relatively among the present invention.

When Fig. 2 flow velocity was big, method, former method and spot sensor flowed to the result relatively among the present invention.

Fig. 3 flow velocity hour, method, former method and the comparison of spot sensor flow rate result among the present invention.

Fig. 4 flow velocity hour, method, former method and spot sensor flow to result's comparison among the present invention.

Fig. 5 technology implementation scheme of the present invention process flow diagram.

Fig. 6 is based on the iteration stream inverting process flow diagram of adaptive threshold selecting technology.

Fig. 7 adaptive threshold selecting technology implementing procedure figure.

When Fig. 8 table 1 flow velocity is big, each methods and results statistical study.

Fig. 9 table 2 flow velocity hour, each methods and results statistical study.

Embodiment

Below in conjunction with Fig. 5,6,7, the present invention is done introduction more in detail:

(5) collection of radar image with D be converted to 32 original radar images.

The X-band navar is operated under the short pulse, and the capillary wave generation Bragg diffraction on the electromagnetic wave of emission and sea, back scattering are received machine and receive, and vision signal converts digital signal into through A/D; Vision signal between two stem signals is formed a width of cloth radar image, in the present invention, needs to use the continuous radar image g (x of 32 width of cloth; Y, t), wherein; X, y are respectively two components in radar image space, and t is the time component.

(6) obtain the image under the Cartesian coordinates and it is obtained image spectrum as Fourier transform.

Rectangle frame need is analyzed in selection in radar image, use the closest approach interpolation realize in the rectangle frame radar image by polar coordinates to the Cartesian coordinates conversion table be shown g1 (x, y, t); (x, y t) become the wave number frequency field to radar image sequence g1 in the rectangle frame after discrete 3 dimension Fourier transforms.

$\mathrm{\&gamma;}(k,w)={\&Integral;}_0^{L_X}{\&Integral;}_0^{L_Y}{\&Integral;}_0^{T}g(x,y,t)\exp \left[i(k_{x}x+k_{y}y-\mathrm{wt})\right]\mathrm{dxdydt}---\left(13\right)$

K=(k wherein _x, k _y), L _x, L _y, T is respectively the 180 degree fuzzy problems of considering actual conditions behind length and width and the seasonal effect in time series total length of time three-dimensional Fourier transform of rectangular area and eliminating spectrum, and we have only kept the part of w＞0 so the energy spectrum of 3 d image is:

$F^{\left(3\right)}(k,w)=\frac{1}{L_X{L}_{Y}T}{\mathrm{\&gamma;}}^2(k,w)---\left(14\right)$

(7) influence of 0 order ripple is only considered in initial valuation, and inverting obtains rough stream.

In the wave number frequency field, choose the bigger point of dozens of energy value, according to the dispersion relation between Ocean-wave Signal wave number and frequency, use the average weighted least square method, simulate stream information, basic mathematical model is following:

In order to draw subsurface runoff u, with formula (15)

${\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&omega;}}_i-S\left(k_i\right))}^2---\left(15\right)$

To u _x, u _yAsk the single order partial derivative, and be zero:

$\begin{array}{c}\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2/\&PartialD;u_x=0\\ \&PartialD;{\mathrm{<figure-callout\; id="2"\; label="Q"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">Q</figure-callout>}}^2/\&PartialD;u_y=0\end{array}---\left(16\right)$

Following formula is write as the form of matrix:

$\left(\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{yx}}& D_{\mathrm{yy}}\end{array}\right)\left(\begin{array}{c}{u}_x\\ u_y\end{array}\right)=\left(\begin{array}{c}{b}_x\\ b_y\end{array}\right)---\left(17\right)$

Be abbreviated as:

Du＝b (18)

The element of matrix D does

$D_{\mathrm{xx}}=\underset{i=1}{\overset{\mathrm{<figure-callout\; id="2"\; label="n\; k\; ix"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">n</figure-callout>}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}^{}{}_2^{},D_{\mathrm{xy}}=D_{\mathrm{yx}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}{k}_{\mathrm{iy}}---\left(19\right)$

$D_{\mathrm{yy}}=\underset{i=1}{\overset{\mathrm{<figure-callout\; id="2"\; label="n\; k\; iy"\; filenames="HDA0000042754790000013.png"\; state="\{\{state\}\}">n</figure-callout>}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{iy}}^{}{}_2^{}$

The element of vector b does

$b_x=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ix}}{w}_{\mathrm{iD}},b_y=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_{\mathrm{iy}}{w}_{\mathrm{iD}}---\left(20\right)$

Wherein, w _ID=w _i-ζ (k _i) in order to calculate the u that effluents, it is D that D is necessary for nonsingular matrix (det (D) ≠ 0) ^-1Exist, at this moment

u＝D ^-1·b (21)

(8) based on the adaptive threshold selecting technology, consider the influence of 0,1 order ripple, use average weighted least square method iterative fitting stream information.Describe in conjunction with accompanying drawing 6.The practical implementation of adaptive threshold selecting technology sees (5) part for details.

The flow velocity (size, direction) that obtains through the initial valuation of (3) step is inaccurate, needs to consider the influence of single order subwave, uses the inversion accuracy of the method raising stream of iteration, and its implementation combines accompanying drawing 6 to do following explanation:

1. use the adaptive threshold selecting technology, obtain threshold value C _It, greater than C _ItEnergy just comprised 0 order and 1 order Wave energy, number is N ₁The stream that utilizes initial valuation to obtain calculates 0 order and 1 order wave frequency by formula (12).

2. judge this N that surveys ₁Individual energy point meets 0 order or 1 order wave dispersion relation, if | w _i-w (K _i) |＜| w _i-w ₁(k _i) |, then this frequency meets 0 order ripple; If | w _i-w (k _i) |＞| w _i-w ₁(k _i) |, then this frequency meets 1 order ripple.The data based different minimal value function that judgement is good is used average least square method suc as formula (15), can match obtain new subsurface runoff.

3. use new subsurface runoff and make up BPF., obtain new wave signal to noise ratio (S/N ratio), use the adaptive threshold selecting technology, obtain new threshold value.

The subsurface runoff that 4. will newly obtain, substitution formula (12) obtains 0 new order and 1 order wave frequency, repeats above-mentioned steps, will obtain constantly accurate subsurface runoff.

(9) in last step, use the adaptive threshold selecting technology, combine accompanying drawing 7 to describe at present,

1. in the short period of time, the variation of wave on time and space is less, and what therefore can be similar to thinks that the signal to noise ratio (S/N ratio) of a preceding sequence chart is the signal to noise ratio snr of this sequence chart;

2. the image spectrum F that obtains by 3 dimension Fourier transforms ⁽³⁾(k w) can obtain the Wave energy P that contains in the image spectrum with signal to noise ratio snr

$P=\frac{\mathrm{SNR}+1}{\mathrm{SNR}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_w}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},w_k)\&CenterDot;{\mathrm{\&Delta;k}}_x{\mathrm{\&Delta;k}}_y\mathrm{\&Delta;w}---\left(22\right)$

3. energy points all in the image spectrum is sorted according to from big to small order, the energy of ordering back i sampled point is P _i,, therefore can make because Wave energy point is generally big than the ground unrest energy point

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i---\left(23\right)$

Can obtain the energy point number m that is used for LSM by following formula, the energy value of m sampled point is P _m, P then _mBe iteration energy threshold C _It

The stream inversion method that the present invention is proposed is applied in the measured data of flow velocity when big (time period be the morning on the 06th October in 2009 10:36:00 to 16:40:00 in afternoon) result like Fig. 1, shown in 2; Statistic analysis result is as shown in table 1, can find out that the more former method of the inventive method result has higher precision and stability.Fig. 3,4 has provided hour result of the inventive method inverting stream (time period be the morning on the 06th October in 2009 08:16:00 to 09:20:00), statistics analysis in table 2, but the more former methods and results precision of knowledge capital inventive method result improves a lot.

Claims (1)

1. iteration ocean current inversion method based on the adaptive threshold selecting technology is characterized in that:

(1) collection of radar image with D be converted to 32 original radar images;

(2) obtain the image under the Cartesian coordinates and it is obtained image spectrum as Fourier transform;

(3) influence of 0 order ripple is only considered in initial valuation, and inverting obtains rough stream;

(4) based on the adaptive threshold selecting technology, consider the influence of 0,1 order ripple, use average weighted least square method iterative fitting stream information;

Said adaptive threshold selecting technology is:

1. in the short period of time, the variation of wave on time and space is less, and the signal to noise ratio (S/N ratio) of a preceding sequence chart is the signal to noise ratio snr of this sequence chart in approximate thinking;

2. the image spectrum F that obtains by 3 dimension Fourier transforms ⁽³⁾(k ω) obtains the Wave energy P that contains in the image spectrum with signal to noise ratio snr

$P=\frac{\mathrm{SNR}+1}{\mathrm{SNR}}\underset{i}{\overset{N_{\mathrm{kx}}}{\mathrm{\&Sigma;}}}\underset{j}{\overset{N_{\mathrm{ky}}}{\mathrm{\&Sigma;}}}\underset{k}{\overset{N_w}{\mathrm{\&Sigma;}}}{F}^{\left(3\right)}(k_{\mathrm{ix}},k_{\mathrm{jy}},{\mathrm{\&omega;}}_k)\&CenterDot;\mathrm{\&Delta;}{k}_x\mathrm{\&Delta;}{k}_y\mathrm{\&Delta;\&omega;}$

Wherein: N _Kx, N _Ky, N _ωScope for spectrum; K=(k _x, k _y) be wave number, k _x, k _yBe respectively the x and the y component of wave number, Δ k _x, Δ k _yBe wavenumber resolution; ω is the frequency of wave, and Δ ω is a frequency resolution;

3. energy points all in the image spectrum is sorted according to from big to small order, the energy of ordering back i sampled point is Pi, because Wave energy point is generally big than the ground unrest energy point, therefore makes

$P=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{P}_i$

Obtain the energy point number m that is used for least square method by following formula, the energy value of m sampled point is P _m, P then _mBe iteration energy threshold C _It

Images