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 PDF

Info

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
Application number
CN202210317866.XA
Other languages
Chinese (zh)
Other versions
CN114675055B (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)
  • 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

Air speed measurement error compensation method for sonde based on inertial system
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:
Figure BDA0003569427230000041
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:
Figure BDA0003569427230000042
continuous condition f by boundary derivativei′(ti+1)＝f′i+1(ti+1) Obtaining:
Figure BDA0003569427230000043
from fi″(ti+1)＝f″i+1(ti+1) Obtaining:
Figure BDA0003569427230000044
let mi＝fi″(ti)＝2ciThen formula (2) is written as mi+6hidi＝mi+1And then:
Figure BDA0003569427230000045
hi,mirespectively intermediate variables;
s3-2-3: a is toi＝ti，
Figure BDA0003569427230000046
Substitution in formula (1) gives:
Figure BDA0003569427230000047
a is toi，bi，ci，diSubstitution in formula (2) gives:
Figure BDA0003569427230000051
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:
Figure BDA0003569427230000052
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；
S3-3: setting a new sampling time interval
Figure BDA0003569427230000053
The interpolated angle θ is then:
θ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
Figure BDA0003569427230000061
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:
Figure BDA0003569427230000062
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:
Figure BDA0003569427230000063
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 b
Figure BDA0003569427230000064
Expressed as:
Figure BDA0003569427230000065
θxis shown as
Figure BDA0003569427230000066
And OXbThe included angle between the two parts:
Figure BDA0003569427230000067
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 b
Figure BDA0003569427230000068
Expressed as:
Figure BDA0003569427230000069
the vector OH is
Figure BDA00035694272300000610
In OYnZnProjection on a plane, n is expressed as:
OHn＝[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
Figure BDA0003569427230000071
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:
Figure BDA0003569427230000072
by forward difference, get tkEast speed of time sonde relative to umbrella
Figure BDA0003569427230000073
Speed in north direction
Figure BDA0003569427230000074
Figure BDA0003569427230000075
wherein ,
Figure BDA0003569427230000076
is tkThe east velocity calculated from the processed three-axis angular velocities at that moment,
Figure BDA0003569427230000077
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 system
Figure BDA0003569427230000078
East speed of sonde relative to geographic coordinate system
Figure BDA0003569427230000079
The wind speed obtained from the GPS/beidou information of the sonde is:
Figure BDA00035694272300000710
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:
Figure BDA00035694272300000711
wherein ,
Figure BDA00035694272300000712
the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,
Figure BDA00035694272300000713
the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to time
Figure BDA00035694272300000714
And
Figure BDA00035694272300000715
the speed of the corresponding moment can be found
Figure BDA00035694272300000716
And
Figure BDA00035694272300000717
the wind speed at this time is:
Figure BDA00035694272300000718
to VwProcessing to obtain the final compensated wind speed Vwind：
Figure BDA0003569427230000081
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:
Figure BDA0003569427230000101
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:
Figure BDA0003569427230000111
continuous condition f by boundary derivativei′(ti+1)＝f′i+1(ti+1) Obtaining:
Figure BDA0003569427230000112
from f to fi″(ti+1)＝f″i+1(ti+1) Obtaining:
Figure BDA0003569427230000113
let m bei＝fi″(ti)＝2ciThen formula (2) is written as mi+6hidi＝mi+1And then:
Figure BDA0003569427230000114
hi,mirespectively intermediate variables;
s3-2-3: a is toi＝ti，
Figure BDA0003569427230000115
Substitution in formula (1) gives:
Figure BDA0003569427230000116
a is toi，bi，ci，diSubstitution in formula (2) gives:
Figure BDA0003569427230000117
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:
Figure BDA0003569427230000118
Figure BDA0003569427230000121
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；
S3-3: setting a new sampling time interval
Figure BDA0003569427230000122
The interpolated angle θ is then:
θ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
Figure BDA0003569427230000123
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:
Figure BDA0003569427230000124
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:
Figure BDA0003569427230000131
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 b
Figure BDA0003569427230000132
Expressed as:
Figure BDA0003569427230000133
as shown in fig. 3, θxIs shown as
Figure BDA00035694272300001312
And OXbThe included angle between:
Figure BDA0003569427230000134
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 b
Figure BDA0003569427230000135
Expressed as:
Figure BDA0003569427230000136
the vector OH is
Figure BDA0003569427230000137
In OYnZnProjection on a plane, n is expressed as:
OHn＝[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
Figure BDA0003569427230000138
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:
Figure BDA0003569427230000139
by forward difference, get tkEast speed of time sonde relative to umbrella
Figure BDA00035694272300001310
Speed in north direction
Figure BDA00035694272300001311
Figure BDA0003569427230000141
wherein ,
Figure BDA0003569427230000142
is tkThe east velocity calculated from the processed three-axis angular velocities at that moment,
Figure BDA0003569427230000143
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 system
Figure BDA0003569427230000144
East speed of sonde relative to geographic coordinate system
Figure BDA0003569427230000145
The wind speed obtained from the GPS/beidou information of the sonde is:
Figure BDA0003569427230000146
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:
Figure BDA0003569427230000147
wherein ,
Figure BDA0003569427230000148
the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,
Figure BDA0003569427230000149
the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to time
Figure BDA00035694272300001410
And
Figure BDA00035694272300001411
the speed of the corresponding moment can be found
Figure BDA00035694272300001412
And with
Figure BDA00035694272300001413
The wind speed at this time is:
Figure BDA00035694272300001414
to VwProcessing to obtain the final compensated wind speed Vwind：
Figure BDA00035694272300001415
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:
Figure FDA0003569427220000021
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:
Figure FDA0003569427220000022
condition f 'is continued by the boundary derivative'i(ti+1)＝f′iv1(ti+1) Obtaining:
Figure FDA0003569427220000023
from f ″)i(ti+1)＝f″i+1(ti+1) Obtaining:
Figure FDA0003569427220000024
let m bei＝f″i(ti)＝2ciThen formula (2) is written as mi+6hidi＝mi+1And then:
Figure FDA0003569427220000025
hi,mirespectively intermediate variables;
s3-2-3: a is toi＝ti，
Figure FDA0003569427220000031
Substitution in formula (1) gives:
Figure FDA0003569427220000032
a is toi，bi，ci，diSubstitution in formula (2) gives:
Figure FDA0003569427220000033
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:
Figure FDA0003569427220000034
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；
S3-3: is provided withSetting a new sampling time interval
Figure FDA0003569427220000035
The interpolated angle θ is then:
θ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
Figure FDA0003569427220000041
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:
Figure FDA0003569427220000042
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:
Figure FDA0003569427220000043
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 b
Figure FDA0003569427220000044
Expressed as:
Figure FDA0003569427220000045
θxis shown as
Figure FDA0003569427220000046
And OXbThe included angle between:
Figure FDA0003569427220000047
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 b
Figure FDA0003569427220000048
Expressed as:
Figure FDA0003569427220000049
the vector OH is
Figure FDA00035694272200000410
In OYnZnProjection on a plane, n is expressed as:
OHn＝[0 -cosγcosψ sinγ]T
according to a geometric relationship, θyIs represented by OYnAngle to OH:
Figure FDA0003569427220000051
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:
Figure FDA0003569427220000052
by forward difference, get tkEast speed of time sonde relative to umbrella
Figure FDA0003569427220000053
Speed in north direction
Figure FDA0003569427220000054
Figure FDA0003569427220000055
wherein ,
Figure FDA0003569427220000056
is tkThe east velocity calculated from the processed three-axis angular velocities at that moment,
Figure FDA0003569427220000057
is tkAnd calculating the north velocity from the processed three-axis angular velocities.
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 system
Figure FDA0003569427220000058
East speed of sonde relative to geographic coordinate system
Figure FDA0003569427220000059
The wind speed obtained from the GPS/beidou information of the sonde is:
Figure FDA00035694272200000510
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:
Figure FDA00035694272200000511
wherein ,
Figure FDA00035694272200000512
the speed of the sonde relative to the parachute, i.e. the n-system, is measured for the inertial system,
Figure FDA00035694272200000513
the wind speed is in an umbrella lower point coordinate system, namely an i system;
tivelocity corresponding to time
Figure FDA00035694272200000514
And
Figure FDA00035694272200000515
the speed of the corresponding moment can be found
Figure FDA00035694272200000516
And with
Figure FDA00035694272200000517
The wind speed at this time is:
Figure FDA00035694272200000518
to VwProcessing to obtain the final compensated wind speed Vwind：
Figure FDA0003569427220000061
Wherein floor (. cndot.) is rounded downward.
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 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)

* Cited by examiner, † Cited by third party
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

Patent Citations (12)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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