JP5394690B2 - Wind prediction apparatus and program - Google Patents

Description

  The present invention relates to a wind prediction device and a program, and more specifically to a device and a computer program for predicting the wind direction and wind speed on the ground, particularly a strong wind region.

  Predicting the wind (wind direction and wind speed) is useful, for example, as a countermeasure against strong winds generated locally, and also contributes to improving the accuracy of rainfall prediction. Prediction of wind requires both a technology or system for precisely observing the current wind direction and speed, and a system that can accurately predict from the observed values.

  As a method for comprehensively predicting temporal changes in weather, there is a method of applying a certain physical equation to meteorological data on meshes or lattice points assumed at predetermined distance intervals in a three-dimensional space and computer-simulating the time conversion. Are known. In Patent Document 1, a weather forecast model is created using the wide-area weather forecast data periodically provided from the weather data server DS as an initial value, and observation results of weather conditions are assimilated into this model, and wind distribution is calculated from the forecast model. Is described, and a system for providing the calculation result as weather prediction information is described.

In addition, as a method of measuring the wind direction and wind speed, there are a high-rise meteorological observation using a sonde or a wind profiler, an observation using a wind direction and anemometer on the ground, and an observation using a Doppler radar.
JP 2007-017316 A

  Observations with a sonde or wind profiler are rough both spatially and temporally, so they cannot be used to predict strong winds, especially for short time predictions such as minutes or tens of minutes (Nauquiast).

  In that respect, observation by Doppler radar can obtain data in a sufficient spatial and temporal range. However, the information obtained from one Doppler radar is the velocity distribution in the direction away from the approaching direction on the radial line, and the wind direction / wind velocity distribution cannot be obtained as it is. Conventional methods such as VAD and VVP estimate the wind direction and wind speed from the observation result of one Doppler radar under the assumption that the wind field is uniform or linear. Average value in a certain area. That is, a phenomenon in which the wind changes suddenly cannot be captured.

  Under such circumstances, there has been no example of wind naust with a high spatial resolution of 500 m so far.

  An object of this invention is to show the wind prediction apparatus and program which can predict the generation | occurrence | production of a local strong wind exactly.

  A wind prediction apparatus according to the present invention includes a storage unit that stores Doppler velocity data observed by a plurality of Doppler radars, a folding correction unit that corrects the Doppler velocity data, and a Doppler corrected by the folding correction unit. Wind direction / wind speed distribution calculation means for calculating wind direction / wind speed distribution data by referring to topographic data from velocity data, and assimilating the wind direction / wind speed distribution data into a numerical weather model to generate wind field data at a predetermined grid interval Data assimilation means, movement vector estimation means for estimating a wind movement vector from a plurality of wind speed field data within a predetermined time, and according to the movement vector, the current wind speed field data is extrapolated in the time direction and the upper wind speed field Predicting a predetermined ground wind region according to the prediction result of the wind speed field prediction means Characterized in that it comprises a surface wind range prediction means that.

  The wind prediction program according to the present invention includes a function of storing Doppler velocity data observed by a plurality of Doppler radars in a storage unit of a computer, a folding correction function for causing the computer to correct the Doppler velocity data, and A wind direction / wind speed distribution calculation function that allows a computer to calculate wind direction / wind speed distribution data with reference to terrain data from Doppler speed data corrected by the aliasing correction function, and the computer calculates the wind direction / wind speed distribution data numerically. A data assimilation function for assimilating data to a weather model and generating wind velocity field data at predetermined grid intervals; a movement vector estimation function for causing the computer to estimate a wind movement vector from a plurality of wind velocity field data within a predetermined time; According to the movement vector, the computer The wind speed field prediction function that extrapolates the wind speed field data in the time direction to predict the sky wind speed field, and the ground wind area prediction that causes the computer to predict a predetermined ground wind area according to the prediction result of the sky wind speed field prediction function And a function.

  According to the present invention, it becomes possible to predict a ground wind after a short time, for example, a strong wind region, with high accuracy locally.

  Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.

  FIG. 1 shows a schematic block diagram of an embodiment of the present invention. In this embodiment, it is assumed that it operates on one computer for easy understanding, but of course, a large amount of data can be handled by operating a plurality of computers in a coordinated manner at high speed. It is clear that results can be obtained.

  The Doppler radar radiates radio waves into the air and measures the moving speed (Doppler speed) in the radial direction of airborne substances such as precipitation particles from the frequency change due to the Doppler effect of the reflected waves. Precipitation particles move in the horizontal direction at the horizontal wind velocity of the spot, and in the vertical direction, the precipitation particles are thought to move in the sum of the vertical wind speed of the spot and the fall velocity relative to the atmosphere. The Doppler velocity is a radial component of the sum of the wind velocity and the falling velocity at the measurement position. In general, the relationship between the precipitation velocity of precipitation particles and the radar reflection factor (received radio wave intensity) is known. The radial component of the precipitation velocity estimated from the radar reflection factor is obtained, and this is calculated from the Doppler velocity. By subtracting, the radial component of the wind speed can be known.

  In general, Doppler radar observation is based on PPI (Plan Position Indicator) observation that rotates in the azimuth direction while being fixed at a certain elevation angle, and performs PPI observation of a plurality of elevation angles within a certain time (volume scan). Thus, three-dimensional data is acquired. The Doppler velocity and radar reflection factor data observed by the Doppler radar are stored together with information on the elevation angle, azimuth angle, and radial distance of the measurement position.

  A single Doppler radar can only measure the velocity component (radial direction component) on the radio wave radiation axis. However, the velocity component measured by each radio wave is obtained by intersecting the radiation beams of multiple Doppler radars at an angle within a predetermined range. The vector can be used to measure the wind speed (wind direction and speed) at the intersection. In this vector synthesis, a significant result can be obtained when the crossing angle is in the range of 90 degrees to plus or minus 60 degrees. It is necessary to install the radar so that the region where the wind speed is to be obtained falls within the range of this intersection angle. For example, when the observation range is 80 km square, a Doppler radar may be arranged at the apex of an equilateral triangle having a side distance of about 40 km.

  In the present embodiment, n Doppler radars 10-1 to 10-n are arranged in a distributed manner so that the wind direction and the wind speed of the measurement target space can be measured. Each of the Doppler radars 10-1 to 10-n analyzes the reflected wave due to the Doppler effect of the radiated radio wave (X band in this embodiment), and outputs a radial Doppler velocity component. By taking into account the delay in reflection, the velocity distribution along the radial line can also be measured. The observation data (Doppler velocity and radar reflection factor data) of each Doppler radar 10-1 to 10-n is recorded via the data network 12 together with information on the elevation angle, azimuth angle and radial distance at the time of each observation. It inputs into the wind direction and wind speed calculation apparatus 14 of an Example.

  The wind direction / wind speed calculation device 14 calculates numerical weather prediction data relating to the wind direction / wind speed from the Doppler velocity data supplied from the Doppler radars 10-1 to 10-n via the data network 12 (in this embodiment, mesoscale). The wind direction / velocity distribution of the target three-dimensional space is calculated with reference to the weather forecast model (MSM output data) and topographic data (altitude / surface roughness data).

The Doppler radar used for precipitation observation is generally a pulse radar, and it is necessary to correct the observed Doppler velocity by turning it back. That is, in the pulse radar, there is a maximum speed (Nyquist speed) determined by the pulse repetition frequency and the wavelength among the Doppler speeds that can be measured. Doppler speeds exceeding this Nyquist speed are expressed in turn. Between the actual Doppler speed V and the detected Doppler speed Vo,
Vnyq = fλ / 4 (1)
V = Vo + 2nVnyq (2)
The relationship is established. However, the Nyquist speed Vnyq is the Nyquist speed, f is the pulse repetition frequency, λ is the wavelength of the radio wave, and n is an integer.

  For example, in the case of a Doppler radar having a pulse repetition frequency f of 2000 Hz and a wavelength λ of 3.2 cm, the expressible Doppler speed is from −16 m / s to +16 m / s. When the measured Doppler velocity Vo is +12 m / s, the actual Doppler velocity V is + 12 + 32 × n (m / s), that is,. . . -52, -20, +12, +44, +76,. . . There is uncertainty. This is because when a signal is sampled, a frequency exceeding 1/2 of the sampling frequency is expressed by being folded.

  Therefore, the wind direction / wind speed calculation device 14 refers to the forecast data of the numerical weather prediction model, corrects the Doppler velocity return correction device 16 for correcting the return of the Doppler velocity, and the Doppler velocity data corrected by the correction device 16 from the topography ( The wind direction / wind speed distribution calculating device 18 calculates the wind direction / wind speed distribution with reference to the altitude / ground surface roughness data). FIG. 2 shows a schematic block diagram of the wind direction / wind speed calculation device 14.

  The hard disk 20a stores Doppler velocity data from the Doppler radars 10-1 to 10-n input via the data network 12, and the hard disk 20b stores radars from the Doppler radars 10-1 to 10-n. Reflection factor data is stored. The hard disk 20c stores forecast data based on a numerical weather model for the target area. This forecast data is represented by a value (Grid Point Value) on a grid point located at a certain distance in a horizontal and vertical three-dimensional space. The hard disk 20d stores altitude / surface roughness data of the target area.

  The Doppler velocity return correction device 16 corrects the Doppler velocity data stored in the hard disk 20a with reference to the radar reflection factor data stored in the hard disk 20b and the numerical weather forecast data stored in the hard disk 20c. FIG. 3 shows a flow of aliasing correction in the Doppler velocity aliasing correction device 16.

  The Doppler velocity folding correction device 16 takes in Doppler velocity data from the hard disk 20a and also takes in radar reflection factor data stored in the hard disk 20b, and excludes apparent abnormal values under a certain quality control standard (S1). . Specifically, since the observed Doppler velocity includes noise due to the reflection of radio waves due to topography and the influence of interference waves, these abnormal data are detected and deleted. Since the terrain does not move, the affected Doppler velocity is almost zero. Therefore, the Doppler velocity data close to zero is first deleted. Next, in order to remove the influence of the interference wave, the moving average deviation of the Doppler velocity data is calculated in the radial direction, and when the value exceeds a set threshold value, the Doppler velocity data is deleted.

  After the abnormal data or noise data is excluded (S1), in this embodiment, aliasing correction is executed by a plurality of methods to improve accuracy.

  First, the aliasing is corrected by comparison with the numerical weather forecast model (S2). A working meteorological agency such as the Japan Meteorological Agency distributes forecast values such as temperature, pressure, wind direction, and wind speed calculated by a numerical meteorological model (for example, Meteorological Agency Meso Numerical Model: MSM) in real time as three-dimensional grid data. The hard disk 20c stores such forecast data of the numerical weather model, and the forecast data of the numerical weather model includes wind speed data at certain time intervals and spatial intervals.

  The Doppler velocity return correction device 16 describes wind speed data of the forecast data as measured by the Doppler radars 10-1 to 10-n (elevation angle, azimuth angle, and radial distance from the radars 10-1 to 10-n). Space-time interpolation. Next, the velocity in the radial direction of the radar 10-1 to 10-n is calculated from the interpolated wind velocity and the precipitation falling velocity estimated from the radar reflection factor observed by the radar. Then, n is determined such that the Doppler velocity V given by Equation (2) is closest to the calculated radial velocity. This operation is performed for all quality-controlled Doppler speed data Vo, and n, that is, Doppler speed V for each is determined.

ｉは動径方向に関する添え字、ｊは方位角方向に関する添え字である。正しく折り返し補正が行われれば、ドップラー速度は空間的になめらかに変化する。そこで、動径方向・方位角方向に隣り合うドップラー速度データの差はナイキスト速度Ｖnyqよりも小さくなると考えられ、Ｄ値は減少する。 In this embodiment, a plurality of simple methods are sequentially used in combination for the aliasing correction. Therefore, the D value of the following equation relating to a certain PPI data is defined according to the following equation as an index indicating how much aliasing correction has been performed by each method. That is,

i is a subscript in the radial direction, and j is a subscript in the azimuth direction. If the aliasing correction is performed correctly, the Doppler speed changes smoothly in space. Therefore, it is considered that the difference between the Doppler velocity data adjacent in the radial direction / azimuth direction is smaller than the Nyquist velocity Vnyq, and the D value decreases.

  The D value is calculated with respect to the Doppler velocity data V corrected by the comparison with the numerical weather forecast data. Let this calculation result be D1.

  Next, the Doppler velocity folding correction device 16 corrects the folding based on the comparison with the horizontal wind by the VAD (Velocity Azimuth Display) method (S3). The VAD method is a method of estimating the altitude distribution of the horizontal wind representing the observation region from one PPI observation data. First, the Doppler velocity folding correction device 16 calculates the altitude distribution of the horizontal wind velocity at the observation time by the VAD method. Next, the altitude of the measurement position by the Doppler radars 10-1 to 10-n is obtained, and the Doppler is calculated from the horizontal wind speed at the altitude obtained by the VAD method and the precipitation falling speed estimated by the radar reflection factor. The radial speed of the radars 10-1 to 10-n is calculated. However, if the elevation angle is up to about 10 degrees, the wind velocity component in the vertical direction can be regarded as zero.

  Similar to the comparison with the numerical weather forecast model in step S2, n that is closest to the radial velocity calculated here is determined by equation (2). That is, the Doppler velocity data corrected in step S2 is substituted as Vo into Equation (2), and the Doppler velocity V obtained by Equation (2) is closest to the radial velocity calculated here. To decide. This operation is executed for all the Doppler speed data corrected in step S2, and n for each, that is, the Doppler speed V is determined.

  With respect to the result of the aliasing correction in step S3, a D value is calculated by the equation (3), and the result is set to D2. When D2 is equal to or less than D1, the aliasing correction in step S3 is validated. If D2 is larger than D1, it is determined that the folding correction is not performed correctly, the Doppler speed data is returned to the one before application of step S3, and D1 is substituted into D2. That is, the loopback correction in step S3 is cancelled.

  Next, the aliasing is corrected based on temporal continuity (S4). The wind speed in the sky changes continuously with respect to time, and when the Doppler radar observation interval is shorter than the fluctuation time scale (for example, about 5 minutes), the comparison with the speed data corrected with the previous aliasing is significant. . Therefore, first, the aliasing corrected Doppler velocity data at the same elevation angle at the immediately preceding time is spatially interpolated at the same position as the measurement position of the aliasing corrected Doppler velocity data from Step S3. Then, the Doppler velocity data corrected in step S3 is substituted for Vo in Equation (2), and the Doppler velocity V obtained in Equation (2) is closest to the spatially interpolated Doppler velocity n. To decide. This operation is executed for all the Doppler speed data corrected in step S3, and n for each, that is, the Doppler speed V is determined.

  With respect to the result of the aliasing correction in step S4, the D value is calculated by the equation (3), and the result is set to D3. When D3 is equal to or smaller than D2, the aliasing correction in step S4 is validated. If D3 is larger than D2, it is determined that the folding correction is not performed correctly, the Doppler speed data is returned to the one before application of step S4, and D2 is substituted into D3. That is, the loopback correction in step S4 is cancelled.

  Finally, aliasing is corrected in consideration of spatial continuity (S5). That is, it is utilized that the wind speeds of adjacent positions at the same time take close values. The wind speed in the sky changes continuously with respect to space, and the data adjacent in the azimuth and radial directions are also continuously distributed with respect to the Doppler velocity.

  First, data points that are more likely to be corrected are extracted from the folding correction results in step S4. The extraction condition is, for example, that the difference between each corrected Doppler velocity and the velocity data used for each comparison is always smaller than a certain threshold value (for example, 1 m / s) in steps S2, S3, and S4. Extract data points that meet the conditions. Next, n that is closest to the Doppler speed of the probable data point is determined with respect to the probable point data and Doppler speed data of the point adjacent in the azimuth / radial direction. Data corrected in this way is also regarded as a probable data point, and data adjacent to this point is corrected. By such recursive processing, the entire Doppler velocity data is corrected.

  With respect to the result of the aliasing correction in step S5, the D value is calculated by the equation (3), and the result is set to D4. When D4 is equal to or smaller than D3, the aliasing correction in step S5 is validated. If D4 is larger than D3, it is determined that the aliasing correction is not performed correctly, and the Doppler velocity data is returned to the one before application of step S5.

  The Doppler speed folding correction device 16 stores the Doppler speed data corrected for folding by a plurality of methods in the hard disk 22 as described above.

  The wind direction / wind speed distribution calculating device 18 refers to the radar reflection factor data of the hard disk 20b and altitude / surface data described later, calculates the wind direction / wind speed distribution of the sky from the aliasing corrected Doppler speed data of the hard disk 22, The calculation result is stored in the hard disk 24.

In general, radar observation makes it difficult to observe the vicinity of the ground surface due to the curvature of the earth and the shielding of radio waves by mountains and structures. However, wind information near the ground surface is very important for disaster prevention. On the other hand, in this embodiment, by assuming a vertical profile model of wind speed called logarithmic rule, the wind direction and wind speed at the specified sea level altitude Z ₁ (for example, 1 km altitude above sea level) in the lower troposphere and the vicinity of the ground surface ( For example, the wind direction and wind speed at a ground altitude of 10 m) are calculated.

Prior to calculations, altitude Z ₁ Oite latitude width and longitude width defining a horizontal grating as lattice spacing of about 500m approximately constant. Next, the altitude / surface roughness data is spatially interpolated on the same latitude / longitude as the defined grid points. However, regarding the spatial interpolation of the ground surface roughness data, the ground surface roughness averaged from the lattice point to a point about 2500 m upwind is used instead of the value of the lattice point. For example, the surface roughness of points separated every 100 m, such as 0 m, 100 m, 200 m,..., 2400 m, 2500 m in the windward direction from the lattice points, is calculated by spatial interpolation, and the average of these is calculated. It is used as the ground surface roughness at the lattice point. This is to consider that the vertical profile of the lower tropospheric wind speed is influenced not only by the ground surface state directly below but also by the ground surface state of the windward side. Further, the wind direction and wind speed used when defining the windward direction are updated every time the calculation is performed immediately before by this embodiment, that is, the data read from the wind direction and wind speed distribution data of the hard disk 24.

と記述される。ここで、ｉは格子点を示す添え字である。ｊは格子点ｉを中心とする、ある大きさの楕円球（たとえば水平軸の半径が１．５ｋｍ、鉛直軸の半径が５００ｍ）の内側に分布するレーダ観測点を示す添え字である。ωは格子点ｉからレーダ観測点までの距離に依存する重みであり、距離が長くなると小さくなる関数で与えられる。α，β，γはそれぞれ、レーダビームがｘ軸（東西）、ｙ軸（南北）及びｚ軸（鉛直）方向となす角度を示す。ｆはレーダ反射因子により推定される降水の落下速度である。Ｋ_ｉは、格子点ｉの周りに分布するレーダ観測点におけるドップラー速度Ｖと格子点ｉにおける風速度の動径方向成分の差についての重み付き二乗平均である。 The wind direction and wind speed on the defined grid points are calculated as follows. Here, the x-axis (east-west), y-axis (north-south), and z-axis (vertical) components of the wind speed on each lattice point are u, v, and w. First, a cost function J that is calculated from the three wind speed components (u, v, w) at each lattice point and the Doppler velocity corrected for aliasing is defined, and (u, v, w) that minimizes the cost function J is defined. Ask. Cost function J is

It is described. Here, i is a subscript indicating a lattice point. j is a subscript indicating radar observation points distributed inside a certain size ellipsoid (for example, the radius of the horizontal axis is 1.5 km and the radius of the vertical axis is 500 m) with the lattice point i as the center. ω is a weight that depends on the distance from the grid point i to the radar observation point, and is given by a function that decreases as the distance increases. α, β, and γ represent angles formed by the radar beam with respect to the x-axis (east-west), y-axis (north-south), and z-axis (vertical) directions, respectively. f is the precipitation falling speed estimated by the radar reflection factor. K _i is a weighted root mean square for the difference between the Doppler velocity V at the radar observation points distributed around the lattice point i and the radial component of the wind velocity at the lattice point i.

If the finally calculated three wind velocity components (u, v, w) are consistent with the observed Doppler velocity, the value of K _i approaches 0 and is the sum of K _i for all grid points. The value of the first item of the cost function J should also approach 0. The second item of the cost function J has a function of a low pass filter related to space, and the wavelength to be attenuated is adjusted by the coefficient C.

  In the case of observation by a Doppler radar installed on the ground, the vertical wind component w is small when calculating (u, v, w) that minimizes the cost function J because the contribution of the vertical wind component is small to the obtained Doppler velocity. Shows unstable behavior. That is, it is difficult to determine the vertical wind speed component w. Therefore, in this embodiment, the assumption of the law of conservation of mass and the logarithm rule of wind speed are applied, and the vertical wind speed component w is described as a function of the x-axis wind speed component u and the y-axis wind speed component v.

と表現できると仮定する。ただし、ｕ＊、ｖ＊は摩擦速度、ｋはカルマン定数、ｚ_ｈは地表面の海抜高度、ｚ_０は地表面粗度である。海抜高度ｚ_１における水平面内の風速を（ｕ、ｖ）とすると、式（１）、（２）から摩擦速度ｕ＊，ｖ＊を消去でき、

と変形できる。 The wind speeds U (z) and V (z) at sea level altitude z near the ground surface are expressed by the logarithmic law of the wind speed,

It can be expressed as However, u * and v * are friction rates, k is the Kalman constant, z _h is the altitude above sea level, and z ₀ is the ground surface roughness. If the wind speed in the horizontal plane at sea level altitude z ₁ is (u, v), the friction velocities u *, v * can be eliminated from equations (1) and (2).

And can be transformed.

となる。ただし、ρは大気密度である。ここでは、地表面付近の現象のみを扱っているので、大気密度ρは、高度によらず一定であると仮定している。Ｌは下記式

のように定義される。 By integrating the equations (4) and (5) from the ground surface height z _h to z ₁ , the mass flux (Fx, Fy) below the altitude z ₁ is

It becomes. Where ρ is the atmospheric density. Here, since only phenomena near the ground surface are handled, it is assumed that the atmospheric density ρ is constant regardless of the altitude. L is the following formula

Is defined as follows.

となる。ただし、ａは地球の半径、λは経度、φは緯度である。水平発散ｄは、

と表される。式（１２）により鉛直風速成分ｗは、水平速度成分（ｕ，ｖ）の関数として表すことができる。これを式（４）に代入することにより、コストファンクションＪは、ｕとｖの関数になる。 From the mass flux (Fx, Fy) of the equations (9) and (10), the vertical wind velocity component w at the altitude z ₁ is expressed by the law of conservation of mass.

It becomes. Where a is the radius of the earth, λ is the longitude, and φ is the latitude. The horizontal divergence d is

It is expressed. The vertical wind velocity component w can be expressed as a function of the horizontal velocity component (u, v) by the equation (12). By substituting this into equation (4), the cost function J becomes a function of u and v.

The cost function J to be minimized is determined by the above calculation. If the wind speed component (u, v) at each lattice point at which the cost function J is minimized is determined, the vertical wind speed component w is obtained by the equation (12), and the altitude z ₁ or less is obtained by the equations (5) and (6). The vertical profile of the wind speed is determined.

  In order to determine the value of the horizontal wind velocity component (u, v) at each grid point that minimizes the cost function J, various solutions handled by nonlinear programming can be used. For example, when the Polak-Ribiere method, which is a kind of conjugate gradient method, is used, each value of (u, v) at each lattice point can be obtained by giving a partial differential for the cost function J and its u, v. The values of the horizontal wind speed components (u, v) on each grid point obtained in this way and the vertical wind speed components (w) obtained from these by the formula (12) are stored in the hard disk 24 and are stored in the data Supplied to the assimilation device 30.

  The data assimilation device 30 predicts the wind direction / wind speed distribution data (wind observation value data) from the wind direction / wind speed calculation device 14 and the cloud resolving model (in this embodiment, Cloud Resolving Storm Simulator (CResS) is adopted). The first estimated value that is the result is assimilated by the three-dimensional variational method, and a wide range of wind speed fields (this is called objective analysis data) that physically matches the observed value is created. By using the three-dimensional variational method as an assimilation method, it is possible to provide objective analysis data of wind speed in real time, unlike the four-dimensional variational method and the ensemble Kalman filter, which have a large calculation cost.

  FIG. 4 shows an operation flow of the data assimilation apparatus 30. The data assimilation device 30 takes in the data of the wind speed field (hard disk) 24) calculated by the wind direction / wind speed calculation device 14 as wind observation values, and takes in the radar reflection factor data as observation values of precipitation particles. The captured observation values are interpolated into the grid points of the cloud resolution numerical model (CReSS model) near the observation location (S21).

  Next, the data corresponding to the observation time is extracted from the prediction data of the cloud resolving numerical model that is operated in real time, and the first estimated value (the wind direction of each grid point, estimated by means other than observation) Wind speed, atmospheric pressure, temperature, water vapor mixing ratio, precipitation particle mixing ratio).

ただし、ρは大気密度である。この合計値を最も小さくする風の推定値が最適な推定値である。総和が一定範囲に収束しない場合（Ｓ２４）、共役勾配法を使って風の推定値を更新し（Ｓ２５）、合計値を再計算する（Ｓ２３）。 Thus, the cloud point of the cloud resolving numerical model in the observation region has two wind data of the observed value and the first estimated value. By repeating the following procedure, an estimated wind value that is as close as possible to both the observed value and the first estimated value and that satisfies the mass conservation law as much as possible is obtained. That is, the difference between the estimated wind value and the observed value, the difference between the estimated wind value and the first estimated value, and the convergence divergence term Dv obtained by Equation (14) are summed for all grid points (S23). ). That is,

Where ρ is the atmospheric density. The estimated wind value that makes the total value the smallest is the optimum estimated value. If the sum does not converge to a certain range (S24), the wind estimation value is updated using the conjugate gradient method (S25), and the total value is recalculated (S23).

  If it is determined that the total value has converged to the minimum value (S24), the data assimilation is terminated and the wind velocity field data (objective analysis data) by the data assimilation is stored in the hard disk 32. The wind speed field data stored in the hard disk 32 is data indicating the wind direction and the wind speed in a space of 1 km or less above sea level in an area assumed by the CReSS model on grid points at predetermined distance intervals. By predicting the moving direction and speed of this wind speed field, it is possible to predict the occurrence of strong winds and their areas.

  The prediction device 34 estimates the movement vector (time change) of the wind speed field from the time change of the wind speed field within a certain time in the present and the past calculated by the data assimilation device 30, and takes 5 minutes in the range up to about one hour ahead Estimate the occurrence of strong wind on the ground and the area. The horizontal resolution was 500 m mesh and the update interval was 5 minutes. FIG. 5 shows an operation flow of the movement vector estimation device 36 and the sky wind velocity field prediction device 38 of the prediction device 34.

  The movement vector estimation device 36 estimates the movement direction and movement speed (hereinafter referred to as a movement vector) of the strong wind region from the current and past wind speed field data of the 500 m mesh of the hard disk 32. Since the purpose of this embodiment is to predict a strong wind region, the current and past three time wind velocity field data are read from the hard disk 32 (S31), and a certain level or more of a strong wind region, here 15 m / s or more is extracted. (S32). A movement vector is calculated using the advection model for the extracted wind speed field data (S33). The values of areas that were not predicted to be strong wind areas are treated as indefinite (wind speed may be weak or strong).

  In this embodiment, an advection model is adopted for estimating the movement vector (S33). The concept of the advection model is briefly explained. When the wind speed is expressed as a function of position (x, y) and time t like Z (x, y, t), Z (x, y, t) represents one curved surface when the time is fixed. The wind speed prediction corresponds to predicting the change of the curved surface.

と表す。式（１５）はいわば質量保存の式であり、ある微小面積内でのＺの時間変化と、その微小面積に周囲から流れ込み又は周囲に流れ出す量との和が、その場所でのＺの増減と等しいことを意味している。ここで、ｕはＺの東西方向の移動速度を示し、ｖはＺの南北方向の移動速度を示す。ｕ，ｖは、Ｚの移動方向と移動速度を表すので、これを移動ベクトルと呼ぶ。また、ｗはＺの発達・衰弱を表す。 Therefore, change in wind speed Z

It expresses. The equation (15) is, so to speak, a mass conservation equation, and the sum of the time change of Z within a certain small area and the amount flowing into or out of the minute area from the periphery is the increase / decrease of Z at that location. Means equal. Here, u represents the moving speed of Z in the east-west direction, and v represents the moving speed of Z in the north-south direction. Since u and v represent the moving direction and moving speed of Z, they are called moving vectors. W represents the development and weakness of Z.

式（１６）により、風速Ｚの平行移動だけでなく、変形も表現できる。 Assuming that (u, v, w) in equation (15) can be expressed by a linear expression of position coordinates using parameters C _{1 to} C ₉ , the following expression is obtained.

Expression (16) can express not only the translation of the wind speed Z but also the deformation.

が得られる。ここで、ｉ，ｊ，ｋは三次元の格子点を示す添え字である。 If the parameters C _{1 to} C ₉ in Expression (16) can be identified, the movement vector of the wind speed Z and the amount of development / weakness can be estimated. Therefore, using the wind speed data (Z) for a plurality of consecutive times and setting the residual (the difference between the left side and the right side) of the equation (15) at point (x, y) and time t _k as V _ijk ,

Is obtained. Here, i, j, and k are subscripts indicating three-dimensional lattice points.

で与えられる評価値Ｊ_ｃを最小にすることに相当する。従って、式（１８）を最小にするように、Ｃ_１〜Ｃ_９を決めれば良く、具体的には、式（１８）のＪ_ｃをＣ_１〜Ｃ_９で偏微分して０とおいて得られる連立一次方程式を解くことにより求めることができる。 If the parameters C _{1 to} C ₉ are correctly identified, the residual V _ijk should be zero. This is the following formula

It corresponds to minimizing the evaluation value J _c given by. Therefore, C _{1 to} C ₉ may be determined so as to minimize Equation (18). Specifically, J _c in Equation (18) is partially differentiated with C _{1 to} C ₉ and set to 0. It can be obtained by solving the simultaneous linear equations.

In this model, it is possible to express even rotation in the wind speed region. However, there are cases where there are problems with accuracy, so the prediction range is divided into several small areas, and after calculating a single movement vector for each small area, they are interpolated and replaced with movement vectors in mesh units. To express the change in the direction of movement of the wind speed region. The movement vector for each small area is set to 0 except for C ₃ and C ₆ among the parameters C _{1 to} C ₉ . This is equivalent to assuming that the movement of the wind speed region is a parallel movement in each small region.

  The upper wind speed field prediction device 38 extrapolates the strong wind region of the current wind speed field data of the hard disk 32 on the assumption that the movement vector estimated by the movement vector estimation device continues for one hour ahead (S34). Thereby, the position or movement of the strong wind area up to 1 hour ahead is predicted.

  The strong wind region predicted by the upper wind speed field prediction device 38 is a strong wind region in the sky (for example, altitude 1000 m). Therefore, the ground strong wind region prediction device 40 predicts the ground strong wind region with respect to the sky strong wind region at each time predicted by the sky wind speed field prediction device 38.

  When observation values of the wind direction and wind speed of the target area exist, the difference between the current wind speed observation value and the upper wind speed value is added to the upper wind speed prediction value obtained by the upper wind speed field prediction device 38, and the addition The result is the ground wind speed prediction value.

  When there is no observation value in the target area, the ground wind speed value is predicted from the predicted upper wind speed value obtained by the upper wind speed field prediction device 38 based on a vertical wind speed profile called a logarithmic rule.

  The ground strong wind region prediction device 40 stores the ground wind direction and wind speed data obtained in this way at each time in the hard disk 42. Since the purpose is to display a strong wind region, the ground wind direction and wind speed data stored in the hard disk 42 is data of an area where the wind is high in the sky.

  In addition, the amount of calculation can be reduced by making the estimation unit of the movement vector estimation device 36 and the prediction unit of the sky wind velocity field prediction device 38 and the ground strong wind region prediction device 40 into small regions obtained by dividing the prediction target region, Calculation time can be shortened.

  The display device 44 displays the wind direction / wind speed distribution information stored in the hard disk 42 as a graphic for displaying the intensity classified by color on the screen of the graphical user interface.

  According to this embodiment, wind, particularly strong wind, can be predicted locally. Since it becomes possible to consider the development and weakness of rainfall based on the wind prediction information obtained, the short-term prediction accuracy of rainfall is improved.

  Although this embodiment is mainly realized by a computer program, it is obvious that the same effect can be obtained even if a part of the function is replaced with dedicated hardware. Further, not only a computer program that operates on a single computer but also a computer program that operates on a large number of computers can be operated in a coordinated manner to achieve the same effects. Any of these configurations belong to the technical scope of the present invention.

  Although the invention has been described with reference to specific illustrative embodiments, various modifications and alterations may be made to the above-described embodiments without departing from the scope of the invention as defined in the claims. This is obvious to an engineer in the field to which the present invention belongs, and such changes and modifications are also included in the technical scope of the present invention.

It is a schematic block diagram of one Example of this invention. 3 is a schematic functional block diagram of a wind direction / wind speed calculation device 14. FIG. A flow of aliasing correction in the Doppler velocity aliasing correction device 16 is shown. 4 is an operation flow of the data assimilation apparatus 30. It is an operation | movement flow of the movement vector estimation apparatus 36 and the sky wind speed field prediction apparatus 38. FIG.

Explanation of symbols

10-1 to 10-n: Doppler radar 12: Data network 14: Wind direction / wind speed calculation device 16: Doppler speed return correction device 18: Wind direction / wind speed distribution calculation device 20a: Hard disk (Doppler speed data)
20b: Hard disk (radar reflection factor data)
20c: Hard disk (Numerical weather model forecast data)
20d: Hard disk (altitude / surface roughness data)
22: Hard disk (Folder corrected Doppler speed data)
24: Hard disk (wind direction / velocity distribution data)
30: Data assimilation device 32: Hard disk (wind velocity field data)
34: Predictor 36: Movement vector estimator 38: Above wind speed field predictor 40: Ground strong wind region predictor 42: Hard disk (wind direction / wind speed distribution prediction data)
44: Display device

Claims (8)

Storage means (20a) for storing Doppler velocity data observed by a plurality of Doppler radars;
A folding correction means (16) for correcting the Doppler speed data by folding;
Wind direction / wind speed distribution calculating means (18) for calculating wind direction / wind speed distribution data by referring to the terrain data from the Doppler speed data corrected by the folding correction means;
Data assimilation means (30) for assimilating the wind direction / velocity distribution data into a numerical weather model and generating wind field data at predetermined grid intervals;
Movement vector estimation means (36) for estimating a wind movement vector from a plurality of wind velocity field data within a predetermined time;
According to the movement vector, an upper wind speed field predicting means (38) for extrapolating the current wind speed field data in the time direction to predict the upper wind speed field,
Ground wind region prediction means (40) for predicting a predetermined ground wind region according to the prediction result of the above-mentioned sky wind velocity field prediction means
The wind prediction apparatus characterized by comprising.   The aliasing correction means considers the first aliasing correction by comparison with the numerical weather forecast data, the second aliasing correction by comparison with the horizontal wind, the third aliasing correction by comparison in the time direction, and spatial continuity. The wind prediction apparatus according to claim 1, wherein two or more of the fourth aliasing corrections are performed. The storage means stores radar reflection factor data from the plurality of Doppler radars;
The wind prediction apparatus according to claim 2, wherein the radar reflection factor is referred to in the first and second folding corrections of the folding correction unit.   4. The ground strong wind region predicting means (40) for predicting a strong wind region on the ground according to a prediction result of the sky wind speed field predicting unit. The wind prediction apparatus described in 1. A function of storing Doppler velocity data observed by a plurality of Doppler radars in a storage means of a computer;
A folding correction function (16) for causing the computer to correct the Doppler velocity data;
A wind direction / wind speed distribution calculation function (18) for causing the computer to calculate wind direction / wind speed distribution data with reference to the terrain data from the Doppler speed data corrected by the turn-back correction function;
A data assimilation function (30) for causing the computer to assimilate the wind direction / velocity distribution data into a numerical meteorological model and generate wind velocity field data of a predetermined grid interval;
A movement vector estimation function (36) for causing the computer to estimate a wind movement vector from a plurality of wind speed field data within a predetermined time;
An upper wind speed field prediction function (38) for causing the computer to extrapolate current wind speed field data in the time direction according to the movement vector and predicting the upper wind speed field;
A ground wind region prediction function (40) for causing the computer to predict a predetermined ground wind region in accordance with the prediction result of the sky wind speed field prediction function.
The wind prediction program characterized by comprising.   The loopback correction function includes a first loopback correction function by comparison with numerical weather forecast data, a second loopback correction function by comparison with horizontal wind, a third loopback correction function by comparison in time direction, and spatial continuity The wind prediction program according to claim 5, wherein the wind prediction program has two or more of the fourth aliasing correction functions in consideration of characteristics. The storage means stores radar reflection factor data from the plurality of Doppler radars;
The wind prediction program according to claim 6, wherein the first and second aliasing correction functions refer to the radar reflection factor.   The ground surface wind region prediction function is a ground strong wind region prediction function (40) for causing the computer to predict a ground strong wind region according to a prediction result of the sky wind speed field prediction function. The wind prediction program according to claim 1.

Images