JP4825183B2 - Surface layer flow velocity estimation method, apparatus and program - Google Patents

Description

  The present invention relates to a surface layer flow velocity estimation method, apparatus, and program. More specifically, the present invention relates to a surface velocity estimation method, apparatus, and program suitable for use in observation of coastal surface current using data obtained by an ocean radar.

  Marine radar is a device that can continuously observe surface currents such as flow velocity, wave height, and wind direction over a long period of time regardless of weather by remote sensing from land. Coastal environment monitoring, port management, coastal disaster prevention, provision of fishery information, etc. Utilization in the field is underway.

  Observation by ocean radar irradiates a radar wave from the transmitting antenna toward the observation area, receives echoes from ocean waves, which are radio waves that are returned by Bragg scattering on the sea surface, and receives the echo from the time series of echo intensity. Obtain the Doppler spectrum and extract the surface flow.

  The echo received by the receiving antenna is a backscattered echo due to a wave component wave that satisfies the Bragg resonance condition, and the Doppler spectrum of this echo is constituted by the spectrum of the primary scattering and the secondary scattering. The flow velocity of the surface layer flow is obtained by analyzing the primary scattering spectrum, and the wave information is obtained by analyzing the primary scattering spectrum and the secondary scattering spectrum.

  The algorithm for detecting the primary scattering spectrum frequency when calculating the flow velocity from the Doppler spectrum of the backscattered echo usually detects the peak with the largest spectral value in the positive and negative frequency regions of the Doppler spectral distribution. Thus, those peak frequencies are obtained (Patent Document 1).

JP-A-11-83992

  However, the echo received by the receiving antenna is not due to the extrinsic noise and the receiver noise, such as receiving a large echo when there is a large ship that crosses the radar wave almost at right angles during radar wave irradiation. It is inevitable to be affected. And if such exogenous noise is mixed, the peak of Doppler spectrum becomes unclear or the peak due to exogenous noise becomes prominent, so it is necessary to detect the legitimate primary scattering spectrum frequency corresponding to the flow velocity. There is a problem that the observation accuracy of the surface layer flow deteriorates.

  In addition, it is necessary to provide information in real time or near real time, especially in coastal environment monitoring, harbor management, and coastal disaster prevention for sudden events, but the conventional method is very useful for detecting the peak frequency of Doppler spectrum. It takes a lot of time, so it is hard to say that information can be provided at the timing that is really needed.

  Therefore, an object of the present invention is to provide a surface layer flow velocity estimation method, apparatus, and program capable of detecting a legitimate peak frequency of a Doppler spectrum in a short time by eliminating the influence of exogenous noise.

In order to achieve this object, the surface layer flow velocity estimation method according to claim 1 includes a step of reading Doppler spectrum data for each beam direction distance and frequency of a processing target region, and a Doppler spectrum data as a spectrum in a positive frequency region. And a negative region spectrum, a positive region spectrum and a negative region spectrum are averaged in a frequency direction, and a positive region spectrum is a frequency direction of a positive region spectrum. Subtracting the average value of the positive region and calculating the deviation of the spectrum of the positive region and subtracting the average value in the frequency direction of the spectrum of the negative region from the spectrum of the negative region to calculate the deviation of the spectrum of the negative region; Calculate the average value of the sub-region for the deviation of the spectrum in the positive region and the spectrum of the negative region Subtracting the average value of the partial area for the deviation of the positive area, and adding the average value of the partial area of the positive area spectral deviation and the average value of the partial area of the negative area spectral deviation for each partial area A binarization matrix is created by binarizing the value of the component corresponding to the partial area depending on whether the total value for each partial area is greater than the binarization threshold, and the step of calculating the total value for each partial area A step of creating an adjacency matrix for the binarized matrix, a step of creating a stationary solution matrix by multiplying the adjacency matrix until it reaches a steady state, and a maximum sum value for each subregion from the stationary solution matrix. Extracting a component corresponding to a connected region composed of partial regions connected to a partial region and having a combined value greater than the binarization threshold for each partial region, and creating a filter base matrix from the binarized matrix; The filter matrix is created by limiting the spread of the filter transmission component in the frequency direction for each row of the filter fundamental matrix, and the spectral deviation in the positive region and the spectral deviation in the negative region are expressed in terms of the beam direction distance and frequency. Summing every time, calculating a spectrum deviation sum, applying a filter matrix to the spectrum deviation sum, creating a filter application spectrum, detecting a peak frequency from the filter application spectrum, And a step of calculating the beam direction flow velocity by multiplying by the conversion coefficient.

According to a second aspect of the present invention, there is provided a surface layer flow velocity estimation apparatus, comprising: means for reading Doppler spectrum data for each beam direction distance and frequency of a processing target region; a Doppler spectrum data for a spectrum having a positive frequency and a negative region; Means for separating the spectrum into a spectrum, means for calculating an average value in the frequency direction for the spectrum in the positive region and the spectrum in the negative region, and subtracting the average value in the frequency direction of the spectrum in the positive region from the spectrum in the positive region. Means for calculating the deviation of the spectrum in the positive region and calculating the deviation of the spectrum in the negative region by subtracting the average value in the frequency direction of the spectrum in the negative region from the spectrum in the negative region; Calculate the average value of the subregions for the deviation of, and calculate the average value of the subregions for the deviation of the spectrum in the negative region And means for calculating the sum value for each partial region by adding the average value of the partial region of the positive region spectrum deviation and the average value of the partial region of the negative region spectrum deviation for each partial region And a means for binarizing a component value corresponding to the partial area depending on whether or not the sum value for each partial area is greater than a binarization threshold, and adjacent to the binarization matrix A means for creating a matrix, a means for creating a stationary solution matrix by multiplying the adjacent matrix until it reaches a steady state, and a subregion connected to the partial region having the maximum sum value for each partial region from the stationary solution matrix and for each partial region Means for extracting a component corresponding to a concatenated region consisting of partial regions whose sum is greater than the binarization threshold and creating a filter base matrix from the binarization matrix, and for each filter transmission component row by row with respect to the filter base matrix Wide in frequency direction Means for creating a filter matrix by restricting the spectrum, means for calculating the spectral deviation sum by adding the spectral deviation of the positive region and the spectral deviation of the negative region for each beam direction distance / frequency, and Means for creating a filter application spectrum by applying a filter matrix to the sum of spectral deviations, means for detecting a peak frequency from the filter application spectrum, and means for calculating a beam direction flow velocity by multiplying the peak frequency and a conversion coefficient And have

Further, the surface layer flow velocity estimation program according to claim 3 , when detecting the peak frequency of the Doppler spectrum for each beam direction distance and frequency, and estimating the surface layer flow velocity, at least for each beam direction distance and frequency of the processing target region. Means for reading another Doppler spectrum data, means for separating the Doppler spectrum data into a spectrum having a positive frequency and a spectrum having a negative frequency, and an average value in the frequency direction for the spectrum in the positive region and the spectrum in the negative region Means for calculating the deviation of the spectrum of the positive region by subtracting the average value of the spectrum of the positive region from the spectrum of the positive region and calculating the frequency direction of the spectrum of the negative region from the spectrum of the negative region Means to calculate the deviation of the spectrum in the negative region by subtracting the average value of Means to calculate the average value of the partial area for deviation of the spectrum and to calculate the average value of the partial area for deviation of the spectrum of the negative area, the average value of the partial area of the deviation of the spectrum of the positive area and the spectrum of the negative area Means for summing the average value of the partial areas of the deviation for each partial area to calculate the total value for each partial area, whether the total value for each partial area is greater than the binarization threshold Means for binarizing corresponding component values to create a binarization matrix, means for creating an adjacency matrix for the binarization matrix, means for creating a stationary solution matrix by multiplying the adjacency matrix until it reaches a steady state, Binary by extracting components corresponding to a connected region consisting of partial regions that are connected to the partial region having the maximum combined value for each partial region and whose combined value is larger than the binarization threshold from the stationary solution matrix From the matrix Means for creating a filter fundamental matrix, means for creating a filter matrix by limiting the spread in the frequency direction of the filter transmission component for each row of the filter fundamental matrix, spectral deviation in the positive region and spectral deviation in the negative region Means for calculating the spectrum deviation sum value for each beam direction distance / frequency, means for applying the filter matrix to the spectrum deviation sum value to create a filter application spectrum, means for detecting the peak frequency from the filter application spectrum The computer is made to function as means for calculating the beam direction flow velocity by multiplying the peak frequency and the conversion coefficient.

  Therefore, according to the surface velocity estimation method, apparatus, and program, the error information verification process is performed to extract a region that satisfies certain conditions from the original Doppler spectrum data, and is detected from the region subjected to the verification process. Since the surface layer velocity is estimated using the measured peak frequency, it is possible to accurately remove the noise from the actual observation data that includes natural noise and exogenous noise such as artificial noise. Surface flow velocity based on flow information is estimated.

  In the present invention, the Doppler spectrum data is used in the sense of frequency Doppler spectrum value data including flow velocity information.

  According to the surface layer flow velocity estimation method, apparatus and program of the present invention, it is possible to appropriately remove noise from actual observation data including exogenous noise such as natural noise and artificial noise. The legitimate peak frequency of the Doppler spectrum can be detected, and the surface layer flow velocity is estimated based on the appropriate flow information, so that the estimation accuracy can be improved. Further, since accurate flow information can be detected, for example, complicated post-processing such as checking of the detection result of the peak frequency is not required, and it is possible to save the labor related to the estimation and to perform the estimation in a short time. It becomes possible to do in. Furthermore, according to the present invention, since flow information contributing to real-time sea area flow monitoring can be provided, information useful for coastal environment monitoring, port management and coastal disaster prevention targeting sudden events is provided. It becomes possible to do.

  Hereinafter, the configuration of the present invention will be described in detail based on the best mode shown in the drawings.

  FIG. 1 to FIG. 12 show an example of embodiments of the surface layer flow velocity estimation method, apparatus, and program of the present invention. This surface velocity estimation method is a method of estimating the surface velocity by detecting the peak frequency of the Doppler spectrum by beam direction distance and frequency, and in the frequency direction average of the spectrum of the positive region from the spectrum of the positive region of the Doppler spectrum. For the deviation calculated by subtracting the value, the average value of the partial region is calculated and the average of the partial region is calculated by subtracting the average value in the frequency direction of the negative region spectrum from the negative region spectrum of the Doppler spectrum. A value is calculated, connected to the partial region where the sum of the average value of the partial region of the positive region spectrum deviation and the average value of the partial region of the negative region spectral deviation is the maximum, and The peak frequency is detected from the concatenated area composed of the partial areas where the total value for each partial area is larger than the predetermined threshold. That.

  As shown in FIG. 1, the surface layer velocity estimation method includes a step (S1) of reading Doppler spectrum data for each beam direction distance and frequency of the processing target region, and the Doppler spectrum data in the positive frequency region. A step of separating the spectrum into a spectrum of a negative region (S2), a step of calculating an average value in the frequency direction for the spectrum of the positive region and the spectrum of the negative region (S3), Subtract the average value in the frequency direction of the spectrum of the negative region to calculate the deviation of the spectrum in the positive region and subtract the average value in the frequency direction of the spectrum of the negative region from the spectrum of the negative region. When calculating the deviation (S4) and calculating the average value of the partial areas for the deviation of the spectrum of the positive area A step (S5) of calculating an average value of the partial areas for the deviation of the spectrum of the negative area, and an average value of the partial areas of the deviation of the spectrum of the positive area and the average value of the partial areas of the deviation of the spectrum of the negative area Are calculated for each partial region to calculate a sum value for each partial region (S6), and whether the sum value for each partial region is larger than a binarization threshold A step of creating a binarization matrix by binarizing values (S7), a step of creating an adjacency matrix for the binarization matrix (S8), and multiplying the adjacency matrix until it reaches a steady state, Corresponding to a step (S9) to be created and a connected region consisting of partial regions connected to the partial region having the maximum sum value for each partial region from the steady solution matrix and having a sum value for each partial region larger than the binarization threshold Ingredients to Generating a filter basic matrix from the binarized matrix (S10), creating a filter matrix by limiting the spread of the filter transmission component in the frequency direction for each row of the filter basic matrix (S11), A step (S12) of adding a spectral deviation of the region and a spectral deviation of the negative region for each beam direction distance / frequency to calculate a combined spectral deviation value, and applying a filter matrix to the combined spectral deviation value. From the step of creating the filter application spectrum (S13), the step of detecting the peak frequency from the filter application spectrum (S14), and the step of calculating the beam direction flow velocity by multiplying the peak frequency and the conversion coefficient (S15). This is realized by the processing configuration.

  Moreover, the surface layer flow velocity estimation method is realized as the surface layer flow velocity estimation apparatus of the present invention. The surface velocity estimation apparatus of the present invention includes means for reading Doppler spectrum data for each beam direction distance and frequency of a processing target region, and separating the Doppler spectrum data into a spectrum in a positive region and a spectrum in a negative region. Means for calculating the average value in the frequency direction for the spectrum in the positive region and the spectrum in the negative region, and subtracting the average value in the frequency direction of the spectrum in the positive region from the spectrum in the positive region. Means for calculating the deviation of the spectrum of the negative region and subtracting the average value in the frequency direction of the spectrum of the negative region from the spectrum of the negative region to calculate the deviation of the spectrum of the negative region; Means for calculating an average value of the area and calculating an average value of the partial area with respect to the deviation of the spectrum of the negative area; Means for adding the average value of the partial area of the positive region spectrum deviation and the average value of the partial region of the negative region spectrum deviation for each partial region; A means for binarizing a component value corresponding to the partial region based on whether or not each sum is greater than a binarization threshold, and creating an adjacency matrix for the binarization matrix Means for multiplying the adjacency matrix until it reaches a steady state, creating a stationary solution matrix, and connecting the partial solution from the stationary solution matrix with the maximum sum value for each partial region, and the sum value for each partial region is Means for extracting a component corresponding to a concatenated region composed of partial regions larger than the binarization threshold and creating a filter base matrix from the binarization matrix, and spreading of the filter transmission component in the frequency direction for each row of the filter base matrix Limit Means for creating a filter matrix, means for summing the deviation of the spectrum in the positive area and the deviation of the spectrum in the negative area for each beam direction distance / frequency, and calculating the spectrum deviation sum, and the filter matrix as the spectrum Means for creating a filter application spectrum by applying to the sum of deviations, means for detecting a peak frequency from the filter application spectrum, and means for calculating the beam direction flow velocity by multiplying the peak frequency and the conversion coefficient.

  The surface layer flow velocity estimation method and the surface layer flow velocity estimation apparatus described above can also be realized by executing the surface layer flow velocity estimation program of the present invention on a computer. In the present embodiment, a case where the surface layer flow velocity estimation program is executed on a computer will be described as an example.

  FIG. 2 shows the overall configuration of the surface layer flow velocity estimation apparatus 10 of the present embodiment for executing the surface layer flow velocity estimation program 17. The surface layer flow velocity estimation apparatus 10 includes a control unit 11, a storage unit 12, an input unit 13, a display unit 14, and a memory 15, and is connected to each other by a signal line such as a bus. In addition, a data server 16 is connected to the surface layer flow velocity estimation apparatus 10 via a communication line or the like, and signals such as data and control commands are transmitted and received (input / output) through the communication line and the like.

  The control unit 11 performs control related to the overall control of the surface layer flow velocity estimation apparatus 10 and the estimation of the surface layer flow velocity by the surface layer flow velocity estimation program 17 stored in the storage unit 12, and is, for example, a CPU (Central Processing Unit). . The storage unit 12 is a device that can store at least data and programs, and is, for example, a hard disk. The memory 15 becomes a memory space that is a work area when the control unit 11 executes various controls and calculations.

  The input unit 13 is an interface for giving at least an operator's command to the control unit 11, and is, for example, a keyboard.

  The display unit 14 displays characters, graphics, and the like under the control of the control unit 11 and is, for example, a display.

  The data server 16 is a server capable of storing at least data.

  The control unit 11 of the surface layer flow velocity estimation apparatus 10 executes a surface layer flow velocity estimation program 17, thereby reading a Doppler spectrum data by a data reading unit 11a for reading the Doppler spectrum data for each beam direction distance and frequency of the processing target region. A spectrum separating unit 11b that separates the spectrum into a spectrum having a positive frequency and a spectrum having a negative frequency, and a spectrum average calculating unit 11c that calculates an average value in the frequency direction for the spectrum in the positive region and the spectrum in the negative region Subtract the average value in the frequency direction of the positive region spectrum from the spectrum in the positive region to calculate the deviation of the positive region spectrum and calculate the average value in the frequency direction of the negative region spectrum from the negative region spectrum. A spectral deviation calculation unit 11d that calculates a deviation of the spectrum of the negative region by pulling; A partial region deviation average value calculation unit 11e that calculates an average value of the partial region for the deviation of the spectrum of the positive region and calculates an average value of the partial region for the deviation of the spectrum of the negative region, and a deviation of the spectrum of the positive region A partial region deviation average value summation unit 11f that calculates a partial region average value and a partial region average value of the negative region spectral deviation for each partial region, and a partial region A binarization matrix creating unit 11g for binarizing a component value corresponding to the partial region to create a binarization matrix depending on whether or not each sum is greater than a binarization threshold; An adjacency matrix creation unit 11h that creates an adjacency matrix, a stationary solution matrix creation unit 11i that creates a steady solution matrix by multiplying the adjacency matrix until it reaches a steady state, and a sum of values for each partial region from the steady solution matrix is maximum. Create a filter base matrix that extracts a component corresponding to a connected area consisting of partial areas that are connected to a partial area and whose combined value is larger than the binarization threshold value. Part 11j, filter matrix creation part 11k for creating a filter matrix by limiting the spread of the filter transmission component in the frequency direction for each row of the filter basic matrix, and the spectral deviation of the positive region and the spectral deviation of the negative region For each beam direction distance and frequency, a spectrum deviation summation unit 11l for calculating a spectrum deviation sum value, a filter application unit 11m for creating a filter application spectrum by applying a filter matrix to the spectrum deviation sum value, A peak frequency detector 11n for detecting a peak frequency from the filter application spectrum, a peak frequency and a conversion coefficient; And a flow velocity calculation unit 11o that calculates the beam direction flow velocity.

  First, the process (S1-S5) which produces partial region spectrum data using the original Doppler spectrum data as a 1st step is demonstrated.

  In executing the surface layer flow velocity estimation method of the present invention, first, the data reading unit 11a of the control unit 11 reads Doppler spectrum data of the processing target region (S1).

  In the present invention, Doppler spectrum data of a backscattered wave from a marine wave, which is a radio wave (that is, an echo) received by the marine radar receiving antenna, is used.

  The radar for collecting data used in the present invention may be of any type as long as it can obtain the Doppler spectrum obtained in the distance direction and the frequency direction. Specifically, for example, the radar transmission / reception method is, for example, an FMCW (abbreviation of frequency modulated continuous wave) method or the FMICW (abbreviation of frequency modulated and interrupted continuous wave) method, and the antenna method is, for example, a mechanical rotation type or DBF (digital). (Beam Forming) or a phased array radar can be used. The center frequency is 24.515 MHz or 41.9 MHz, the frequency sweep width is about 300 kHz, the transmission output is about 50 W to 100 W, the beam width is about 13 degrees to 20 degrees, and the distance resolution is about 0.5 km, for example. Can be used. The method of obtaining the Doppler spectrum from the echo itself is a well-known technique, so details are omitted here (for example, Noro Nozaki and Toshihiko Umehara: Observation of ocean currents and waves using a short-wave ocean radar using the 24 MHz band (Part 1) ), HAM Journal, 63, pp. 77-85, 1989; Nozaki Norio and Umehara Toshihiko: Observation of ocean currents and waves with a short-wave ocean radar using the 24 MHz band (Part 2), HAM Journal, 64, pp. 103- 108, 1989; Japan Society of Civil Engineers, Coastal Engineering Committee, Research Status Review Subcommittee: Coastal Ocean Observation with Land-based Radar, 2001; Satoshi Fujii: Technology and History of Marine Radar, Coastal Ocean Research, 41, pp.97 -108, 2004).

  In the present embodiment, description will be made using Doppler spectra of echoes actually received by two radar stations installed on the ocean coast of the observation target region so that the same region can be observed from different directions. To do. A Doppler spectrum used in this embodiment is shown in FIG. In FIG. 3, the horizontal axis represents the frequency (unit: Hz), the vertical axis represents the distance in the beam direction (unit: km), and the contour value represents the spectrum intensity (unit: dB × 100). The radar station has at least a transmission antenna, a reception antenna, and a control device for controlling both antennas.

  In the present embodiment, Doppler spectrum data of up to 32 km (that is, 64 ranges) with a beam direction distance of 500 m is collected for 8 directions at intervals of 8.2 degrees in one observation of each radar station. As shown in FIG. 3, the Doppler spectrum used in this embodiment has noise around the beam direction distance of 8 km, and has a plurality of peaks even at a position where the spectrum is large as seen around the beam direction distance of 13 km. It is what you have.

  In the present embodiment, the Doppler spectrum data is stored in advance in the data server 16 as the spectrum database 18. The Doppler spectrum data is stored in the spectrum database 18 for each observation region, each observation date, each radar, each beam direction, each beam direction distance, and each frequency.

  Here, although not explicitly shown in the description of each process, all the subsequent processes such as various calculations, storage in the memory 15 and reading from the memory 15 are performed for each radar and each beam direction.

  The data reading unit 11 a reads the Doppler spectrum data for each beam direction distance and frequency for the processing target date and time of the processing target region from the spectrum database 18 and stores them in the memory 15. The Doppler spectrum used in this embodiment shown in FIG. 3 is obtained by displaying the Doppler spectrum of echoes of a specific beam direction received by a specific radar station in a certain processing target region / date and time by frequency and beam direction distance. It is.

  Next, the spectrum separation unit 11b of the control unit 11 separates the Doppler spectrum data read in the process of S1 into spectrum data in the positive frequency domain and spectrum data in the negative frequency domain (S2).

  Specifically, the spectrum separation unit 11b reads from the memory 15 Doppler spectrum data for each beam direction distance and frequency at the processing target date and time of the processing target region stored in the memory 15 in the processing of S1, and the frequency is positive. The spectrum is divided into the spectrum Ai, j of the region and the spectrum Bi, j of the negative region. Here, A and B are spectrum values (unit is dB), i is a beam direction distance number, and j is a frequency number.

  In this embodiment, since the width of the beam direction distance is 500 m, a number 1 is assigned to a distance 0.5 km, a number 2 is assigned to a distance 1 km, and a serial number is assigned in order every 500 m thereafter, and a number 64 is assigned to a distance 32 km. (For example, Shinichi Sakai et al .: Development of a high-performance coastal ocean radar for wide-area flow observation, Research Report U02056, Central Research Laboratory of Electric Power, 2003). Therefore, in this embodiment, the beam direction distance number i is a natural number from 1 to 64.

  In the present embodiment, the center frequency is 41.9 MHz, and a frequency band of plus or minus 1.755375 Hz is received across this center frequency. Here, in the description of the present embodiment, the center frequency is described as 0 Hz unless otherwise specified. The frequency pitch is 0.00596 Hz. Furthermore, in the present embodiment, processing is performed except between −0.09536 Hz and 0.09536 Hz. In this embodiment, for the positive region, number 1 is assigned to the frequency 0.10132 Hz, number 2 is assigned to the frequency 0.10728 Hz, and thereafter, consecutive numbers are assigned in order every 0.00596 Hz, and the frequency is 1.22776 Hz. Is given the number 191. In addition, for the negative region, number 1 is assigned to the frequency -1.22776 Hz, number 2 is assigned to the frequency -1.22180 Hz, and consecutive numbers are assigned sequentially every 0.00596 Hz in the same manner, and numbers are assigned to the frequency -0.1032 Hz. 191 is given. Therefore, in this embodiment, the frequency number j is a natural number from 1 to 191.

  In the Doppler spectrum data, the data in the positive frequency region is the spectrum data of the scattering response from the sea surface wave approaching the radar station, and the data in the negative frequency region is the scattering from the sea surface wave moving away from the radar station. It is spectral data of response.

  Then, the spectrum separation unit 11b includes a positive region spectrum Ai, j and a negative region spectrum Bi, j (where i = 1, 2, 3,... 64, and j = 1, 2, 3,. 191) is stored in the memory 15.

  Next, the spectrum average calculation unit 11c of the control unit 11 calculates an average value in the frequency direction for the positive region spectrum and the negative region spectrum set in the process of S2 (S3).

Specifically, the spectrum average calculation unit 11c reads the positive region spectrum Ai, j and the negative region spectrum Bi, j stored in the memory 15 in the process of S2 from the memory 15, and determines the beam direction for each positive and negative region. The average value A ^- i of the positive region spectrum A i, j, which is the average value of the spectrum for each distance number i (that is, the average value in the frequency direction), and the average value B ^- i, of the negative region spectrum Bi, j・ Calculate

Then, the spectrum average calculation unit 11c stores the positive region spectrum average value A ^- i,... And the negative region spectrum average value B ^- i, (where i = 1, 2, 3,..., 64) in the memory 15. Remember me.

  Next, the spectrum deviation calculating unit 11d of the control unit 11 includes the positive region spectrum and the negative region spectrum set in the process of S2, and the positive region spectrum average value and the negative region spectrum average calculated in the process of S3. The deviation from the average value is calculated using the value (S4).

Specifically, the spectrum deviation calculation unit 11d reads the positive region spectrum Ai, j and the negative region spectrum Bi, j stored in the memory 15 in the process of S2 from the memory 15, and also stores the memory 15 in the process of S3. stored positive region average spectral value a in ^- i, - and negative area average spectral value B ^- i, a-read from the memory 15, positive region spectrum deviation a ¹ i according to equation 1a, calculates the j The negative region spectral deviation B ¹ i, j is calculated by Equation 1b.

(Number ^{1a) A 1 i, j =} Ai, j-Ai, ·
(Number ^{1b) B 1 i, j =} Bi, j-Bi, ·

  The process of S4 is to calculate the difference between the value for each beam direction distance number i and the frequency number j for the Doppler spectrum and the average value for each beam direction distance number i. Since it is often negative, it is intended to make the peak frequency and the spectrum value in the surrounding area positive in consideration of the ease of subsequent processing.

FIG. 4 shows a positive region spectral deviation A ¹ i, j and a negative region spectral deviation B ¹ i calculated by performing the processing from S2 to S4 on the Doppler spectral data in the present embodiment shown in FIG. , j. FIG. 4A shows data in a negative frequency region, and is spectral data of a scattering response from a sea wave moving away from the radar station. FIG. 4B is data in a positive frequency region, and is spectral data of a scattering response from a sea surface wave approaching the radar station.

  In FIG. 4, the vertical axis represents the beam direction distance (km) as in FIG. 3, and the horizontal axis represents the beam direction flow velocity (unit: cm / s) calculated by multiplying the frequency and the conversion coefficient. Note that the method of calculating the beam velocity by multiplying the frequency and the conversion coefficient is a well-known technique, so details are omitted here (for example, Shinichi Sakai et al .: High-performance coastal ocean for wide-area flow observation) Radar development, research report U02056, Central Research Laboratory of Electric Power, 2003).

Then, the spectral deviation calculation unit 11d performs positive domain spectral deviation A ¹ i, j and negative domain spectral deviation B ¹ i, j (where i = 1, 2, 3,..., 64, and j = 1. , 2, 3,..., 191) are stored in the memory 15.

  Next, the partial region deviation average value calculation unit 11e of the control unit 11 uses the positive region spectral deviation and the negative region spectral deviation calculated in the process of S4, and the average value of the partial region of the positive region spectral deviation and The average value of the partial regions having the negative region spectral deviation is calculated (S5).

Specifically, the partial region deviation average value calculation unit 11e stores the positive region spectral deviation A ¹ i, j and the negative region spectral deviation B ¹ i, j stored in the memory 15 in the process of S4 from the memory 15. Read.

Then, according to Equations 2 and 3, the beam direction distance of the original Doppler spectrum data is bundled in order from the smallest one, and the distance width is wider than the original distance width and the frequency of the original Doppler spectrum data is smaller. The partial area average value A ² k, l of the positive area spectral deviation A ¹ i, j and the negative area spectral deviation B for each partial area having a frequency width wider than the original frequency pitch. ¹ The partial area average value B ² k, l of i, j is calculated. At this time, the distance range numbers k (= 1, 2, 3,...) Of the partial areas are sequentially assigned in ascending order of the distance with respect to the distance range of the partial area formed by collectively combining some of the beam direction distances of the original data. And the frequency range number l (= 1, 2, 3,...) Of the partial region in order from the one with the smallest original frequency number with respect to the frequency range of the partial region created by enclosing some of the frequencies of the original data. ).

  In Equation 2 and Equation 3, ist = im × (k−1) +1 and jst = jm × (l−1) +1. im is the number of bundles in the beam direction distance of the original Doppler spectrum data when forming a partial region of the Doppler spectrum data, and jm is the number of bundles of frequencies in the original Doppler spectrum data.

The calculation of the partial region average values A ² k, l and B ² k, l of the positive region and the negative region spectral deviation is to reduce the resolution of the original Doppler spectral data and smooth the change of the spectral value. This is the process to be performed. Therefore, the number im of beam direction distances and the number jm of frequency bundles are set as natural numbers of 2 or more in principle. Specifically, the numbers im and jm of the binding are the original Doppler because the width of the beam direction distance and the width of the frequency (that is, the flow velocity) of the partial region determined by these are resolutions when the peak frequency is captured. It is set as appropriate in consideration of the estimation accuracy required based on the purpose of the flow velocity estimation together with the beam direction distance pitch and the frequency pitch of the spectrum data. For example, the width of the distance in the beam direction of the partial region may be about 1 to 4 km, and the width of the flow velocity obtained by converting the frequency may be about 10 to 30 cm / s. However, the width of the distance in the beam direction and the width of the flow velocity in the partial region are not limited to these, and both may be smaller or larger. For example, when high estimation accuracy is required or when the pitch of the beam direction distance of the original Doppler spectrum data is large, the beam direction distance is not bundled, and the number of beam direction distance bundles im = 1. Also good.

  In this embodiment, when the beam pitch distance of the original data is 500 m, the number of distance links im = 4 is set, and when the frequency pitch of the original data is 0.00596 Hz, the number of frequency links jm = 7 is set. Thereby, in the present embodiment, the width of the distance in the beam direction of the partial region is 2 km, and the width of the flow velocity converted from the frequency is about 14 cm / s.

  In this embodiment, since the beam direction distance number i = 1, 2, 3,..., 64, the distance range number k of the partial region is from 1 to 16 when the number of distance links im = 4. In the present embodiment, jst in Equation 2 and Equation 3 is jst = jm × (l−1) +2 and the frequency number j = 1, 2, 3,. Assuming 7, the frequency range number l of the partial region is 1 to 27.

FIG. 5 shows a partial region of the positive region spectral deviation calculated by performing the process of S5 on the positive region spectral deviation A ¹ i, j and the negative region spectral deviation B ¹ i, j shown in FIG. The average value A ² k, l (FIG. 5B) and the partial region average value B ² k, l of the negative region spectral deviation (FIG. 5A). The horizontal and vertical axes in FIG. 5 are the same as those in FIG.

  As shown in FIG. 5, by performing the process of S5, the change in value becomes smooth compared with the original spectrum data, and the spectrum noise is leveled.

  The values of the beam direction distance bundle im and the frequency bundle jm are preliminarily defined on the surface layer flow velocity estimation program 17, and these values are calculated at the stage of executing the process of S5. The unit 11e may read from the program, or the operator specifies the values of the number im of the beam direction distance and the number jm of the frequency band, and the specified values are the partial region deviation average value calculation unit 11e. May be read.

  When the operator designates values of the beam direction distance bundle number im and the frequency bundle number jm, the partial area deviation average value calculation unit 11e performs the beam direction distance bundle at the stage of executing the process of S5. A message requesting the specification of the value of the number im and the frequency jm of the frequency is displayed on the display unit 14 and the specified value of the operator input via the input unit 13 is read.

Then, the partial region deviation average value calculation unit 11e generates a partial region average value A ² k, l of the positive region spectral deviation and a partial region average value B ² k, l of the negative region spectral deviation (where k = 1, 2, 3, ..., 16 and l = 1, 2, 3, ..., 27) are stored in the memory 15.

  Next, the process (S6-S11) which produces a filter matrix using partial region spectrum data as a 2nd step is demonstrated.

  As a filter for detecting the peak frequency from the Doppler spectrum data developed for the frequency and the beam direction distance, a region including a peak in the Doppler spectrum data and connected to the peak and having a spectrum larger than a predetermined threshold (concatenation) (Referred to as a region) is a feature of the present invention, and the second stage process is a process for creating this filter.

  In the present embodiment, an adjacency matrix is used as a method of extracting a connected region including a peak and including a spectrum larger than a predetermined threshold from Doppler spectrum data.

  The partial region deviation average value summation unit 11f of the control unit 11 sums up the partial region average value of the positive region spectral deviation and the partial region average value of the negative region spectral deviation calculated in the process of S5 for each partial region. (S6).

Specifically, the partial region deviation average value summing unit 11f performs the partial region average value A ² k, l of the positive region spectral deviation and the partial region average of the negative region spectral deviation stored in the memory 15 in the process of S5. The value B ² k, l is read from the memory 15, and the partial region average value A ² k, l of the positive region spectral deviation and the partial region average value B ² k, l of the negative region spectral deviation for each partial region according to Equation 4. Is added to calculate the total value C ² k, l of the partial area average value of the spectral deviation between the positive area and the negative area (hereinafter referred to as the partial area deviation average combined value C ² k, l as appropriate). To do.

(Equation 4) C ² k, l = A ² k, l + B ² k, l

FIG. 6 shows the sum of the partial region deviation averages obtained by adding the partial region average value A ² k, l of the positive region spectral deviation and the partial region average value B ² k, l of the negative region spectral deviation shown in FIG. The value C ² k, l. The horizontal and vertical axes in FIG. 6 are the same as those in FIG.

  As shown in FIG. 6, by performing the process of S6, the degree of protrusion around the peak frequency with respect to other regions becomes more prominent.

Then, the partial area deviation average summation unit 11f is configured to generate the partial area deviation average summation value C ² k, l (where k = 1, 2, 3,..., 16 and l = 1, 2, 3,... 27) is stored in the memory 15.

  Next, the binarization matrix creating unit 11g of the control unit 11 binarizes the component value corresponding to the partial region for each partial region using the partial region deviation average sum calculated in the process of S6. A binarization matrix is created (S7).

Specifically, the binarization matrix creation unit 11g reads the partial region deviation average sum C ² k, l stored in the memory 15 in the process of S6 from the memory 15, and subregion deviation average sum for each partial region. When the value C ² k, l is greater than the binarization threshold Tb, the value of the component fk, l corresponding to the partial area is set to 1, and when the value C ² k, l is equal to or less than the binarization threshold Tb, the partial area is handled. The binarization matrix F (component fk, l) is created by setting the value of the component fk, l to be set to zero.

  The binarization for each partial area is a process performed to extract spectrum information including flow velocity information in Doppler spectrum data obtained via a radar station. Therefore, the binarization threshold value Tb is set to a value that can eliminate the noise-affected portion that does not include the flow velocity information in the original Doppler spectrum data. Specifically, electromagnetic waves are attenuated in the distance band near the limit of the beam direction distance, and it is considered that the influence of noise is superior to the flow velocity information. Therefore, the maximum value of the spectrum in the distance band is set to the binarization threshold. It can be considered to be Tb.

In this embodiment, the maximum value of the partial area deviation average sum value C ² k, l in all frequency directions at a position in the beam direction distance 30 km where the spectrum is considered to be attenuated in the distance direction is defined as the binarization threshold Tb. Set.

Further, the binarization matrix creation unit 11g further determines a partial region where the partial region deviation average sum C ² k, l is maximum, and sets the distance range number k of the partial region to the maximum partial region distance range number k. In addition, the frequency range number l of the partial region is set as the maximum partial region frequency range number l ′.

FIG. 7 is a binarization matrix F composed of components fk, l corresponding to the partial regions created by performing the processing of S7 on the partial region deviation average sum C ² k, l shown in FIG. . In FIG. 7, the horizontal axis is the frequency range number l, and the vertical axis is the distance range number k. The squares in the figure represent partial areas corresponding to combinations of the frequency range number l and the distance range number k.

The gray squares indicate that the partial area deviation average sum C ² k, l of the corresponding partial area is equal to or less than the binarization threshold Tb and the component fk, l = 0 of the binarization matrix F. The white squares indicate that the sum C ² k, l of the corresponding partial area is larger than the binarization threshold Tb and the component fk, l = 1. Further, the black squares indicate a partial region where the total value C ² k, l is the maximum. Therefore, the distance range number k of the partial area corresponding to the black cell becomes the maximum partial area distance range number k ′, and the frequency range number l becomes the maximum partial area frequency range number l ′.

  In the present embodiment, as shown in FIG. 7, the maximum partial area distance range number k ′ = 3 and the maximum partial area frequency range number l ′ = 22.

  Then, the binarization matrix creation unit 11g has components fk, l (where k = 1, 2, 3,..., 16 and l = 1, 2, 3,..., 27) of the binarization matrix F. And the maximum partial region distance range number k ′ and the maximum partial region frequency range number l ′ are stored in the memory 15.

  Next, the adjacency matrix creating unit 11h of the control unit 11 creates an adjacency matrix using the binarized matrix created in the process of S7 (S8).

  Specifically, the adjacency matrix creation unit 11h first assigns the partial region numbers m to the partial regions in the order of the frequency range number l in the order of the distance range number k. For example, the number m = 1 is assigned to the partial region where the distance range number k = 1 and the frequency range number l = 1, and the number m is assigned to the partial region where the distance range number k = 1 and the frequency range number l = 2. = 2.

  In this embodiment, since the maximum value of the frequency range number l is 27, the number m = 28 is assigned to the partial area where the distance range number k = 2 and the frequency range number l = 1, and the distance range number The number m = 406 is assigned to the partial region where the frequency range number l = 1 and k = 16 (which is the maximum value of the distance range number k of the present embodiment).

  The adjacency matrix creation unit 11h also reads the maximum partial area distance range number k ′ and the maximum partial area frequency range number l ′ stored in the memory 15 in the process of S7 from the memory 15, and the distance range number k is the maximum partial area. The number m of the partial area that is the distance range number k ′ and the frequency range number l is the maximum partial area frequency range number l ′ is the maximum partial area number m ′.

  Then, the adjacency matrix creation unit 11h generates a partial area number m corresponding to a combination of the distance range number k and the frequency range number l (where k = 1, 2, 3,..., 16, l = 1, 2, 3). ,..., 27 and m = 1, 2, 3,..., 432) and the maximum partial area number m ′ are stored in the memory 15.

  Next, the adjacency matrix creation unit 11h includes the binarization matrix F stored in the memory 15 in the process of S7, and the combination of the distance range number k and the frequency range number l stored in the memory 15 in the above process. The correspondence relationship with the partial area number m is read from the memory 15.

  The adjacency matrix creation unit 11h first sets all the components gm, m of the adjacency matrix G to zero. Then, when the value of the component fk, l of the binarization matrix F of the partial region is 1, the adjacency matrix creating unit 11h corresponds to the partial region corresponding to the combination of the distance range number k and the frequency range number 1 of the partial region. Let component gm, m = 1 of adjacency matrix G for number m.

  Also, the adjacency matrix creation unit 11h performs vertical, horizontal, and left-right movements with a distance range number of k ± 1 or a frequency range number of l ± 1 for a partial region where the value of the component fk, l of the binarized matrix F is 1. The values of the components fk ± 1, l and fk, l ± 1 of the binarization matrix F of the partial region of the position are confirmed.

  Then, when the value of the component fk ± 1, l or fk, l ± 1 of the binarized matrix F is 1, the adjacency matrix creation unit 11h sets the partial region number m and the distance range number k ± 1 as The component gm, n = 1 of the adjacency matrix G for the number n corresponding to the combination with the frequency range number m or the combination of the distance range number m and the frequency range number l ± 1.

  As described above, the adjacency matrix creation unit 11h sets the value of the component gm, n of the adjacency matrix G. The method for setting the number n is the same as the method for setting the partial area number m.

  FIG. 8 is an adjacency matrix G created by performing the process of S8 on the binarized matrix F shown in FIG. As in FIG. 7, the horizontal axis in FIG. 8 is the frequency range number l, and the vertical axis is the distance range number k. When partial areas are adjacent vertically and horizontally, that is, when the component gm, n = 1 of the adjacency matrix G, bar lines are displayed between the numbers in the cells.

  Then, the adjacency matrix creation unit 11h stores the values of the components gm, n (where m = 1, 2, 3,..., 432 and n = 1, 2, 3,..., 432) of the adjacency matrix G. 15 is stored.

  Next, the steady solution matrix creation unit 11i of the control unit 11 creates a steady solution matrix using the adjacency matrix created in the process of S8 (S9).

Specifically, the stationary solution matrix creation unit 11i reads the adjacency matrix G stored in the memory 15 in the process of S8 from the memory 15, and each component g ^p m of the matrix G ^{p obtained} by multiplying the adjacency matrix G by p times. , n and 1 are compared, and a small value is adopted as a new component g ^pm , n (Equation 5).

(Equation 5) g ^p m, n = min (1, matrix G ^p = G ^p-1 component g ^p m, n of G)

The steady-state solution matrix creating unit 11i performs a process of multiplying the adjacency matrix G to matrix G ^p which an adjacency matrix G and multiply p is equal to the matrix G ^p-1, is a G ^p = G ^{p-1 Let the} matrix G ^p be a stationary solution matrix Gs.

Then, the stationary solution matrix creation unit 11i has components g ^p m, n (where m = 1, 2, 3,..., 432 and n = 1, 2, 3,..., 432) of the stationary solution matrix Gs. Is stored in the memory 15.

  Next, the filter basic matrix creating unit 11j of the control unit 11 creates a filter basic matrix using the binarized matrix created in the process of S7 and the steady solution matrix created in the process of S9 (S10). .

Specifically, the filter basic matrix creation unit 11j reads from the memory 15 the binarized matrix F stored in the memory 15 in the process of S7 and the steady solution matrix Gs stored in the memory 15 in the process of S9. Then, the filter basic matrix creation unit 11j includes the distance range number k and the frequency range number corresponding to the column number n in which the component g ^p m ′, n is 1 with the maximum partial region number m ′ in the steady solution matrix Gs. A filter basic matrix H (component hk, l) is created by setting all components fk, l of the binarized matrix F excluding l to zero.

FIG. 9 is a filter base matrix H created by performing the processes of S9 and S10 on the adjacency matrix G shown in FIG. As in FIG. 7, the horizontal axis in FIG. 9 is the frequency range number l and the vertical axis is the distance range number k, and the squares in the figure represent partial areas. A gray cell represents that the value of the component hk, l corresponding to the partial region is zero, and a white cell has a value of 1 for the component hk, l corresponding to the partial region (hereinafter, , Called a filter transmission component). Further, the black squares indicate that the partial area has the maximum partial area deviation average combined value C ² k, l.

As shown in FIG. 9, by performing the processes of S9 and S10, the partial region where the partial region deviation average combined value C ² k, l is maximum with respect to the binarization matrix F shown in FIG. Only the partial area connected to the black square) (ie, the connected area, white square in the figure) remains with the component hk, l having a value of 1 (ie, as a filter transmission component), while the sum C ² k, The value of the component hk, l of the partial region not connected to the partial region where l is maximum becomes zero, and the noise of the spectrum is removed.

  The filter basic matrix creation unit 11j then calculates the values of the components hk, l (where k = 1, 2, 3,..., 16 and l = 1, 2, 3,..., 27) of the filter basic matrix H. Is stored in the memory 15.

  Next, the filter matrix creation unit 11k of the control unit 11 creates a filter matrix using the filter basic matrix created in the process of S10 (S11).

Specifically, the filter matrix creation unit 11k reads the filter basic matrix H stored in the memory 15 in the process of S10 from the memory 15. Then, the filter matrix creation unit 11k has a partial region for a row k in which the component having a value of 1 in the filter basic matrix H (that is, the filter transmission component) hk, l is connected more than the maximum transmission component frequency width lmax. A filter matrix X (component xk, l) is created by setting the component values excluding the left and right (lmax-1) / 2 components from the position of the partial region where the deviation average sum C ² k, l is maximum to zero.

  The transmission component maximum frequency width lmax defines the width in the frequency direction when extracting the region including the peak frequency from the original Doppler spectrum data, that is, the extent of the frequency direction of the region corresponding to the filter transmission component, It is defined as the number of partial areas. The value of the transmission component maximum frequency width lmax is not limited to a specific value, and is appropriately set in consideration of the width of the flow velocity when the frequency is converted into the flow velocity, for example.

  In the present embodiment, the objective is to obtain a filter having a spread of about 50 cm / s on both sides of the peak frequency when converted to a flow rate, and the maximum transmission component is about 100 cm / s when converted to a flow rate. The frequency width lmax = 7 is set. As in this embodiment, when the transmission component maximum frequency width lmax is defined as the entire width centered on the peak, the value of lmax is set to an odd number.

FIG. 10 is a filter matrix X created by performing the process of S11 on the filter basic matrix H shown in FIG. As in FIG. 7, the horizontal axis in FIG. 10 is the frequency range number l, the vertical axis is the distance range number k, and the squares in the figure represent partial areas. A gray cell indicates that the value of the component xk, l corresponding to the partial region is zero, and a white cell indicates that the value of the component xk, l corresponding to the partial region is 1. To express. Further, the black squares indicate that the partial area has the maximum partial area deviation average combined value C ² k, l.

As shown in FIG. 10, by performing the process of S11, the partial region (the white cell in the figure) in which the value of the component xk, l is 1 compared to the filter basic matrix H shown in FIG. It is limited to a certain range from the partial region (black square in the figure) where the deviation average sum C ² k, l is the maximum, and noise in the spectrum is further removed.

  The filter matrix creation unit 11k stores the values of the components xk, l (where k = 1, 2, 3,..., 16 and l = 1, 2, 3,..., 27) of the filter matrix X in the memory. 15 is stored.

  Then, the process (S12-S14) which detects a peak frequency using a filter matrix as a 3rd step, and the process (S15) which estimates a flow velocity are demonstrated.

  The spectral deviation summation unit 11l of the control unit 11 sums the positive region spectral deviation and the negative region spectral deviation calculated in the process of S4 for each beam direction distance and for each frequency (S12).

Specifically, the spectral deviation summation unit 11l reads from the memory 15 the positive area spectral deviation A ¹ i, j and the negative area spectral deviation B ¹ i, j stored in the memory 15 in the process of S4, The total spectral deviation value C ¹ i, j is calculated by adding the positive region spectral deviation A ¹ i, j and the negative region spectral deviation B ¹ i, j.

Figure 11 is a positive region spectrum deviation A ¹ i, j and negative regions spectrum deviation B ¹ i, j and the sum combined spectral deviation total value C ¹ i, j shown in FIG. The horizontal and vertical axes in FIG. 11 are the same as those in FIG.

Then, the spectrum deviation summation unit 11l stores the spectrum deviation summation value C ¹ i, j in the memory 15.

  Next, the filter application unit 11m of the control unit 11 applies the filter matrix created in the process of S11 to the spectrum deviation sum calculated in the process of S12 (S13).

Specifically, the filter application unit 11m reads the filter matrix X stored in the memory 15 in the process of S11 from the memory 15, and uses the spectrum deviation sum C ¹ i, j stored in the memory 15 in the process of S12. The filter matrix X is read from the memory 15 and applied to the spectrum deviation sum value C ¹ i, j.

That is, the filter application unit 11m performs spectral deviation only for the spectrum of the beam direction distance number i and frequency number j corresponding to the distance range number k and frequency range number l of the component xk, l whose value is 1 in the filter matrix X. A filter-applied spectrum Di, j having a sum value C ¹ i, j is created. In this embodiment, for i = 1 and 191 of the filter application spectrum Di, j, the values of the components xk, l of k = 1 and 27 of the filter matrix X are used.

FIG. 12 shows a filter application spectrum Di, j obtained as a result of applying the filter matrix X shown in FIG. 10 to the spectrum deviation sum C ¹ i, j shown in FIG. The horizontal and vertical axes in FIG. 12 are the same as those in FIG.

  As shown in FIG. 12, by performing the process of S13, the extrinsic noise such as the intensity of the spectrum becoming extremely strong or the intensity varies locally is removed, and a legitimate peak is included from the original spectrum data. Only regions with contiguous spectra remain.

  The filter application unit 11m stores the filter application spectrum Di, j in the memory 15.

  Next, the peak frequency detection unit 11n of the control unit 11 detects the peak frequency using the filter application spectrum created in the process of S13 (S14).

Specifically, the peak frequency detector 11n is filter application spectrum Di stored in the memory 15 in the processing of S13, reads a j, filter application spectrum Di, moment D in the frequency direction j ^m i, · and the average value D ^- i, · is calculated for each beam direction distance number i, and the center of gravity in the frequency direction is calculated using these. Note that the calculation of the center of gravity itself is a well-known technique as the method of moments, so details are omitted here (for example, Hiroyuki Kameda et al .: Probability / statistical analysis, Gihodo Publishing, 1981).

  Then, the peak frequency detector 11n stores the center of gravity for each beam direction distance number i in the memory 15 as the peak frequency.

  Next, the flow velocity calculation unit 11o of the control unit 11 calculates the beam direction flow velocity using the peak frequency detected in the process of S14 (S15).

  Specifically, the flow velocity calculation unit 11o reads the peak frequency for each beam direction distance number i stored in the memory 15 in the process of S14, and multiplies this peak frequency by the conversion coefficient for each beam direction distance number i. The beam direction velocity is calculated. The method of calculating the beam direction velocity by multiplying the frequency and the conversion coefficient is a well-known technique, so details are omitted here (for example, Shinichi Sakai et al .: High-performance coastal ocean radar for wide-area flow observation) Development, Research Report U02056, Central Research Laboratory of Electric Power, 2003). The conversion coefficient is defined in advance on the surface layer flow velocity estimation program 17.

  Then, the flow velocity calculation unit 11o stores the beam direction flow velocity for each beam direction distance number i in the memory 15.

  Through the above processing, the beam direction flow velocity, that is, the line-of-sight flow velocity is estimated for each beam direction distance number for each radar and each beam direction.

  The surface layer flow velocity estimation apparatus 10 accumulates the estimated flow velocity as estimation result data in the data server 16 or superimposes it on the topographic map of the processing target area according to the purpose of the flow velocity estimation and the purpose of utilizing the estimation result. 14 or the like. At this time, since it is the beam direction flow velocity for each radar, that is, the line-of-sight flow velocity, which is estimated by the processing up to S15, the surface layer flow velocity estimation apparatus 10 performs vector synthesis of the line-of-sight flow velocity for each radar for each beam azimuth and beam direction distance. Then, the horizontal flow velocity for each point may be calculated.

  According to the surface velocity estimation method, apparatus, and program of the present invention having the above configuration, the error (noise) information verification process is superior to the conventional technology, and the flow information is accurately detected by the effect. can do. Therefore, the present invention can avoid a decrease in estimation accuracy and complexity of post-processing with respect to actual observation data including external noise such as natural noise and artificial noise. As a result, this method makes it possible to provide flow information that contributes to real-time ocean flow monitoring in a short time.

  In addition, although the above-mentioned form is an example of the suitable form of this invention, it is not limited to this, A various deformation | transformation implementation is possible in the range which does not deviate from the summary of this invention. For example, in the present embodiment, the Doppler spectrum data of echoes received by two radar stations installed on the ocean coast in the observation target region are spectrumd in the data server 16 connected to the surface layer flow velocity estimation device 10 by a communication line or the like. However, the present invention is not limited to this, and data from a radar station may be directly input to the surface layer flow velocity estimation apparatus 10 via a wireless communication line using PHS, for example. In this case, it is possible to estimate and provide information in almost real time.

  In the present embodiment, Doppler spectrum data of echoes received by two radar stations are aggregated into one as the spectrum database 18, and the surface layer flow velocity estimation apparatus 10 estimates the flow velocity using the data, and if necessary. However, the present invention is not limited to this, and the surface layer flow velocity estimation device 10 is installed in each radar station to estimate the flow velocity at each radar station and use flow velocity data. A device for providing information may be further provided, and the flow rate data estimated by the device may be transferred and aggregated via a wired / wireless communication line so that the device provides the information.

  In this embodiment, an adjacency matrix is used when creating a filter matrix for detecting a peak frequency (S8 to S11). However, in the Doppler spectrum data developed for the frequency and the beam direction distance, The method of creating the filter matrix is not limited to the method of this embodiment as long as it is a method that can extract a region that includes a peak and that is connected to the peak and has a spectrum that is larger than a predetermined threshold. Specifically, in the process of S8 of the present embodiment, the adjacency matrix G is created for the binarized matrix F shown in FIG. 7, but the present invention is not limited to this. For example, a labeling process for image processing is performed. May be used.

In this case, first, a white cell of the binarization matrix F shown in FIG. 7, that is, a connection region where partial regions whose partial region deviation average sum C ² k, l is larger than the binarization threshold Tb is connected. Extract by image processing. Subsequently, a different label for each extracted connected area is given to a partial area included in the connected area. Then, the partial region deviation average sum C ² k, l is assigned to the maximum partial region (that is, the partial region having the maximum partial region distance range number k ′ and the maximum partial region frequency range number 1 ′ in this embodiment). The filter matrix H is created by setting the component corresponding to the partial region of the same label as the label to 1 and setting the component corresponding to the partial region of a different label to 0.

  An embodiment in which the surface layer flow velocity estimation method, apparatus and program of the present invention are applied to the estimation of the flow velocity using the data of the actual sea area observation in the vicinity of the Osaka Bay Akashi Strait will be described with reference to FIGS.

  In this example, the sea area observation result obtained by using two digital beam forming radars in the Akashi Strait area of Osaka Bay was used. The installation points and beam orientations of the two radars are as shown in FIG. 13, and the symbols ● marked with St. A and St. B in the figure are emitted radially from each of the radar installation points and the symbols ●. The straight line is the beam direction. Sea area observation was carried out for about a month from January 18, 2003 to February 20, 2003.

  Furthermore, in order to compare with the estimation result of the flow velocity by the surface layer flow velocity estimation method, apparatus and program of the present invention, on January 31 and February 17 during the observation period, ADCP (Acoustic Doppler Current Profiler: Using the ultrasonic Doppler velocimeter, the sea area flow was directly observed. The ADCP used for observation was a 600 kHz work hose and broadband type, and recorded an average flow velocity of 15 pings (15 times) at 10-second intervals with a measurement layer thickness of 50 cm. The standard deviation on the specification was 2.3 cm.

  The flow velocity observation using ADCP was performed by installing ADCP on the side of three ships. The measurement points of flow rate using ADCP are as shown in FIG. The measurement time of flow velocity using ADCP is 12 stations and 3 tides from 7:00 to 13:00 on January 31st, and 15 stations and 2 tides from 13:30 to 16:30 on February 17th. Met.

  The flow velocity observation using ADCP is taken on line AA ′, BB ′, and CC ′ in FIG. 14A on January 31st, and DD in FIG. 14B on February 17th. The ', EE' and FF 'lines were moved by different ships, and the test was performed at the point where the beam azimuths of the two stations of St. A and St. B overlap (symbol ◯ in the figure). And in the measurement of one point, the flow velocity result which continued for several tens of minutes was recorded. In the comparison of the estimation result according to the present invention and the measurement result using ADCP in this example, the result at the water depth of about 2 m, which is the uppermost layer in the specifications of the device, is used for the flow velocity obtained by observation using ADCP. It was.

  The flow velocity was estimated by the present invention using actual spectrum data obtained by the radar, and the result shown in FIG. 15 was obtained. FIG. 15 shows the average flow velocity from January 18th to 31st. In FIG. 15, the direction of the arrow in the figure represents the direction of flow at the position of the arrow, and the length of the arrow represents the speed of flow at the position of the arrow.

The results shown in Figure 16 to create a scatter plot of both to compare match specific between the measured flow velocity ^{A V} using the estimated beam direction velocity ^{R V} and ADCP the present invention was obtained . In FIG. 16, the horizontal axis is the value of the estimation result according to the present invention, and the vertical axis is the value of the measurement result using ADCP. A straight line with a value = value on the horizontal axis (a broken line with y = x) is shown. Note that the flow velocity measured using ADCP is the result of converting into the gaze direction velocity of the radar. The estimation results according to the present invention and the measurement results using ADCP are obtained by combining the data of January 31 and February 17 using an average flow velocity of 15 minutes.

  FIG. 16 confirms that the estimation result according to the present invention and the measurement result using ADCP are distributed in the vicinity of y = x, and the estimation result according to the present invention shows the same result as the measurement result using ADCP. It was. Further, the correlation coefficient between the estimation result according to the present invention and the measurement result using ADCP was 0.98, indicating a very high correlation.

  From the result of this example, the estimation result according to the present invention is very close to the actual measurement result using ADCP, and it was confirmed that the estimation according to the present invention is very reliable. Therefore, according to the method of the present invention, accurate flow information can be detected, so that complicated post-processing such as checking of the detection result of the peak frequency is not required, and the labor for estimation can be saved. It was confirmed that the estimation can be performed in a short time.

  In addition, regarding the estimation according to the present invention in this example, the processing time for estimating the temporary point was about 30 seconds (PC used for the processing is 1 CPU: Pentium 4 2.4 GHz; Pentium is a registered trademark). From this result, the present invention can estimate the flow velocity very quickly, and is suitable for remote flow monitoring, for example, and can provide information on the flow in almost real time. It was confirmed that.

It is a flowchart explaining an example of embodiment of the surface layer flow velocity estimation method of this invention. It is a functional block diagram of the surface layer flow velocity estimation apparatus in the case of implementing the surface layer flow velocity estimation method of embodiment using a program. It is a figure which shows the Doppler spectrum of embodiment. It is a figure which shows the spectrum deviation of embodiment. (A) is a figure which shows a negative area | region spectrum deviation. (B) is a figure which shows a positive area | region spectrum deviation. It is a figure which shows the partial area average value of the spectrum deviation of embodiment. (A) is a figure which shows the partial area average value of a negative area | region spectral deviation. (B) is a figure which shows the partial area average value of a positive area | region spectral deviation. It is a figure which shows the total value which added together the partial area average value of the positive area | region spectral deviation of embodiment, and the partial area average value of the negative area | region spectral deviation. It is a figure which shows the binarization matrix corresponding to the partial area | region of embodiment. It is a figure which shows the adjacency matrix of embodiment. It is a figure which shows the filter basic matrix of embodiment. It is a figure which shows the filter matrix of embodiment. It is a figure which shows the spectrum deviation total value of embodiment. It is a figure which shows the filter application spectrum of embodiment. It is a figure which shows the installation point and beam direction of the radar of an Example. It is a figure which shows the measuring point of the flow velocity using ADCP of an Example. (A) is a figure which shows the measurement point in the flow velocity observation on January 31. FIG. (B) is a figure which shows the measurement point in the flow velocity observation on February 17. It is a figure which shows the estimation result of the flow velocity by this invention of an Example. It is a scatter diagram for the comparison with the flow velocity estimated by this invention of the Example, and the flow velocity measured using ADCP.

Explanation of symbols

DESCRIPTION OF SYMBOLS 10 Surface layer flow velocity estimation apparatus 11 Control part 12 Storage part 13 Input part 14 Display part 15 Memory 16 Data server 17 Surface layer flow rate estimation program 18 Spectrum database

Claims (3)

  Reading the Doppler spectrum data for each beam direction distance and frequency of the processing target region, separating the Doppler spectrum data into a spectrum of a positive region and a spectrum of a negative region, and the positive region And calculating a mean value in the frequency direction for the spectrum of the positive region and the spectrum of the negative region, and subtracting a mean value in the frequency direction of the spectrum of the positive region from the spectrum of the positive region to obtain a deviation of the spectrum of the positive region And calculating a deviation of the spectrum of the negative region by subtracting an average value in the frequency direction of the spectrum of the negative region from the spectrum of the negative region, and a partial region for the deviation of the spectrum of the positive region And calculating the average value of the subregion with respect to the spectral deviation of the negative region Calculating the average value, and adding the average value of the partial areas of the positive region spectrum deviation and the average value of the partial areas of the negative region spectrum deviation for each partial region. A step of calculating a summation value of each of the partial regions, and a step of creating a binarization matrix by binarizing the value of the component corresponding to the partial region depending on whether or not the summation value for each partial region is larger than a binarization threshold A step of creating an adjacency matrix for the binarization matrix, a step of creating a stationary solution matrix by multiplying the adjacency matrix until a steady state is obtained, and a sum value for each partial region from the stationary solution matrix A filter base matrix is created from the binarization matrix by extracting a component corresponding to a concatenation area consisting of partial areas that are connected to the largest partial area and whose sum value for each partial area is larger than the binarization threshold. Creating a filter matrix by limiting the spread of the filter transmission component in the frequency direction for each row of the filter basic matrix, and the spectral deviation of the positive region and the spectral deviation of the negative region, For each beam direction distance and for each frequency, to calculate a spectrum deviation sum, to apply the filter matrix to the spectrum deviation sum, create a filter application spectrum, and to peak from the filter application spectrum A method for estimating a surface flow velocity, comprising: detecting a frequency; and calculating a beam direction flow velocity by multiplying the peak frequency and a conversion coefficient.   Means for reading Doppler spectrum data for each beam direction distance and frequency of the region to be processed; means for separating the Doppler spectrum data into a spectrum of a positive region and a spectrum of a negative region; and the positive region And means for calculating an average value in the frequency direction for the spectrum of the positive region and the spectrum of the negative region, and subtracting the average value in the frequency direction of the spectrum of the positive region from the spectrum of the positive region, thereby deviating the spectrum of the positive region And calculating a deviation of the spectrum of the negative region by subtracting an average value in the frequency direction of the spectrum of the negative region from the spectrum of the negative region, and a partial region for the deviation of the spectrum of the positive region And calculate the average value of the partial area with respect to the deviation of the spectrum of the negative area. And the average value of the partial region of the deviation of the spectrum of the positive region and the average value of the partial region of the spectrum of the negative region are summed for each partial region to calculate a sum value for each partial region Means for binarizing a component value corresponding to the partial region depending on whether or not a sum value for each partial region is larger than a binarization threshold, and the binary Means for creating an adjacency matrix for the quantization matrix, means for creating a stationary solution matrix by multiplying the adjacency matrix until it reaches a steady state, and a partial region having a maximum sum value for each partial region from the stationary solution matrix Means for extracting a component corresponding to a concatenated region composed of partial regions that are concatenated with each other and whose sum value for each partial region is larger than the binarization threshold, and creating a filter base matrix from the binarization matrix, Filter base matrix A filter matrix by limiting the spread in the frequency direction for each row of the filter transmission component, and the deviation of the spectrum in the positive region and the deviation of the spectrum in the negative region for each distance in the beam direction and for each frequency. Means for calculating a spectrum deviation sum value by summing, a means for applying the filter matrix to the spectrum deviation sum value to create a filter application spectrum, a means for detecting a peak frequency from the filter application spectrum, and A surface layer flow velocity estimation apparatus comprising means for calculating a beam flow velocity by multiplying a peak frequency and a conversion coefficient.   At least means for reading Doppler spectrum data for each beam direction distance and frequency of the processing target region when detecting the peak frequency of the Doppler spectrum for each beam direction distance and frequency and estimating the surface layer flow velocity, the Doppler spectrum Means for separating data into a positive frequency spectrum and a negative frequency spectrum; means for calculating an average value in a frequency direction for the positive frequency spectrum and the negative frequency spectrum; Subtract the average value in the frequency direction of the spectrum in the positive region from the spectrum to calculate the deviation of the spectrum in the positive region and subtract the average value in the frequency direction of the spectrum in the negative region from the spectrum in the negative region Means for calculating the deviation of the spectrum of the negative region, the deviation of the spectrum of the positive region Means for calculating the average value of the partial area and calculating the average value of the partial area for the deviation of the spectrum of the negative area, the average value of the partial area of the deviation of the spectrum of the positive area and the spectrum of the negative area Means for summing the average value of the partial areas of the deviation for each partial area to calculate a sum value for each partial area, whether the sum value for each partial area is greater than a binarization threshold Means for binarizing a component value corresponding to a region to create a binarization matrix, means for creating an adjacency matrix for the binarization matrix, and multiplying the adjacency matrix until it reaches a steady state to obtain a stationary solution matrix Means for creating, corresponding to a connected region consisting of partial regions connected from the steady solution matrix to a partial region having the maximum combined value for each partial region and having a combined value for each partial region greater than the binarization threshold You Means for extracting a component to create a filter base matrix from the binarization matrix, means for limiting the spread of the filter transmission component in the frequency direction for each row of the filter base matrix, and creating the filter matrix, the positive region Means for calculating a spectral deviation sum by adding the deviation of the spectrum of the spectrum and the deviation of the spectrum of the negative region for each beam direction distance / frequency, and applying the filter by applying the filter matrix to the sum of the spectral deviations A surface layer flow velocity estimation program for causing a computer to function as a means for creating a spectrum, a means for detecting a peak frequency from the filter application spectrum, and a means for calculating a beam direction flow velocity by multiplying the peak frequency and a conversion coefficient.