CN104698462A - Sea surface wind field fusion method for SAR (Synthetic Aperture Radar) based on variation - Google Patents

Abstract

The invention discloses a sea surface wind field fusion method for SAR (Synthetic Aperture Radar) based on variation. The method comprises the following steps: firstly, selecting an area with relatively obvious wind stripe characteristics in an SAR image, and analyzing to obtain the wind direction of ocean surface with high precision by adopting a WT (Wavelet Transform) method, so as to obtain a sea surface wind speed by inverting a geophysical model function; secondly, inverting a sea surface wind field in an SAR observation area by using a WRF (Weather Research Forecast) mode wind direction; and finally, optimizing and adjusting a background filed through an observation field by using a variation method by taking the sea surface wind field as an observation field based on wind stripe inversion and taking the sea surface wind field as a background wind field based on wind stripe inversion, so as to obtain a fused sea surface wind field. According to the sea surface wind field fusion method for the SAR based on the variation provided by the invention, the problem that sea surface wind direction inversion depends on the wind stripes can be solved, the inverted sea surface wind direction can be optimized and adjusted by using a variational method, and the integral precision of the inverted regional sea surface wind field can be improved.

Description

Based on the synthetic-aperture radar Ocean Wind-field fusion method of variation

Technical field

The invention belongs to remote sensing technology field, particularly a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation.

Background technology

SAR (synthetic-aperture radar, synthetic aperture radar, SAR) has the feature of round-the-clock, round-the-clock and high spatial resolution, can obtain the Ocean Wind-field of high-spatial and temporal resolution.But the single incident angle mechanism of SAR makes it first will obtain wind direction when WIND FIELDS, then calculation of wind speed.The obtaining means of wind direction has two classes, one class is based on SAR image wind striped inverting wind direction, the wind direction precision that the method obtains is high, but statistics shows not have in the image of nearly 60% obvious wind striped can for inverting, so be difficult to utilize wind striped to carry out the Wind-field Retrieval of large area region; Another kind of be adopt scatterometer or Numerical Prediction Models wind direction as wind outside to, and realize Wind-field Retrieval in conjunction with SAR image information Wind Speed Inversion, the method is reliable and stable, be easy to carry out the work of regional SAR ocean surface wind retrieving, but the wind direction data spatial resolution that it provides and precision all lower, the changing condition of sea surface wind vector in small scale cannot be obtained, also have impact on the precision of wind speed retrieval.

In air and Marine Sciences, fusion method is a lot, be applied to wind field fusion method to mainly contain: objective analysis method, successive correction analysis, B-spline method, the variational method, variation are in conjunction with regularization method, Loess low pass filtering method etc., wherein variational method has become the mainstream technology of current Data Assimilation, is the dominant direction of Meteorological Data Assimilation.Because wind field merges the complex nature of the problem, the theoretical research of merging about Ocean Wind-field is relative with operational use less, and SAR is subject to the restriction of wind striped in ocean surface wind retrieving, the detection datas such as fusion data many employings microwave altimeter, scatterometer and radiometer, the high resolving power feature of SAR data is not fully utilized.SAR image and numerical model data work in coordination with the problem that inverting Ocean Wind-field solves wind stripe loss, high-resolution wind striped WIND FIELDS is utilized to be optimized adjustment, for the operational use realizing SAR high precision ocean surface wind retrieving provides important support to pattern wind direction as the WIND FIELDS of initial wind direction.At present, the method merged for SAR Ocean Wind-field have not been reported, and therefore, needs to find a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation.

Summary of the invention

The object of the present invention is to provide a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation, utilize variational method to be optimized adjustment with the overall precision improving regional ocean surface wind retrieving to the Ocean Wind-field after SAR inverting.

The technical solution realizing the object of the invention is: a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation, said method comprising the steps of:

Based on a synthetic-aperture radar Ocean Wind-field fusion method for variation, comprise the following steps:

Step 1, carries out experience differentiation to SAR image after calibration, SAR image is divided into the obvious region of feature and feature Fuzzy region; (class is the obvious region of wind streak feature, and a class is that above-mentioned classification foundation experience differentiates not containing wind streak feature or the fuzzyyer region of wind streak feature; )

Step 2, the obvious region of the feature described in selecting step 1, utilizes D continue wavelet transform inversion method wind direction, and in conjunction with physical geography module function calculation of wind speed, obtain the obvious region of described feature based on wind striped inverting Ocean Wind-field;

Step 3, choose NCEP (the Environmental forecasting centre matched with SAR observed image space-time, NationalCenters for Environmental Prediction) data are input to WRF pattern and carry out analog computation WRF (weather forecast pattern, Weather Research and Forecasting Mode) wind direction, then WRF wind direction and corresponding SAR data are input to physical geography module function calculation of wind speed jointly, obtain the Ocean Wind-field of whole SAR observation area based on the inverting of WRF wind direction;

Step 4, with in step 2 obtain based on wind striped inverting Ocean Wind-field for observation field, with the Ocean Wind-field based on the inverting of WRF wind direction obtained in step 3 for ambient field, utilize variational method to carry out fusion adjustment by observation field to ambient field, obtain the Ocean Wind-field after merging.

Feature Fuzzy region of the present invention does not need process, because the inventive method is exactly the wind field adjusting the SAR image entirety utilizing WRF pattern simulation to go out with the wind field that the obvious region of feature in SAR image and wind fringe area are finally inversed by.

More preferably, step 2 specifically comprises the following steps,

201, D continue wavelet transform is carried out to the SAR image in the obvious region of wind streak feature, obtains the Wavelet Energy Spectrum image under different scale,

Described D continue wavelet transform method is two-dimentional Mexican-hat wavelet transformation,

Described two-dimentional Mexican-hat wavelet transformation formula is formula (1):

${\mathrm{\Ψ}}_H\left(a\right)=(a\·a)e^{\left(\frac{1}{2(a\·a)}\right)}---\left(1\right)$

In formula: a represents the variable of two-dimensional space-frequency field; Represent inner product of vectors, Ψ _ha () represents wavelet transform function;

202, two-dimentional fast Fourier (FFT) conversion is carried out to Wavelet Energy Spectrum image, calculates the wavenumber spectrum of SAR image apoplexy striped:

Two-dimensional fast fourier transform formula is formula (2):

$Y_{l,m}=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{b=1}{\overset{N}{\mathrm{\Σ}}}{X}_{j,b}{e}^{-2\mathrm{\πi}(\mathrm{jl}+\mathrm{bm})/N}---\left(2\right)$

Wherein, Y is the wavenumber spectrum of SAR image apoplexy striped, and X is image intensity value, and l, m represent that SAR image pixel is horizontal, lengthwise position; L, m=1,2 ..., N, j, b represent SAR image pixel location variable.

203, the line of two-dimentional wavenumber spectrum peak value is done vertical line, after wind direction deblurring is carried out to described vertical line, obtains wind direction of ocean surface.

More preferably, geophysical model is C-band geophysical model CMOD (C-band models), and described C-band physical geography module function is:

σ _C(θ,φ,u)＝10 ^A(u,θ)(1+B(u,θ)cos(φ)+C(u,θ)cos2φ) (3)

In formula, σ _crepresent VV polarimetric radar backscattering coefficient, φ represent wind direction relative to radar look-directions apparent wind to, θ represents radar antenna incident angle, and u represents ocean surface wind speed; A (u, θ), B (u, θ) and C (u, θ) represent by ten meters, sea At The Height wind speed, apparent wind to, polarization mode, radar frequency and incident angle determine coefficient.

More preferably, utilize variational method to carry out fusion adjustment by observation field to ambient field in step 4, obtain the Ocean Wind-field after merging, specifically comprise the following steps:

401, adopt the method for two-dimentional variation (2D-Var) to carry out fusion calculation, the objective function defined by formula (4) is minimum to be retrained:

J＝J _q+J _m(4)

The objective function J of its apoplexy striped wind field _qwith the objective function J of WRF pattern _mbe defined as formula (5) and formula (6):

$J_q=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^{-1}(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)---\left(5\right)$

$J_m=\frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)---\left(6\right)$

Wherein, U and V merges the zonal wind and meridional wind (on region mode net point) that obtain, q and m represents wind striped wind field and region WRF pattern, i.e. U _q, V _qrepresent zonal wind and the meridional wind of wind striped wind field respectively, U _m, V _mrepresent the zonal wind under the WRF pattern of region and meridional wind respectively, described Q _u, Q _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of wind striped wind field respectively, M _u, M _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of WRF pattern respectively, q and m represents wind striped wind field and region WRF pattern, H _qfor space interpolation operator, analysis field is projected to observation space, () ^trepresent matrix transpose, () ^-1represent inverse of a matrix;

402, utilize Kriging interpolation method to construct space interpolation operator H _q, assuming that V _ifor the wind speed and direction value of i-th point in outside wind field, i=1,2 ..., n (n is the sum of observation station in outside wind field), for outside wind field is interpolated into the wind speed and direction value of a kth point (observation station) on wind striped wind field region, k be a kth point on wind striped wind field region (k=1,2 ..., K, K are observation station sum on wind striped wind field region), be the contribution weights of i-th data to a kth impact point, order

$\stackrel{\‾}{Y}={({\stackrel{\‾}{V}}_1,{\stackrel{\‾}{V}}_2,...,{\stackrel{\‾}{V}}_k)}^T,$

X＝(V ₁,V ₂,…V _N) ^T,

Then:

$\stackrel{\‾}{Y}=H_{q}X$

Observation Operators H _qbe expressed as formula (7),

$H_q={\left(h_i^k\right)}_{K\×N}---\left(7\right)$

Contribution weight obtained by Kriging system of equations:

Wherein, Lagrange multiplier when μ is minimization process, c (x _i, x _j) be net point x in outside wind field _iwith x _jbetween covariance function, x ^kfor the interpolation point on wind striped wind field region;

403 utilize the remaining difference method structural setting error co-variance matrix of observation, and background error covariance matrix is formula (9) and (10):

＜U _m-U _q＞ ²＝dQ _u+dM _u(9)

＜V _m-V _q＞ ²＝dQ _v+dM _v(10)

Wherein < > represents time average, dQ _u, dM _urepresent the zonal wind error variance of wind striped wind field and the zonal wind error variance of WRF pattern respectively, dQ _v, dM _vrepresent the meridional wind error variance of wind striped wind field and the meridional wind error variance of WRF pattern respectively, dQ _u, dQ _v, dM _u, dM _vrespectively with Q _u, Q _v, M _u, M _vdiagonal element consistent;

404, based target function is minimum, adopts variation analysis method wind striped wind field and WRF pattern wind field to be merged to the zonal wind u and meridional wind v that obtain optimized analysis wind field.

More preferably, step 404 specifically comprises the following steps:

The zonal wind u of optimized analysis wind field and meridional wind v is the zonal wind U and meridional wind V that after meeting fusion, the fusion of objective function least commitment condition (value of formula (4) is minimum) obtains, and namely meets formula (11):

$\begin{array}{c}J=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^{-1}(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)+\\ \frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)=\min !\end{array}---\left(11\right)$

Based on the variational method, can obtain:

$\begin{array}{c}\mathrm{\&delta;J}=H_q^T{Q}_u^{-1}(H_{q}U-U_q)\mathrm{\&delta;U}+H_q^T{Q_v}^{-1}(H_{q}V-V_q)\mathrm{\&delta;V}+\\ {M_u}^{-1}(U-U_m)\mathrm{\&delta;U}+{M_v}^{-1}(V-V_m)\mathrm{\&delta;V}=0!\end{array}---\left(12\right)$

Utilize the arbitrariness of δ U, δ V, U, V meet following condition:

$H_q^T{Q}_u^{-1}(U_q-H_{q}U)+{M_u}^{-1}(U-U_m)=0$

$H_q^T{Q_v}^{-1}(V_q-H_{q}V)+{M_v}^{-1}(V_m-V)=0---\left(13\right)$

Thus show that the zonal wind u of optimized analysis wind field and meridional wind v is:

$u=U={[H_q^T{Q}_u^{-1}{H}_q+{M_u}^{-1}]}^{-1}[H_q^T{Q}_u^{-1}{U}_q+{M_u}^{-1}{U}_m]$

$v=V={[H_q^T{Q_v}^{-1}{H}_q+{M_v}^{-1}]}^{-1}[H_q^T{{Q_v}^{-1}V}_q+{M_v}^{-1}{V}_m].$

Compared with prior art, the present invention includes following beneficial effect:

The invention provides a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation, for the feature of the high-spatial and temporal resolution of SAR data, SAR image and numerical model data is adopted to work in coordination with inverting Ocean Wind-field, solve the problem of wind stripe loss, and utilize variational method to be optimized adjustment to the Ocean Wind-field after inverting, improve the overall precision of regional ocean surface wind retrieving, for the operational use realizing SAR high precision ocean surface wind retrieving provides important support.

Accompanying drawing explanation

Fig. 1 is the synthetic-aperture radar Ocean Wind-field fusion method process flow diagram based on variation of the present invention;

Fig. 2 is embodiments of the invention data used;

Fig. 3 (a), 3 (b) are SAR data Wind-field Retrieval images, and wherein 4 (a) is based on wind striped WIND FIELDS image, and 4 (b) is based on WRF pattern wind direction WIND FIELDS image.

Embodiment

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

Below with reference to accompanying drawing of the present invention; clear, complete description and discussion are carried out to the technical scheme in the embodiment of the present invention; obviously; as described herein is only a part of example of the present invention; it is not whole examples; based on the embodiment in the present invention, the every other embodiment that those of ordinary skill in the art obtain under the prerequisite not making creative work, all belongs to protection scope of the present invention.

For the ease of the understanding to the embodiment of the present invention, be further explained for specific embodiment below in conjunction with accompanying drawing, and each embodiment does not form the restriction to the embodiment of the present invention.

A kind of SAR Ocean Wind-field fusion method based on variation of the present invention.Composition graphs 1, comprises the following steps:

1, Feature Selection: experience differentiation is carried out to SAR image after calibration, SAR image is divided into the obvious region of feature and feature Fuzzy region; (class is the obvious region of wind streak feature, and a class is not containing wind streak feature or the fuzzyyer region of wind streak feature; )

The obvious region of wind streak feature is chosen in SAR image after calibration.

2, based on the ocean surface wind retrieving of SAR image wind striped: the obvious region of the feature described in selecting step 1, utilize D continue wavelet transform inversion method wind direction, and in conjunction with physical geography module function calculation of wind speed, obtain the obvious region of described feature based on wind striped inverting Ocean Wind-field;

Step 2 specifically comprises the following steps,

201, D continue wavelet transform is carried out to the SAR image in the obvious region of wind streak feature, obtains the Wavelet Energy Spectrum image under different scale,

Described D continue wavelet transform method is two-dimentional Mexican-hat wavelet transformation,

Described two-dimentional Mexican-hat wavelet transformation formula is formula (1):

${\mathrm{\&Psi;}}_H\left(a\right)=(a\&CenterDot;a)e^{\left(\frac{1}{2(a\&CenterDot;a)}\right)}---\left(1\right)$

In formula: a represents the variable of two-dimensional space-frequency field; Represent inner product of vectors, Ψ _ha () represents wavelet transform function;

202, two-dimentional fast Fourier (FFT) conversion is carried out to Wavelet Energy Spectrum image, calculates the wavenumber spectrum of SAR image apoplexy striped:

Two-dimensional fast fourier transform formula is formula (2):

$Y_{l,m}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{b=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_{j,b}{e}^{-2\mathrm{\&pi;i}(\mathrm{jl}+\mathrm{bm})/N}---\left(2\right)$

Wherein, Y: be the wavenumber spectrum of SAR image apoplexy striped, X is image intensity value, l, m represent that SAR image pixel is horizontal, lengthwise position; J, b represent SAR image pixel location variable.

203, the line of two-dimentional wavenumber spectrum peak value is done vertical line, after wind direction deblurring is carried out to described vertical line, obtains wind direction of ocean surface.

More preferably, geophysical model is C-band geophysical model CMOD (C-band models), and described C-band physical geography module function is:

σ _C(θ,φ,u)＝10 ^A(u,θ)(1+B(u,θ)cos(φ)+C(u,θ)cos2φ) (3)

In formula, σ _crepresent VV polarimetric radar backscattering coefficient, φ represent wind direction relative to radar look-directions apparent wind to, θ represents radar antenna incident angle, and u represents ocean surface wind speed; A (u, θ), B (u, θ) and C (u, θ) represent by ten meters, sea At The Height wind speed, apparent wind to, polarization mode, radar frequency and incident angle determine coefficient.

3, the ocean surface wind retrieving based on WRF wind direction:

Choose the NCEP data matched with SAR observed image space-time to be input to WRF pattern and to carry out analog computation WRF wind direction, then WRF wind direction and corresponding SAR data are input to physical geography module function calculation of wind speed jointly, obtain the Ocean Wind-field of whole SAR observation area based on the inverting of WRF wind direction.

4, utilize variation fusion method to optimize and revise Ocean Wind-field: with in step 2 obtain based on wind striped inverting Ocean Wind-field for observation field, with the Ocean Wind-field based on the inverting of WRF wind direction obtained in step 3 for ambient field, utilize variational method to carry out fusion adjustment by observation field to ambient field, obtain the Ocean Wind-field after merging.

Utilize variational method to carry out fusion adjustment by observation field to ambient field in step 4, obtain the Ocean Wind-field after merging, specifically comprise the following steps:

4.1 adopt the method for two-dimentional variation (2D-Var) to carry out fusion calculation, and the objective function defined by formula (4) is minimum to be retrained:

J＝J _q+J _m(4)

The objective function J of its apoplexy striped wind field _qwith the objective function J of WRF pattern _mbe defined as formula (5) and formula (6):

$J_q=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^{-1}(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)---\left(5\right)$

$J_m=\frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)---\left(6\right)$

Wherein, U and V merges the zonal wind and meridional wind (on region mode net point) that obtain, q and m represents wind striped wind field and region WRF pattern, i.e. U _q, V _qrepresent zonal wind and the meridional wind of wind striped wind field respectively, U _m, V _mrepresent the zonal wind under the WRF pattern of region and meridional wind respectively, described Q _u, Q _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of wind striped wind field respectively, M _u, M _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of WRF pattern respectively, q and m represents wind striped wind field and region WRF pattern, H _qfor space interpolation operator, analysis field is projected to observation space, () ^trepresent matrix transpose, () ^-1represent inverse of a matrix;

4.2 utilize Kriging interpolation method to construct space interpolation operator H _q, assuming that V _ifor the wind speed and direction value of i-th point in outside wind field, i=1,2 ..., n, for outside wind field is interpolated into the wind speed and direction value of a kth point (observation station) on wind striped wind field region, k be a kth point on wind striped wind field region (k=1,2 ..., K, K represent observation station numbers all on wind striped wind field region), be the contribution weights of i-th data to a kth impact point, order:

$\stackrel{\&OverBar;}{Y}={({\stackrel{\&OverBar;}{V}}_1,{\stackrel{\&OverBar;}{V}}_2,...,{\stackrel{\&OverBar;}{V}}_k)}^T,$

X＝(V ₁,V ₂,…V _N) ^T,

Then

$\stackrel{\&OverBar;}{Y}=H_{q}X$

Observation Operators H _qbe expressed as formula (7),

$H_q={\left(h_i^k\right)}_{K\&times;N}---\left(7\right)$

Contribution weight obtained by Kriging system of equations:

Wherein, Lagrange multiplier when μ is minimization process, c (x _i, x _j) be net point x in outside wind field _iwith x _jbetween covariance function, x ^kfor the interpolation point on wind striped wind field region;

4.3 utilize the remaining difference method structural setting error co-variance matrix of observation, and the statistical value that the error of the present embodiment observational data is got is 1.7m/s, and background error covariance matrix is formula (9) and (10):

＜U _m-U _q＞ ²＝dQ _u+dM _u(9)

＜V _m-V _q＞ ²＝dQ _v+dM _v(10)

Wherein < > represents time average, dQ _u, dM _urepresent the zonal wind error variance of wind striped wind field and the zonal wind error variance of WRF pattern respectively, dQ _v, dM _vrepresent the meridional wind error variance of wind striped wind field and the meridional wind error variance of WRF pattern respectively, dQ _u, dQ _v, dM _u, dM _vrespectively with Q _u, Q _v, M _u, M _vdiagonal element consistent;

4.4 based target functions are minimum, adopt variation analysis method wind striped wind field and WRF pattern wind field to be merged to the zonal wind u and meridional wind v that obtain optimized analysis wind field.

The zonal wind u of optimized analysis wind field and meridional wind v is the zonal wind U and meridional wind V that after meeting fusion, the fusion of objective function least commitment condition (value of formula (4) is minimum) obtains, and namely meets formula (11):

$\begin{array}{c}J=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^{-1}(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)+\\ \frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)=\min !\end{array}---\left(11\right)$

Based on variational method, can obtain:

$\begin{array}{c}\mathrm{\&delta;J}=H_q^T{Q}_u^{-1}(H_{q}U-U_q)\mathrm{\&delta;U}+H_q^T{Q_v}^{-1}(H_{q}V-V_q)\mathrm{\&delta;V}+\\ {M_u}^{-1}(U-U_m)\mathrm{\&delta;U}+{M_v}^{-1}(V-V_m)\mathrm{\&delta;V}=0!\end{array}---\left(12\right)$

Utilize the arbitrariness of δ U, δ V, U, V meet following condition:

$H_q^T{Q}_u^{-1}(U_q-H_{q}U)+{M_u}^{-1}(U-U_m)=0$

$H_q^T{Q_v}^{-1}(V_q-H_{q}V)+{M_v}^{-1}(V_m-V)=0---\left(13\right)$

Thus show that the zonal wind u of optimized analysis wind field and meridional wind v is:

$u=U={[H_q^T{Q}_u^{-1}{H}_q+{M_u}^{-1}]}^{-1}[H_q^T{Q}_u^{-1}{U}_q+{M_u}^{-1}{U}_m]$

$v=V={[H_q^T{Q_v}^{-1}{H}_q+{M_v}^{-1}]}^{-1}[H_q^T{{Q_v}^{-1}V}_q+{M_v}^{-1}{V}_m].$ Below in conjunction with embodiment, further detailed description is done to the present invention:

As shown in Figure 1, the WSM imaging pattern VV polarization data of SAR data ENVISAT/ASAR used, when detection time is 24 days 8 May in 2011 to the flow process of the present embodiment, spatial resolution is 75m, and image size is 3 ° × 3 °, and image center geographic coordinate is 56.5N, 152.5W, as shown in Figure 2; The buoy observational data that the buoy data acquisition NDBC that compare of analysis uses provides, buoy number is 40678, and position is 56.07N, 152.57W, and the time is 7:50AM, differs 10min with SAR image detection time.Concrete steps are:

The first step, carries out experience differentiation to SAR image after calibration and chooses the obvious region of wind streak feature.As shown in Figure 2, white edge region is the obvious region of wind streak feature, and white point is buoy position.

Second step, adopts two-dimension continuous wavelet transform to carry out ocean surface wind retrieving to the obvious region of wind striped, first carries out two-dimentional continuous N exican-hat wavelet transformation to SAR image, obtain the Wavelet Energy Spectrum image under different scale, extract wind stripe information; Then Two-dimensional FFT conversion is carried out to energy spectrum image, calculate the wavenumber spectrum of SAR image apoplexy striped; Finally the line of two-dimentional wavenumber spectrum peak value is done vertical line, after wind direction deblurring is carried out to it, obtain wind direction of ocean surface, and in conjunction with physical geography module function calculation of wind speed.The Wind-field Retrieval result that Fig. 3 (a) obtains for utilizing the inverting of Mexican-hat wavelet method.

3rd step, wind outside is to the 20Km resolution WRF wind field data selecting 8h, and identical with SAR data detection time, the time mates the most.First the segmentation of extra large land is carried out to SAR image, and by 20Km resolution, gridding is carried out to image and divide, finally WRF wind direction is matched each wind field unit as initial wind direction, and utilize physical geography module function to calculate corresponding ocean surface wind speed.Fig. 3 (b) utilizes WRF wind direction as the Wind-field Retrieval result of initial wind direction for SAR image.

4th step, adopts the method for two-dimentional variation (2D-Var) to carry out Ocean Wind-field fusion:

(1) Kriging interpolation method structure space interpolation operator H is adopted _q;

(2) error of this example observational data gets [Choisnard J, Laroche S.Properties of variational dataassimilation for synthetic aperture radar wind retrieval [J] .Journal of Geophysical Research:Oceans (1978 – 2012), 2008,113 (C5) .] statistical value 1.7m/s, then utilize the remaining difference method structural setting error co-variance matrix of observation;

(3) finally adopt variation analysis method to carry out fusion to wind striped wind field and pattern wind field and obtain optimized analysis wind field.

5th step, observation field, ambient field and analysis field and buoy measured data are compared, verify validity of the present invention, the result is as shown in table 1.

The root-mean-square error statistical nature of table 1 observation field, ambient field and analysis field

As known from Table 1, be optimized adjustment by observation field (Ocean Wind-field based on the inverting of wind striped) to ambient field (Ocean Wind-field based on the inverting of WRF pattern wind direction) through the present invention to obtain analyzing wind field, the root-mean-square error analyzing wind field is significantly less than the root-mean-square error of Background Winds, inversion accuracy obtains effective raising, demonstrates validity of the present invention.

Below be only the preferred embodiment of the present invention; be noted that for those skilled in the art; under the premise without departing from the principles of the invention, can also make some improvements and modifications, these improvements and modifications also should be considered as protection scope of the present invention.

Claims (5)

1., based on a synthetic-aperture radar Ocean Wind-field fusion method for variation, it is characterized in that, comprise the following steps:

Step 1, carries out experience differentiation to SAR image after calibration, SAR image is divided into the obvious region of feature and feature Fuzzy region;

Step 2, the obvious region of the feature described in selecting step 1, utilizes D continue wavelet transform inversion method wind direction, and in conjunction with physical geography module function calculation of wind speed, obtain the obvious region of described feature based on wind striped inverting Ocean Wind-field;

Step 3, choose the NCEP data matched with SAR observed image space-time to be input to WRF pattern and to carry out analog computation WRF wind direction, then WRF wind direction and corresponding SAR data are input to physical geography module function calculation of wind speed jointly, obtain the Ocean Wind-field of whole SAR observation area based on the inverting of WRF wind direction;

Step 4, with in step 2 obtain high-precision based on wind striped inverting Ocean Wind-field for observation field, with the Ocean Wind-field based on the inverting of WRF wind direction obtained in step 3 for ambient field, utilize variational method to carry out fusion adjustment by observation field to ambient field, obtain the Ocean Wind-field after merging.

2. a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation according to claim 1, it is characterized in that, step 2 specifically comprises the following steps,

201, D continue wavelet transform is carried out to the SAR image in the obvious region of wind streak feature, obtains the Wavelet Energy Spectrum image under different scale,

Described D continue wavelet transform method is two-dimentional Mexican-hat wavelet transformation,

Described two-dimentional Mexican-hat wavelet transformation formula is formula (1):

${\mathrm{\&psi;}}_H\left(a\right)=(a\&CenterDot;a)e^{\left(\frac{1}{2(a\&CenterDot;a)}\right)}---\left(1\right)$

In formula: a represents the variable of two-dimensional space-frequency field; Represent inner product of vectors, Ψ _ha () represents wavelet transform function;

202, two-dimensional fast fourier transform is carried out to Wavelet Energy Spectrum image, calculates the wavenumber spectrum of SAR image apoplexy striped:

Two-dimensional fast fourier transform formula is formula (2):

$Y_{l,m}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{b=1}{\overset{N}{\mathrm{\&Sigma;}}}{X}_{j,b}{e}^{-2\mathrm{\&pi;i}(\mathrm{jl}+\mathrm{bm})/N}---\left(2\right)$

Wherein, Y is the wavenumber spectrum of SAR image apoplexy striped, and X is image intensity value, l, m=1,2 ..., N; L, m represent that SAR image pixel is horizontal, lengthwise position; J, b represent that SAR image pixel is horizontal, lengthwise position variable;

203, the line of two-dimentional wavenumber spectrum peak value is done vertical line, after wind direction deblurring is carried out to described vertical line, obtains wind direction of ocean surface.

3. a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation according to claim 1, is characterized in that,

Described geophysical model is C-band geophysical model CMOD5, and described C-band physical geography module function is:

σ _C(θ,φ,u)＝10 ^A(u,θ)(1+B(u,θ)cos(φ)+C(u,θ)cos2φ) (3)

In formula, σ _crepresent VV polarimetric radar backscattering coefficient, φ represent wind direction relative to radar look-directions apparent wind to, θ represents radar antenna incident angle, and u represents ocean surface wind speed; A (u, θ), B (u, θ) and C (u, θ) represent by ten meters, sea At The Height wind speed, apparent wind to, polarization mode, radar frequency and incident angle determine coefficient.

4. a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation according to claim 1, it is characterized in that, utilize variational method to carry out fusion adjustment by observation field to ambient field in step 4, obtain the Ocean Wind-field after merging, specifically comprise the following steps:

401, adopt two-dimentional variational approach to carry out fusion calculation, the objective function defined by formula (4) is minimum to be retrained:

J＝J _q+J _m(4)

The objective function J of its apoplexy striped wind field _qwith the objective function J of WRF pattern _mbe defined as formula (5) and formula (6):

$J_q=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^{-1}(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)---\left(5\right)$

$J_m=\frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)---\left(6\right)$

Wherein, U and V merges the zonal wind and meridional wind that obtain, q and m represents wind striped wind field and region WRF pattern respectively, i.e. U _q, V _qrepresent zonal wind and the meridional wind of wind striped wind field respectively, U _m, V _mrepresent the zonal wind under the WRF pattern of region and meridional wind respectively, described Q _u, Q _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of wind striped wind field respectively, M _u, M _vrepresent zonal wind error co-variance matrix and the meridional wind error co-variance matrix of WRF pattern respectively, H _qfor space interpolation operator, analysis field is projected to observation space;

402, utilize Kriging interpolation method to construct space interpolation operator H _q, assuming that V _ifor the wind speed and direction value of i-th point in outside wind field, i=1,2, n, for outside wind field is interpolated into the wind speed and direction value of a kth observation station on wind striped wind field region, k is a kth point on wind striped wind field region, k=1,2, K, be the contribution weights of i-th data to a kth impact point, order:

$\stackrel{\&OverBar;}{Y}={({\stackrel{\&OverBar;}{V}}_1,{\stackrel{\&OverBar;}{V}}_2,...,{\stackrel{\&OverBar;}{V}}_k)}^T,$

$X={(V_1,V_2,...V_N)}^T,$

Then:

$\stackrel{\&OverBar;}{Y}=H_{q}X$

Observation Operators H _qbe expressed as formula (7),

$H_q={\left(h_i^k\right)}_{K\&times;N}---\left(7\right)$

Contribution weight obtained by Kriging system of equations:

Wherein, Lagrange multiplier when μ is minimization process, c (x _i, x _j) be net point x in outside wind field _iwith x _jbetween covariance function, x ^kfor the interpolation point on wind striped wind field region;

403, utilize the remaining difference method structural setting error co-variance matrix of observation, background error covariance matrix is formula (9) and (10):

＜U _m-U _q＞ ²＝dQ _u+dM _u(9)

＜V _m-V _q＞ ²＝dQ _v+dM _v(10)

Wherein < > represents time average, dQ _u, dM _urepresent the zonal wind error variance of wind striped wind field and the zonal wind error variance of WRF pattern respectively, dQ _v, dM _vrepresent the meridional wind error variance of wind striped wind field and the meridional wind error variance of WRF pattern respectively, dQ _u, dQ _v, dM _u, dM _vrespectively with Q _u, Q _v, M _u, M _vdiagonal element consistent;

404, based target function is minimum, adopts variation analysis method wind striped wind field and WRF pattern wind field to be merged to the zonal wind u and meridional wind v that obtain optimized analysis wind field.

5. a kind of synthetic-aperture radar Ocean Wind-field fusion method based on variation according to claim 4, it is characterized in that, step 404 specifically comprises the following steps:

The zonal wind u of optimized analysis wind field and meridional wind v is the zonal wind U and meridional wind V that after meeting fusion, the fusion of objective function least commitment condition obtains, and meets formula (11):

$\begin{array}{c}J=\frac{1}{2}{(H_{q}U-U_q)}^T{Q_u}^1(H_{q}U-U_q)+\frac{1}{2}{(H_{q}V-V_q)}^T{Q_v}^{-1}(H_{q}V-V_q)+\\ \frac{1}{2}{(U-U_m)}^T{M_u}^{-1}(U-U_m)+\frac{1}{2}{(V-V_m)}^T{M_v}^{-1}(V-V_m)=\min !\end{array}---\left(11\right)$ Based on variational method, can obtain:

$\begin{array}{c}\mathrm{\&delta;J}=H_q^T{Q}_u^{-1}(H_{q}U-U_q)\mathrm{\&delta;U}+H_q^T{Q_v}^{-1}(H_{q}V-V_q)\mathrm{\&delta;V}+\\ {M_u}^{-1}(U-U_m)\mathrm{\&delta;U}+{m_v}^{-1}\left({V-V}_m\right)\mathrm{\&delta;V}=0!\end{array}---\left(12\right)$

Utilize the arbitrariness of δ U, δ V, U, V meet following condition:

$H_q^T{Q}_u^{-1}(U_q-H_{q}U)+{M_u}^{-1}(U-U_m)=0$

$H_q^T{Q_v}^{-1}(V_q-H_{q}V)+{M_v}^{-1}(V_m-V)=0---\left(13\right)$

Thus show that the zonal wind u of optimized analysis wind field and meridional wind v is:

$u=U={[H_q^T{Q}_u^{-1}{H}_q+{M_u}^{-1}]}^{-1}\left[H_q^T{{Q_u^{-1}{U}_q+M}_u}^{-1}{U}_m\right]$

$v=V={[H_q^T{Q_v}^{-1}{H}_q+{H_v}^{-1}]}^{-1}[H_q^T{Q_v}^{-1}{V}_q+{M_V}^{-1}{V}_m].$