CN114675055B - Air speed measurement error compensation method of sonde based on inertial system - Google Patents

Air speed measurement error compensation method of sonde based on inertial system Download PDF

Info

Publication number
CN114675055B
CN114675055B CN202210317866.XA CN202210317866A CN114675055B CN 114675055 B CN114675055 B CN 114675055B CN 202210317866 A CN202210317866 A CN 202210317866A CN 114675055 B CN114675055 B CN 114675055B
Authority
CN
China
Prior art keywords
sonde
wind speed
axis
umbrella
angle
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
CN202210317866.XA
Other languages
Chinese (zh)
Other versions
CN114675055A (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202210317866.XA priority Critical patent/CN114675055B/en
Publication of CN114675055A publication Critical patent/CN114675055A/en
Application granted granted Critical
Publication of CN114675055B publication Critical patent/CN114675055B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P21/00—Testing or calibrating of apparatus or devices covered by the preceding groups
    • G01P21/02—Testing or calibrating of apparatus or devices covered by the preceding groups of speedometers
    • G01P21/025—Testing or calibrating of apparatus or devices covered by the preceding groups of speedometers for measuring speed of fluids; for measuring speed of bodies relative to fluids
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/18—Testing or calibrating meteorological apparatus
    • 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

  • Environmental & Geological Engineering (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Atmospheric Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Navigation (AREA)

Abstract

The invention discloses a sonde wind speed measurement error compensation method based on an inertial system. The error compensation is carried out by utilizing the measured data of the inertial system, so that the point-by-point high-precision compensation can be carried out according to the actual condition of each test of the sonde, and the wind speed detection precision is improved.

Description

Air speed measurement error compensation method of sonde based on inertial system
Technical Field
The invention belongs to the technical field of meteorological detection, and particularly relates to a sonde wind speed measurement error compensation method based on an inertial system.
Background
Typhoons are one of the most serious natural disasters in the world. In order to realize the more essential scientific cognition of the natural phenomenon of typhoons, the key technical bottleneck is to directly detect the typhoons in a multi-element, long-process and fine manner to obtain key meteorological data such as temperature, humidity, pressure, wind speed, wind direction and the like of typhoons in an inner core area, wherein the measurement of the wind speed of the typhoons in the inner core area is particularly important for predicting typhoons paths and strength. At present, the measurement of typhoon wind speed is mainly performed by a wind profile radar, a laser wind-finding radar, a sonde and the like, and in order to meet the in-situ detection of typhoon nuclear area wind speed, the sonde can only be used for measurement. However, the traditional sonde equipment can only calculate the wind speed by means of position change data detected by a sonde GPS or a Beidou system, and the error is large. When the sonde is in extreme environments such as typhoons, the sonde can roll around a parachute rapidly in the falling process, the severe change of the attitude of the sonde can cause great error influence on wind speed measurement, at the moment, the speed of the sonde cannot directly represent the flow speed of an actual wind field, and the wind speed measurement error caused by the change of the motion attitude of the sonde must be eliminated.
At present, a great deal of study is carried out on a navigation type sonde wind speed measurement algorithm by a plurality of students researching high altitude detection at home and abroad. Jaakko study VIASALA RS, 92 sonde indicates that wind speed measurements are to be data correlation analyzed and smoothed. Jauhiaain et al in summary of the international sonde comparison test indicate that the wind speed measurement process should take into account sonde oscillation factors, but do not mention the method of processing. Zhang Engong according to the calculation method of the compound pendulum period, the pendulum period of the sonde is researched, and the influence of the sonde pendulum on wind speed measurement is removed by using a filtering method. However, the above methods idealize the motion of the sonde into simple periodic cone-swing motion, so as to filter all detection data, and the accuracy is low. When the sonde is in the typhoon-like environment with rapid change and extremely strong wind speed, the motion is complex and unknown, and the problem of under-compensation or over-compensation and the like can be caused by simply compensating the detection data by using the periodic cone pendulum motion, so that the wind speed and direction information of the sonde in the current space time can not be accurately obtained.
Patent application CN114019583A discloses a high-precision safe wind measuring system and a high-precision safe wind measuring method based on inertial compensation, wherein the speed is obtained by integrating the acceleration of a sonde measured by an inertial module, and the wind speed measured by a positioning module is compensated in a weighted average mode; according to the attitude information of the sonde measured by the inertia module, the swing period of the sonde is calculated, and the influence of the swing of the sonde is filtered by multi-point filtering aiming at the data of different periods. First, it is incorrect and meaningless to integrate the speed of the sonde acceleration measured by the inertial module for the first point. In the sonde system, two components of a parachute and a sonde mainly exist, the meaning of wind speed compensation is mainly to compensate complex motion of the sonde around the parachute, if the wind speed is calculated by simply utilizing acceleration information obtained by an inertia module, the measurement is equivalent to the instantaneous wind speed of the whole sonde system, the effect of the measurement is the same as that of GPS, the measurement error caused by the sonde moving around the parachute is not solved, and the wind speed compensation effect cannot be achieved. For the second point, in order to solve the error, the measured attitude information in the inertial system must be utilized, but the method only solves the swing period of the sonde, which has the same problem as the research article mentioned above, namely, the problem that the wind measurement data at all moments are overcompensated or undercompensated by utilizing the solved period, and not each moment needs to be subjected to the uniform swing period influence, and each moment has unique characteristics, so that the wind speed error measured by the sonde is compensated one by one according to the actual situation.
In summary, it is desirable to design a method for compensating for the wind speed measurement error of a sonde, which compensates for the wind speed measurement error generated by the attitude change of the sonde.
Disclosure of Invention
Aiming at the wind speed measurement error of the current sonde caused by the self-posture change under the severe working condition, the invention provides a sonde wind speed measurement error compensation method based on an inertial system. And solving the true motion gesture track of the sonde by combining the inertial navigation information on the sonde, and compensating the wind speed measurement error caused by the motion of the sonde. The specific technical scheme of the invention is as follows:
the utility model provides a sonde wind speed measurement error compensation method based on inertial system, the structure of throwing down in the wind speed measurement process includes umbrella, rope and sonde, and the sonde passes through the rope and links to each other with the umbrella bottom, and the umbrella provides the time of floating for the sonde, and the rope keeps taut state in throwing down the in-process, the compensation method includes following steps:
s1: integrating an inertial system on the sonde, performing high-altitude downward casting detection to obtain three-dimensional space meteorological data, and recording and storing in-situ detection data of the sonde and triaxial inertial navigation information by a ground receiver in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received sonde GPS/Beidou information;
s3: calculating attitude information in the falling process of the sonde by utilizing the received sonde triaxial inertial navigation information;
s4: the result of the step S3 is transformed and resolved through a space coordinate system to obtain real-time relative space motion attitude information of the sonde;
s5: calculating a wind speed measurement error of the sonde due to self posture change according to the relative space movement posture information of the sonde obtained in the step S4, and substituting the wind speed measurement error into the three-dimensional space wind speed information calculated in the step S2 to finish wind speed compensation;
s6: and (3) performing rolling time window smoothing on the wind speed data compensated in the step (S5), and further eliminating the influence of environmental noise.
Further, the triaxial inertial navigation information in the step S3 comprises triaxial acceleration, triaxial angular velocity and triaxial angle information, the triaxial angle information with the angle range of (-180 degrees, 180 degrees) is reserved, the triaxial angle original data is interpolated through cubic spline interpolation, and differential solution is carried out on displacement to obtain the movement velocity;
the inertial system uses Euler angles to represent the coordinate system rotation sequence when the attitude is represented by the Euler angles, namely, the inertial system rotates around the Z axis, then rotates around the Y axis and finally rotates around the X axis, so that the angle range of rotation around the Y axis is only +/-90 degrees;
the barycenter of the sonde is taken as an origin O, and the direction of the connecting point of the barycenter pointing rope of the sonde and the sonde is taken as OX b Axis with the direction from the mass center of the sonde to the transverse center line of the sonde as OZ b Shaft to and OX b Shaft and OZ b The vertical axis direction is OY b Shaft, establishing coordinate system OX fixedly connected with sonde b Y b Z b Namely b;
with the umbrella bottom as the origin O and the east direction of the geographic coordinate system pointing from the umbrella bottom as OY n Axis with north direction from umbrella bottom to geographic coordinate system as OZ n Axis with the direction from the umbrella bottom to the geographic coordinate system as the zenith direction OX n Shaft, establishing coordinate system OX fixedly connected with umbrella bottom n Y n Z n Namely n is; when the sonde is static, i.e. the sonde and the umbrella do not swing relatively, the b tie can be obtained by translating the n tie along the rope;
the bottom of the umbrella at the following throwing time is taken as an origin O, OX i Axis-directed upward, OY i Axis-pointing east, OZ i The axis points to north direction to establish a space rectangular coordinate system OX i Y i Z i I is;
definition of rope and OX i The included angle of the axes is theta x Rope and OY i The included angle of the axes is theta y ；
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for wrap OX b The axis rotates by an angle theta, and the original data has n+1 data nodes: (θ) 0 ,t 0 )， (θ 1 ,t 1 )，…，(θ i ,t i )，…，(θ n ,t n )，i＝0, 1..n, wherein t i For the current data recording time, θ i At t i Time-of-day corresponding OX b The shaft rotation angle value, in the cubic spline interpolation process, satisfies: on each segment interval [ t ] i ,t i+1 ]On f (t) =f i (t) is a cubic equation; interpolation condition f (t i )＝θ i The method comprises the steps of carrying out a first treatment on the surface of the The curve is smooth, i.e., f (t), f' (t) are continuous; let each coefficient of the third equation be a i 、b i 、c i 、d i The following steps are:
Figure BDA0003569427230000041
wherein f (t) is the definition domain after cubic spline interpolation in (t) 0 ,t n ) The upper expression is a function of the angle θ=f (t) at time t, f i (t) is a segment interval [ t ] i ,t i+1 ]The upper represents a function corresponding to the angle of the moment, f i ′(t)、f i "(t) is f respectively i A first and second derivative function of (t);
s3-2: obtaining a system of linear equations for the unknowns:
s3-2-1: from f i (t i )＝a i +b i (t i -t i )+c i (t i -t i ) 2 +d i (t i -t i ) 3 ＝θ i Obtaining a i ＝θ i, wherein , ai ,b i ,c i ,d i Are all coefficients;
s3-2-2: let h i ＝t i+1 -t i By boundary continuous condition f i (t i+1 )＝θ i+1 Pushing out:
Figure BDA0003569427230000042
from the boundary derivative continuous condition f i ′(t i+1 )＝f′ i+1 (t i+1 ) The method comprises the following steps of:
Figure BDA0003569427230000043
from f i ″(t i+1 )＝f″ i+1 (t i+1 ) The method comprises the following steps of:
Figure BDA0003569427230000044
let m i ＝f i ″(t i )＝2c i Then formula (2) is written as m i +6h i d i ＝m i+1 Then:
Figure BDA0003569427230000045
h i ,m i respectively intermediate variables;
s3-2-3: will a i ＝t i ，
Figure BDA0003569427230000046
Substituting into formula (1) to obtain:
Figure BDA0003569427230000047
will a i ，b i ，c i ，d i Substituting into formula (2), obtain:
Figure BDA0003569427230000051
s3-2-4: expanding equation (5) to all time points gives m 0 ,m 1 ,m 2 ,...,m n For a linear system of equations of unknown variables, a natural boundary condition, i.e. m, is selected during interpolation 0 ＝0，m n =0, then the system of equations is:
Figure BDA0003569427230000052
to this end, in each segment interval [ t ] i ,t i+1 ]Above, a can be determined by solving the system of equations i ，b i ，c i ，d i ；
S3-3: setting a new sampling time interval
Figure BDA0003569427230000053
The interpolated angle θ is:
θ i+j ＝f i (t i +jΔt i )，j＝1,2,...,N-1
the corresponding time precision information is obtained by controlling the N value;
s3-4: for the winding OY b Rotation angle of axis gamma, around OZ b The rotation angle psi data of the shaft is subjected to data interpolation by the same cubic spline interpolation method, the total number of data before interpolation is n+1, and the number of data after interpolation is n+1.
Further, in the step S4, at a time t i Coordinate system OX of sonde b Y b Z b Relative to the coordinate system OX in which the umbrella is located n Y n Z n The displacement and angle transformation relationship occurs under the limitation of the rope only under tension, and the motion speed of the sonde relative to the umbrella is calculated through the attitude angle
Figure BDA0003569427230000061
S4-1: according to the rotation sequence of the z-y-x coordinate axes, the coordinate transformation matrix from the n system to the b system is as follows:
Figure BDA0003569427230000062
wherein ,Cθ Is around OX n Rotation matrix of rotation θ, C γ Is wound around OY n Rotation matrix of rotation gamma, C ψ Is wound around OZ n The rotation matrix of rotation ψ is expressed as:
Figure BDA0003569427230000063
/>
s4-2: when only the rotation relationship of the b-series to the n-series is considered and the displacement relationship is not considered, OX n ＝[1 0 0] T The representation under b
Figure BDA0003569427230000064
Expressed as:
Figure BDA0003569427230000065
θ x represented as
Figure BDA0003569427230000066
With OX b Included angle between:
Figure BDA0003569427230000067
s4-3: when only the rotation relation of the b-series to the n-series is considered and the displacement relation is not considered, OX b ＝[1 0 0] T The representation under b
Figure BDA0003569427230000068
Expressed as:
Figure BDA0003569427230000069
vector OH is
Figure BDA00035694272300000610
In OY n Z n Projection onto a plane, n is expressed below:
OH n ＝[0 -cosγcosψ sinγ] T
according to the geometric relationship, θ y Represented as OY n Included angle with OH:
Figure BDA0003569427230000071
s4-4: let the rope length be l, then t k At time, k=0, 1, nN, and θ is obtained through angle calculation x (k) And theta y (k) The position of the sonde relative to the umbrella bottom is:
Figure BDA0003569427230000072
by forward difference, t is obtained k East speed of time sonde relative to umbrella
Figure BDA0003569427230000073
North speed->
Figure BDA0003569427230000074
Figure BDA0003569427230000075
wherein ,
Figure BDA0003569427230000076
at t k East speed calculated from the processed triaxial angular speed at time +.>
Figure BDA0003569427230000077
At t k The north velocity is calculated from the processed triaxial angular velocity.
Further, in the step S5, the GPS/beidou information of the sonde includes the north speed of the sonde relative to the geographic coordinate system
Figure BDA0003569427230000078
East speed of sonde relative to geographic coordinate system>
Figure BDA0003569427230000079
The wind speed obtained by the GPS/Beidou information of the sonde is as follows:
Figure BDA00035694272300000710
because the sonde is connected with the umbrella by virtue of the rope, the wind speed obtained by GPS measurement is actually the superposition of the motion of the sonde relative to the umbrella and the wind speed, so the wind speed obtained by GPS measurement is decomposed into:
Figure BDA00035694272300000711
wherein ,
Figure BDA00035694272300000712
measuring the speed of the sonde relative to the umbrella, n, for the inertial system,/for the inertial system>
Figure BDA00035694272300000713
The wind speed is the wind speed under the coordinate system of the umbrella lower point, namely the i system;
t i speed corresponding to time
Figure BDA00035694272300000714
And->
Figure BDA00035694272300000715
A speed at which the corresponding moment can be found +.>
Figure BDA00035694272300000716
And->
Figure BDA00035694272300000717
The wind speed at this time is:
Figure BDA00035694272300000718
for V w Processing to obtain the final compensated wind speed V wind ：
Figure BDA0003569427230000081
Wherein floor (·) is rounded down.
The invention has the beneficial effects that:
1. according to the air speed measurement error compensation method of the sonde based on the inertial system, which is provided by the invention, the real attitude information of the sonde motion is obtained by integrating the inertial component on the traditional sonde system, so that the air speed measurement error caused by the sonde motion is compensated.
2. Compared with the attitude errors estimated by other methods, the method has the advantages of authenticity and reliability, and can compensate the real-time wind speed information measured by the sonde point by point with high precision according to the actual condition of each test of the sonde, rather than compensating the wind speed by estimating the approximate swing period of the sonde, thereby avoiding the problems of overcompensation and undercompensation of the method.
Drawings
For a clearer description of an embodiment of the invention or of the solutions of the prior art, reference will be made to the accompanying drawings, which are used in the embodiments and which are intended to illustrate, but not to limit the invention in any way, the features and advantages of which can be obtained according to these drawings without inventive labour for a person skilled in the art. Wherein:
FIG. 1 is a flow chart of the inertial system compensated sonde wind speed measurement error of the present invention;
FIG. 2 is a coordinate system definition of the sonde and umbrella;
FIG. 3 is an attitude definition;
FIG. 4 is a spatial coordinate system positional relationship;
FIG. 5 is a measured motion trajectory of the sonde relative to the parachute;
FIG. 6 is a physical view of a sonde with an inertial navigation system;
FIG. 7 is an inertial navigation module mounting location and triaxial definition on a sonde;
FIG. 8 shows the result of the sonde east and north wind speed compensation, wherein (a) is the sonde east wind speed compensation result, the dotted line is the east wind speed obtained by GPS calculation, the solid line is the east wind speed compensated based on the inertial system, (b) is the sonde north wind speed compensation result, the dotted line is the north wind speed obtained by GPS calculation, and the solid line is the north wind speed compensated based on the inertial system;
fig. 9 shows the compensation result of the sonde after east and north wind speed vector synthesis, the dotted line is the wind speed obtained by the GPS calculation, and the solid line is the wind speed after inertial system compensation.
Detailed Description
In order that the above-recited objects, features and advantages of the present invention will be more clearly understood, a more particular description of the invention will be rendered by reference to the appended drawings and appended detailed description. It should be noted that, without conflict, the embodiments of the present invention and features in the embodiments may be combined with each other.
In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention, however, the present invention may be practiced in other ways than those described herein, and therefore the scope of the present invention is not limited to the specific embodiments disclosed below.
The invention integrates an inertial system on the sonde and analyzes the captured true motion attitude information of the sonde so as to compensate the wind speed measurement error caused by the motion of the sonde under severe working conditions.
As shown in FIG. 1, a method for compensating wind speed measurement error of a sonde based on an inertial system, wherein a down-casting structure in the wind speed measurement process comprises an umbrella, a rope and a sonde, the sonde is connected with an umbrella bottom through the rope, the umbrella provides a floating time for the sonde, and the rope is kept in a tensioned state in the down-casting process, and the method comprises the following steps:
s1: integrating an inertial system on the sonde, performing high-altitude downward casting detection to obtain three-dimensional space meteorological data, and recording and storing in-situ detection data of the sonde and triaxial inertial navigation information by a ground receiver in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received sonde GPS/Beidou information;
s3: calculating attitude information in the falling process of the sonde by utilizing the received sonde triaxial inertial navigation information;
the three-axis inertial navigation information in the step S3 comprises three-axis acceleration, three-axis angular velocity and three-axis angle information, the three-axis angle information with the angle range of (-180 degrees, 180 degrees) is reserved, three-time spline interpolation is used for interpolating original data of the three-axis angle, and differential solution is carried out on displacement to obtain the motion speed;
as shown in fig. 2, the inertial system uses the euler angle to represent the coordinate system rotation sequence when the attitude is represented by the rotation sequence of the inertial system around the Z axis, then around the Y axis and finally around the X axis, so that the rotation angle range around the Y axis is only +/-90 degrees;
the barycenter of the sonde is taken as an origin O, and the direction of the connecting point of the barycenter pointing rope of the sonde and the sonde is taken as OX b Axis with the direction from the mass center of the sonde to the transverse center line of the sonde as OZ b Shaft to and OX b Shaft and OZ b The vertical axis direction is OY b Shaft, establishing coordinate system OX fixedly connected with sonde b Y b Z b Namely b;
with the umbrella bottom as the origin O and the east direction of the geographic coordinate system pointing from the umbrella bottom as OY n Axis with north direction from umbrella bottom to geographic coordinate system as OZ n Axis with the direction from the umbrella bottom to the geographic coordinate system as the zenith direction OX n Shaft, establishing coordinate system OX fixedly connected with umbrella bottom n Y n Z n Namely n is; when the sonde is static, i.e. the sonde and the umbrella do not swing relatively, the b tie can be obtained by translating the n tie along the rope;
the bottom of the umbrella at the following throwing time is taken as an origin O, OX i Axis-directed upward, OY i Axis-pointing east, OZ i The axis points to north direction to establish a space rectangular coordinate system OX i Y i Z i I is;
definition of rope and OX i The included angle of the axes is theta x Rope and OY i The included angle of the axes is theta y ；
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for wrap OX b The axis rotates by an angle theta, and the original data has n+1 data nodes: (θ) 0 ,t 0 )， (θ 1 ,t 1 )，…，(θ i ,t i )，…，(θ n ,t n ) I=0, 1,..n, where t i For the current data recording time, θ i At t i Time-of-day corresponding OX b The shaft rotation angle value, in the cubic spline interpolation process, satisfies: on each segment interval [ t ] i ,t i+1 ]On f (t) =f i (t) is a cubic equation; interpolation condition f (t i )＝θ i The method comprises the steps of carrying out a first treatment on the surface of the The curve is smooth, i.e., f (t), f' (t) are continuous; let each coefficient of the third equation be a i 、b i 、c i 、d i The following steps are:
Figure BDA0003569427230000101
wherein f (t) is the definition domain after cubic spline interpolation in (t) 0 ,t n ) The upper expression is a function of the angle θ=f (t) at time t, f i (t) is a segment interval [ t ] i ,t i+1 ]The upper represents a function corresponding to the angle of the moment, f i ′(t)、f i "(t) is f respectively i A first and second derivative function of (t);
s3-2: obtaining a system of linear equations for the unknowns:
s3-2-1: from f i (t i )＝a i +b i (t i -t i )+c i (t i -t i ) 2 +d i (t i -t i ) 3 ＝θ i Obtaining a i ＝θ i, wherein , ai ,b i ,c i ,d i Are all coefficients;
s3-2-2: let h i ＝t i+1 -t i By boundary continuous condition f i (t i+1 )＝θ i+1 Pushing out:
Figure BDA0003569427230000111
from the boundary derivative continuous condition f i ′(t i+1 )＝f′ i+1 (t i+1 ) The method comprises the following steps of:
Figure BDA0003569427230000112
from f i ″(t i+1 )＝f″ i+1 (t i+1 ) The method comprises the following steps of:
Figure BDA0003569427230000113
let m i ＝f i ″(t i )＝2c i Then formula (2) is written as m i +6h i d i ＝m i+1 Then:
Figure BDA0003569427230000114
h i ,m i respectively intermediate variables;
s3-2-3: will a i ＝t i ，
Figure BDA0003569427230000115
Substituting into formula (1) to obtain:
Figure BDA0003569427230000116
will a i ，b i ，c i ，d i Substituting into formula (2), obtain:
Figure BDA0003569427230000117
s3-2-4: expanding equation (5) to all time points gives m 0 ,m 1 ,m 2 ,...,m n For a linear system of equations of unknown variables, a natural boundary condition, i.e. m, is selected during interpolation 0 ＝0，m n =0, then the system of equations is:
Figure BDA0003569427230000118
Figure BDA0003569427230000121
to this end, in each segment interval [ t ] i ,t i+1 ]Above, a can be determined by solving the system of equations i ，b i ，c i ，d i ；
S3-3: setting a new sampling time interval
Figure BDA0003569427230000122
The interpolated angle θ is:
θ i+j ＝f i (t i +jΔt i )，j＝1,2,...,N-1
the corresponding time precision information is obtained by controlling the N value;
s3-4: for the winding OY b Rotation angle of axis gamma, around OZ b The rotation angle psi data of the shaft is subjected to data interpolation by the same cubic spline interpolation method, the total number of data before interpolation is n+1, and the number of data after interpolation is n+1.
S4: the result of the step S3 is transformed and resolved through a space coordinate system to obtain real-time relative space motion attitude information of the sonde; in step S4, at time t i Coordinate system OX of sonde b Y b Z b Relative to the coordinate system OX in which the umbrella is located n Y n Z n The displacement being related to angular transformation under the limitation of a tensioned rope only, see figure 4As shown, the motion speed of the sonde relative to the umbrella is calculated through the attitude angle
Figure BDA0003569427230000123
S4-1: according to the rotation sequence of the z-y-x coordinate axes, the coordinate transformation matrix from the n system to the b system is as follows:
Figure BDA0003569427230000124
wherein ,Cθ Is around OX n Rotation matrix of rotation θ, C γ Is wound around OY n Rotation matrix of rotation gamma, C ψ Is wound around OZ n The rotation matrix of rotation ψ is expressed as:
Figure BDA0003569427230000131
s4-2: when only the rotation relationship of the b-series to the n-series is considered and the displacement relationship is not considered, OX n ＝[1 0 0] T The representation under b
Figure BDA0003569427230000132
Expressed as:
Figure BDA0003569427230000133
as shown in fig. 3, θ x Represented as
Figure BDA00035694272300001312
With OX b Included angle between:
Figure BDA0003569427230000134
s4-3: when only the rotation relation of the b-series to the n-series is considered and the displacement relation is not considered, OX b ＝[1 0 0] T In b seriesThe following expression
Figure BDA0003569427230000135
Expressed as:
Figure BDA0003569427230000136
vector OH is
Figure BDA0003569427230000137
In OY n Z n Projection onto a plane, n is expressed below:
OH n ＝[0 -cosγcosψ sinγ] T
according to the geometric relationship, θ y Represented as OY n Included angle with OH:
Figure BDA0003569427230000138
s4-4: let the rope length be l, then t k At time, k=0, 1, nN, and θ is obtained through angle calculation x (k) And theta y (k) The position of the sonde relative to the umbrella bottom is:
Figure BDA0003569427230000139
by forward difference, t is obtained k East speed of time sonde relative to umbrella
Figure BDA00035694272300001310
North speed->
Figure BDA00035694272300001311
Figure BDA0003569427230000141
wherein ,
Figure BDA0003569427230000142
at t k East speed calculated from the processed triaxial angular speed at time +.>
Figure BDA0003569427230000143
At t k The north velocity is calculated from the processed triaxial angular velocity.
S5: calculating a wind speed measurement error of the sonde due to self posture change according to the relative space movement posture information of the sonde obtained in the step S4, and substituting the wind speed measurement error into the three-dimensional space wind speed information calculated in the step S2 to finish wind speed compensation; in step S5, GPS/Beidou information of the sonde comprises the north speed of the sonde relative to a geographic coordinate system
Figure BDA0003569427230000144
East speed of sonde relative to geographic coordinate system>
Figure BDA0003569427230000145
The wind speed obtained by the GPS/Beidou information of the sonde is as follows:
Figure BDA0003569427230000146
because the sonde is connected with the umbrella by virtue of the rope, the wind speed obtained by GPS measurement is actually the superposition of the motion of the sonde relative to the umbrella and the wind speed, so the wind speed obtained by GPS measurement is decomposed into:
Figure BDA0003569427230000147
wherein ,
Figure BDA0003569427230000148
measuring the speed of the sonde relative to the umbrella, n, for the inertial system,/for the inertial system>
Figure BDA0003569427230000149
The wind speed is the wind speed under the coordinate system of the umbrella lower point, namely the i system;
t i speed corresponding to time
Figure BDA00035694272300001410
And->
Figure BDA00035694272300001411
A speed at which the corresponding moment can be found +.>
Figure BDA00035694272300001412
And->
Figure BDA00035694272300001413
The wind speed at this time is:
Figure BDA00035694272300001414
for V w Processing to obtain the final compensated wind speed V wind ：
Figure BDA00035694272300001415
Wherein floor (·) is rounded down.
S6: and (3) performing rolling time window smoothing on the wind speed data compensated in the step (S5), and further eliminating the influence of environmental noise.
In order to facilitate understanding of the above technical solutions of the present invention, the following detailed description of the above technical solutions of the present invention is provided by specific embodiments.
Example 1
As shown in fig. 6-7, a sounding device with a self-grinding inertial navigation system is used for performing a 20km drop-in experiment in Xinjiang to detect three-dimensional space meteorological data. The actual measurement motion track of the sonde relative to the parachute is shown in fig. 5, the wind speed measurement drift caused by random swing of the sonde is obtained through calculation, and the east speed and the north speed obtained through GPS measurement are respectively compensated. The compensation effect is shown in fig. 8, wherein (a) is the result of compensating the eastern wind speed of the sonde, the dotted line is the eastern wind speed obtained by the calculation of the GPS, the solid line is the eastern wind speed compensated based on the inertial system, (b) is the result of compensating the northward wind speed of the sonde, the dotted line is the northward wind speed obtained by the calculation of the GPS, and the solid line is the northward wind speed compensated based on the inertial system. It can be seen that the signal-to-noise ratio of the east wind speed and the north wind speed after inertial navigation data compensation is obviously improved. Vector synthesis is carried out on the wind speed obtained through GPS measurement and the wind speed subjected to inertial navigation compensation, a final wind speed measurement result is shown in fig. 9, a dotted line is the wind speed obtained through GPS calculation, a solid line is the wind speed subjected to inertial system compensation, and it can be seen that the signal to noise ratio of the wind speed subjected to inertial navigation data compensation is obviously improved.
The above description is only of the preferred embodiments of the present invention and is not intended to limit the present invention, but various modifications and variations can be made to the present invention by those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (3)

1. The utility model provides a sonde wind speed measurement error compensation method based on inertial system, the structure of throwing down in the wind speed measurement process includes umbrella, rope and sonde, and the sonde passes through the rope and links to each other with the umbrella bottom, and the umbrella provides the time of floating for the sonde, and the rope keeps taut state in throwing down the in-process, its characterized in that, the compensation method includes following steps:
s1: integrating an inertial system on the sonde, performing high-altitude downward casting detection to obtain three-dimensional space meteorological data, and recording and storing in-situ detection data of the sonde and triaxial inertial navigation information by a ground receiver in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received sonde GPS/Beidou information;
s3: calculating attitude information in the falling process of the sonde by utilizing the received sonde triaxial inertial navigation information;
s4: the result of the step S3 is transformed and resolved through a space coordinate system to obtain real-time relative space motion attitude information of the sonde;
s5: calculating a wind speed measurement error of the sonde due to self posture change according to the relative space movement posture information of the sonde obtained in the step S4, and substituting the wind speed measurement error into the three-dimensional space wind speed information calculated in the step S2 to finish wind speed compensation;
s6: performing rolling time window smoothing on the wind speed data compensated in the step S5, and further eliminating the influence of environmental noise;
the three-axis inertial navigation information in the step S3 comprises three-axis acceleration, three-axis angular velocity and three-axis angle information, the three-axis angle information with the angle range of (-180 degrees, 180 degrees) is reserved, three-axis angle original data are interpolated through cubic spline interpolation, and differential solution is carried out on displacement to obtain the motion speed;
the inertial system uses Euler angles to represent the coordinate system rotation sequence when the attitude is represented by the Euler angles, namely, the inertial system rotates around the Z axis, then rotates around the Y axis and finally rotates around the X axis, so that the angle range of rotation around the Y axis is only +/-90 degrees;
the barycenter of the sonde is taken as an origin O, and the direction of the connecting point of the barycenter pointing rope of the sonde and the sonde is taken as OX b Axis with the direction from the mass center of the sonde to the transverse center line of the sonde as OZ b Shaft to and OX b Shaft and OZ b The vertical axis direction is OY b Shaft, establishing coordinate system OX fixedly connected with sonde b Y b Z b Namely b;
with the umbrella bottom as the origin O and the east direction of the geographic coordinate system pointing from the umbrella bottom as OY n Axis with north direction from umbrella bottom to geographic coordinate system as OZ n Axis with the direction from the umbrella bottom to the geographic coordinate system as the zenith direction OX n Shaft, establishing coordinate system OX fixedly connected with umbrella bottom n Y n Z n Namely n is; when the sonde is static, i.e. the sonde and the umbrella do not swing relatively, the b tie can be obtained by translating the n tie along the rope;
the bottom of the umbrella at the following throwing time is taken as an origin O, OX i Axis-directed upward, OY i Axis-pointing east, OZ i The axis points to north direction to establish a space rectangular coordinate systemOX i Y i Z i I is;
definition of rope and OX i The included angle of the axes is theta x Rope and OY i The included angle of the axes is theta y ；
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for wrap OX b The axis rotates by an angle theta, and the original data has n+1 data nodes: (θ) 0 ,t 0 )，(θ 1 ,t 1 )，…，(θ i ,t i )，…，(θ n ,t n ) I=0, 1,..n, where t i For the current data recording time, θ i At t i Time-of-day corresponding OX b The shaft rotation angle value, in the cubic spline interpolation process, satisfies: on each segment interval [ t ] i ,t i+1 ]On f (t) =f i (t) is a cubic equation; interpolation condition f (t i )＝θ i The method comprises the steps of carrying out a first treatment on the surface of the The curve is smooth, i.e., f (t), f' (t) are continuous; let each coefficient of the third equation be a i 、b i 、c i 、d i The following steps are:
Figure FDA0004087535690000021
wherein f (t) is the definition domain after cubic spline interpolation in (t) 0 ,t n ) The upper expression is a function of the angle θ=f (t) at time t, f i (t) is a segment interval [ t ] i ,t i+1 ]The upper represents a function corresponding to the angle of the moment, f i ′(t)、f i "(t) is f respectively i A first and second derivative function of (t);
s3-2: obtaining a system of linear equations for the unknowns:
s3-2-1: from f i (t i )＝a i +b i (t i -t i )+c i (t i -t i ) 2 +d i (t i -t i ) 3 ＝θ i Obtaining a i ＝θ i, wherein ,ai ,b i ,c i ,d i Are all coefficients;
s3-2-2: let h i ＝t i+1 -t i By boundary continuous condition f i (t i+1 )＝θ i+1 Pushing out:
Figure FDA0004087535690000022
from the boundary derivative continuous condition f i ′(t i+1 )＝f′ i+1 (t i+1 ) The method comprises the following steps of:
Figure FDA0004087535690000023
from f i ″(t i+1 )＝f″ i+1 (t i+1 ) The method comprises the following steps of:
Figure FDA0004087535690000024
let m i ＝f i ″(t i )＝2c i Then formula (2) is written as m i +6h i d i ＝m i+1 Then:
Figure FDA0004087535690000025
h i ,m i respectively intermediate variables;
s3-2-3: will a i ＝t i ，
Figure FDA0004087535690000031
Substituting into formula (1) to obtain:
Figure FDA0004087535690000032
will a i ，b i ，c i ，d i Substituted into (2) to obtain：
Figure FDA0004087535690000033
S3-2-4: expanding equation (5) to all time points gives m 0 ,m 1 ,m 2 ,...,m n For a linear system of equations of unknown variables, a natural boundary condition, i.e. m, is selected during interpolation 0 ＝0，m n =0, then the system of equations is:
Figure FDA0004087535690000034
to this end, in each segment interval [ t ] i ,t i+1 ]Above, a can be determined by solving the system of equations i ，b i ，c i ，d i ；
S3-3: setting a new sampling time interval
Figure FDA0004087535690000035
The interpolated angle θ is:
θ i+j ＝f i (t i +jΔt i )，j＝1,2,...,N-1
the corresponding time precision information is obtained by controlling the N value;
s3-4: for the winding OY b Rotation angle of axis gamma, around OZ b The rotation angle psi data of the shaft is subjected to data interpolation by the same cubic spline interpolation method, the total number of data before interpolation is n+1, and the number of data after interpolation is n+1.
2. The compensation method according to claim 1, wherein in step S4, at time t i Coordinate system OX of sonde b Y b Z b Relative to the coordinate system OX in which the umbrella is located n Y n Z n The displacement and angle transformation relation occurs under the limitation of the rope only under tension, and the sonde phase is calculated through the attitude angleFor speed of movement of umbrella
Figure FDA0004087535690000041
S4-1: according to the rotation sequence of the z-y-x coordinate axes, the coordinate transformation matrix from the n system to the b system is as follows:
Figure FDA0004087535690000042
wherein ,Cθ Is around OX n Rotation matrix of rotation θ, C γ Is wound around OY n Rotation matrix of rotation gamma, C ψ Is wound around OZ n The rotation matrix of rotation ψ is expressed as:
Figure FDA0004087535690000043
/>
s4-2: when only the rotation relationship of the b-series to the n-series is considered and the displacement relationship is not considered, OX n ＝[100] T The representation under b
Figure FDA0004087535690000044
Expressed as:
Figure FDA0004087535690000045
θ x represented as
Figure FDA0004087535690000046
With OX b Included angle between:
Figure FDA0004087535690000047
s4-3: when only the rotation relation of the b-series to the n-series is considered and the displacement relation is not considered, OX b ＝[100] T The representation under b
Figure FDA0004087535690000048
Expressed as:
Figure FDA0004087535690000049
vector OH is
Figure FDA00040875356900000410
In OY n Z n Projection onto a plane, n is expressed below:
OH n ＝[0 -cosγcosψ sinγ] T
according to the geometric relationship, θ y Represented as OY n Included angle with OH:
Figure FDA0004087535690000051
s4-4: let the rope length be l, then t k At moment, k=0, 1, k, nn, and θ is obtained through angle calculation x (k) And theta y (k) The position of the sonde relative to the umbrella bottom is:
Figure FDA0004087535690000052
by forward difference, t is obtained k East speed of time sonde relative to umbrella
Figure FDA0004087535690000053
North speed->
Figure FDA0004087535690000054
Figure FDA0004087535690000055
wherein ,
Figure FDA0004087535690000056
at t k East speed calculated from the processed triaxial angular speed at time +.>
Figure FDA0004087535690000057
At t k The north velocity is calculated from the processed triaxial angular velocity.
3. The compensation method according to claim 2, wherein in step S5, the GPS/beidou information of the sonde includes a north speed of the sonde with respect to a geographic coordinate system
Figure FDA0004087535690000058
East speed of sonde relative to geographic coordinate system>
Figure FDA0004087535690000059
The wind speed obtained by the GPS/Beidou information of the sonde is as follows:
Figure FDA00040875356900000510
/>
because the sonde is connected with the umbrella by virtue of the rope, the wind speed obtained by GPS measurement is actually the superposition of the motion of the sonde relative to the umbrella and the wind speed, so the wind speed obtained by GPS measurement is decomposed into:
Figure FDA00040875356900000511
wherein ,
Figure FDA00040875356900000512
measuring and calculating for inertial systemThe speed of the sonde relative to the umbrella, i.e. the n-system,>
Figure FDA00040875356900000513
the wind speed is the wind speed under the coordinate system of the umbrella lower point, namely the i system;
t i speed corresponding to time
Figure FDA00040875356900000514
And->
Figure FDA00040875356900000515
A speed at which the corresponding moment can be found +.>
Figure FDA00040875356900000516
And->
Figure FDA00040875356900000517
The wind speed at this time is:
Figure FDA00040875356900000518
for V w Processing to obtain the final compensated wind speed V wind ：
Figure FDA0004087535690000061
Wherein floor (g) is rounded down.
CN202210317866.XA 2022-03-29 2022-03-29 Air speed measurement error compensation method of sonde based on inertial system Active CN114675055B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210317866.XA CN114675055B (en) 2022-03-29 2022-03-29 Air speed measurement error compensation method of sonde based on inertial system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210317866.XA CN114675055B (en) 2022-03-29 2022-03-29 Air speed measurement error compensation method of sonde based on inertial system

Publications (2)

Publication Number Publication Date
CN114675055A CN114675055A (en) 2022-06-28
CN114675055B true CN114675055B (en) 2023-05-05

Family

ID=82075227

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210317866.XA Active CN114675055B (en) 2022-03-29 2022-03-29 Air speed measurement error compensation method of sonde based on inertial system

Country Status (1)

Country Link
CN (1) CN114675055B (en)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110221333A (en) * 2019-04-11 2019-09-10 同济大学 A kind of error in measurement compensation method of vehicle-mounted INS/OD integrated navigation system
CN113701747A (en) * 2021-07-20 2021-11-26 北京航天控制仪器研究所 Inertial measurement system attitude angle error separation method based on centrifuge excitation

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104252010B (en) * 2013-06-27 2019-05-31 深圳航天东方红海特卫星有限公司 A kind of radiosonde and its meteorological data measurement method
CN103472503B (en) * 2013-07-24 2016-08-10 中国人民解放军理工大学 Sonde and upper air wind finding method based on INS
CN106290969B (en) * 2015-05-12 2019-01-18 湖北航天飞行器研究所 A kind of wind speed and direction detection method considering drag parachute aerodynamic influence
CN107358006B (en) * 2017-07-25 2021-10-22 华北电力大学（保定） Lorenz disturbance wind speed prediction method based on principal component analysis
US10775532B2 (en) * 2018-01-26 2020-09-15 Innovative Automation Technologies, LLC Adaptive guided wind sonde
US11415592B2 (en) * 2019-06-20 2022-08-16 James Eagleman Computing device and related methods for determining wind speed values from local atmospheric events
CN110673228B (en) * 2019-08-30 2020-09-01 北京航空航天大学 Formula of throwing sonde under imitative dandelion structure
CN113311436B (en) * 2021-04-30 2022-07-12 中国人民解放军国防科技大学 Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform
CN113671598B (en) * 2021-08-13 2024-03-19 中国人民解放军国防科技大学 Combined high-altitude wind detection method
CN114019583B (en) * 2021-10-29 2024-04-30 北京理工大学 High-precision wind measuring system and method based on inertial compensation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110221333A (en) * 2019-04-11 2019-09-10 同济大学 A kind of error in measurement compensation method of vehicle-mounted INS/OD integrated navigation system
CN113701747A (en) * 2021-07-20 2021-11-26 北京航天控制仪器研究所 Inertial measurement system attitude angle error separation method based on centrifuge excitation

Also Published As

Publication number Publication date
CN114675055A (en) 2022-06-28

Similar Documents

Publication Publication Date Title
CN111878064B (en) Gesture measurement method
CN108225308B (en) Quaternion-based attitude calculation method for extended Kalman filtering algorithm
CN1740746B (en) Micro-dynamic carrier attitude measuring apparatus and measuring method thereof
CN108362282B (en) Inertial pedestrian positioning method based on self-adaptive zero-speed interval adjustment
CN107270893B (en) Lever arm and time asynchronous error estimation and compensation method for real estate measurement
CN111024064B (en) SINS/DVL combined navigation method for improving Sage-Husa adaptive filtering
CN109425339B (en) Ship heave error compensation method considering lever arm effect based on inertia technology
CN103175530B (en) Method for estimating and compensating coupling torque of aerial remote sensing inertially stabilized platform
CN108036785A (en) A kind of aircraft position and orientation estimation method based on direct method and inertial navigation fusion
CN102393201B (en) Dynamic lever arm compensating method of position and posture measuring system (POS) for aerial remote sensing
CN106289246A (en) A kind of rods arm measure method based on position and orientation measurement system
CN103245360A (en) Autocollimation method of carrier aircraft rotating type strapdown inertial navigation system under shaking base
CN109682377A (en) A kind of Attitude estimation method based on the decline of dynamic step length gradient
CN112923924B (en) Method and system for monitoring posture and position of anchoring ship
CN109709628B (en) Calibration method for gravity gradiometer of rotating accelerometer
CN106403952A (en) Method for measuring combined attitudes of Satcom on the move with low cost
CN104698485A (en) BD, GPS and MEMS based integrated navigation system and method
CN104182614A (en) System and method for monitoring attitude of mechanical arm with six degrees of freedom
CN109443342A (en) NEW ADAPTIVE Kalman's UAV Attitude calculation method
CN109931957A (en) SINS self-alignment method for strapdown inertial navigation system based on LGMKF
CN103017787A (en) Initial alignment method suitable for rocking base
CN108007477A (en) A kind of inertia pedestrian's Positioning System Error suppressing method based on forward and reverse filtering
CN112556724A (en) Initial coarse alignment method for low-cost navigation system of micro aircraft in dynamic environment
CN105115505A (en) Two-rank dynamic disturbance torque compensation method of four-axis inertial stabilization platform system
CN111766397A (en) Meteorological wind measurement method based on inertia/satellite/atmosphere combination

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