CN109100723A - Upper-level winds inversion method based on Doppler radar data - Google Patents
Upper-level winds inversion method based on Doppler radar data Download PDFInfo
- Publication number
- CN109100723A CN109100723A CN201810826385.5A CN201810826385A CN109100723A CN 109100723 A CN109100723 A CN 109100723A CN 201810826385 A CN201810826385 A CN 201810826385A CN 109100723 A CN109100723 A CN 109100723A
- Authority
- CN
- China
- Prior art keywords
- radar
- data
- alt
- sounding
- point
- 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.)
- Granted
Links
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 present invention proposes a kind of upper-level winds inversion method based on Doppler radar data, and this method carries out multi-level high-altitude wind field using the extrapolation product and radial velocity product and sounding report data of Doppler radar and merges.The present invention utilizes the extrapolating results of radar return image, converts the velocity vector field with U, V component for the average radial velocity field in each elevation angle face, improves the readability of high-altitude wind field subjective analysis.Radar wind field is suitably corrected using Sounding Data, modified dynamics is successively decreased with the increase of distance between radar mesh point and sounding station, so that wind field has preferable globality and continuity, and is more in line with the situation of real atmosphere wind field.The method that the present invention compares some multivariate data fusions, the computational complexity of inventive algorithm is moderate, and the timeliness of the secondary wind field output of radar can be improved.
Description
Technical field
The invention belongs to field of earth sciences more particularly to it is a kind of using the detection informations such as Doppler radar carry out it is more
The method of level upper-level winds inverting.
Background technique
High-altitude wind field is characterized in carrying out one of the important object of Synoptic Diagnostic, in meteorological, military and aerospace etc.
Play the role of vital in field.Currently, obtain earth overhead atmospheric wind approach have very much, such as utilize sounding balloon,
The equipment such as Doppler radar, wind profile radar and meteorological satellite are observed or detect.These equipment and technological means
It is each advantageous, wherein Doppler radar (hereinafter referred to as radar) has space exploration high resolution, data renewal frequency fast
It the features such as strong with actual effect, is widely studied and is answered in terms of local weather information monitors and closes in short-term
With.
Average radial velocity is one of basic product of radar, and it is substantially former that the acquisition of the speed, which is with Doppler effect,
Reason, the wind field information characterized is the point centered on radar, a certain velocity component radially swept along radar body, and non-real
Real wind vector.In order to obtain in radar detection coverage area different height layer practical wind field information, it is many for many years to learn
Person proposes respective method, has plenty of the multiple product using radar itself, such as reflectivity factor and average radial velocity, passes through
Mathematical method constructs constraint equation, is finally inversed by wind field information；Also have plenty of using and combine a variety of detection informations, such as Wind outline thunder
Wind field information is calculated up to the multivariate datas such as wind, satellite derived wind are surveyed.Due to the complexity and prior art hand of weather system
The limitation of section, in practical business application, for the inverting of high-altitude fact wind field, there are still certain defect and deficiencies.Example
Such as, Wind-field Retrieval is carried out merely with the reflectivity factor and average radial velocity of single portion's radar, method is relatively easy, but by radar
Clutter, electromagnetic wave beam with detection range increase become deviation, bending etc. adverse effect so that the wind field accuracy being finally inversed by compared with
It is low.In another example fusion radar radial velocity and the multivariate datas such as Cloud-motion wind or steam wind-guiding, such methods are primarily limited to difference
The problems such as spatial resolution of data is inconsistent, and space lattice point correspondence is inconsistent and observation time is inconsistent, so that
Low when high when calculated wind field accuracy rate, stability is poor.
Summary of the invention
In order to overcome the shortcomings of the above method, the present invention proposes a kind of extrapolation product and diameter using Doppler radar
The method that multi-level high-altitude wind field merges is carried out to Velocity products and sounding report data.
Specific technical solution of the present invention is as follows:
A kind of upper-level winds inversion method based on Doppler radar data, method includes the following steps:
Step 1: assuming that preparing the upper-level winds of a certain moment t of inverting, a radar is selected first, when using portion radar t
It carves and before closest to the n times radar data of t moment, storm body that radar data is characterized is obtained in each spy using extrapolation algorithm
The displacement characteristic vector of the displacement characteristic vector of measuring point, any sensing point is denoted asWherein, N takes the whole of 2 < N≤10
Number；Indicate a certain elevation angle face that radar body is swept, γ indicates the distance between radar detection point and radar center point library, and ω is indicated
The azimuth that radar body is swept, t indicate that handled radar data corresponds to the time；
Step 2: by displacement characteristic vector in upper stepThe period of time T swept according to radar bodyit, unit away from
Length L corresponding from libraryu, it is converted into moving velocity vector as the following formula
As vector, the numerical value of speed and the information in direction are characterized, is definedSignified side
It is θ to the angle with direct north, whenWhen being directed toward direct north, θ=0 is rotated in a clockwise direction, and θ is gradually
Increase；
Step 3: by moving velocity vector in upper stepBe decomposed into two components of U, V, wherein U component with work as
Preceding radial direction is vertical, and V component is overlapped with current radial direction, it is assumed thatIt indicatesIn the component of U axis；It indicatesIn the component of V axis；
Step 4: assuming that any sensing point in radar dataThe radar average radial velocity of position isWhenWhen, indicate that the direction of the speed is along being radially directed towards the center of circle；Conversely, working asWhen, the direction for indicating the speed is along radial away from the center of circle；It willAs benchmark speed
Degree calculatesWithMultiplying power, and zoomed in or out with this multiplying powerNumerical value,
It obtainsCalculation formula is as follows:
Step 5: by the component in the direction UWith the component in the direction VIt merges as the following formula,
Revised moving velocity vector is calculated, is denoted as
Step 6: according to step 2~step 5, to the displacement characteristic vector of each sensing point of radarBy
One is calculated, and obtains t moment in each elevation angle face that radar body is swept, each range bin, each azimuthal revised shifting
Dynamic velocity vector, i.e.,
Step 7: to sensing point any in radar dataIt searches and the most similar sounding in its geographical location
It calls out the stops a little, reads the sounding sounding of detection time closest to t moment in a little of calling out the stops and call out the stops point data, the data record earth
A certain longitude and latitude position is gone up, vertically upward the wind direction and wind speed of multiple height layers or multiple barosphere standard levels；Obtain with often
The corresponding sounding of a sensing point is called out the stops point data, and V is recorded ash(lot,lat,alt)；
Step 8: the above step sounding is called out the stops point data VhHeight on the basis of each alt in (lot, lat, alt), by radar
Moving velocity vector after the amendment in each elevation angle faceUsing spatial interpolation algorithm, it is right in each alt institute to be converted to
The moving velocity vector for answering height, is denoted asγ ' expression sensing point and radar center point in the horizontal direction
The distance between, the value of γ ' by γ,And LuIt is acquired by triangle method of geometry；
Step 9: defining an a radius of influence r and biggest impact radius Rif, meet 0 < r < Rif, r expression Vip(alt,
γ ', ω, t) in the space length called out the stops between a little of any point and sounding, the distance can be called out the stops by sounding longitude a little, latitude and
VipWarp, the latitude of point in (alt, γ ', ω, t) are calculated by method of geometry；Vip(alt, γ ', ω, t) each point
Longitude and longitude by portion's radar geographical location and γ ' be calculated；
Construct a weighing factor function Wif(r), which is the decreasing function that must meet the following conditions: 0 < Wif(r)<
1；It is called out the stops point data V according to soundingh(lot, lat, alt), to the V in step 8ip(alt, γ ', ω, t) is carried out inclined by following formula
Difference is corrected:
Vbs(alt, γ ', ω, t)=(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,
alt)。
Further design of the invention is:
In step 1, the data mode of the n times radar data is baseis reflectivity data and average radial velocity data.
In step 1, the extrapolation algorithm uses Maximum correlation method or cross-correlation technique.
In step 7, if the sounding is called out the stops, the record of point data only has atmospheric pressure information, the height letter not apart from ground
Breath, converts to obtain by pressure-height formula meteorology method.
In step 7, point data is called out the stops using non-standard level for sounding, using interpolation related algorithm, is inside inserted into mark
Quasi- level.
In step 8, moving velocity vector after amendmentBe interpolated into alt height by with the locational space phase
The moving velocity vector in adjacent upper and lower each 1 elevation angle face is calculated using anti-distance weighting algorithm.
In step 9, Vip(alt, γ ', ω, t) parallelogram can be used or proper orthogonal decomposition solves.
The present invention has the advantage that as follows compared with prior art:
1, the present invention utilizes the extrapolating results of radar return image, converts the average radial velocity field in each elevation angle face to
Velocity vector field with U, V component improves the readability of high-altitude wind field subjective analysis.
2, the present invention suitably corrects radar wind field using Sounding Data, and modified dynamics is with radar mesh point and visits
The increase of distance between sky station and successively decrease so that wind field has preferable globality and continuity, and be more in line with real atmosphere
The situation of wind field.
3, the method that the present invention compares some multivariate data fusions, the computational complexity of inventive algorithm is moderate, can be improved
The timeliness of radar afterproduct (wind field) output.
Detailed description of the invention
Fig. 1 is the main control flow chart of the embodiment of the present invention one.
Fig. 2 is moving velocity vector in embodiment oneBe decomposed into two components of U, V, wherein U component with work as
Preceding radial direction is vertical, and V component is overlapped with current radial direction.
Fig. 3 is in one step 8 of application example, to " velocity vector in each elevation angle face of radar being interpolated into an isobaris surface
Calculating process " auxiliary description.
Fig. 4 is the W of construction in one step 9 of application exampleif(r) functional image of function.
Specific embodiment
Technical solution of the present invention is further illustrated with reference to the accompanying drawing.
Embodiment one:
As shown in Figure 1, the present invention is based on the upper-level winds inversion method of Doppler radar data, this method includes following
Step:
Step 1: assuming that preparing the upper-level winds of a certain moment t of inverting, a radar is selected first, when using portion radar t
It carves and (data mode is baseis reflectivity data and average radial velocity number closest to the n times radar data of t moment before
According to), using extrapolation algorithm, such as Maximum correlation method, cross-correlation technique, storm body that radar data is characterized is calculated in each spy
The displacement characteristic vector of the displacement characteristic vector of measuring point, any sensing point is denoted asWherein, N takes the whole of 2 < N≤10
Number；The elevation scale in a certain elevation angle face that expression radar body is swept, the number in elevation angle face and each elevation angle face depends on radar skill
The mode that art parameter and radar body are swept；γ indicates the distance between radar detection point and radar center point library, and range bin is radar
The basic unit of radially upper detection resolution, the length in unit distance library also depend on Radar Technology ginseng in echo signal processing
The mode that several and radar body is swept；ω indicates the azimuth that radar body is swept, usually with the azimuth of direct north for 0 °, along clockwise
Direction be incremented by, body sweep one week be 360 °；T indicates that handled radar data corresponds to the time.Physical significance are as follows:
On the basis of any time t, using t, t-1, t-2 ..., the radar observation data at t-N+1 moment, calculated by extrapolation algorithm
To the period of time T swept in radar bodyit, by any elevation angleDetermined by any range bin γ and any azimuth ω
On spatial position, the mobile speed of radar return, the unit of the speed are as follows: range bin number/body flyback time period.
Step 2: by displacement characteristic vector in upper stepThe period of time T swept according to radar bodyit, unit away from
Length L corresponding from libraryu, it is converted into moving velocity vectorWherein, TitUnit be minute (min), LuList
Position is kilometer (km),Unit be meter per second (m/s), calculation method are as follows:
As vector,Other than having the numerical values recited of characterization speed, there are also the information for characterizing direction, there is no harm in
DefinitionThe angle in signified direction and direct north is θ, whenWhen being directed toward direct north, θ=0,
It is rotated in a clockwise direction, θ is incrementally increased.
Step 3: by moving velocity vector in upper stepBe decomposed into two components of U, V, wherein U component with work as
Preceding radial direction is vertical, and V component is overlapped with current radial direction, as shown in Figure 2.U, two components of V can be calculated using the methods of Orthogonal Decomposition
It obtains, it is assumed thatIt indicatesIn the component of U axis；It indicatesIn V
The component of axis.
Step 4: assuming that any sensing point in radar dataThe radar average radial velocity of position isThe value is one of basic product of radar, can be directly obtained.WhenWhen, indicate the speed
The direction of degree is the radar center point in the center of circle, that is, as shown in Figure 2 along being radially directed towards the center of circle；Conversely, working as
When, the direction for indicating the speed is along radial away from the center of circle.It willAs reference speed, calculateWithMultiplying power, and zoomed in or out with this multiplying powerNumerical value, obtainCalculation method are as follows:
Step 5: by the component in the direction UWith the component in the direction VIt merges, calculates
To revised moving velocity vector, it is denoted asCalculation method are as follows:
Step 6: according to step 2~step 5, to the displacement characteristic vector of each sensing point of radarOne by one
It is calculated, obtains t moment in each elevation angle face that radar body is swept, each range bin, each azimuthal revised movement
Velocity vector, i.e.,
Step 7: to sensing point any in radar dataIt searches and the most similar sounding in its geographical location
It calls out the stops a little, reads the sounding sounding of detection time closest to t moment in a little of calling out the stops and call out the stops point data, the data record earth
A certain longitude and latitude position is gone up, vertically upward the wind direction and wind speed of multiple height layers or multiple barosphere standard levels；Obtain with often
The corresponding sounding of a sensing point is called out the stops point data, and V is recorded ash(lot,lat,alt)；
It, can be by the way that " pressure is high public without the elevation information apart from ground if the data record only has atmospheric pressure information
The meteorology method such as formula " converts to obtain.For Vh(lot, lat, alt) wherein, lot indicates that longitude, lat indicate longitude, alt table
Show height above sea level, VhIndicate the wind vector of additional space position；The height layer for a little detecting and recording every time in general, sounding is called out the stops
Or barosphere, quantity are not fixed, every part of sounding point data of calling out the stops has recorded how many groups of probe values, with regard to how many alt value, no
Harm is defined as alt1, alt2 ..., altn, alt ∈ { alt1, alt2 ..., altn }；Every part of sounding is called out the stops point data, longitude and latitude
Degree coordinate position is fixed and invariable；Different soundings are called out the stops a little, and latitude and longitude coordinates position is not often identical.
It should be added that the sounding report used in meteorological department of China business, the height layer recorded or air pressure
Interpolation related algorithm can be used for non-standard level in several standard levels that layer is usually fixed, and is inside inserted into index bed.
Following steps are the calculating for standard level, the case where for non-standard level, are also suitable after interpolation calculates.
Step 8: the above step sounding is called out the stops point data VhHeight on the basis of each alt in (lot, lat, alt), by radar
Moving velocity vector after the amendment in each elevation angle faceUsing spatial interpolation algorithm, it is right in each alt institute to be converted to
The moving velocity vector for answering height, is denoted as Vip(alt,γ′,ω,t).Wherein, spatial interpolation algorithm can be used anti-distance weighting and calculate
Method, γ ' expression the distance between sensing point and radar center point in the horizontal direction, unit are kilometer (km), the value of γ ' by
γ、And LuIt is acquired by triangle method of geometry；
Step 9: sounding calls out the stops longitude a little and latitude is known, VipThe longitude and warp of (alt, γ ', ω, t) each point
Degree by portion's radar geographical location and γ ' be calculated；
Define a radius of influence r and a biggest impact radius Rif, meet 0 < r < Rif, r expression Vip(alt,γ′,ω,
T) any point and sounding are called out the stops the space length between a little in, longitude, latitude and the V which can be called out the stops a little by soundingip(alt,
γ ', ω, t) in warp, the latitude of point be calculated by mathematical method；
Construct a weighing factor function Wif(r), which is an empirical function, must meet condition: it is one and passs
Subtraction function and 0 < Wif(r)<1；Since sounding report record is the true wind field information observed, it is finally inversed by compared to Radar Products
Various wind field information, have higher confidence level.It is called out the stops point data V according to soundingh(lot, lat, alt), in step 8
Vip(alt, γ ', ω, t) carries out deviation and corrects, calculation method are as follows:
Vbs(alt, γ ', ω, t)=(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)
Due to Vbs(alt,γ′,ω,t)、Vip(alt, γ ', ω, t0 and Vh(lot, lat, alt) is all vector, therefore,
Calculating process should meet the computing rule of vector plus-minus, such as use parallelogram or proper orthogonal decomposition.
Thus a sounding is completed to call out the stops point data VhThe deviation of each of (lot, lat, alt) height alt is corrected,
It obtains through the revised calibrated altitude layer winds retrieval V of deviationbs(alt,γ′,ω,t)；All soundings are further completed to call out the stops a little
The deviation of data is corrected, and the winds retrieval information of multiple calibrated altitude layers is obtained.
Theoretically, in radar coverage, effective sounding website and effective sounding wind field record more, Vbs(alt,
γ ', ω, t) it more can truly reflect the Characteristics of Wind Field of each height layer in high-altitude；, whereas if sounding website is excessively sparse or sounding
Wind field record it is very few, due to biggest impact radius RifLimitation, some regions of radar coverage then may be without step 9
Calculating process, this also allows for discontinuity of the wind field in large scale, is avoided " exaggerated correction ".
The winds retrieval information of multiple calibrated altitude layers in radar coverage can be obtained in calculating through the above steps.
These information can be engaged in Synoptic Diagnostic for professional and provide aid decision after graphical representation；These can also be believed
Breath is input in numerical model as initial fields, to improve synoptic model accuracy of the forecast.
Application example one:
The Doppler radar data conduct of the present embodiment selection t=" 45 divide when 14 days 8 May in 2018 " In Nanjing
Objective for implementation, the radar are the radars of a S-band, sweep mode using VCP21 body, and body is swept the period and is generally 6 minutes.This implementation
The final purpose of example is the atmospheric wind of multiple calibrated altitude layers to be calculated by following step, and it is preferable that wind field is had
Globality and continuity and the situation for being more in line with real atmosphere wind field.
Step 1: firstly, reading radar detection time t moment and before closest to N=3 radar detection number of t moment
According to that is, (practical body sweeps the mistake of period presence ± 1 minute to the GPR Detection Data at 2018-5-148:32,8:39 and 8:45 moment
Difference).Then, according to the baseis reflectivity value in data, extrapolation calculating is carried out using cross-correlation technique, calculates radar data institute
In the displacement characteristic vector of each sensing point, the displacement characteristic vector of any sensing point is denoted as the storm body of characterization
Wherein,Indicate a certain elevation angle face that radar body is swept, effective elevation angle face quantity of the present embodiment radar is 9, standard elevation scale
It is respectively as follows: 0.5 °, 1.45 °, 2.4 °, 3.35 °, 4.3 °, 6.0 °, 9.9 °, 14.6 °, 19.5 °；γ indicates radar detection point and thunder
Up to the distance between central point library, range bin is the basic unit of radially upper detection resolution in radar echo signal processing,
The mode that the length in unit distance library also depends on Radar Technology parameter and radar body is swept, the unit distance library of the present embodiment radar
Length is 1km；ω indicates the azimuth that radar body is swept, to be usually incremented by along clockwise direction with the azimuth of direct north for 0 °,
Body sweeps one week as 360 °；T indicates that handled radar data corresponds to the time.Physical significance are as follows: be with 8:45
Benchmark is calculated by extrapolation algorithm in radar body using the radar observation data at this 3 moment of 8:32,8:39 and 8:45
The period of time T sweptit=6min, by any elevation angleSpace determined by any range bin γ and any azimuth ω
On position, the mobile speed of the storm body that radar data is characterized, the unit of the speed are as follows: range bin number/6min.
Step 2: by displacement characteristic vectorThe period of time T swept according to radar bodyit=6min, unit away from
Length L corresponding from libraryu=1km, is converted into moving velocity vectorCalculation method are as follows:
As vector,Other than having the numerical values recited of characterization speed, there are also the information for characterizing direction, no
Harm definitionThe angle in signified direction and direct north is θ, whenWhen being directed toward direct north, θ=
0, it is rotated in a clockwise direction, θ is incrementally increased.
Step 3: by moving velocity vectorIt is decomposed into two components of U, V, wherein U component and current radial direction
Vertically, V component is overlapped with current radial direction.U, the calculation method of two components of V are as follows:
Wherein,It indicatesIn the component of U axis；It indicates
In the component of V axis.
Step 4: assuming thatThe radar average radial velocity of position isThe value is radar
One of basic product can directly obtain.WhenWhen, indicate that the direction of the speed is along being radially directed towards circle
The heart；Conversely, working asWhen, the direction for indicating the speed is along radial away from the center of circle.It willMake
For benchmark speed, calculateWithMultiplying power, and zoomed in or out with this multiplying power
Numerical value, obtainCalculation method are as follows:
Step 5: by the component in the direction UWith the component in the direction VIt merges, is repaired
Moving velocity vector after just, is denoted asCalculation method are as follows:
Step 6: according to step 2 to step 5,3 GPR Detection Datas corresponding to t moment can be calculated in radar body
On each elevation angle face swept, each range bin, each azimuth, the moving velocity vector after being corrected, i.e.,
Step 7: according to the geographical location of current radar, longitude and latitude is (118.698 °, 32.191 °) and radar detection have
Radius is imitated, to sensing point any in radar dataEffective sounding in the coverage area is searched to call out the stops a little.This implementation
In example, only a sounding is called out the stops a little, and website number is 25238, and longitude and latitude is (118.8 °, 32.0 °).When with current radar data
Between on most similar sounding report be on the same day 08 when data, data record 1000hPa, 925hPa, 850hPa ...,
The wind direction and wind speed of 100hPa totally 11 atmosphere pressurized layers.Current sounding report data might as well be defined is recorded as Vh(lot,lat,
Alt), wherein lot indicates that longitude, lat indicate that longitude, alt indicate height above sea level, VhIndicate the wind speed arrow of additional space position
Amount.
In meteorology, analysis for upper atmosphere situation, ordinary practice use atmospheric pressure replace height above sea level as pair
The measurement of different height, and can be calculated by the meteorology method such as " pressure-height formula " between atmospheric pressure and height above sea level.
Since current sounding report record only has atmospheric pressure information, not apart from the elevation information on ground, also, in each step of the invention
In rapid, no matter using atmospheric pressure or using height above sea level as the measurement to different height, essential distinction is had no.Therefore,
Step 7VhAlt in (lot, lat, alt), unit hPa, and the atmospheric pressure that unit is hPa will be used in following steps
Indicate different height layers.
Step 8: V is recorded to sounding reporthEach of (lot, lat, alt) alt will be each using spatial interpolation algorithm
Moving velocity vector after the amendment in a elevation angle faceIt is interpolated into alt height, is as a result denoted as Vip(alt,γ′,ω,
t0.Wherein, in the horizontal direction, the distance between sensing point and radar center point, unit are kilometer (km), γ's ' for γ ' expression
Value can by γ,And LuIt is acquired by mathematics method of geometry.
The method of space interpolation has very much, and the method that the present embodiment calculates is somebody's turn to do as shown in figure 3, by taking alt=925hPa as an example
The a certain radially V as shown in pitch black color dot of height layeripVelocity vector at (925hPa, γ ', ω, t) is by empty with the position
Between adjacent upper and lower each 1 elevation angle face Vw(0.5 °, γ, ω, t) and Vw(1.45 °, γ, ω, t use anti-distance weighting algorithm
Calculating acquires.
Step 9: defining an a radius of influence r and biggest impact radius Rif, meet 0 < r < Rif, r expression Vip(alt,
γ, ω, t) on any one mesh point website corresponding with sounding report space length, the distance can by sounding report warp, latitude and
VipThe warp, latitude of (alt, γ, ω, t) mesh point are calculated by method of geometry.In the present embodiment, Rif=50km.
Construct a weighing factor function Wif(r), which is an empirical function, must meet condition: it is one and passs
Subtraction function and 0 < Wif(r)<1.In the present embodiment, Wif(r) function are as follows:
Wif(r) functional image is as shown in Figure 4.
Since sounding report record is the true wind field information observed, the various wind fields being finally inversed by compared to Radar Products
Information has higher confidence level.V is recorded according to sounding reporthThe velocity vector information of (lot, lat, alt), in step 8
Vip(alt, γ ', ω, t) carries out deviation and corrects, calculation method are as follows:
Vbs(alt, γ ', ω, t)=(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)
Due to Vbs(alt,γ′,ω,t)、Vip(alt, γ ', ω, t) and Vh(lot, lat, alt) is all vector, therefore,
Calculating process should meet the computing rule of vector plus-minus, such as use parallelogram or proper orthogonal decomposition.
Thus a sounding is completed to call out the stops point data VhThe deviation of each of (lot, lat, alt) height alt is corrected,
It obtains through the revised calibrated altitude layer winds retrieval V of deviationbs(alt,γ′,ω,t)；All soundings are further completed to call out the stops a little
The deviation of data is corrected, and the winds retrieval information of multiple calibrated altitude layers is obtained.
Theoretically, in radar coverage, effective sounding website and effective sounding wind field record more, Vbs(alt,
γ ', ω, t) it more can truly reflect the Characteristics of Wind Field of each height layer in high-altitude；, whereas if sounding website is excessively sparse or sounding
Wind field record it is very few, due to biggest impact radius RifLimitation, some regions of radar coverage then may be without step 9
Calculating process, this also allows for discontinuity of the wind field in large scale, is avoided " exaggerated correction ".
The winds retrieval information of multiple calibrated altitude layers in radar coverage can be obtained in calculating through the above steps.
These information can be engaged in Synoptic Diagnostic for professional and provide aid decision after graphical representation；These can also be believed
Breath is input in numerical model as initial fields, to improve synoptic model accuracy of the forecast.
Claims (7)
1. a kind of upper-level winds inversion method based on Doppler radar data, method includes the following steps:
Step 1: assuming that preparing the upper-level winds of a certain moment t of inverting, select a radar first, using portion's radar t moment and
Before closest to the n times radar data of t moment, storm body that radar data is characterized is obtained in each sensing point using extrapolation algorithm
Displacement characteristic vector, the displacement characteristic vector of any sensing point is denoted asWherein, N takes the integer of 2 < N≤10；Indicate a certain elevation angle face that radar body is swept, γ indicates the distance between radar detection point and radar center point library, and ω indicates thunder
The azimuth swept up to body, t indicate that handled radar data corresponds to the time；
Step 2: by displacement characteristic vector in upper stepThe period of time T swept according to radar bodyit, unit distance library pair
The length L answeredu, it is converted into moving velocity vector as the following formula
As vector, the numerical value of speed and the information in direction are characterized, is definedSignified direction with
The angle of direct north is θ, whenWhen being directed toward direct north, θ=0 is rotated in a clockwise direction, and θ is incrementally increased；
Step 3: by moving velocity vector in upper stepIt is decomposed into two components of U, V, wherein U component and current diameter
To vertical, V component radial is overlapped with current, it is assumed thatIt indicatesIn the component of U axis；It indicatesIn the component of V axis；
Step 4: assuming that any sensing point in radar dataThe radar average radial velocity of position isWhenWhen, indicate that the direction of the speed is along being radially directed towards the center of circle；Conversely, working asWhen, the direction for indicating the speed is along radial away from the center of circle；It willAs benchmark speed
Degree calculatesWithMultiplying power, and zoomed in or out with this multiplying powerNumerical value,
It obtainsCalculation formula is as follows:
Step 5: by the component in the direction UWith the component in the direction VIt merges, calculates as the following formula
Revised moving velocity vector is obtained, is denoted as
Step 6: according to step 2~step 5, to the displacement characteristic vector of each sensing point of radarIt carries out one by one
It calculates, obtains t moment in each elevation angle face that radar body is swept, each range bin, each azimuthal revised movement speed
Vector, i.e.,
Step 7: to sensing point any in radar dataLookup is called out the stops with the most similar sounding in its geographical location
Point reads the sounding sounding of detection time closest to t moment in a little of calling out the stops and calls out the stops point data, on the data record earth certain
One longitude and latitude position, the vertically upward wind direction and wind speed of multiple height layers or multiple barosphere standard levels；It obtains and each spy
The corresponding sounding of measuring point is called out the stops point data, and V is recorded ash(lot,lat,alt)；
Step 8: the above step sounding is called out the stops point data VhHeight, each by radar on the basis of each alt in (lot, lat, alt)
Moving velocity vector after the amendment in elevation angle faceUsing spatial interpolation algorithm, the height corresponding to each alt is converted to
The moving velocity vector of degree, is denoted as Vip(alt,γ′,ω,t)；γ ' expression in the horizontal direction sensing point and radar center point it
Between distance, the value of γ ' by γ,And LuIt is acquired by triangle method of geometry；
Step 9: defining an a radius of influence r and biggest impact radius Rif, meet 0 < r < Rif, r expression Vip(alt,γ′,
ω, t) in the space length called out the stops between a little of any point and sounding, longitude, latitude and the V which can be called out the stops a little by soundingip
Warp, the latitude of point in (alt, γ ', ω, t) are calculated by method of geometry；VipThe warp of (alt, γ ', ω, t) each point
Degree and longitude by portion's radar geographical location and γ ' be calculated；
Construct a weighing factor function Wif(r), which is the decreasing function that must meet the following conditions: 0 < Wif(r)<1；According to
Sounding is called out the stops point data Vh(lot, lat, alt), to the V in step 8ip(alt, γ ', ω, t) carries out deviation by following formula and orders
Just:
Vbs(alt, γ ', ω, t)=(1-Wif(r))×Vip(alt,γ′,ω,t)+Wif(r)×Vh(lot,lat,alt)。
2. according to claim 1 based on the upper-level winds inversion method of Doppler radar data, in step 1, the n times
The data mode of radar data is baseis reflectivity data and average radial velocity data.
3. described outer in step 1 according to claim 2 based on the upper-level winds inversion method of Doppler radar data
Predication method uses Maximum correlation method or cross-correlation technique.
4. according to claim 3 based on the upper-level winds inversion method of Doppler radar data, in step 7, if should
The call out the stops record of point data of sounding only has atmospheric pressure information, not apart from the elevation information on ground, passes through pressure-height formula meteorology
Method converts to obtain.
5. -4 any upper-level winds inversion method based on Doppler radar data according to claim 1, in step 7,
Point data is called out the stops using non-standard level for sounding, using interpolation related algorithm, is inside inserted into index bed.
6. according to claim 5 based on the upper-level winds inversion method of Doppler radar data, in step 8, after amendment
Moving velocity vectorIt is interpolated into alt height and passes through upper and lower each 1 elevation angle face adjacent with the locational space
Moving velocity vector is calculated using anti-distance weighting algorithm.
7. according to claim 6 based on the upper-level winds inversion method of Doppler radar data, in step 9, Vip(alt,
γ ', ω, t) parallelogram or proper orthogonal decomposition solution can be used.
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 true CN109100723A (en) | 2018-12-28 |
| CN109100723B 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) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109782019A (en) * | 2019-03-15 | 2019-05-21 | 中国科学技术大学 | Atmosphere pollution two dimensional motion speed measurement method and device |
| CN110288856A (en) * | 2019-06-21 | 2019-09-27 | 中国民用航空总局第二研究所 | The Scheduled Flight monitoring system and method for fine forecast based on wind |
| CN112946657A (en) * | 2021-02-03 | 2021-06-11 | 南京信息工程大学 | Method for identifying ground wind field in strong convection weather |
| CN113009490A (en) * | 2021-02-20 | 2021-06-22 | 江苏省气象台 | Radar three-dimensional wind field inversion method based on high-resolution mode dynamic constraint |
| CN113311436A (en) * | 2021-04-30 | 2021-08-27 | 中国人民解放军国防科技大学 | Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform |
| CN114063465A (en) * | 2021-09-22 | 2022-02-18 | 中国航空工业集团公司西安飞机设计研究所 | Distributed countermeasure simulation system visual jitter elimination method and visual nodes |
| CN115825894A (en) * | 2022-11-17 | 2023-03-21 | 中国能源建设集团广东省电力设计研究院有限公司 | Method and device for determining wind energy capture position, terminal equipment and medium |
| CN116381692A (en) * | 2023-04-18 | 2023-07-04 | 南京市气象台 | Precipitation phase state identification QPE algorithm based on X-band dual-polarization radar |
Citations (6)
| 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 |
| CN104035096A (en) * | 2014-06-06 | 2014-09-10 | 南京大学 | Vertical wind profile nonlinear inversion method based on Doppler weather radar |
| CN104503000A (en) * | 2014-12-15 | 2015-04-08 | 深圳航天东方红海特卫星有限公司 | Sonde wind measuring system and wind measuring method |
| CN106199605A (en) * | 2016-07-06 | 2016-12-07 | 西南技术物理研究所 | Wind field error correcting method |
| CN107741587A (en) * | 2017-09-15 | 2018-02-27 | 北京无线电测量研究所 | A kind of meteorological detection method and system of dimensional wind |
| CN108107434A (en) * | 2017-12-26 | 2018-06-01 | 厦门市气象灾害防御技术中心(海峡气象开放实验室厦门市避雷监测技术中心) | Area three-dimensional wind field pattern splicing method based on Dual-Doppler weather radar inverting |
-
2018
- 2018-07-25 CN CN201810826385.5A patent/CN109100723B/en active Active
Patent Citations (6)
| 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 |
| CN104035096A (en) * | 2014-06-06 | 2014-09-10 | 南京大学 | Vertical wind profile nonlinear inversion method based on Doppler weather radar |
| CN104503000A (en) * | 2014-12-15 | 2015-04-08 | 深圳航天东方红海特卫星有限公司 | Sonde wind measuring system and wind measuring method |
| CN106199605A (en) * | 2016-07-06 | 2016-12-07 | 西南技术物理研究所 | Wind field error correcting method |
| CN107741587A (en) * | 2017-09-15 | 2018-02-27 | 北京无线电测量研究所 | A kind of meteorological detection method and system of dimensional wind |
| CN108107434A (en) * | 2017-12-26 | 2018-06-01 | 厦门市气象灾害防御技术中心(海峡气象开放实验室厦门市避雷监测技术中心) | Area three-dimensional wind field pattern splicing method based on Dual-Doppler weather radar inverting |
Non-Patent Citations (2)
| Title |
|---|
| 周海光: "多普勒天气雷达风场反演技术研究进展", 《地球科学进展》 * |
| 周鑫等: "多普勒雷达高空风场反演方法研究", 《现代雷达》 * |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109782019A (en) * | 2019-03-15 | 2019-05-21 | 中国科学技术大学 | Atmosphere pollution two dimensional motion speed measurement method and device |
| CN110288856A (en) * | 2019-06-21 | 2019-09-27 | 中国民用航空总局第二研究所 | The Scheduled Flight monitoring system and method for fine forecast based on wind |
| CN112946657A (en) * | 2021-02-03 | 2021-06-11 | 南京信息工程大学 | Method for identifying ground wind field in strong convection weather |
| CN112946657B (en) * | 2021-02-03 | 2023-10-27 | 南京信息工程大学 | Method for identifying ground wind field in strong convection weather |
| CN113009490A (en) * | 2021-02-20 | 2021-06-22 | 江苏省气象台 | 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 |
| CN113311436A (en) * | 2021-04-30 | 2021-08-27 | 中国人民解放军国防科技大学 | Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform |
| CN114063465A (en) * | 2021-09-22 | 2022-02-18 | 中国航空工业集团公司西安飞机设计研究所 | Distributed countermeasure simulation system visual jitter elimination method and visual nodes |
| CN114063465B (en) * | 2021-09-22 | 2024-05-03 | 中国航空工业集团公司西安飞机设计研究所 | View jitter elimination method and view node of distributed countermeasure simulation system |
| CN115825894A (en) * | 2022-11-17 | 2023-03-21 | 中国能源建设集团广东省电力设计研究院有限公司 | Method and device for determining wind energy capture position, terminal equipment and medium |
| CN115825894B (en) * | 2022-11-17 | 2023-08-18 | 中国能源建设集团广东省电力设计研究院有限公司 | Method, device, terminal equipment and medium for determining wind energy capturing position |
| CN116381692A (en) * | 2023-04-18 | 2023-07-04 | 南京市气象台 | Precipitation phase state identification QPE algorithm based on X-band dual-polarization radar |
| CN116381692B (en) * | 2023-04-18 | 2024-02-23 | 南京市气象台 | Precipitation phase state identification QPE algorithm based on X-band dual-polarization radar |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109100723B (en) | 2022-05-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109100723A (en) | Upper-level winds inversion method based on Doppler radar data | |
| CN111583214B (en) | Sea surface wind speed inversion method based on RBF neural network and based on marine radar image | |
| Kurtz et al. | Sea ice thickness, freeboard, and snow depth products from Operation IceBridge airborne data | |
| CN110609287B (en) | Double-frequency radar scatterometer and method for simultaneously measuring sea surface wind field and flow field | |
| Liang et al. | A composite approach of radar echo extrapolation based on TREC vectors in combination with model-predicted winds | |
| CN107843895B (en) | A kind of Dual-Doppler weather radar dimensional wind inversion method | |
| CN110879428B (en) | Method for measuring rain through networking in overlapped area between radars in same wave band under complex terrain condition | |
| JP6689396B2 (en) | Weather forecasting device, weather forecasting method, and program | |
| Crétaux et al. | Absolute calibration of Jason radar altimeters from GPS kinematic campaigns over Lake Issykkul | |
| Zhao et al. | Theoretical analysis and numerical experiments of variational adjoint approach for refractivity estimation | |
| CN114186484B (en) | Surface soil humidity remote sensing inversion method suitable for high-altitude area | |
| Mohammadloo et al. | Correcting multibeam echosounder bathymetric measurements for errors induced by inaccurate water column sound speeds | |
| CN113009531A (en) | Small-scale high-precision low-altitude troposphere atmosphere refractive index model | |
| CN113960625B (en) | Water depth inversion method based on satellite-borne single-photon laser active and passive remote sensing fusion | |
| CN106872979B (en) | Sea wave parameter acquisition method based on sea surface fluctuation moving target SAR image refocusing | |
| CN114252875B (en) | High-precision meshing method for imaging altitude data | |
| CN116609859A (en) | Weather disaster high-resolution regional mode forecasting system and method | |
| CN109444892B (en) | Calibration method for radar quantitative measurement of precipitation | |
| CN116500648A (en) | Wind profile inversion method for foundation laser radar target area | |
| CN104360335A (en) | Ionized layer parametric inversion method based on AIS target indication | |
| US20240027575A1 (en) | Vertical situation display using cell-based architecture | |
| Yang et al. | VDatum for eastern Louisiana and Mississippi coastal waters: tidal datums, marine grids, and sea surface topography | |
| CN113670951B (en) | Microwave radiometer self-adaptive inversion algorithm based on shipborne mobile platform | |
| CN116203591B (en) | Method for jointly estimating high-precision ionosphere VTEC of China area based on multiple stations | |
| CN114020725B (en) | Window sliding GPM data correction method considering spatial distribution |
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 |