CN109100723B - High-altitude wind inversion method based on Doppler weather radar data - Google Patents

High-altitude wind inversion method based on Doppler weather radar data Download PDF

Info

Publication number
CN109100723B
CN109100723B CN201810826385.5A CN201810826385A CN109100723B CN 109100723 B CN109100723 B CN 109100723B CN 201810826385 A CN201810826385 A CN 201810826385A CN 109100723 B CN109100723 B CN 109100723B
Authority
CN
China
Prior art keywords
radar
data
alt
sounding
gamma
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810826385.5A
Other languages
Chinese (zh)
Other versions
CN109100723A (en
Inventor
王兴
卞浩瑄
朱彬
杨胜朋
苗春生
王军
徐经纬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN201810826385.5A priority Critical patent/CN109100723B/en
Publication of CN109100723A publication Critical patent/CN109100723A/en
Application granted granted Critical
Publication of CN109100723B publication Critical patent/CN109100723B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/95—Radar or analogous systems specially adapted for specific applications for meteorological use
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention provides a high-altitude wind inversion method based on Doppler weather radar data. The method utilizes the extrapolation result of the radar echo image to convert the average radial velocity field of each elevation surface into a velocity vector field with U, V components, and improves the readability of the high-altitude wind field subjective analysis. The radar wind field is properly corrected by utilizing the sounding data, and the correction strength is gradually reduced along with the increase of the distance between the radar grid point and the sounding station, so that the wind field has better integrity and continuity, and better conforms to the condition of the actual atmospheric wind field. Compared with a plurality of multivariate data fusion methods, the algorithm disclosed by the invention is moderate in operation complexity, and can improve the timeliness of radar secondary wind field output.

Description

High-altitude wind inversion method based on Doppler weather radar data
Technical Field
The invention belongs to the field of geoscience, and particularly relates to a method for performing multilevel high-altitude wind inversion by using detection data such as Doppler weather radar and the like.
Background
The characteristics of the high-altitude wind field are one of important objects for weather diagnosis and analysis, and play a vital role in the fields of meteorology, military, aerospace and the like. Currently, there are many ways to acquire the atmospheric wind field over the earth, such as observing or detecting by using sounding balloons, doppler weather radar, wind profile radar, meteorological satellites, and other devices. The devices and the technical means have the advantages that the Doppler weather radar (hereinafter referred to as radar) has the characteristics of high detection spatial resolution, high data updating frequency, high effectiveness and the like, and is widely researched and applied in the aspects of local weather condition monitoring, short-time near weather forecast and the like.
The average radial velocity is one of basic products of the radar, the velocity is obtained by taking the Doppler effect as a basic principle, and the wind field information represented by the average radial velocity is a velocity component in a certain radial direction swept along the radar body by taking the radar as a central point, but not a true wind velocity vector. In order to obtain actual wind field information of different height layers in a radar detection coverage range, a plurality of scholars propose respective methods for many years, and some methods utilize various products of the radar, such as a reflectivity factor and an average radial speed, construct a constraint equation through a mathematical method and reversely generate the wind field information; and in some methods, wind field information is calculated by using and combining various detection data, such as wind profile radar wind measurement, satellite cloud wind and other metadata. Due to the complexity of a weather system and the limitation of the prior art means, in practical business application, certain defects and shortcomings still exist for inversion of a high-altitude live wind field. For example, the wind field inversion is performed only by using the reflectivity factor and the average radial velocity of a single radar, the method is relatively simple, but the wind field inversion accuracy is low due to the adverse effects that radar clutter and electromagnetic wave beams become deflected and bent as the detection distance increases. For another example, the method is mainly limited by the problems of inconsistent spatial resolution, inconsistent spatial grid point correspondence, inconsistent observation time and the like of different data by fusing multiple data such as radar radial velocity and cloud wind or water vapor wind guide, so that the calculated wind field accuracy is high and low in time and poor in stability.
Disclosure of Invention
In order to overcome the defects of the method, the invention provides a method for carrying out multi-level high-altitude wind field fusion by using an extrapolation product and a radial velocity product of a Doppler weather radar and sounding report data.
The invention adopts the following specific technical scheme:
a high altitude wind inversion method based on Doppler weather radar data comprises the following steps:
step 1: supposing that high-altitude wind at a certain time t is to be inverted, firstly, a radar is selected, the radar data of the radar at the time t and N times of radar data closest to the time t before are utilized, an extrapolation algorithm is adopted to obtain displacement characteristic vectors of a storm body represented by the radar data at each detection point, and the displacement characteristic vector of any detection point is recorded as
Figure GDA0003592351250000021
Wherein N is 2<N is an integer not more than 10;
Figure GDA0003592351250000022
representing a certain elevation surface of radar body scanning, gamma representing a distance library between a radar detection point and a radar central point, omega representing the azimuth angle of radar body scanning, and t representing the corresponding time of processed radar data;
step 2: shifting the feature vector in the previous step
Figure GDA0003592351250000023
According to the time period T of the radar body scanitLength L corresponding to Unit distance libraryuConverted into a moving velocity vector according to the following formula
Figure GDA0003592351250000024
Figure GDA0003592351250000025
Figure GDA0003592351250000026
Information characterizing the value and direction of the velocity as a vector, definition
Figure GDA0003592351250000027
The included angle between the pointed direction and the north direction is theta, when
Figure GDA0003592351250000028
When pointing to the north direction, theta is 0, and rotates clockwise, and theta increases gradually;
and step 3: moving velocity vector in the last step
Figure GDA0003592351250000029
Decomposed into U, V two components, where the U component is perpendicular to the current radial direction and the V component coincides with the current radial direction, assuming
Figure GDA00035923512500000210
To represent
Figure GDA00035923512500000211
A component in the U axis;
Figure GDA00035923512500000212
to represent
Figure GDA00035923512500000213
A component at the V axis;
and 4, step 4: assuming any probe point in the radar data
Figure GDA00035923512500000214
Radar average radial velocity of position of
Figure GDA00035923512500000215
When in use
Figure GDA00035923512500000216
Then, the direction indicating the speed is directed to the center of the circle along the radial direction; on the contrary, when
Figure GDA00035923512500000217
Then, the direction representing the velocity is radially away from the center of the circle; will be provided with
Figure GDA00035923512500000218
As a reference speed, calculating
Figure GDA00035923512500000219
And
Figure GDA00035923512500000220
and enlarging or reducing at the magnification
Figure GDA00035923512500000221
Is given a numerical value of
Figure GDA00035923512500000222
The calculation formula is as follows:
Figure GDA00035923512500000223
and 5: component of U direction
Figure GDA00035923512500000224
Component with V direction
Figure GDA00035923512500000225
The motion vectors are combined according to the following formula, and the corrected motion velocity vector is obtained through calculation and is recorded as
Figure GDA00035923512500000226
Figure GDA00035923512500000227
Step 6: according to the steps 2 to 5, displacement characteristic vectors of each detection point of the radar are subjected to
Figure GDA00035923512500000228
Calculating one by one to obtain the corrected moving velocity vectors of each elevation angle surface, each distance base and each azimuth angle swept by the radar body at the time t, namely
Figure GDA00035923512500000229
And 7: for any detection point in radar data
Figure GDA00035923512500000230
Searching a sounding report station closest to the geographical position of the sounding report station, and reading sounding report station data of which the detection time is closest to the time t in the sounding report station, wherein the data records the wind direction and the wind speed of a certain latitude position on the earth and vertically upwards at a plurality of height layers or a plurality of atmospheric pressure layers in a standard level; obtaining the sounding station reporting data corresponding to each detection pointIs recorded as Vh(lot,lat,alt)；
And 8: the station data V is reported in the above stepsh(lot, lat, alt) each alt is a reference height, and the corrected moving velocity vector of each elevation surface of the radar is used
Figure GDA0003592351250000031
Converting the altitude into a moving velocity vector V at the corresponding altitude of each alt by using a spatial interpolation algorithmip(alt, γ', ω, t); gamma 'represents the distance between the detection point and the radar central point in the horizontal direction, and the value of gamma' is formed by gamma,
Figure GDA0003592351250000032
And LuObtaining the target through a triangular geometric method;
and step 9: defining an influence radius R and a maximum influence radius RifSatisfy 0<r<RifR represents VipThe spatial distance between any point (alt, gamma', omega, t) and the sounding site can be determined by the longitude, latitude and V of the sounding siteipThe longitude and latitude of the point in (alt, gamma', omega, t) are obtained by geometric calculation; vipThe longitude and the longitude of each point (alt, gamma ', omega, t) are calculated according to the geographic position of the radar and gamma';
constructing an influence weight function Wif(r) which is a decreasing function subject to the following conditions: 0<Wif(r)<1; reporting station data V according to soundingh(lot, lat, alt), for V in step 8ip(alt, γ', ω, t) deviation correction is performed by the following formula:
Vbs(alt,γ′,ω,t)＝(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)。
the invention is further designed in that:
in step 1, the data form of the N times of radar data is basic reflectivity data and average radial velocity data.
In step 1, the extrapolation algorithm uses a maximum correlation method or a cross correlation method.
In step 7, if the record of the sounding station data only has atmospheric pressure information and no height information from the ground, the recording is obtained by conversion through a high-pressure formula meteorological method.
In step 7, the data of the sounding report station is interpolated to a standard level by adopting a non-standard level and an interpolation correlation algorithm.
In step 8, the corrected moving velocity vector
Figure GDA0003592351250000033
The interpolated alt height is calculated by the inverse distance weighting algorithm through the moving speed vectors of the upper and lower 1 elevation surfaces adjacent to the position space.
In step 9, Vip(alt, γ', ω, t) can be solved using parallelogram or orthonormal decomposition methods.
Compared with the prior art, the invention has the following advantages:
1. the method utilizes the extrapolation result of the radar echo image to convert the average radial velocity field of each elevation surface into a velocity vector field with U, V components, and improves the readability of the high-altitude wind field subjective analysis.
2. The radar wind field is properly corrected by utilizing the sounding data, and the correction strength is gradually reduced along with the increase of the distance between the radar grid point and the sounding station, so that the wind field has better integrity and continuity and better conforms to the condition of the actual atmospheric wind field.
3. Compared with a plurality of methods for fusing multivariate data, the algorithm disclosed by the invention is moderate in operation complexity, and can improve the timeliness of the output of a radar secondary product (wind field).
Drawings
FIG. 1 is a main control flow chart of the first embodiment of the present invention.
FIG. 2 is a moving velocity vector in the first embodiment
Figure GDA0003592351250000041
Is decomposed into U, V two components, where the U component is perpendicular to the current radial direction and the V component coincides with the current radial direction.
Fig. 3 is an auxiliary description of a calculation process for interpolating velocity vectors of each elevation plane of the radar to an isobaric plane in step 8 of an application example.
FIG. 4 shows W constructed in step 9 of the application exampleif(r) a function image of the function.
Detailed Description
The technical scheme of the invention is further explained by combining the attached drawings.
The first embodiment is as follows:
as shown in fig. 1, the high altitude wind inversion method based on doppler weather radar data of the present invention includes the following steps:
step 1: assuming that high-altitude wind at a certain time t is to be inverted, firstly, a radar is selected, the radar data (in the form of basic reflectivity data and average radial velocity data) of the radar at the time t and the N times before the time t and closest to the time t are utilized, an extrapolation algorithm, such as a maximum correlation method, a cross correlation method and the like, is adopted to calculate the displacement characteristic vector of a storm body represented by the radar data at each detection point, and the displacement characteristic vector of any detection point is recorded as
Figure GDA0003592351250000042
Wherein N is 2<N is an integer not more than 10;
Figure GDA0003592351250000043
the number of elevation surfaces and the elevation angle degree of each elevation surface are determined by the technical parameters of the radar and the mode of the radar body sweep; gamma represents a distance library between a radar detection point and a radar central point, the distance library is a basic unit of detection resolution along the radial direction in radar echo signal processing, and the length of the unit distance library also depends on radar technical parameters and a radar body scanning mode; omega represents the azimuth angle of the radar body sweep, the azimuth angle in the north direction is usually 0 degrees, the azimuth angle is increased in the clockwise direction, and the body sweep is 360 degrees for one circle; t represents the time corresponding to the processed radar data.
Figure GDA0003592351250000044
The physical meaning of (A) is as follows: using any time T as a reference, and calculating a time period T of the radar body scanning by using radar detection data at the time T, T-1, T-2, … and T-N +1 through an extrapolation algorithmitFrom any elevation angle
Figure GDA0003592351250000045
The speed of the radar echo moving at the space position determined by any distance library gamma and any azimuth angle omega is given by the following unit: distance bin number/volume sweep time period.
Step 2: shifting the feature vector in the previous step
Figure GDA0003592351250000046
According to the time period T of the radar body scanitLength L corresponding to unit distance libraryuConverted into a moving velocity vector
Figure GDA0003592351250000047
Wherein, TitIn units of minutes (min), LuIn units of kilometers (km),
Figure GDA0003592351250000048
the unit of (a) is meter/second (m/s), and the calculation method is as follows:
Figure GDA0003592351250000051
as a vector, the vector is a vector,
Figure GDA0003592351250000052
besides the value of the speed, the direction information is also provided, and the definition is not limited
Figure GDA0003592351250000053
The included angle between the pointed direction and the north direction is theta
Figure GDA0003592351250000054
When pointing in the north direction, theta is 0, androtating clockwise, θ increases gradually.
And step 3: moving velocity vector in the last step
Figure GDA0003592351250000055
Is decomposed into U, V two components, where the U component is perpendicular to the current radial direction and the V component coincides with the current radial direction, as shown in fig. 2. U, V the two components can be calculated by orthogonal decomposition, etc., assuming that
Figure GDA0003592351250000056
To represent
Figure GDA0003592351250000057
A component in the U axis;
Figure GDA0003592351250000058
represent
Figure GDA0003592351250000059
The component at the V axis.
And 4, step 4: assuming any probe point in the radar data
Figure GDA00035923512500000510
Radar average radial velocity of position of
Figure GDA00035923512500000511
This value is one of the basic radar products and is directly available. When in use
Figure GDA00035923512500000512
Then, the direction of the velocity is directed to the center of a circle along the radial direction, wherein the center of the circle is the center point of the radar shown in fig. 2; on the contrary, when
Figure GDA00035923512500000513
The direction representing the velocity is radially away from the center of the circle. Will be provided with
Figure GDA00035923512500000514
As a reference speedDegree, calculation
Figure GDA00035923512500000515
And
Figure GDA00035923512500000516
and enlarging or reducing at the magnification
Figure GDA00035923512500000517
To obtain a value of
Figure GDA00035923512500000518
The calculation method comprises the following steps:
Figure GDA00035923512500000519
and 5: component of U direction
Figure GDA00035923512500000520
Component with V direction
Figure GDA00035923512500000521
Merging, calculating to obtain the corrected moving speed vector, and recording as
Figure GDA00035923512500000522
The calculation method comprises the following steps:
Figure GDA00035923512500000523
step 6: according to the steps 2 to 5, displacement characteristic vectors of each detection point of the radar are subjected to
Figure GDA00035923512500000524
Calculating one by one to obtain the corrected moving velocity vectors of each elevation angle surface, each distance base and each azimuth angle swept by the radar body at the time t, namely
Figure GDA00035923512500000525
And 7: for any detection point in radar data
Figure GDA00035923512500000526
Searching a sounding report station closest to the geographical position of the sounding report station, and reading sounding report station data of which the detection time is closest to the time t in the sounding report station, wherein the data records the wind direction and the wind speed of a certain latitude position on the earth and vertically upwards at a plurality of height layers or a plurality of atmospheric pressure layers in a standard level; obtaining sounding station reporting data corresponding to each detection point, and recording the data as Vh(lot,lat,alt)；
If the data record only has atmospheric pressure information and no height information from the ground, the data record can be converted by a meteorological method such as a high pressure formula. For Vh(lot, lat, alt) wherein lot represents longitude, lat represents latitude, alt represents altitude, VhA wind speed vector representing a respective spatial position; generally, the number of height layers or air pressure layers detected and recorded by each sounding station is not fixed, and how many groups of detection values are recorded by each sounding station data, how many alt values are, which are not defined as alt1, alt2, …, altn, alt ∈ { alt1, alt2, …, altn }; the longitude and latitude coordinate position of each piece of station data is fixed and unchanged; the longitude and latitude coordinate positions of different sounding report stations are often different.
It should be noted that, in the sounding report used in the meteorological department service in China, the recorded height layer or air pressure layer is usually fixed in several standard layers, and for the nonstandard layer, an interpolation correlation algorithm can be adopted to interpolate to the standard layer. The following steps are for the calculation of the standard hierarchy, and are also applicable to the case of the non-standard hierarchy after the interpolation calculation.
And 8: the station data V is reported in the above stepsh(lot, lat, alt) each alt is a reference height, and the corrected moving velocity vector of each elevation surface of the radar is used
Figure GDA0003592351250000061
Using air spaceInterpolation algorithm, converting into moving velocity vector of height corresponding to each alt, and recording as Vip(alt, γ', ω, t). Wherein, the spatial interpolation algorithm can adopt an inverse distance weight algorithm, gamma 'represents the distance between the detection point and the radar central point in the horizontal direction, the unit is kilometer (km), and the value of gamma' is formed by gamma,
Figure GDA0003592351250000062
And LuObtaining the target through a triangular geometric method;
and step 9: the longitude and latitude of the sounding site are known, VipThe longitude and the longitude of each point (alt, gamma ', omega, t) are calculated according to the geographic position of the radar and gamma';
defining an influence radius R and a maximum influence radius RifSatisfy 0<r<RifR represents VipThe spatial distance between any point (alt, gamma', omega, t) and the sounding site can be determined by the longitude, latitude and V of the sounding siteipThe longitude and latitude of the point in (alt, gamma', omega, t) are obtained by calculation through a mathematical method;
constructing an influence weight function Wif(r), which is an empirical function, subject to the condition: it is a decreasing function and 0<Wif(r)<1; because the sounding report is the observed real wind field information, the sounding report has higher reliability compared with various wind field information of radar product inversion. Reporting station data V according to soundingh(lot, lat, alt), for V in step 8ip(alt, gamma', omega, t) to make deviation correction, the calculation method is:
Vbs(alt,γ′,ω,t)＝(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)
due to Vbs(alt,γ′,ω,t)、Vip(alt, γ', ω, t) and VhThe (lot, lat, alt) are all vectors, so the calculation process should conform to the calculation rule of vector addition and subtraction, such as parallelogram or orthogonal decomposition.
Thereby completing a probing report station data VhCorrecting deviation of each height alt in (lot, lat, alt) to obtain a standard height layer wind field vector V after deviation correctionbs(alt, γ', ω, t); and further completing deviation correction of all the data of the station detection data, and obtaining wind field vector information of a plurality of standard height layers.
Theoretically, the more effective sounding sites and effective sounding wind field records are within the radar detection range, the more Vbs(alt, gamma', omega, t) can reflect the wind field characteristics of each high-altitude layer more truly; on the contrary, if the sounding site is too sparse or the sounding wind field records are too few, the radius R is influenced to the maximum extentifIn this case, the calculation process of step 9 may not be performed in some regions of the radar detection range, which also takes into account the discontinuity of the wind field on a large scale and avoids "overcorrection".
Through the calculation of the steps, the wind field vector information of a plurality of standard height layers in the radar detection range can be obtained. After the information is graphically displayed, auxiliary decision can be provided for professional personnel to carry out weather diagnosis and analysis; this information can also be entered into the numerical model as an initial field to improve the accuracy of the weather pattern forecast.
Application example one:
in the embodiment, doppler weather radar data in the Nanjing area, where t is "5, 14, 8, 45 min in 2018", is selected as a target to be implemented, the radar is a radar in an S-band, a VCP21 volume sweep mode is adopted, and the volume sweep period is generally 6 minutes. The final purpose of the embodiment is to calculate the atmospheric wind field of a plurality of standard height layers through the following steps, so that the wind field has better integrity and continuity and better conforms to the conditions of the actual atmospheric wind field.
Step 1: first, the radar detection data at the time t and the nearest time N before the time t are read, namely, the radar detection data at the times 2018-5-148:32, 8:39 and 8:45 (the actual volume sweep period has an error of ± 1 minute). Then, according to the basic reflectivity value in the data, adopting a cross correlation method to perform extrapolation calculation, calculating the displacement characteristic vector of the storm body represented by the radar data at each detection point, and recording the displacement characteristic vector of any detection point as
Figure GDA0003592351250000071
Wherein,
Figure GDA0003592351250000072
the effective elevation angle surfaces of the radar of the embodiment are 9, and the standard elevation angle degrees are respectively as follows: 0.5 °, 1.45 °, 2.4 °, 3.35 °, 4.3 °, 6.0 °, 9.9 °, 14.6 °, 19.5 °; gamma represents a distance library between a radar detection point and a radar central point, the distance library is a basic unit of detection resolution along the radial direction in radar echo signal processing, the length of the unit distance library also depends on radar technical parameters and a radar body scanning mode, and the length of the unit distance library of the radar in the embodiment is 1 km; omega represents the azimuth angle of the radar body sweep, the azimuth angle in the north direction is usually 0 degrees, the azimuth angle is increased in the clockwise direction, and the body sweep is 360 degrees for one circle; t represents the time corresponding to the processed radar data.
Figure GDA0003592351250000073
The physical meaning of (A) is as follows: using the radar detection data at 3 times of 8:32, 8:39 and 8:45 by taking 8:45 as a reference, and calculating a time period T of the radar volume sweep through an extrapolation algorithmit6min, from any elevation angle
Figure GDA0003592351250000074
The radar data represents the velocity of the storm body movement at any spatial position determined by any range bin γ and any azimuth ω, and the velocity is given by: distance bin number/6 min.
Step 2: shift feature vectors
Figure GDA0003592351250000075
According to the time period T of radar body scanningitLength L corresponding to unit distance library of 6minu1km, converted into a moving velocity vector
Figure GDA0003592351250000076
The calculation method comprises the following steps:
Figure GDA0003592351250000077
as a vector, the vector is a vector,
Figure GDA0003592351250000078
besides the value of the speed, the direction information is also provided, and the definition is not limited
Figure GDA0003592351250000079
The included angle between the pointed direction and the north direction is theta
Figure GDA00035923512500000710
When pointing to the north direction, θ is 0, and rotates clockwise, θ increases gradually.
And step 3: will move velocity vector
Figure GDA0003592351250000081
Is decomposed into U, V two components, where the U component is perpendicular to the current radial direction and the V component coincides with the current radial direction. U, V the calculation method of the two components is as follows:
Figure GDA0003592351250000082
wherein,
Figure GDA0003592351250000083
to represent
Figure GDA0003592351250000084
A component in the U axis;
Figure GDA0003592351250000085
to represent
Figure GDA0003592351250000086
The component at the V axis.
And 4, step 4: suppose in
Figure GDA0003592351250000087
Radar average radial velocity of position of
Figure GDA0003592351250000088
This value is one of the basic radar products and is directly available. When the temperature is higher than the set temperature
Figure GDA0003592351250000089
Then, the direction indicating the speed is directed to the center of the circle along the radial direction; on the contrary, when
Figure GDA00035923512500000810
The direction representing the velocity is radially away from the center of the circle. Will be provided with
Figure GDA00035923512500000811
As a reference speed, calculating
Figure GDA00035923512500000812
And
Figure GDA00035923512500000813
and enlarging or reducing at the magnification
Figure GDA00035923512500000814
To obtain a value of
Figure GDA00035923512500000815
The calculation method comprises the following steps:
Figure GDA00035923512500000816
and 5: component of U direction
Figure GDA00035923512500000817
Component with V direction
Figure GDA00035923512500000818
The motion vectors are combined to obtain the corrected motion velocity vector, which is recorded as
Figure GDA00035923512500000819
The calculation method comprises the following steps:
Figure GDA00035923512500000820
step 6: according to the steps 2 to 5, the movement velocity vector, namely the corrected movement velocity vector of the 3 times radar detection data corresponding to the time t on each elevation surface, each distance bank, each azimuth angle of the radar body scan can be calculated
Figure GDA00035923512500000821
And 7: according to the geographic position of the current radar, the longitude and latitude are (118.698 degrees, 32.191 degrees) and the effective radius of radar detection, any detection point in radar data is detected
Figure GDA00035923512500000822
And searching for the effective air-detecting report station in the coverage range. In this embodiment, there is only one sounding station, the station number is 25238, and the longitude and latitude are (118.8 °,32.0 °). The sounding report closest in time to the current radar data is data at 08 days, and the data records wind directions and wind speeds of 11 atmospheric pressure layers including 1000hPa, 925hPa, 850hPa, … and 100 hPa. Record of current sounding report is not defined as Vj(lot, lat, alt), where lot represents longitude, lat represents latitude, alt represents altitude, VhRepresenting the wind velocity vector at the corresponding spatial location.
In meteorology, for the analysis of high-altitude atmospheric conditions, it is common practice to use atmospheric pressure instead of altitude as a measure for different heights, and the atmospheric pressure and the altitude can be calculated by a meteorology method such as a 'pressure raising formula'. Because the current sounding report only has atmospheric pressure information and no altitude information from the ground, and no matter the atmospheric pressure is adopted in each step of the inventionThere is no essential difference in adopting altitude as a measure for different heights. Thus, step 7VhAlt in (lot, lat, alt) in units hPa, and the different height layers will be represented in the following steps using atmospheric pressure in units hPa.
And 8: for sounding empty report record Vh(lot, lat, alt) using a spatial interpolation algorithm to correct the respective elevation surface to obtain a corrected motion velocity vector
Figure GDA0003592351250000091
Interpolate to alt height and record the result as Vip(alt, γ', ω, t). Wherein, gamma 'represents the distance between the detection point and the radar central point in the horizontal direction, the unit is kilometer (km), and the value of gamma' can be selected from gamma,
Figure GDA0003592351250000092
And LuIs obtained by a mathematical geometric method.
There are many methods for spatial interpolation, and the calculation method of this embodiment is shown in fig. 3, taking alt equal to 925hPa as an example, V shown as dark black dot in a certain radial direction of the height layeripThe velocity vector at (925hPa, γ', ω, t) is defined by V at 1 elevation plane each of the upper and lower adjacent to the position spacew(0.5 °, γ, ω, t) and VwAnd (1.45 degrees, gamma, omega, t) is calculated by adopting an inverse distance weight algorithm.
And step 9: defining an influence radius R and a maximum influence radius RifSatisfy 0<r<RifR represents VipThe space distance between any grid point on (alt, gamma, omega, t) and the corresponding station of the sounding report can be determined by the longitude, latitude and V of the sounding reportipAnd (alt, gamma, omega, t) grid points are obtained by calculating the longitude and latitude by a geometric method. In this example, Rif＝50km。
Constructing an influence weight function Wif(r), which is an empirical function, subject to the condition: it is a decreasing function and 0<Wif(r)<1. In this example, WifThe (r) function is:
Figure GDA0003592351250000093
Wifthe functional image of (r) is shown in fig. 4.
Because the sounding report is the observed real wind field information, the sounding report has higher reliability compared with various wind field information of radar product inversion. According to the sounding report record Vh(lot, lat, alt) velocity vector information for V in step 8ip(alt, gamma', omega, t) to make deviation correction, the calculation method is:
Vbs(alt,γ′,ω,t)＝(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)
due to Vbs(alt,γ′,ω,t)、Vip(alt, γ', ω, t) and VhThe (lot, lat, alt) are all vectors, so the calculation process should conform to the calculation rule of vector addition and subtraction, such as parallelogram or orthogonal decomposition.
Thereby completing a probing report station data VhCorrecting deviation of each height alt in (lot, lat, alt) to obtain a standard height layer wind field vector V after deviation correctionbs(alt, γ', ω, t); and further completing deviation correction of all the data of the station detection data, and obtaining wind field vector information of a plurality of standard height layers.
Theoretically, the more effective sounding sites and effective sounding wind field records are within the radar detection range, the more Vbs(alt, gamma', omega, t) can reflect the wind field characteristics of each high-altitude layer more truly; on the contrary, if the sounding site is too sparse or the sounding wind field records are too few, the radius R is influenced to the maximum extentifIn this case, the calculation process of step 9 may not be performed in some regions of the radar detection range, which also takes into account the discontinuity of the wind field on a large scale and avoids "overcorrection".
And obtaining wind field vector information of a plurality of standard height layers in the radar detection range through the calculation of the steps. After the information is graphically displayed, auxiliary decision can be provided for professional personnel to carry out weather diagnosis and analysis; this information can also be entered into the numerical model as an initial field to improve the accuracy of the weather pattern forecast.

Claims (7)

1. A high altitude wind inversion method based on Doppler weather radar data comprises the following steps:
step 1: supposing that high-altitude wind at a certain time t is to be inverted, firstly, a radar is selected, the radar data of the radar at the time t and N times of radar data closest to the time t before are utilized, an extrapolation algorithm is adopted to obtain displacement characteristic vectors of a storm body represented by the radar data at each detection point, and the displacement characteristic vector of any detection point is recorded as
Figure FDA0003592351240000011
Wherein N is an integer more than 2 and less than or equal to 10;
Figure FDA0003592351240000012
the method comprises the following steps of representing a certain elevation surface of radar body scanning, representing a distance library between a radar detection point and a radar central point, representing an azimuth angle of the radar body scanning by omega, and representing time corresponding to processed radar data by t;
step 2: shifting the feature vector in the previous step
Figure FDA0003592351240000013
According to the time period T of the radar body scanitLength L corresponding to unit distance libraryuConverted into a moving velocity vector according to the following formula
Figure FDA0003592351240000014
Figure FDA0003592351240000015
Figure FDA0003592351240000016
As a vector, a value representing the speedAnd information of direction, definition
Figure FDA0003592351240000017
The included angle between the pointed direction and the north direction is theta
Figure FDA0003592351240000018
When pointing to the north direction, theta is 0, and rotates clockwise, and theta increases gradually;
and step 3: moving velocity vector in the last step
Figure FDA0003592351240000019
Decomposed into U, V two components, where the U component is perpendicular to the current radial direction and the V component coincides with the current radial direction, assuming
Figure FDA00035923512400000110
To represent
Figure FDA00035923512400000111
A component in the U axis;
Figure FDA00035923512400000112
to represent
Figure FDA00035923512400000113
A component at the V axis;
and 4, step 4: assuming any probe point in the radar data
Figure FDA00035923512400000114
Radar average radial velocity of position of
Figure FDA00035923512400000115
When in use
Figure FDA00035923512400000116
Then, the direction indicating the speed is directed to the center of the circle along the radial direction; on the contrary, when
Figure FDA00035923512400000117
Then, the direction representing the velocity is radially away from the center of the circle; will be provided with
Figure FDA00035923512400000118
As a reference speed, calculating
Figure FDA00035923512400000119
And
Figure FDA00035923512400000120
and enlarging or reducing at the magnification
Figure FDA00035923512400000121
To obtain a value of
Figure FDA00035923512400000122
The calculation formula is as follows:
Figure FDA00035923512400000123
and 5: component of U direction
Figure FDA00035923512400000124
Component with V direction
Figure FDA00035923512400000125
The motion vectors are combined according to the following formula, and the corrected motion velocity vector is obtained through calculation and is recorded as
Figure FDA00035923512400000126
Figure FDA00035923512400000127
Step 6: according to the steps 2 to 5, displacement characteristic vectors of each detection point of the radar are subjected to
Figure FDA00035923512400000128
Calculating one by one to obtain the corrected moving velocity vectors of each elevation angle surface, each distance base and each azimuth angle swept by the radar body at the time t, namely
Figure FDA00035923512400000129
And 7: for any detection point in radar data
Figure FDA00035923512400000130
Searching a sounding report station closest to the geographical position of the sounding report station, and reading sounding report station data of which the detection time is closest to the time t in the sounding report station, wherein the data records the wind direction and the wind speed of a certain latitude position on the earth and vertically upwards at a plurality of height layers or a plurality of atmospheric pressure layers in a standard level; obtaining sounding station reporting data corresponding to each detection point, and recording the data as Vh(lot, lat, alt); wherein, lot represents longitude, lat represents latitude, alt represents altitude;
and 8: the station data V is reported in the above stepsh(lot, lat, alt) each alt is a reference height, and the corrected moving velocity vector of each elevation surface of the radar is used
Figure FDA0003592351240000021
Converting into moving velocity vector at height corresponding to each alt by using spatial interpolation algorithm, and recording as
Figure FDA0003592351240000022
Gamma 'represents the distance between the detection point and the radar central point in the horizontal direction, and the value of gamma' is formed by gamma,
Figure FDA0003592351240000023
And LuBy passingSolving by a triangular geometric method;
and step 9: defining an influence radius R and a maximum influence radius RifR is more than 0ifR represents VipThe spatial distance between any point (alt, gamma', omega, t) and the sounding site can be determined by the longitude, latitude and V of the sounding siteipThe longitude and latitude of the point in (alt, gamma', omega, t) are obtained by calculation through a geometric method; vipThe longitude and the longitude of each point (alt, gamma ', omega, t) are calculated according to the geographic position of the radar and gamma';
constructing an influence weight function Wif(r) which is a decreasing function subject to the following conditions: w is more than 0if(r) < 1; reporting station data V according to soundingh(lot, lat, alt), for V in step 8ip(alt, γ', ω, t) deviation correction is performed by the following formula:
Vbs(alt，γ′，ω，t)＝(1-Wif(r))×Vip(alt，γ′，ω，t)+Wif(r)×Vh(lot，lat，alt)。
2. the Doppler weather radar data-based high altitude wind inversion method according to claim 1, wherein in step 1, the data of the N times of radar data are in the form of basic reflectivity data and average radial velocity data.
3. The Doppler weather radar data-based high altitude wind inversion method according to claim 2, wherein in the step 1, the extrapolation algorithm adopts a maximum correlation method or a cross correlation method.
4. The high altitude wind inversion method based on Doppler weather radar data as claimed in claim 3, wherein in step 7, if the record of the sounding station data only has atmospheric pressure information and no altitude information from the ground, the recording is obtained by conversion through a high altitude formula meteorology method.
5. The high altitude wind inversion method based on Doppler weather radar data as recited in any one of claims 1 to 4, wherein in step 7, the sounding station data is interpolated to a standard level by using an interpolation correlation algorithm and a non-standard level.
6. The Doppler weather radar data-based high altitude wind inversion method as claimed in claim 5, wherein in step 8, the corrected moving velocity vector
Figure FDA0003592351240000024
The interpolated alt height is calculated by the inverse distance weighting algorithm through the moving speed vectors of the upper and lower 1 elevation surfaces adjacent to the position space.
7. The Doppler weather radar data-based high altitude wind inversion method as claimed in claim 6, wherein V is obtained in step 9ip(alt, γ', ω, t) can be solved using parallelogram or orthonormal decomposition methods.
CN201810826385.5A 2018-07-25 2018-07-25 High-altitude wind inversion method based on Doppler weather radar data Active CN109100723B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810826385.5A CN109100723B (en) 2018-07-25 2018-07-25 High-altitude wind inversion method based on Doppler weather radar data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810826385.5A CN109100723B (en) 2018-07-25 2018-07-25 High-altitude wind inversion method based on Doppler weather radar data

Publications (2)

Publication Number Publication Date
CN109100723A CN109100723A (en) 2018-12-28
CN109100723B true CN109100723B (en) 2022-05-27

Family

ID=64847361

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810826385.5A Active CN109100723B (en) 2018-07-25 2018-07-25 High-altitude wind inversion method based on Doppler weather radar data

Country Status (1)

Country Link
CN (1) CN109100723B (en)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109782019B (en) * 2019-03-15 2020-04-24 中国科学技术大学 Method and device for measuring two-dimensional movement speed of atmospheric pollutants
CN110288856A (en) * 2019-06-21 2019-09-27 中国民用航空总局第二研究所 The Scheduled Flight monitoring system and method for fine forecast based on wind
CN112946657B (en) * 2021-02-03 2023-10-27 南京信息工程大学 Method for identifying ground wind field in strong convection weather
CN113009490B (en) * 2021-02-20 2022-10-21 江苏省气象台 Radar three-dimensional wind field inversion method based on high-resolution mode dynamic constraint
CN113311436B (en) * 2021-04-30 2022-07-12 中国人民解放军国防科技大学 Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform
CN114063465B (en) * 2021-09-22 2024-05-03 中国航空工业集团公司西安飞机设计研究所 View jitter elimination method and view node of distributed countermeasure simulation system
CN115825894B (en) * 2022-11-17 2023-08-18 中国能源建设集团广东省电力设计研究院有限公司 Method, device, terminal equipment and medium for determining wind energy capturing position
CN116381692B (en) * 2023-04-18 2024-02-23 南京市气象台 Precipitation phase state identification QPE algorithm based on X-band dual-polarization radar

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080169975A1 (en) * 2007-01-12 2008-07-17 Young Paul Yee Process for generating spatially continuous wind profiles from wind profiler measurements
CN104035096B (en) * 2014-06-06 2017-05-03 南京大学 Vertical wind profile nonlinear inversion method based on Doppler weather radar
CN104503000B (en) * 2014-12-15 2017-10-27 深圳航天东方红海特卫星有限公司 A kind of sonde wind measuring system and wind detection method
CN106199605B (en) * 2016-07-06 2018-08-24 西南技术物理研究所 Wind field error correcting method
CN107741587A (en) * 2017-09-15 2018-02-27 北京无线电测量研究所 A kind of meteorological detection method and system of dimensional wind
CN108107434B (en) * 2017-12-26 2020-06-09 厦门市气象灾害防御技术中心(海峡气象开放实验室厦门市避雷监测技术中心) Regional three-dimensional wind field picture splicing method based on double-Doppler radar inversion

Also Published As

Publication number Publication date
CN109100723A (en) 2018-12-28

Similar Documents

Publication Publication Date Title
CN109100723B (en) High-altitude wind inversion method based on Doppler weather radar data
CN110058236B (en) InSAR and GNSS weighting method oriented to three-dimensional surface deformation estimation
CN107843895B (en) A kind of Dual-Doppler weather radar dimensional wind inversion method
US7200491B1 (en) System for producing high-resolution, real-time synthetic meteorological conditions for a specified location
CN105182339A (en) Method for correcting environmental influences at slope deformation monitoring on the basis of corner reflector
CN112731564B (en) Intelligent thunder forecasting method based on Doppler weather radar data
CN110879428B (en) Method for measuring rain through networking in overlapped area between radars in same wave band under complex terrain condition
CN113325422B (en) Space-based rain radar target positioning and rainfall information three-dimensional processing method and system
JPWO2018168165A1 (en) Weather forecasting device, weather forecasting method, and program
CN113009531A (en) Small-scale high-precision low-altitude troposphere atmosphere refractive index model
CN111487621A (en) Sea surface flow field inversion method based on radar image and electronic equipment
Feng et al. Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation
CN116485857B (en) High-time-resolution glacier thickness inversion method based on multi-source remote sensing data
CN116500648A (en) Wind profile inversion method for foundation laser radar target area
CN114252875B (en) High-precision meshing method for imaging altitude data
CN111505626B (en) Method for measuring two-dimensional terrain gradient by using bottom view differential interference
Idris et al. Comparison of retracked coastal altimetry sea levels against high frequency radar on the continental shelf of the Great Barrier Reef, Australia
CN113568008A (en) Sunflower-8 satellite-based minute-level precipitation real-time inversion estimation method
Yang et al. Can CINRAD Radar with VCP-21 Mode Capture the Accumulated Rainfall Pattern and Intensity of Fast-moving Storms?
CN117348116B (en) Local strong convection weather forecast method based on non-sharp inclined pressure belt
CN117054984A (en) Method for compensating detection capability of airborne radar by utilizing lightning data
CN116203591B (en) Method for jointly estimating high-precision ionosphere VTEC of China area based on multiple stations
Liu et al. Verification of Three-Dimensional Gridding Algorithm for Weather Radar Reflectivity
Xing et al. Geometric Correction of Spaceborne SAR Image Based on DEM and ICESat Database
GCR InSAR Analysis over the Eastern Coastline of Florida

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant