CN102073037B

CN102073037B - Iterative Ocean Current Inversion Method Based on Adaptive Threshold Selection Technique

Info

Publication number: CN102073037B
Application number: CN2011100008584A
Authority: CN
Inventors: 卢志忠; 戴运桃; 刘利强; 贾瑞才; 王淑娟; 李焱; 李英; 张立娜; 王凯; 胡浩
Original assignee: Harbin Engineering University
Current assignee: Harbin Hatran Navigation Technology Co ltd
Priority date: 2011-01-05
Filing date: 2011-01-05
Publication date: 2012-07-11
Anticipated expiration: 2031-01-05
Also published as: CN102073037A

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

Iterative Ocean Current Inversion Method Based on Adaptive Threshold Selection Technique

技术领域 technical field

本发明涉及一种遥感技术，特别涉及基于X波段导航雷达的海洋遥感技术。The invention relates to a remote sensing technology, in particular to an ocean remote sensing technology based on X-band navigation radar.

背景技术 Background technique

测流仪(海流计、声学海流剖面仪)能够测量海表面的流信息，但其容易受人为或自然灾害的破坏，所以目前仅仅用于海洋科考调查中，目前我国还没有对海面流场形成实时有效的业务化监测能力。Current meters (current meters, acoustic current profilers) can measure current information on the sea surface, but they are easily damaged by man-made or natural disasters, so they are currently only used in marine scientific investigations. Form a real-time and effective operational monitoring capability.

X波段雷达不仅可以监视移动的目标，其工作在短脉冲模式下时，其可以测量海流和海浪信息，流速(大小、方向)的测量是X波段雷达海浪监测的关键技术。与其它海洋传感器(波浪浮标，测流仪)相比，X波段雷达探测范围更广，安全性高，维护简单等优点，可以作为波浪浮标和测流仪的替代品。国内已有多家单位就X波段雷达测流问题开展初步研究。目前，常用的海流反演技术实现如下：X-band radar can not only monitor moving targets, but also measure ocean current and wave information when it works in short pulse mode. The measurement of current velocity (size, direction) is the key technology for X-band radar ocean wave monitoring. Compared with other marine sensors (wave buoys, current meters), the X-band radar has a wider detection range, high safety, and simple maintenance. It can be used as a substitute for wave buoys and current meters. Many units in China have carried out preliminary research on the X-band radar flow measurement problem. At present, the commonly used ocean current inversion techniques are implemented as follows:

(1)雷达图像的时间和空间序列g(x，y，t)经3维离散傅里叶变换到波数频率域F⁽³⁾(k_x，k_y，ω)，其中，x，y分别表示雷达图像空间的两个分量，t表示雷达图像的时间轴，k_x，k_y分别为波数的x与y分量，ω为波浪的频率。(1) The time and space sequence g(x, y, t) of the radar image is transformed into the wavenumber frequency domain F ⁽³⁾ (k _x , _ky , ω) by 3D discrete Fourier transform, where x, y are respectively Represents the two components of the radar image space, t represents the time axis of the radar image, k _x , k _y are the x and y components of the wave number, and ω is the frequency of the wave.

(2)由于海浪本身的波数和频率符合色散关系，如下所示：(2) Since the wave number and frequency of the ocean wave itself conform to the dispersion relation, as shown below:

$S S ((k k)) = = \sqrt{g g | | k k | | tanh tanh ((| | k k | | h h))} + + k k \cdot &Center Dot; u u - - - - - - ((11))$

其中g为重力加速度，k＝(k_x，k_y)为波数，h为水深，u为流速。Where g is the acceleration of gravity, k=(k _x , _ky ) is the wave number, h is the water depth, and u is the flow velocity.

在图像谱中选一个能量阈值C，要求能量值大于C的离散点的数目一般为数十个，使用这些的离散能量点，根据海浪本身波数和频率的特性，使用能量加权的最小二乘法拟合得到粗略的流速。Select an energy threshold C in the image spectrum, requiring that the number of discrete points whose energy value is greater than C is generally tens of points. Using these discrete energy points, according to the wave number and frequency characteristics of the wave itself, use the energy-weighted least squares method to fit to get a rough flow rate.

(3)首先根据得到的流速，确定0、1阶色散关系。然后根据经验选取能量阈值C₁，要求C₁远小于C，选取大于C₁的能量离散点，然后分别判断这些点属于0阶还是1阶次波，根据能量加权最小二乘法拟合出流速，将新的流速确定新的0、1阶色散关系，重复本步骤，将得到不断准确的流速。(3) First, determine the 0- and 1-order dispersion relations based on the obtained flow velocity. Then select the energy threshold C ₁ based on experience, requiring C ₁ to be much smaller than C, select energy discrete points greater than C ₁ , and then judge whether these points belong to 0-order or 1-order waves, and fit the flow velocity according to the energy-weighted least square method, Determine the new 0-1 order dispersion relationship with the new flow velocity, and repeat this step to obtain continuously accurate flow velocity.

${w w}_{p p} = = ((p p + + 11)) \sqrt{\frac{gk gk}{p p + + 11} tanh tanh ((\frac{kh kh}{p p + + 11}))} + + k k \cdot &Center Dot; u u - - - - - - ((22))$

其中，p＝0，1时分别对应0、1阶色散关系。Wherein, p=0 and 1 correspond to 0- and 1-order dispersion relations respectively.

这种方法主要存在以下不足：①X波段雷达测量得到的海表层流应该为分析的区域内平均流信息，而能量加权的最小二乘方法拟合得到的流，倾向于能量较大波的流；②依靠经验选取阈值的方法适应性不好，不同海域、海况的经验值不同，需要经过大量实验得到，通用性不好。This method mainly has the following deficiencies: ①The surface laminar current measured by X-band radar should be the average flow information in the analyzed area, while the flow obtained by fitting the energy-weighted least squares method tends to flow with relatively large energy waves; ② The method of relying on experience to select the threshold value has poor adaptability, and the experience values of different sea areas and sea conditions are different, and it needs to be obtained through a large number of experiments, and the versatility is not good.

发明内容 Contents of the invention

本发明的目的在于提供一种对流速测量准确，通用性好，能提高X波段雷达测量海表层流的精度的基于自适应阈值选取技术的迭代海流反演方法。The purpose of the present invention is to provide an iterative ocean current inversion method based on an adaptive threshold selection technology that is accurate in flow velocity measurement, has good versatility, and can improve the accuracy of X-band radar in measuring sea surface laminar flow.

本发明的目的是这样实现的：The purpose of the present invention is achieved like this:

(1)雷达图像的采集和A\D转换得到32幅原始雷达图像，(1) Radar image collection and A\D conversion to get 32 original radar images,

X波段导航雷达工作在短脉冲下，发射的电磁波与海面的毛细波发生布拉格散射，后向散射被接收机接收，视频信号经A/D转换为数字信号，两个船艏信号间的视频信号组成一幅雷达图像，选择32幅连续的雷达图像g(x，y，t)，其中，x，y分别为雷达图像空间的两分量，t为时间分量；The X-band navigation radar works under short pulses. Bragg scattering occurs between the emitted electromagnetic waves and the capillary waves on the sea surface, and the backscattering is received by the receiver. The video signal is converted into a digital signal by A/D, and the video signal between the two bow signals To form a radar image, select 32 consecutive radar images g(x, y, t), where x, y are the two components of the radar image space, and t is the time component;

(2)得到笛卡尔坐标下的图像并将其作傅里叶变换得到图像谱；(2) Obtain the image under the Cartesian coordinates and do Fourier transform to obtain the image spectrum;

在雷达图像中选择需要分析矩形框，使用最近点插值实现矩形框内雷达图像由极坐标向笛卡尔坐标转换表示为g1(x，y，t)；矩形框内雷达图像序列g1(x，y，t)经离散3维傅里叶变换后成为波数频率域；Select the rectangular frame to be analyzed in the radar image, and use the nearest point interpolation to realize the transformation of the radar image in the rectangular frame from polar coordinates to Cartesian coordinates and express it as g1(x, y, t); the radar image sequence in the rectangular frame g1(x, y , t) becomes wave number frequency domain after discrete 3-dimensional Fourier transform;

$γ γ ((k k,, w w)) = = {&Integral; &Integral;}_{00}^{{L L}_{X x}} {&Integral; &Integral;}_{00}^{{L L}_{Y Y}} {&Integral; &Integral;}_{00}^{T T} g g ((x x,, y the y,, t t)) exp exp [[i i (({k k}_{x x} x x + + {k k}_{y the y} y the y - - wt wt))]] dxdydt dxdydt$

其中k＝(k_x，k_y)，L_x，L_y，T分别为矩形区域的长、宽和时间序列的时间总长度三维傅里叶变换后考虑都实际情况，和消除谱的180度模糊问题，只保留w＞0的部分因此3维图像的能量谱为：Among them k=(k _x , k _y ), L _x , Ly _y , T are respectively the length and width of the rectangular area and the time total length of the time series. After the three-dimensional Fourier transform, consider the actual situation, and eliminate the 180 degrees of the spectrum Fuzzy problem, only keep the part of w>0, so the energy spectrum of the 3D image is:

${F f}^{((33))} ((k k,, w w)) = = \frac{11}{{L L}_{X x} {L L}_{Y Y} T T} {γ γ}^{22} ((k k,, w w));;$

(3)初始估值，仅考虑0阶次波的影响，反演得到粗略的流；(3) Initial estimation, only considering the influence of the 0th order wave, the inversion obtains a rough flow;

在波数频率域中，选取数十个能量值较大的点，根据海浪信号波数和频率间的色散关系，使用平均加权的最小二乘法，拟合出流信息，基本的数学模型如下：In the wave number frequency domain, select dozens of points with large energy values, and use the average weighted least square method to fit the outflow information according to the dispersion relationship between the wave number and frequency of the ocean wave signal. The basic mathematical model is as follows:

为了得出表层流u，将式In order to obtain the surface flow u, the formula

${Q Q}^{22} = = {Σ Σ}_{i i = = 11}^{n no} {(({ω ω}_{i i} - - S S (({k k}_{i i}))))}^{22}$

对u_x、u_y求一阶偏导数，并使其为零：Find the first partial derivative of u _x , u _y and make it zero:

$&PartialD; &PartialD; {Q Q}^{22} / / &PartialD; &PartialD; {u u}_{x x} = = 00$

$&PartialD; &PartialD; {Q Q}^{22} / / &PartialD; &PartialD; {u u}_{y the y} = = 00$

将上式写成矩阵的形式：Write the above formula in matrix form:

$(\begin{matrix} {D D.}_{xx xxx} & {D D.}_{xy xy} \\ {D D.}_{yx yx} & {D D.}_{yy yy} \end{matrix}) (\begin{matrix} {u u}_{x x} \\ {u u}_{y the y} \end{matrix}) = = (\begin{matrix} {b b}_{x x} \\ {b b}_{y the y} \end{matrix})$

简写为：Abbreviated as:

Du＝bDu = b

矩阵D的元素为The elements of matrix D are

${D D.}_{xx xxx} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{}^{ix ix},, {D D.}_{xy xy} = = {D D.}_{yx yx} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{ix ix} {k k}_{iy iy}$

${D D.}_{yy yy} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{}^{iy iy}$

矢量b的元素为The elements of the vector b are

${b b}_{x x} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{ix ix} {w w}_{iD iD},,$ ${b b}_{y the y} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{iy iy} {w w}_{iD iD}$

其中，w_iD＝w_i-ζ(k_i)为了计算出流u，D必须为非奇异矩阵(det(D)≠0)即D^-1是存在的，此时Among them, w _iD =w _i -ζ(k _i ) In order to calculate the flow u, D must be a non-singular matrix (det(D)≠0), that is, D ^-1 exists, at this time

u＝D^-1·b；u=D ^-1 ·b;

(4)基于自适应阈值选取技术，考虑0、1阶次波的影响，使用平均加权最小二乘方法迭代拟合流信息；(4) Based on the adaptive threshold selection technology, considering the influence of 0 and 1 order waves, the average weighted least squares method is used to iteratively fit the flow information;

①使用自适应阈值选取技术，得到阈值C_it，大于C_it的能量就包含了0阶次和1阶次海浪能量，数目为N₁，利用初始估值得到的流，由式

计算0阶次和1阶次波频率；①Using the adaptive threshold selection technique, the threshold C _it is obtained. The energy greater than C _it includes the wave energy of the 0th order and the 1st order, and the number is N ₁ . The flow obtained by the initial estimate is expressed by

Calculate the 0-order and 1-order wave frequencies;

应用平均最小二乘法，拟合得到新的表层流；② Determine whether the measured N ₁ energy points conform to the 0-order or 1-order wave dispersion relation, if |w _i -w(k _i )|<|w _i -w ₁ (k _i )|, then the The frequency conforms to the 0th order wave; if |w _i -w(k _i )|＞|w _i -w ₁ (k _i )|, then the frequency conforms to the 1st order wave, and the judged data will be judged according to different minimum value function as

Apply the average least squares method to fit the new surface flow;

③应用新的表层流构建带通滤波器，得到新的海浪信噪比，使用自适应阈值选取技术，得到新的阈值；③Use the new surface flow to construct a band-pass filter to obtain a new wave signal-to-noise ratio, and use the adaptive threshold selection technology to obtain a new threshold;

④将新得到的表层流，代入式

得到新的0阶次和1阶次波频率，重复上述步骤，将得到不断精确的表层流。④ Substituting the newly obtained surface flow into the formula

Obtain the new 0-order and 1-order wave frequencies, and repeat the above steps to obtain continuously accurate surface flow.

所述自适应阈值选取技术为：The adaptive threshold selection technique is:

①在较短的时间内，海浪在时间和空间上的变化较小，近似的认为前个序列图的信噪比为本序列图的信噪比SNR；①In a relatively short period of time, the changes in time and space of the ocean waves are small, and the signal-to-noise ratio of the previous sequence diagram is approximately considered to be the signal-to-noise ratio SNR of the sequence diagram;

②由3维傅里叶变换得到的图像谱F⁽³⁾(k，w)和信噪比SNR得到图像谱中含有的海浪能量P② The wave energy P contained in the image spectrum is obtained from the image spectrum F ⁽³⁾ (k, w) obtained by 3-dimensional Fourier transform and the signal-to-noise ratio SNR

$P P = = \frac{SNR SNR + + 11}{SNR SNR} {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {Σ Σ}_{k k}^{{N N}_{w w}} {F f}^{((33))} (({k k}_{ix ix},, {k k}_{jy jy},, {w w}_{k k})) \cdot &Center Dot; {Δk Δk}_{x x} {Δk Δk}_{y the y} Δw Δw$

③对图像谱中所有的能量点按照从大到小的顺序进行排序，排序后第i采样点的能量为P_i，因为海浪能量点一般较背景噪声能量点大，因此可令③Sort all the energy points in the image spectrum in order from large to small. After sorting, the energy of the i-th sampling point is P _i , because the wave energy point is generally larger than the background noise energy point, so it can be set

$P P = = {Σ Σ}_{i i = = 11}^{m m} {P P}_{i i}$

由上式求出用于LSM的能量点数目m，第m个采样点的能量值为P_m，则P_m为迭代能量阈值C_it。The number m of energy points used for LSM is obtained from the above formula, and the energy value of the mth sampling point is P _m , then P _m is the iteration energy threshold C _it .

本发明的主要技术要点为：The main technical points of the present invention are:

(1)初始估值，获得粗略的流速(大小、方向)。假定矩形框内海浪空间均匀，时间稳定，则海浪模型符合高斯分布。当没有流速的情况下，海表面张力波和重力波的色散关系为(1) Initial estimation to obtain a rough flow velocity (size, direction). Assuming that the waves in the rectangular box are uniform in space and stable in time, the wave model conforms to the Gaussian distribution. When there is no current velocity, the dispersion relationship between sea surface tension waves and gravity waves is

$ζ ζ ((k k)) \approx \approx \sqrt{g g | | k k | | tanh tanh ((| | k k | | h h))} - - - - - - ((33))$

其中，ζ(k)是波频率，k是波数矢量，h是水深，g是重力加速度。如果有一个相对于雷达的表层流u，频率中就引进了多普勒频移项，色散关系变为：where ζ(k) is the wave frequency, k is the wave number vector, h is the water depth, and g is the gravitational acceleration. If there is a surface flow u relative to the radar, a Doppler shift term is introduced in frequency, and the dispersion relation becomes:

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

其中，S(k)为理论的海浪频率，k·u为多普勒频移，u＝(u_x，u_y)是海浪场与雷达天线平台间的相对速度，包括雷达平台运动速度(比如船速)和表层流速度矢量，θ为波数k与u之间夹角，当θ＝90时，多普勒频移项为0，只有流u在波数k的方向上有分量时才能影响多普勒频移，当雷达平台静止时，u指的是海表层流。Among them, S(k) is the theoretical wave frequency, k u is the Doppler frequency shift, u=(u _x , u _y ) is the relative velocity between the wave field and the radar antenna platform, including the movement speed of the radar platform (such as Vessel speed) and surface flow velocity vector, θ is the angle between wave number k and u, when θ = 90, the Doppler frequency shift item is 0, only when the flow u has a component in the direction of wave number k, it can affect the Puller shift, when the radar platform is stationary, u refers to the sea surface laminar current.

通过平均加权最小二乘法可得到粗略的流信息。平均加权的最小二乘法的极小值函数定义如式(5)所示：Rough flow information can be obtained by mean weighted least squares method. The definition of the minimum value function of the average weighted least squares method is shown in formula (5):

${Q Q}^{22} = = {Σ Σ}_{i i = = 11}^{n no} {(({ω ω}_{i i} - - S S (({k k}_{i i}))))}^{22} - - - - - - ((55))$

其中，ω_i为雷达观测的海浪频率，n为用于最小二乘法的能量点个数。Among them, ω _i is the wave frequency observed by the radar, and n is the number of energy points used for the least square method.

(2)阈值自动选取技术获得阈值C_it。(2) Threshold automatic selection technology obtains threshold C _it .

当风速较大、波高较高时，雷达图像的回波较强，海浪信号的信噪比较大，此时的数据的质量较好。反之，海浪信号信噪比较小，数据质量差。通过海浪信号的信噪比判断数据质量的好坏，依此选取阈值C_it。When the wind speed is high and the wave height is high, the echo of the radar image is strong, the signal-to-noise ratio of the wave signal is high, and the quality of the data at this time is better. On the contrary, the signal-to-noise ratio of the ocean wave signal is small, and the data quality is poor. The quality of the data is judged by the signal-to-noise ratio of the ocean wave signal, and the threshold C _it is selected accordingly.

首先给出海浪图像能量信噪比的定义：Firstly, the definition of energy signal-to-noise ratio of ocean wave image is given:

$SNR SNR = = \frac{SIG SIG}{BGN BGN} - - - - - - ((66))$

$SIG SIG = = {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {F f}^{((22))} (({k k}_{ix ix},, {k k}_{jy jy})) {Δk Δk}_{x x} {Δk Δk}_{y the y} - - - - - - ((77))$

$BGN BGN = = {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {Σ Σ}_{k k}^{{N N}_{ω ω}} {F f}^{((33))} (({k k}_{ix ix},, {k k}_{jy jy},, {ω ω}_{k k})) {Δk Δk}_{x x} {Δk Δk}_{y the y} Δω Δω - - - - - - ((88))$

$- - {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {F f}^{((22))} (({k k}_{ix ix},, {k k}_{jy jy})) {Δk Δk}_{x x} {Δk Δk}_{y the y}$

${F f}^{((22))} ((k k)) = = 22 {&Integral; &Integral;}_{ω ω > > 00} {F f}^{((33))} (({k k}_{x x},, {k k}_{y the y},, ω ω)) δ δ ((ω ω - - {ω ω}_{00})) dω dω - - - - - - ((99))$

其中，SIG为波浪谱能量，BGN为背景噪声，F⁽³⁾是3维图像谱，F⁽²⁾(k)是2维图像谱，N_kx、N_ky、N_ω为谱的范围，Δk_x、Δk_y为波数分辨率，Δω为频率分辨率，δ(ω-ω₀)是带通滤波器。Among them, SIG is wave spectrum energy, BGN is background noise, F ⁽³⁾ is 3D image spectrum, F ⁽²⁾ (k) is 2D image spectrum, N _kx , N _ky , N _ω are spectrum ranges, Δk _x and Δk _y are wave number resolution, Δω is frequency resolution, and δ(ω-ω ₀ ) is a bandpass filter.

阈值自动选取方法的基本实现，结合附图7说明。The basic realization of the threshold automatic selection method is described in conjunction with accompanying drawing 7.

①在较短的时间内，海浪的信噪比变化较小即在时间上是稳定的，因此假设之前得到的信噪比为本序列图的信噪比SNR。①In a short period of time, the SNR of ocean waves changes little, that is, it is stable in time. Therefore, it is assumed that the SNR obtained before is the SNR of this sequence diagram.

②由3维傅里叶变换得到的图像谱F⁽³⁾(k，w)和信噪比SNR可以得到图像谱中含有的海浪能量P。② The wave energy P contained in the image spectrum can be obtained from the image spectrum F ⁽³⁾ (k, w) obtained by 3-dimensional Fourier transform and the signal-to-noise ratio SNR.

$P P = = \frac{SNR SNR + + 11}{SNR SNR} {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {Σ Σ}_{k k}^{{N N}_{w w}} {F f}^{((33))} (({k k}_{ix ix},, {k k}_{jy jy},, {w w}_{k k})) \cdot &Center Dot; {Δk Δk}_{x x} {Δk Δk}_{y the y} Δw Δw - - - - - - ((1010))$

$P P = = {Σ Σ}_{i i = = 11}^{m m} {P P}_{i i} - - - - - - ((1111))$

由上式即可求出用于LSM的能量点数目m，第m个采样点的能量值为P_m，则P_m为迭代能量阈值C_it。The number m of energy points used for LSM can be obtained from the above formula, and the energy value of the mth sampling point is P _m , then P _m is the iteration energy threshold C _it .

(3)基于阈值自动选取技术的迭代估值。(3) Iterative estimation based on automatic threshold selection technology.

在此，我们只考虑0、1阶次谐波的情况。Here, we only consider the case of 0 and 1 order harmonics.

①应用自适应阈值选取技术，选取阈值C_it，可以认为大于C_it的能量部分包含了0阶次和1阶次海浪能量，数目为N₁。利用初始估值得到的流，分别由式(12)计算0阶次和1阶次波频率。①Applying the adaptive threshold selection technology to select the threshold C _it t, it can be considered that the energy part greater than C _it includes the wave energy of the 0th order and the 1st order, and the number is N ₁ . Using the flow obtained from the initial estimate, the 0-order and 1-order wave frequencies are calculated by Equation (12), respectively.

${w w}_{p p} = = ((p p + + 11)) \sqrt{\frac{gk gk}{p p + + 11} tanh tanh ((\frac{kh kh}{p p + + 11}))} + + k k \cdot &Center Dot; u u - - - - - - ((1212))$

p＝0、1分别对应0、1阶次波的色散关系。p=0, 1 correspond to the dispersion relations of 0 and 1 order waves respectively.

③应用新的表层流构建带通滤波器，得到新的海浪信号信噪比，应用自适应阈值选取技术，得到新的阈值。③Use the new surface flow to construct a band-pass filter to obtain a new wave signal-to-noise ratio, and apply the adaptive threshold selection technology to obtain a new threshold.

④将新得到的表层流，代入式(12)中，得到新的0阶次和1阶次波色散关系，重复上述步骤，将得到不断精确的表层流。④ Substitute the newly obtained surface flow into Equation (12) to obtain the new 0-order and 1-order wave dispersion relations, and repeat the above steps to obtain continuously accurate surface flow.

发明的基于自动阈值选取的流反演新方法与现有技术相比的优点：The advantages of the invented new flow inversion method based on automatic threshold selection compared with the existing technology:

(1)不同风速和区域下雷达序列图中海浪信号的是不同的，主要是通过信噪比体现。由信噪比和图像的总能量可以计算出图像序列中海浪的能量，由此可以得出阈值。根据以上思路本发明提出了阈值自动选取技术，它能够自动识别海浪信号好坏，并通过迭代技术，精确的从图像谱中分离出海浪能量，弥补了经验值的不足。(1) The wave signals in the radar sequence diagrams are different under different wind speeds and regions, mainly reflected by the signal-to-noise ratio. The energy of the waves in the image sequence can be calculated from the signal-to-noise ratio and the total energy of the image, and thus the threshold can be obtained. According to the above ideas, the present invention proposes an automatic threshold value selection technology, which can automatically identify whether the wave signal is good or bad, and accurately separate the wave energy from the image spectrum through iterative technology, making up for the lack of empirical value.

(2)从图1、2、表1可知，当流速较大时(＞0.3m/s)，本发明方法的流反演结果较原方法结果具有两方面优势：①能更好反映流的趋势；②具有更小的误差和更好的稳定性。(2) From Figures 1, 2, and Table 1, it can be seen that when the flow velocity is large (>0.3m/s), the flow inversion results of the method of the present invention have two advantages compared with the results of the original method: ①It can better reflect the flow trend; ② smaller error and better stability.

(3)流速较小时(＜0.2m/s)，信号较弱，从弱信号中提取流信息是一个难题。本发明的方法的反演结果较原方法已有很大的改进如图3、4、表1所示。(3) When the flow velocity is small (<0.2m/s), the signal is weak, and it is a difficult problem to extract flow information from the weak signal. Compared with the original method, the inversion result of the method of the present invention has been greatly improved as shown in Fig. 3, 4 and Table 1.

附图说明Description of drawings

图1流速较大时，本发明中方法、原方法与现场传感器流速结果比较。When Fig. 1 flow rate is bigger, method among the present invention, original method and field sensor flow rate result comparison.

图2流速较大时，本发明中方法、原方法与现场传感器流向结果比较。When the flow velocity in Fig. 2 is relatively large, the method in the present invention, the original method and the field sensor flow direction results are compared.

图3流速较小时，本发明中方法、原方法与现场传感器流速结果比较。When the flow velocity of Fig. 3 is small, method among the present invention, original method and field sensor flow velocity result comparison.

图4流速较小时，本发明中方法、原方法与现场传感器流向结果比较。When the flow rate in Fig. 4 is small, the method in the present invention, the original method and the field sensor flow direction results are compared.

图5本发明的技术实施方案流程图。Fig. 5 is a flow chart of the technical implementation of the present invention.

图6基于自适应阈值选取技术的迭代流反演流程图。Figure 6 is a flowchart of iterative flow inversion based on adaptive threshold selection technology.

图7自适应阈值选取技术实施流程图。Figure 7 is a flow chart of the implementation of adaptive threshold selection technology.

图8表1流速较大时，各方法结果统计分析。Figure 8 Table 1 Statistical analysis of the results of each method when the flow rate is relatively large.

图9表2流速较小时，各方法结果统计分析。Statistical analysis of the results of each method when the flow rate in Table 2 of Fig. 9 is small.

具体实施方式 Detailed ways

下面结合图5、6、7，对本发明做更详细介绍：Below in conjunction with Fig. 5,6,7, the present invention is described in more detail:

(5)雷达图像的采集和A\D转换得到32幅原始雷达图像。(5) Radar image collection and A\D conversion to get 32 original radar images.

X波段导航雷达工作在短脉冲下，发射的电磁波与海面的毛细波发生布拉格散射，后向散射被接收机接收，视频信号经A/D转换为数字信号，两个船艏信号间的视频信号组成一幅雷达图像，在本发明中，需要使用32幅连续的雷达图像g(x，y，t)，其中，x，y分别为雷达图像空间的两分量，t为时间分量。The X-band navigation radar works under short pulses. Bragg scattering occurs between the emitted electromagnetic waves and the capillary waves on the sea surface, and the backscattering is received by the receiver. The video signal is converted into a digital signal by A/D, and the video signal between the two bow signals To form a radar image, in the present invention, 32 consecutive radar images g(x, y, t) need to be used, where x, y are respectively two components of the radar image space, and t is a time component.

(6)得到笛卡尔坐标下的图像并将其作傅里叶变换得到图像谱。(6) Obtain the image under Cartesian coordinates and perform Fourier transform to obtain the image spectrum.

在雷达图像中选择需要分析矩形框，使用最近点插值实现矩形框内雷达图像由极坐标向笛卡尔坐标转换表示为g1(x，y，t)；矩形框内雷达图像序列g1(x，y，t)经离散3维傅里叶变换后成为波数频率域。Select the rectangular frame to be analyzed in the radar image, and use the nearest point interpolation to realize the transformation of the radar image in the rectangular frame from polar coordinates to Cartesian coordinates and express it as g1(x, y, t); the radar image sequence in the rectangular frame g1(x, y , t) becomes wave number frequency domain after discrete 3-dimensional Fourier transform.

$γ γ ((k k,, w w)) = = {&Integral; &Integral;}_{00}^{{L L}_{X x}} {&Integral; &Integral;}_{00}^{{L L}_{Y Y}} {&Integral; &Integral;}_{00}^{T T} g g ((x x,, y the y,, t t)) exp exp [[i i (({k k}_{x x} x x + + {k k}_{y the y} y the y - - wt wt))]] dxdydt dxdydt - - - - - - ((1313))$

其中k＝(k_x，k_y)，L_x，L_y，T分别为矩形区域的长、宽和时间序列的时间总长度三维傅里叶变换后考虑都实际情况，和消除谱的180度模糊问题，我们只保留了w＞0的部分因此3维图像的能量谱为：Among them k=(k _x , k _y ), L _x , Ly _y , T are respectively the length and width of the rectangular area and the time total length of the time series. After the three-dimensional Fourier transform, consider the actual situation, and eliminate the 180 degrees of the spectrum For the blur problem, we only keep the part where w>0, so the energy spectrum of the 3D image is:

${F f}^{((33))} ((k k,, w w)) = = \frac{11}{{L L}_{X x} {L L}_{Y Y} T T} {γ γ}^{22} ((k k,, w w)) - - - - - - ((1414))$

(7)初始估值，仅考虑0阶次波的影响，反演得到粗略的流。(7) Initial estimation, only consider the influence of 0th order wave, and invert to obtain a rough flow.

为了得出表层流u，将式(15)In order to obtain the surface flow u, the formula (15)

${Q Q}^{22} = = {Σ Σ}_{i i = = 11}^{n no} {(({ω ω}_{i i} - - S S (({k k}_{i i}))))}^{22} - - - - - - ((1515))$

$\begin{matrix} &PartialD; &PartialD; {Q Q}^{22} / / &PartialD; &PartialD; {u u}_{x x} = = 00 \\ &PartialD; &PartialD; {Q Q}^{22} / / &PartialD; &PartialD; {u u}_{y the y} = = 00 \end{matrix} - - - - - - ((1616))$

将上式写成矩阵的形式：Write the above formula in matrix form:

$(\begin{matrix} {D D.}_{xx xxx} & {D D.}_{xy xy} \\ {D D.}_{yx yx} & {D D.}_{yy yy} \end{matrix}) (\begin{matrix} {u u}_{x x} \\ {u u}_{y the y} \end{matrix}) = = (\begin{matrix} {b b}_{x x} \\ {b b}_{y the y} \end{matrix}) - - - - - - ((1717))$

简写为：Abbreviated as:

Du＝b (18)Du=b (18)

矩阵D的元素为The elements of matrix D are

${D D.}_{xx xxx} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{}^{ix ix},, {D D.}_{xy xy} = = {D D.}_{yx yx} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{ix ix} {k k}_{iy iy} - - - - - - ((1919))$

${D D.}_{yy yy} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{}^{iy iy}$

矢量b的元素为The elements of the vector b are

${b b}_{x x} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{ix ix} {w w}_{iD iD},, {b b}_{y the y} = = {Σ Σ}_{i i = = 11}^{n no} {k k}_{iy iy} {w w}_{iD iD} - - - - - - ((2020))$

u＝D^-1·b (21)u＝D ^-1 ·b (21)

(8)基于自适应阈值选取技术，考虑0、1阶次波的影响，使用平均加权最小二乘方法迭代拟合流信息。结合附图6进行说明。自适应阈值选取技术具体实施详见第(5)部分。(8) Based on the adaptive threshold selection technology, considering the influence of 0 and 1 order waves, the average weighted least squares method is used to iteratively fit the flow information. Description will be made in conjunction with accompanying drawing 6. For the specific implementation of adaptive threshold selection technology, see section (5).

通过第(3)步骤初始估值得到的流速(大小、方向)是不准确的，需要考虑一阶次波的影响，使用迭代的方法提高流的反演精度，其实现方法结合附图6做如下说明：The flow velocity (size, direction) obtained through the initial estimation in step (3) is inaccurate. It is necessary to consider the influence of the first-order wave, and use an iterative method to improve the inversion accuracy of the flow. The implementation method is combined with Figure 6. As explained below:

①使用自适应阈值选取技术，得到阈值C_it，大于C_it的能量就包含了0阶次和1阶次海浪能量，数目为N₁。利用初始估值得到的流，由式(12)计算0阶次和1阶次波频率。①Use the adaptive threshold selection technology to obtain the threshold C _it , the energy greater than C _it includes the wave energy of the 0th order and the 1st order, and the number is N ₁ . Using the flow obtained from the initial estimate, the 0-order and 1-order wave frequencies are calculated by Equation (12).

③应用新的表层流构建带通滤波器，得到新的海浪信噪比，使用自适应阈值选取技术，得到新的阈值。③Use the new surface flow to construct a band-pass filter to obtain a new wave signal-to-noise ratio, and use the adaptive threshold selection technology to obtain a new threshold.

④将新得到的表层流，代入式(12)，得到新的0阶次和1阶次波频率，重复上述步骤，将得到不断精确的表层流。④ Substitute the newly obtained surface flow into equation (12) to obtain the new 0-order and 1-order wave frequencies, and repeat the above steps to obtain continuously accurate surface flow.

(9)在上步骤中使用到自适应阈值选取技术，现结合附图7进行说明，(9) In the last step, the adaptive threshold selection technology is used, which is now described in conjunction with accompanying drawing 7,

①在较短的时间内，海浪在时间和空间上的变化较小，因此可近似的认为前个序列图的信噪比为本序列图的信噪比SNR；①In a relatively short period of time, the changes in time and space of the ocean waves are small, so it can be approximately considered that the signal-to-noise ratio of the previous sequence diagram is the signal-to-noise ratio SNR of this sequence diagram;

②由3维傅里叶变换得到的图像谱F⁽³⁾(k，w)和信噪比SNR可以得到图像谱中含有的海浪能量P② The wave energy P contained in the image spectrum can be obtained from the image spectrum F ⁽³⁾ (k, w) obtained by 3-dimensional Fourier transform and the signal-to-noise ratio SNR

$P P = = \frac{SNR SNR + + 11}{SNR SNR} {Σ Σ}_{i i}^{{N N}_{kx x}} {Σ Σ}_{j j}^{{N N}_{ky ky}} {Σ Σ}_{k k}^{{N N}_{w w}} {F f}^{((33))} (({k k}_{ix ix},, {k k}_{jy jy},, {w w}_{k k})) \cdot &Center Dot; {Δk Δk}_{x x} {Δk Δk}_{y the y} Δw Δw - - - - - - ((22 twenty two))$

$P P = = {Σ Σ}_{i i = = 11}^{m m} {P P}_{i i} - - - - - - ((23 twenty three))$

将本发明提出的流反演方法应用到流速较大时的实测数据中(时间段为2009年10月06日上午10:36:00到下午16:40:00)结果如图1、2所示，统计分析结果如表1所示，可以看出本发明方法结果较原方法有更高的精度和稳定性。图3、4给出了本发明方法反演流较小时的结果(时间段为2009年10月06日上午08:16:00到09:20:00)，统计结果分析见表2，可知本发明方法结果较原方法结果精度有很大的提高。Apply the flow inversion method proposed by the present invention to the measured data when the flow velocity is relatively large (the time period is from 10:36:00 am to 16:40:00 pm on October 6, 2009) and the results are shown in Figures 1 and 2 Show, statistical analysis result is as shown in table 1, can find out that the method result of the present invention has higher precision and stability than former method. Fig. 3, 4 has provided the result (time section is 08:16:00 am to 09:20:00 am on October 6, 2009) when the inversion flow of the method of the present invention is less, and statistical result analysis is shown in Table 2, it can be seen that this Compared with the original method, the precision of the result of the invented method is greatly improved.

Claims

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{SNR + 1}{SNR} Σ_{i}^{N_{kx}} Σ_{j}^{N_{ky}} Σ_{k}^{N_{w}} F^{(3)} (k_{ix}, k_{jy}, ω_{k}) \cdot Δ k_{x} Δ k_{y} Δω

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 = Σ_{i = 1}^{m} 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