CN114675055A - Air speed measurement error compensation method for sonde based on inertial system - Google Patents
Air speed measurement error compensation method for sonde based on inertial system Download PDFInfo
- Publication number
- CN114675055A CN114675055A CN202210317866.XA CN202210317866A CN114675055A CN 114675055 A CN114675055 A CN 114675055A CN 202210317866 A CN202210317866 A CN 202210317866A CN 114675055 A CN114675055 A CN 114675055A
- Authority
- CN
- China
- Prior art keywords
- sonde
- axis
- wind speed
- angle
- umbrella
- 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
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000005259 measurement Methods 0.000 title claims abstract description 40
- 238000001514 detection method Methods 0.000 claims abstract description 11
- 230000008569 process Effects 0.000 claims description 14
- 230000006870 function Effects 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000006073 displacement reaction Methods 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 230000001133 acceleration Effects 0.000 claims description 6
- 238000006467 substitution reaction Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 238000004804 winding Methods 0.000 claims description 6
- 238000011065 in-situ storage Methods 0.000 claims description 4
- 230000007613 environmental effect Effects 0.000 claims description 3
- 238000009499 grossing Methods 0.000 claims description 3
- 238000005096 rolling process Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000011218 segmentation Effects 0.000 claims description 3
- 230000003068 static effect Effects 0.000 claims description 3
- HLLICFJUWSZHRJ-UHFFFAOYSA-N tioxidazole Chemical compound CCCOC1=CC=C2N=C(NC(=O)OC)SC2=C1 HLLICFJUWSZHRJ-UHFFFAOYSA-N 0.000 claims description 3
- 238000005266 casting Methods 0.000 claims description 2
- 230000000717 retained effect Effects 0.000 claims description 2
- 238000012360 testing method Methods 0.000 abstract description 4
- 230000000694 effects Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 206010034719 Personality change Diseases 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
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)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Atmospheric Sciences (AREA)
- Biodiversity & Conservation Biology (AREA)
- Ecology (AREA)
- Environmental Sciences (AREA)
- Navigation (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
The invention discloses a method for compensating a wind speed measurement error of a sonde based on an inertial system. The error compensation is carried out by utilizing the measured data of the inertial system, the high-precision compensation can be carried out point by point according to the actual situation of each test of the sonde, and the wind speed detection precision is improved.
Description
Technical Field
The invention belongs to the technical field of meteorological detection, and particularly relates to a method for compensating wind speed measurement errors of a sonde based on an inertial system.
Background
Typhoon is one of the most serious natural disasters in the world. In order to realize more essential scientific cognition on the natural phenomenon of the typhoon, the key technical bottleneck is to directly detect the typhoon 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 a typhoon kernel area, wherein the measurement of the wind speed of the typhoon kernel area is particularly important for predicting the path and the strength of the typhoon. At present, wind profile radar, laser wind measuring radar, sonde and the like are mainly used for measuring the typhoon wind speed, and the sonde can only be used for measuring the typhoon core area wind speed in order to meet the requirement of in-situ detection. However, the traditional sonde equipment can only solve the wind speed by depending on the position change data detected by the sonde GPS or the Beidou system, and the error is large. When the sonde is in an extreme environment such as typhoon, the sonde can rapidly roll and rotate around a parachute in the falling process, the severe change of the posture of the sonde can cause large error influence on wind speed measurement, 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 moving posture of the sonde must be eliminated.
At present, a great number of scholars studying high-altitude exploration at home and abroad make a great deal of research on wind speed measurement algorithms of navigation sondes. The Jaakko study VIASALA RS92 sonde indicates that wind speed measurements are to be data-dependent and smoothed. Jauhiaine et al in summarizing the international sonde alignment test indicate that the sonde oscillation factor should be taken into account in the wind speed measurement process, but the method of processing is not mentioned. According to a calculation method of a compound pendulum period, Zhanhong researches a swing period of the sonde, and influences of the swing of the sonde on wind speed measurement are removed by using a filtering method. However, the methods are used for optimizing the motion of the sonde into simple and periodic conical pendular motion so as to filter all detection data, and the precision is low. When the sonde is in an environment with rapid change of typhoon and extremely strong wind speed, the movement of the sonde is complex and unknown, the problem of under-compensation or over-compensation and the like can be caused by only using periodic conical pendulum movement to compensate detection data, and the wind speed and direction information of the sonde in the current space-time can not be accurately obtained.
The patent application CN114019583A discloses a high-precision safe wind measuring system and method based on inertia compensation, wherein the speed is obtained by integrating the acceleration of a sonde measured by an inertia module, and the wind speed measured by a positioning module is compensated in a weighted average mode; and resolving the swing period of the sonde according to the attitude information of the sonde measured by the inertia module, and carrying out multi-point filtering on data in different periods to filter the influence of the swing of the sonde. First, it is incorrect and meaningless to use the integration of the sonde acceleration measured by the inertial module to obtain the velocity for the first point. In the sonde system, two components, namely a parachute and a sonde, mainly exist, the wind speed compensation is carried out, the wind speed compensation meaning is mainly to compensate the complex movement of the sonde around the parachute, if the wind speed is solved by simply utilizing the acceleration information obtained by the inertia module, the instantaneous wind speed of the whole sonde system is equivalent to the measurement, the effect is the same as that of the GPS measurement, the measurement error caused by the movement of the sonde around the parachute is not solved, and the wind speed compensation effect cannot be achieved. For the second point, attitude information measured in an inertial system must be used for solving errors, but the method only solves the swing period of the sonde, and the method has the same problem with the previously mentioned research article, namely the problem that overcompensation or undercompensation exists when the calculation period is used for compensating wind measurement data at all moments, the influence of consistent swing period is not required to be compensated at every moment, each moment has own unique characteristic, and wind speed errors measured by the sonde need to be compensated one by one according to actual conditions.
In summary, it is necessary to design a method for compensating the wind speed measurement error of the sonde, so as to compensate the wind speed measurement error caused by the attitude change of the sonde.
Disclosure of Invention
The invention provides a method for compensating a wind speed measurement error of a sonde based on an inertial system, aiming at the wind speed measurement error of the current sonde caused by self attitude change under severe working conditions. And (4) solving the true motion attitude track of the sonde by combining 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:
a method for compensating a wind speed measurement error of a sonde based on an inertial system is characterized in that a lower dropping structure in a wind speed measurement process comprises an umbrella, a rope and the sonde, the sonde is connected with the bottom of the umbrella through the rope, the umbrella provides floating time for the sonde, and the rope keeps a tensioned state in the lower dropping process, and the compensation method comprises the following steps:
s1: an inertial system is integrated on the sonde to perform high-altitude downward projection detection to obtain meteorological data of a three-dimensional space, and a ground receiver records and stores in-situ detection data of the sonde and triaxial inertial navigation information in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received GPS/Beidou information of the sonde;
s3: calculating attitude information of the sonde in the falling process by utilizing the received sonde triaxial inertial navigation information;
s4: converting the result of the step S3 through a space coordinate system and calculating the angle to obtain the real-time relative space motion attitude information of the sonde;
s5: calculating the wind speed measurement error of the sonde caused by the posture change of the sonde through the relative space motion 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 complete wind speed compensation;
s6: the wind speed data compensated in step S5 is subjected to rolling time window smoothing processing to further eliminate the influence of environmental noise.
Further, the triaxial inertial navigation information in the step S3 includes triaxial acceleration, triaxial angular velocity and triaxial angle information, the triaxial angle information with an angle range of (-180 °,180 °) is retained, the triaxial angle original data is interpolated by cubic spline interpolation, and then the differential solution is performed on the shift to obtain the motion velocity;
the rotation sequence of the coordinate system when the inertial system uses the Euler angle to express the posture is that the inertial system firstly rotates around a Z axis, then rotates around a Y axis and finally rotates around an X axis, so that the angle range of the rotation around the Y axis is only +/-90 degrees;
using the center of mass of the sonde as an origin O and using the direction from the center of mass of the sonde to the connecting point of the rope and the sonde as OXbThe axis is OZ from the center of mass of the sonde to the transverse center line of the sondebShaft to OXbAxis and OZbAxis perpendicular direction is OYbAxis establishing a coordinate system OX fixed to the sondebYbZbI.e., system b;
the umbrella bottom is taken as an origin O, and the east direction of the geographic coordinate system pointed from the umbrella bottom is taken as OYnShaft to point from the bottom of the umbrella to the groundNorth direction of physical coordinate system is OZnAxis, perpendicular to the umbrella base to the geographic coordinate system as OXnAxis establishing a coordinate system OX solidly connected to the base of the umbrellanYnZnI.e. n is; when the sonde is static, namely the sonde and the umbrella do not swing relatively, the b series can be obtained by translating the n series along the rope;
the bottom of the umbrella is the original point O, OX at the time of castingiAxis pointing to sky, OYiAxis pointing east, OZiThe axis points to the north direction, and a space rectangular coordinate system OX is establishediYiZiI.e., the i series;
definition of cords and OXiAngle of axis thetaxRope and OYiAngle of axis thetay;
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for winding OXbThe axis rotation angle theta is, n +1 data nodes exist in the original data: (theta)0,t0), (θ1,t1),…,(θi,ti),…,(θn,tn) I is 0, 1.. times.n, wherein t isiRecording time, θ, for the current dataiIs tiOX corresponding to timebThe shaft rotation angle value meets the following requirements in the cubic spline interpolation process: on each segment interval ti,ti+1]F (t) fi(t) are all cubic equations; interpolation condition f (t)i)=θi(ii) a The curve is smooth, i.e. f (t), f' (t) are continuous; let each coefficient of cubic equation be ai、bi、ci、diThen, there are:
wherein f (t) is the definition domain after cubic spline interpolation in (t)0,tn) The above represents the function of the angle θ at time t ═ f (t), fi(t) is a segmentation interval [ t ]i,ti+1]A function representing the angle of the corresponding time instant, fi′(t)、fi"(t) is each fi(t) first and second derivative functions;
s3-2: a system of linear equations for the unknowns is obtained:
s3-2-1: from fi(ti)=ai+bi(ti-ti)+ci(ti-ti)2+di(ti-ti)3=θiTo obtain ai=θi, wherein , ai,bi,ci,diAre all coefficients;
s3-2-2: let hi=ti+1-tiFrom a boundary continuation condition fi(ti+1)=θi+1And (3) pushing out:
continuous condition f by boundary derivativei′(ti+1)=f′i+1(ti+1) Obtaining:
from fi″(ti+1)=f″i+1(ti+1) Obtaining:
let mi=fi″(ti)=2ciThen formula (2) is written as mi+6hidi=mi+1And then:hi,mirespectively intermediate variables;
a is toi,bi,ci,diSubstitution in formula (2) gives:
s3-2-4: extending equation (5) to all time points to give m0,m1,m2,...,mnFor a linear system of equations of unknown variables, a natural boundary condition, m, is selected during interpolation0=0,mnWhen 0, the system of equations is:
to this end, in each segment interval [ t ]i,ti+1]Above, a can be determined by solving a system of equationsi,bi,ci,di;
θi+j=fi(ti+jΔti),j=1,2,...,N-1
obtaining corresponding time precision information by controlling the N value;
s3-4: for winding OYbRotation angle of axis gamma, about OZbAnd carrying out data interpolation on the data of the rotation angle psi of the axis by the same cubic spline interpolation method, wherein n +1 data are obtained before the interpolation, and n +1 data are obtained after the interpolation.
Further, the stepsIn step S4, at time tiCoordinate system OX where sonde is locatedbYbZbRelative to the coordinate system OX in which the umbrella is locatednYnZnThe displacement and angle transformation relation is generated under the limitation of a tensioned rope, and the movement speed of the sonde relative to the umbrella is calculated through the attitude angle
S4-1: according to the rotation sequence of the z-y-x coordinate axis, the coordinate transformation matrix from the n system to the b system is as follows:
wherein ,CθIs wound around OXnRotation matrix of rotation theta, CγIs wound around OYnRotation matrix of rotation gamma, CψIs wound around OZnThe rotation matrices of the rotation psi are respectively expressed as:
s4-2: OX when only the rotational relationship of the b system with respect to the n system is considered and the displacement relationship is not consideredn=[1 0 0]TRepresentation under bExpressed as:
s4-3: OX when only the rotational relationship of the b system with respect to the n system is considered and the displacement relationship is not consideredb=[1 0 0]TRepresentation under bExpressed as:
OHn=[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
s4-4: let the length of the rope be l, then tkThe time, k is 0,1,.. cne, and θ is obtained through angle calculationx(k) And thetay(k) Then, the position of the sonde relative to the umbrella bottom is as follows:
wherein ,is tkThe east velocity calculated from the processed three-axis angular velocities at that moment,is tkAnd calculating the north direction velocity from the processed three-axis angular velocity.
Further, in step S5, the GPS/beidou information of the sonde includes the north speed of the sonde relative to the geographic coordinate systemEast speed of sonde relative to geographic coordinate systemThe wind speed obtained from the GPS/beidou information of the sonde is:
because the sonde is connected with the umbrella by a rope, the wind speed measured and calculated by the GPS is actually the superposition of the movement of the sonde relative to the umbrella and the wind speed, and therefore the wind speed measured and calculated by the GPS is decomposed into:
wherein ,the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to timeAndthe speed of the corresponding moment can be foundAndthe wind speed at this time is:
to VwProcessing to obtain the final compensated wind speed Vwind:
Wherein floor (. cndot.) is rounded downward.
The invention has the beneficial effects that:
1. according to the method for compensating the wind speed measurement error of the sonde based on the inertial system, the inertial component is integrated on the traditional sonde system, so that the real attitude information of the movement of the sonde is obtained, and the wind speed measurement error caused by the movement of the sonde is further compensated.
2. The method utilizes the measured data of the inertial system to carry out error compensation, has more authenticity and reliability compared with attitude errors estimated by other methods, can compensate point by point with high precision according to the actual situation of each test of the sonde, compensates the real-time wind speed information measured by the sonde, and avoids the problems of overcompensation and undercompensation by the method instead of only estimating the approximate swing period compensation wind speed of the sonde.
Drawings
In order to illustrate embodiments of the present invention or technical solutions in the prior art more clearly, the drawings which are needed in the embodiments will be briefly described below, so that the features and advantages of the present invention can be understood more clearly by referring to the drawings, which are schematic and should not be construed as limiting the present invention in any way, and for a person skilled in the art, other drawings can be obtained on the basis of these drawings without any inventive effort. Wherein:
FIG. 1 is a flow chart of the inertial system compensation sonde wind speed measurement error of the present invention;
FIG. 2 is a coordinate system definition of the sonde and the parachute;
FIG. 3 is an attitude definition;
FIG. 4 is a spatial coordinate system position relationship;
fig. 5 is a measured movement trajectory of the sonde relative to the parachute;
FIG. 6 is a diagram of a sonde in a self-developed inertial navigation system;
FIG. 7 shows the mounting position of the inertial navigation module on the sonde and the three-axis definition;
fig. 8 shows compensation results of the wind speeds in the east and north directions of the sonde, where (a) is the compensation result of the wind speed in the east direction of the sonde, the dotted line is the wind speed in the east direction calculated by the GPS, the solid line is the wind speed in the east direction compensated by the inertial system, and (b) is the compensation result of the wind speed in the north direction of the sonde, the dotted line is the wind speed in the north direction calculated by the GPS, and the solid line is the wind speed in the north direction compensated by the inertial system;
fig. 9 shows compensation results of the sonde after east and north wind speed vector synthesis, the dotted line represents wind speed obtained by GPS calculation, and the solid line represents wind speed compensated by an inertial system.
Detailed Description
In order that the above objects, features and advantages of the present invention can be more clearly understood, a more particular description of the invention will be rendered by reference to the appended drawings. It should be noted that the embodiments and features of the embodiments of the present invention may be combined with each other without conflict.
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 specifically described herein, and therefore the scope of the present invention is not limited by 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, thereby compensating the wind speed measurement error of the sonde caused by self motion under severe working conditions.
As shown in fig. 1, a method for compensating a wind speed measurement error of a sonde based on an inertial system, wherein a dropsonde structure in a wind speed measurement process comprises an umbrella, a rope and the sonde, the sonde is connected with an umbrella bottom through the rope, the umbrella provides floating time for the sonde, and the rope keeps a tensioned state in the dropsonde process, and the compensation method comprises the following steps:
s1: an inertial system is integrated on the sonde to perform high-altitude downward projection detection to obtain meteorological data of a three-dimensional space, and a ground receiver records and stores in-situ detection data of the sonde and triaxial inertial navigation information in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received GPS/Beidou information of the sonde;
s3: calculating attitude information of the sonde in the falling process by utilizing the received sonde triaxial inertial navigation information;
the triaxial inertial navigation information in the step S3 comprises triaxial acceleration, triaxial angular velocity and triaxial angle information, the triaxial angle information with an angle range of (-180 degrees and 180 degrees) is reserved, the triaxial angle original data is interpolated through cubic spline interpolation, and the motion velocity is obtained by performing differential solution on the interpolated value;
as shown in fig. 2, the rotation sequence of the coordinate system when the inertial system uses the euler angle to represent the posture is that the inertial system firstly rotates around the Z axis, then rotates around the Y axis, and finally rotates around the X axis, so that the rotation angle around the Y axis is only +/-90 degrees;
using the center of mass of the sonde as the origin O to perform secondary explorationThe direction of the connecting point of the mass center pointing rope and the sonde of the sonde is OXbThe axis is OZ from the center of mass of the sonde to the transverse center line of the sondebShaft to OXbAxis and OZbAxis perpendicular direction is OYbShaft establishing a coordinate system OX fixed to the sondebYbZbI.e., system b;
the umbrella bottom is taken as an origin O, and the east direction of the geographic coordinate system pointed from the umbrella bottom is taken as OYnAxis, north direction from umbrella bottom to geographic coordinate system as OZnAxis, perpendicular to the umbrella base to the geographic coordinate system as OXnAxis establishing a coordinate system OX solidly connected to the base of the umbrellanYnZnI.e. n is; when the sonde is static, namely the sonde and the umbrella do not swing relatively, the b series can be obtained by translating the n series along the rope;
the bottom of the umbrella is the original point O, OXiAxis pointing to sky, OYiAxis pointing east, OZiThe axis points to the north direction, and a space rectangular coordinate system OX is establishediYiZiI.e., the i series;
definition of cords and OXiAngle of axis thetaxRope and OYiAngle of axis thetay;
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for winding OXbThe axis rotation angle theta is, n +1 data nodes exist in the original data: (theta)0,t0), (θ1,t1),…,(θi,ti),…,(θn,tn) I is 0, 1.. times.n, wherein t isiRecording time, θ, for the current dataiIs tiOX corresponding to timebThe shaft rotation angle value meets the following requirements in the cubic spline interpolation process: on each segment interval ti,ti+1]F (t) fi(t) are all cubic equations; interpolation condition f (t)i)=θi(ii) a The curve is smooth, i.e. f (t), f' (t) are continuous; let each coefficient of cubic equation be ai、bi、ci、diThen, there are:
wherein f (t) is the definition domain after cubic spline interpolation in (t)0,tn) The above represents the function of the angle θ at time t ═ f (t), fi(t) is a segmentation interval [ t ]i,ti+1]A function representing the angle of the corresponding time instant, fi′(t)、fi"(t) is each fi(t) first and second derivative functions;
s3-2: a system of linear equations for the unknowns is obtained:
s3-2-1: from fi(ti)=ai+bi(ti-ti)+ci(ti-ti)2+di(ti-ti)3=θiTo obtain ai=θi, wherein , ai,bi,ci,diAre all coefficients;
s3-2-2: let hi=ti+1-tiFrom a boundary continuation condition fi(ti+1)=θi+1And (3) pushing out:
continuous condition f by boundary derivativei′(ti+1)=f′i+1(ti+1) Obtaining:
from f to fi″(ti+1)=f″i+1(ti+1) Obtaining:
let m bei=fi″(ti)=2ciThen formula (2) is written as mi+6hidi=mi+1And then:hi,mirespectively intermediate variables;
a is toi,bi,ci,diSubstitution in formula (2) gives:
s3-2-4: extending equation (5) to all time points to give m0,m1,m2,...,mnSelecting a natural boundary condition, m, for a system of linear equations of unknown variables during interpolation0=0,mnWhen 0, the system of equations is:
to this end, in each segment interval [ t ]i,ti+1]Above, a can be determined by solving a system of equationsi,bi,ci,di;
θi+j=fi(ti+jΔti),j=1,2,...,N-1
obtaining corresponding time precision information by controlling the N value;
s3-4: for winding OYbRotation angle of axis gamma, about OZbAnd carrying out data interpolation on the data of the rotation angle psi of the axis by the same cubic spline interpolation method, wherein n +1 data are obtained before the interpolation, and n +1 data are obtained after the interpolation.
S4: converting the result of the step S3 through a space coordinate system and calculating the angle to obtain the real-time relative space motion attitude information of the sonde; in step S4, at time tiThe coordinate system OX of the sondebYbZbRelative to the coordinate system OX in which the umbrella is locatednYnZnThe displacement and angle transformation relation occurs under the limitation of the rope which is only tensioned, and the movement speed of the sonde relative to the umbrella is calculated through the attitude angle as shown in figure 4
S4-1: according to the rotation sequence of the z-y-x coordinate axis, the coordinate transformation matrix from the n system to the b system is as follows:
wherein ,CθIs wound around OXnRotation matrix of rotation theta, CγIs wound around OYnRotation matrix of rotation gamma, CψIs wound around OZnThe rotation matrices of the rotation psi are respectively expressed as:
s4-2: considering only the rotation relationship of b system relative to n system and not considering bitWhen moving in relation to each other, OXn=[1 0 0]TRepresentation under bExpressed as:
s4-3: OX when only the rotational relationship of the b system with respect to the n system is considered and the displacement relationship is not consideredb=[1 0 0]TRepresentation under bExpressed as:
OHn=[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
s4-4: the length of the rope is set asl, then tkThe time, k is 0,1,.. cne, and θ is obtained through angle calculationx(k) And thetay(k) Then, the position of the sonde relative to the umbrella bottom is as follows:
wherein ,is tkThe east velocity calculated from the processed three-axis angular velocities at that moment,is tkAnd calculating the north direction velocity from the processed three-axis angular velocity.
S5: calculating the wind speed measurement error of the sonde caused by the posture change of the sonde through the relative space motion 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 complete wind speed compensation; in step S5, the GPS/beidou information of the sonde includes the north speed of the sonde relative to the geographic coordinate systemEast speed of sonde relative to geographic coordinate systemThe wind speed obtained from the GPS/beidou information of the sonde is:
because the sonde is connected with the umbrella by a rope, the wind speed measured and calculated by the GPS is actually the superposition of the movement of the sonde relative to the umbrella and the wind speed, and therefore the wind speed measured and calculated by the GPS is decomposed into:
wherein ,the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to timeAndthe speed of the corresponding moment can be foundAnd withThe wind speed at this time is:
to VwProcessing to obtain the final compensated wind speed Vwind:
Wherein floor (. cndot.) is rounded downward.
S6: the wind speed data compensated in step S5 is subjected to rolling time window smoothing processing, thereby further eliminating the influence of environmental noise.
For the convenience of understanding the above technical aspects of the present invention, the following detailed description will be given of the above technical aspects of the present invention by way of specific examples.
Example 1
As shown in FIGS. 6-7, a 20km drop test was performed in Xinjiang using a sonde with a self-developed inertial navigation system to detect meteorological data in a three-dimensional space. The actual measurement movement track of the sonde relative to the parachute is shown in fig. 5, wind speed measurement drift caused by random swing of the sonde is obtained through calculation, and east speed and north speed obtained through GPS measurement and calculation are compensated respectively. The compensation effect is shown in fig. 8, where (a) is the sonde east wind speed compensation result, the dotted line is the east wind speed calculated by the GPS, the solid line is the east wind speed compensated by the inertial system, (b) is the sonde north wind speed compensation result, the dotted line is the north wind speed calculated by the GPS, and the solid line is the north wind speed compensated by the inertial system. It can be seen that the signal-to-noise ratios of the east and north wind speeds compensated by the inertial navigation data are obviously improved. The wind speed measured by the GPS and the wind speed compensated by the inertial navigation are subjected to vector synthesis to obtain a final wind speed measurement result as shown in FIG. 9, wherein the dotted line is the wind speed calculated by the GPS, and the solid line is the wind speed compensated by the inertial system, so that the signal-to-noise ratio of the wind speed compensated by the inertial navigation data is obviously improved.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (4)
1. A method for compensating a wind speed measurement error of a sonde based on an inertial system is characterized in that a lower dropping structure in a wind speed measurement process comprises an umbrella, a rope and the sonde, the sonde is connected with the bottom of the umbrella through the rope, the umbrella provides floating time for the sonde, and the rope keeps a tensioned state in the lower dropping process, and the method comprises the following steps:
s1: an inertial system is integrated on the sonde to perform high-altitude downward projection detection to obtain meteorological data of a three-dimensional space, and a ground receiver records and stores in-situ detection data of the sonde and triaxial inertial navigation information in real time;
s2: calculating the three-dimensional space wind speed information detected by the sonde by utilizing the received GPS/Beidou information of the sonde;
s3: calculating attitude information of the sonde in the falling process by utilizing the received sonde triaxial inertial navigation information;
s4: converting the result of the step S3 through a space coordinate system and calculating the angle to obtain the real-time relative space motion attitude information of the sonde;
s5: calculating the wind speed measurement error of the sonde caused by the posture change of the sonde through the relative space motion 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 complete wind speed compensation;
s6: the wind speed data compensated in step S5 is subjected to rolling time window smoothing processing, thereby further eliminating the influence of environmental noise.
2. The compensation method of claim 1, wherein the triaxial inertial navigation information in step S3 includes triaxial acceleration, triaxial angular velocity and triaxial angle information, the triaxial angle information with an angle range of (-180 °,180 °) is retained, the triaxial angle raw data is interpolated by cubic spline interpolation, and the motion velocity is obtained by performing differential solution on the shift;
the rotation sequence of the coordinate system when the inertial system uses the Euler angle to express the posture is that the inertial system firstly rotates around a Z axis, then rotates around a Y axis and finally rotates around an X axis, so that the angle range of the rotation around the Y axis is only +/-90 degrees;
using the center of mass of the sonde as an origin O and using the direction pointing from the center of mass of the sonde to the connecting point of the rope and the sonde as OXbThe axis is OZ from the center of mass of the sonde to the transverse center line of the sondebShaft to OXbAxis and OZbAxis perpendicular direction is OYbShaft establishing a coordinate system OX fixed to the sondebYbZbI.e., system b;
the umbrella bottom is taken as an origin O, and the east direction pointing from the umbrella bottom to the geographic coordinate system is taken as OYnAxis, north direction from umbrella bottom to geographic coordinate system as OZnAxis, perpendicular to the umbrella base to the geographic coordinate system as OXnAxis establishing a coordinate system OX solidly connected to the base of the umbrellanYnZnI.e. n is; when the sonde is static, namely the sonde and the umbrella do not swing relatively, the system b can be obtained by translating the system n along the rope;
the bottom of the umbrella is the original point O, OX at the time of castingiAxis pointing to sky, OYiAxis pointing east, OZiThe axis points to the north direction, and a space rectangular coordinate system OX is establishediYiZiI.e., the i series;
definition of cords and OXiAngle of axis thetaxRope and OYiAngle of axis thetay;
The cubic spline interpolation method for the original data comprises the following steps:
s3-1: for winding OXbThe axis rotation angle theta is, n +1 data nodes exist in the original data: (theta)0,t0),(θ1,t1),…,(θi,ti),…,(θn,tn) I is 0, 1.. times.n, wherein t isiRecording time, θ, for the current dataiIs tiOX corresponding to timebThe shaft rotation angle value meets the following requirements in the cubic spline interpolation process: on each segment interval ti,ti+1]F (t) fi(t) are all cubic equations; interpolation condition f (t)i)=θi(ii) a The curve is smooth, i.e. f (t), f' (t) are continuous; let the terms of cubic equationCoefficient ai、bi、ci、diThen, there are:
wherein f (t) is the definition domain after cubic spline interpolation in (t)0,tn) The above represents the function of the angle θ at time t ═ f (t), fi(t) is a segmentation interval [ t ]i,ti+1]Is represented by a function of the angle of the corresponding instant, f'i(t)、f″i(t) are each fi(t) first and second derivative functions;
s3-2: a system of linear equations for the unknowns is obtained:
s3-2-1: from fi(ti)=ai+bi(ti-ti)+ci(ti-ti)2+di(ti-ti)3=θiTo obtain ai=θi, wherein ,ai,bi,ci,diAre all coefficients;
s3-2-2: let hi=ti+1-tiFrom the boundary condition fi(ti+1)=θi+1And (3) pushing out:
condition f 'is continued by the boundary derivative'i(ti+1)=f′iv1(ti+1) Obtaining:
from f ″)i(ti+1)=f″i+1(ti+1) Obtaining:
let m bei=f″i(ti)=2ciThen formula (2) is written as mi+6hidi=mi+1And then:hi,mirespectively intermediate variables;
a is toi,bi,ci,diSubstitution in formula (2) gives:
s3-2-4: extending equation (5) to all time points to give m0,m1,m2,...,mnSelecting a natural boundary condition, m, for a system of linear equations of unknown variables during interpolation0=0,mnWhen 0, the system of equations is:
to this end, in each segment interval [ t ]i,ti+1]Above, a can be determined by solving a system of equationsi,bi,ci,di;
θi+j=fi(ti+jΔti),j=1,2,...,N-1
obtaining corresponding time precision information by controlling the N value;
s3-4: for winding OYbRotation angle of axis gamma, about OZbAnd carrying out data interpolation on the data of the rotation angle psi of the axis by the same cubic spline interpolation method, wherein n +1 data are obtained before the interpolation, and n +1 data are obtained after the interpolation.
3. The compensation method as claimed in claim 2, wherein in step S4, at time tiThe coordinate system OX of the sondebYbZbRelative to the coordinate system OX in which the umbrella is locatednYnZnThe displacement and angle transformation relation is generated under the limitation of a tensioned rope, and the movement speed of the sonde relative to the umbrella is calculated through the attitude angle
S4-1: according to the rotation sequence of the z-y-x coordinate axis, the coordinate transformation matrix from the n system to the b system is as follows:
wherein ,CθIs wound around OXnRotation matrix of rotation theta, CγIs wound around OYnRotation matrix of rotation gamma, CψIs wound around OZnThe rotation matrices of the rotation psi are respectively expressed as:
s4-2: OX when only the rotational relationship of the b system with respect to the n system is considered and the displacement relationship is not consideredn=[1 0 0]TRepresentation under bExpressed as:
s4-3: OX when only the rotational relationship of the b system with respect to the n system is considered and the displacement relationship is not consideredb=[1 0 0]TRepresentation under bExpressed as:
OHn=[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
s4-4: let the length of the rope be l, then tkThe time, k is 0,1,.. cne, and θ is obtained through angle calculationx(k) And thetay(k) Then, the position of the sonde relative to the umbrella bottom is as follows:
4. The compensation method of claim 3, wherein in step S5, the GPS/Beidou information of the sonde comprises the north velocity of the sonde relative to the geographic coordinate systemEast speed of sonde relative to geographic coordinate systemThe wind speed obtained from the GPS/beidou information of the sonde is:
because the sonde is connected with the umbrella by a rope, the wind speed measured and calculated by the GPS is actually the superposition of the movement of the sonde relative to the umbrella and the wind speed, and therefore the wind speed measured and calculated by the GPS is decomposed into:
wherein ,the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to timeAndthe speed of the corresponding moment can be foundAnd withThe wind speed at this time is:
to VwProcessing to obtain the final compensated wind speed Vwind:
Wherein floor (. cndot.) is rounded downward.
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 true CN114675055A (en) | 2022-06-28 |
CN114675055B 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 (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103472503A (en) * | 2013-07-24 | 2013-12-25 | 中国人民解放军理工大学 | Sonde and upper-air-wind detecting method based on INS |
CN104252010A (en) * | 2013-06-27 | 2014-12-31 | 深圳航天东方红海特卫星有限公司 | Radiosonde and weather data measuring method thereof |
CN106290969A (en) * | 2015-05-12 | 2017-01-04 | 湖北航天飞行器研究所 | A kind of wind speed and direction detection method considering drag parachute aerodynamic influence |
CN107358006A (en) * | 2017-07-25 | 2017-11-17 | 华北电力大学(保定) | A kind of Lorenz disturbance wind speed forecasting methods based on principal component analysis |
CN110221333A (en) * | 2019-04-11 | 2019-09-10 | 同济大学 | A kind of error in measurement compensation method of vehicle-mounted INS/OD integrated navigation system |
CN110673228A (en) * | 2019-08-30 | 2020-01-10 | 北京航空航天大学 | Formula of throwing sonde under imitative dandelion structure |
US20200025972A1 (en) * | 2018-01-26 | 2020-01-23 | Innovative Automation Technologies, LLC d/b/a IA Tech | Adaptive guided wind sonde |
US20200400706A1 (en) * | 2019-06-20 | 2020-12-24 | James Eagleman | Computing device and related methods for determining wind speed values from local atmospheric events |
CN113311436A (en) * | 2021-04-30 | 2021-08-27 | 中国人民解放军国防科技大学 | Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform |
CN113671598A (en) * | 2021-08-13 | 2021-11-19 | 中国人民解放军国防科技大学 | Combined high-altitude wind detection method |
CN113701747A (en) * | 2021-07-20 | 2021-11-26 | 北京航天控制仪器研究所 | Inertial measurement system attitude angle error separation method based on centrifuge excitation |
CN114019583A (en) * | 2021-10-29 | 2022-02-08 | 北京理工大学 | High-precision wind measuring system and method based on inertia compensation |
-
2022
- 2022-03-29 CN CN202210317866.XA patent/CN114675055B/en active Active
Patent Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104252010A (en) * | 2013-06-27 | 2014-12-31 | 深圳航天东方红海特卫星有限公司 | Radiosonde and weather data measuring method thereof |
CN103472503A (en) * | 2013-07-24 | 2013-12-25 | 中国人民解放军理工大学 | Sonde and upper-air-wind detecting method based on INS |
CN106290969A (en) * | 2015-05-12 | 2017-01-04 | 湖北航天飞行器研究所 | A kind of wind speed and direction detection method considering drag parachute aerodynamic influence |
CN107358006A (en) * | 2017-07-25 | 2017-11-17 | 华北电力大学(保定) | A kind of Lorenz disturbance wind speed forecasting methods based on principal component analysis |
US20200025972A1 (en) * | 2018-01-26 | 2020-01-23 | Innovative Automation Technologies, LLC d/b/a IA Tech | Adaptive guided wind sonde |
CN110221333A (en) * | 2019-04-11 | 2019-09-10 | 同济大学 | A kind of error in measurement compensation method of vehicle-mounted INS/OD integrated navigation system |
US20200400706A1 (en) * | 2019-06-20 | 2020-12-24 | James Eagleman | Computing device and related methods for determining wind speed values from local atmospheric events |
CN110673228A (en) * | 2019-08-30 | 2020-01-10 | 北京航空航天大学 | Formula of throwing sonde under imitative dandelion structure |
CN113311436A (en) * | 2021-04-30 | 2021-08-27 | 中国人民解放军国防科技大学 | Method for correcting wind measurement of motion attitude of laser wind measuring radar on mobile platform |
CN113701747A (en) * | 2021-07-20 | 2021-11-26 | 北京航天控制仪器研究所 | Inertial measurement system attitude angle error separation method based on centrifuge excitation |
CN113671598A (en) * | 2021-08-13 | 2021-11-19 | 中国人民解放军国防科技大学 | Combined high-altitude wind detection method |
CN114019583A (en) * | 2021-10-29 | 2022-02-08 | 北京理工大学 | High-precision wind measuring system and method based on inertia compensation |
Non-Patent Citations (2)
Title |
---|
李大鹏等: "基于临空飞艇的大气参数原位探测***及试验", 《测控技术》 * |
赵世军: "基于INS的高空风测量技术研究", 《第八届全国优秀青年气象科技工作者学术研讨会论文汇编中国气象学会会议论文集》 * |
Also Published As
Publication number | Publication date |
---|---|
CN114675055B (en) | 2023-05-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111878064B (en) | Gesture measurement method | |
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 | |
CN111024064B (en) | SINS/DVL combined navigation method for improving Sage-Husa adaptive filtering | |
CN104698485B (en) | Integrated navigation system and air navigation aid based on BD, GPS and MEMS | |
CN109425339B (en) | Ship heave error compensation method considering lever arm effect based on inertia technology | |
CN108759845A (en) | A kind of optimization method based on inexpensive multi-sensor combined navigation | |
CN103196445B (en) | Based on the carrier posture measuring method of the earth magnetism supplementary inertial of matching technique | |
CN108036785A (en) | A kind of aircraft position and orientation estimation method based on direct method and inertial navigation fusion | |
CN106482734A (en) | A kind of filtering method for IMU Fusion | |
CN109682377A (en) | A kind of Attitude estimation method based on the decline of dynamic step length gradient | |
CN109709628B (en) | Calibration method for gravity gradiometer of rotating accelerometer | |
CN106017452B (en) | Double tops disturbance rejection north finding method | |
CN103076026B (en) | A kind of method determining Doppler log range rate error in SINS | |
CN106441372B (en) | A kind of quiet pedestal coarse alignment method based on polarization with gravitation information | |
CN107255474A (en) | A kind of PDR course angles of fusion electronic compass and gyroscope determine method | |
CN106403952A (en) | Method for measuring combined attitudes of Satcom on the move with low cost | |
CN115574816B (en) | Bionic vision multi-source information intelligent perception unmanned platform | |
CN108444468B (en) | Directional compass integrating downward vision and inertial navigation information | |
CN105735969A (en) | Oil well bore track plotting device and method | |
CN116448145A (en) | Navigation attitude determination method based on polarization vector space difference | |
CN113418499B (en) | Method and system for resolving roll angle of rotary aircraft | |
CN113175926B (en) | Self-adaptive horizontal attitude measurement method based on motion state monitoring | |
CN114675055A (en) | Air speed measurement error compensation method for sonde based on inertial system | |
Zhe et al. | Adaptive complementary filtering algorithm for imu based on mems |
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 |