• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Drag derived altitude aided navigation method

    2017-11-17 08:31:59HuSONGMeilingWANGJingYANGDntongSUNLongtoKONGChunlingWEI
    CHINESE JOURNAL OF AERONAUTICS 2017年5期

    Hu SONG,Meiling WANG,Jing YANG,*,Dntong SUN,Longto KONG,Chunling WEI

    aSchool of Automation Science and Electrical Engineering,Beihang University,Beijing 100083,China

    bCollaborative Innovation Center of Advanced Aero-Engine,Beihang University,Beijing 100083,China

    cNational Key Laboratory of Space Intelligent Control Technology,Beijing 100094,China

    Drag derived altitude aided navigation method

    Hua SONGa,b,Meiling WANGa,Jing YANGa,*,Dantong SUNa,Longtao KONGa,Chunling WEIc

    aSchool of Automation Science and Electrical Engineering,Beihang University,Beijing 100083,China

    bCollaborative Innovation Center of Advanced Aero-Engine,Beihang University,Beijing 100083,China

    cNational Key Laboratory of Space Intelligent Control Technology,Beijing 100094,China

    The navigation problem of the lifting reentry vehicles has attracted much research interest in the past decade.This paper researches the navigation in the blackout zone during the reentry phase of the aircraft,when the communication signals are attenuated and even interrupted by the blackout zone.However,when calculating altitude,a pure classic inertial navigation algorithm appears imprecise and divergent.In order to obtain a more precise aircraft altitude,this paper applies an integrated navigation method based on inertial navigation algorithms,which uses drag derived altitude to aid the inertial navigation during the blackout zone.This method can overcome the shortcomings of the inertial navigation system and improve the navigation accuracy.To furthe r improve the navigation accuracy,the applicable condition and the main error factors,such as the atmospheric coefficient error and drag coefficient error are analyzed in detail.Then the damping circuit design of the navigation control system and the damping coefficients determination is introduced.The feasibility of the method is verified by the typical reentry trajectory simulation,and the influence of the iterative times on the accuracy is analyzed.Simulation results show that iterative three times achieves the best effect.

    1.Introduction

    Owing to the characteristics of the rocket,spacecraft,reentry vehicle and aircraft,the lifting reentry vehicle has been widely used in the field of space technology,and the navigation problem of the reentry vehicle has been getting more and more attention.1,2The reentry flight is characterized as high speed,large attack angle,large aerodynamic interference,and the existence of the phenomenon of blackout zone.3,4According to the characteristics of the reentry flight and autonomous navigation,the main reentry vehicle navigation is Inertial Navigation System(INS)/Global Positioning System(GPS)/Celestial Navigation System(CNS)integrated navigation at present.5,6

    When the reentry vehicle returns to the atmosphere at a high speed,the ground communication will be interrupted in a certain altitude area,and the area is called blackout zone.7,8The blackout zone occurred between 35 km and 80 km over the Earth’s atmosphere,and when the reentry vehicle is in the blackout zone,the temperature of the vehicle hulls will reach 2000°C,hence the GPS,radio navigation and communication will be interrupted,and the n the reentry vehicle losses external contact(high temperature ionizes the air around the aircraft into plasma,and the n shields the electromagnetic wave).Thus the groundings cannot obtain the information of the shuttle’s real-time condition.9–12At this point,the navigation of reentry vehicle can only rely on inertial navigation,a kind ofindependent navigation technology that does not rely on external signals.But due to the defect of INS,relying on pure inertial navigation willlead to large navigation errors.So when the vehicle returns,the blackout phenomenon makes real time communication and measurement rathe r hard.

    In view of the problem of blackout phenomenon,Ref.8introduced several methods,such as changing the shape design of reentry vehicle,adding an external magnetic field to weaken the influence of the blackout zone.But with the se methods,the aircraft structure or flight environment needs changing,the implementation is difficult and of high cost.13Some researches use only INS for the aircraft navigation during the blackout zone,14,15which employs inertial navigation system directly to calculate the altitude of the aircraft and proves to have a large error.In Refs.5,6,the GPS’s prediction information is used for the vehicle navigation during the blackout zone,but the uncertainty of the prediction information will affect the navigation accuracy.So in order to ensure the navigation accuracy when the aircraft is out of the blackout zone,it is necessary to use auxiliary method to restrain the divergence of INS.The Refs.16,17proposed a combined navigation method which uses Drag Derived Altitude(DDA)aided inertial navigation system to the reentry vehicle navigation in the blackout zone. The nongravitational acceleration of the reentry vehicle,and the aerodynamic and atmospheric models are used to estimate the altitude of the craft in DDA.And using DDA is able to get a higher accuracy than a pure classic inertial navigation algorithm.But the main error factors of DDA and the damping circuit design of the integrated navigation control system has not been introduced in detail in the two articles.

    This paper mainly discusses the navigation method during the blackout zone.Firstly it introduces the principle of DDA measurement and the solution of the parameters used,the n a two order damping loop is established to combine the altitude of DDA with the inertial navigation information.Next,this paper analyses the applicable condition,the main error factors and the effect of the iterative times in detail.And at last,the integrated navigation method is applied to the navigation of the reentry vehicle during the blackout zone and the feasibility of this method is verified by simulation results.

    2.Drag derived altitude

    2.1.Definition of coordinate systems

    (1)Earth-Centered Inertial(ECI)coordinate system

    The inertial coordinate system of Earth defines a coordinate which is fixed with the rotating Earth and whose origin Oiis at Earth center.The principal direction xiis through the equatorial at RAAN(Right Ascension of the plane Ascending Node of the orbit),ziis through the Earth’s rotation axis to north and yicompletes the right-handed system.

    (2)Earth-Centered Earth-Fixed(ECEF)coordinate system

    A coordinate system whose origin Oeis at the center of the Earth and axis xeis through the equatorial at Greenwich,yeis through the equatorial and is perpendicular to xe,and zecompletes the right-handed system.

    (3)Body coordinate system

    A coordinate system whose origin Obis at the aircraft center and axes move with the aircraft.The xbis typically aligned with body spin axis and points to the head of the aircraft,and the ybis perpendicular to xband points to the top of the aircraft,and zbcompletes the right-handed system.

    (4)Navigation coordinate system

    Navigation coordinate system(Onxnynzn)is the same as the Earth-centered inertial coordinate system in this paper.

    (5)Speed coordinate system

    A coordinate system whose origin Obis at the aircraft center and whose axes move with the aircraft.The xvis along the direction of the speed,the yvis in the longitudinal symmetry plane of the aircraft and perpendicular to xvand points to the top of the aircraft,and the zvcompletes the right-handed system.

    The transformation between different coordinate systems are usually defined as the form of C21,which means the transformation matrix from 1 coordinate system to 2 coordinate system.And the transformation between different coordinate systems are given as follows:

    (A)Navigation coordinate system and body coordinate system

    The transformation matrix can be acquired by rotating the navigation coordinate system three times to the body coordinate system z(θ)→ y(ψ)→ x(γ).The transformation matrix is where θ,ψ and γ are respectively pitch angle,heading angle and roll angle of the aircraft.

    (B)Speed coordinate system and body coordinate system

    According to the definition,the Obyvaxis of the speed coordinate system is in longitudinal symmetry plane of the aircraft,that is the xbObybplane of the body coordinate system.Therefore,the conversion between the two coordinate systems only exists two Euler angles,and the body coordinate system can be acquired by the body coordinate system rotating two times:y(β)→ z(α),and the transformation matrix is

    where α is the angle of the attack,the angle between the velocity vector projection in the main symmetry plane of the aircraft and the longitudinal axis of the body,which selects the downward direction as positive direction.β is the sideslip angle,which is the angle between the velocity vector and the main symmetry plane of the aircraft,which chooses the right direction of the main symmetry plane as positive direction.

    2.2.DDA principle

    DDA is a pseudo-measurement of the altitude obtained from aerodynamic model,atmospheric model and the estimated non-gravitational acceleration.When flying without other external forces,the aircraft is subjected to the gravity and aerodynamic force.The aerodynamic force can be decomposed into drag and lift.And the force analysis of the aircraft is shown in Fig.1,where P is the push along the longitudinal axis of the aircraft and pointing to the front,G is the Earth gravity,pointing to the center of the Earth center,D is the drag along the opposite direction of the movement,L is the lift,perpendicular the motion direction and pointing upward.Mris the rolling moment,and Vris the aircraft velocity.

    By the the ory of aerodynamics,the air drag is

    where ρ is the local atmospheric density,Srefis the atmospheric reentry aircraft referential area,and CDis the drag coefficient.

    According to Newton’s second law:

    Fig.1 Force analysis of aircraft.

    where m is the aircraft mass,Γ is the acceleration produced by air drag,which can be obtained by accelerometers.

    According to Eqs.(3)and(4),the non-gravitational acceleration along the velocity in the opposite direction is:

    Assume that the atmospheric density evolution is locally exponential,the n the atmospheric density model is given as:

    where h is the altitude,H(h)is the local atmospheric coefficient,and ρ0is the sea level atmospheric density.According to Eqs.(5)and(6),the ‘pseudo-measurement” of the altitude h can be derived by

    2.3.Altitude calculation method using DDA combined with INS

    It is clear that the local atmospheric coefficient H(h)and the drag coefficient CDare the key to solving Eq.(7),and other parameters can be obtained from aircraft parameters and accelerometers.

    (1)Atmospheric coefficient H(h)calculation

    Velocity and position information provided by INS can help to calculate the local atmospheric coefficient,and provide the relative velocity.In order to obtain the local atmospheric coefficient,it is necessary to convert the position information Xp(t)of INS output to the altitude h.

    Then the local atmospheric coefficient H(h)can be acquired by using the estimated altitude h to look up the atmosphere model table.

    (2)Relative velocity Vrcalculation

    Relative velocity,the modulus of the aircraft velocity,can be calculated according to the three axis velocity in navigation coordinate,

    where Vn= [Vnx,Vny,Vnz]is the triaxial velocity component of aircraft obtained by INS in the navigation frame.

    (3)Drag coefficient CDcalculation

    Drag coefficient is determined by the velocity and angle of attack.This paper uses Mach number Ma=Vr/Vsoundto represent speed,where Vsound=340 m/s.

    The blackout zone occurred between 35 km and 80 km over the Earth’s atmosphere,and this altitude ranges from stratosphere to mesosphere.In stratosphere,wind speed is generally tens of meters per second,18and this velocity is very small relative to reentry vehicle velocity.Therefore,wind speed can be ignored when calculating the angle of attack during the blackout zone.

    When the re is no wind,the ground velocity is equal to the relative air velocity,the n the angle of attack and sideslip angle can be calculated by using the ground velocity.Converting the velocity vector Vninto the body system, we get Vb=CbnVn= [Vbx,Vby,Vbz]T.Thus α and β can be assessed by

    Then,using Mach number and attack angle we get drag coefficient according to the aerodynamic model.Since it is hard to establish the aerodynamic model analytic equation,the drag coefficient is obtained by looking-up table method in this paper.

    (4)Non-gravitational acceleration calculation

    According to the force analysis in Fig.1,the output of aircraft accelerometer is determined by the resultant force except gravity.Since the re is only air drag along the velocity direction except gravity,the acceleration Γ produced by the air drag is the projection of resultant acceleration in the opposite direction of velocity.Converting the accelerometer output from the body system to the speed system Obxvyvzv,the component along the opposite direction of x axis is Γ.

    According to the above analysis,a ‘pseudo-measureme nt”of the altitude is available,and the process is shown in Fig.2,whereis the angular velocity of the b coordinate system relative to the i coordinate system projected in the b coordinate system,is the specific force of the b coordinate system relative to the i coordinate system projected in the b coordinate system,is the specific force of the v coordinate system relative to the i coordinate system projected in the v coordinate system,Xp(t)and Vn(t)are the position and velocity information calculated by the INS respectively.

    Fig.2 Calculation process of DDA.

    3.Referential height aided INS navigation method

    3.1.INS divergence mechanism

    In the inertial coordinate system,the three channels of INS are divergent.So without auxiliary information,the altitude calculated by INS would have a great error during the blackout zone.The reasons are analyzed as follows.

    3.1.1.INS solution process

    In the inertial system,the inertial navigation system is mainly divided into three steps:attitude update,position update and speed update.

    (1)Attitude update

    The update equation of attitude matrix is

    (2)Speed update

    The velocity differential equation in the inertial system is

    wherefbis the specific force measured by accelerometer;gnis gravitational acceleration in the inertial coordinate system,and its calculation formula is

    where rx,ry,rzare the location coordinates of aircraft in inertialsystem,and ris the geocentric distance,thatisμ is the product of the Earth’s mass and gravitational constant,Reis the radius of the Earth’s equator,and J2is the perturbation term.Different referential ellipsoid models use different Earth parameters,and the constants used in this scheme are μ =3.986035 × 1014m3/s2,Re=6378140 m and J2=1.08230×10-3.

    (3)Position update

    The position differential equation in the inertial coordinate system is

    where rn= [rx,ry,rz]Tis the position vector of aircraft in inertial coordinate system.

    3.1.2.Error propagation equation of INS

    The error propagation equation is derived based on the INS mechanics.

    (1)Propagation equation of attitude misalignment angle

    The error angle between the ma the matical platform^n simulated by INS and the ideal navigation coordinate system n is defined as the attitude misalignment angle, that is φ = [φx,φy,φz]T,where φ is small.Then n can be obtained by^n through the following equation:

    where [φ×]is anti-symmetric matrix of φ,that is

    Bring Eqs.(10)–(16),the n

    According to Eq.(10),the partial differential equation of

    Compare Eqs.(17)and(18),it is obvious that

    According to [(AB)×]=A[B×]AT,where A is an orthogonal matrix,B is a column vector.Then

    This is the propagation equation of attitude misalignment angle,whereis the error of angular velocity.It can be seen that the platform drift of inertial coordinate system is caused by the measurement error of the gyroscope totally.

    (2)Propagation equation of velocity error

    According to Eq.(11),the partial differential equation of

    Eq.(22)is the speed error propagation equation,where δVn= [δVnx,δVny,δVnz]T, [fn×] isanti-symmetric matrix formed by the projection component of accelerometer output in navigation system,that is fn=Cnbfb.δfbis the measurement error of acceleration,δgnis the projection of gravitational acceleration error in the navigation coordinate system caused by the position error and can be obtained by the differential calculation of gn.

    (3)Propagation equation of position error

    The propagation characteristics of position error is relatively simple,with the expression as

    where δrn= [δrx,δry,δrz]T.

    According to the error propagation equation,the reason why the three channels are divergent is explained as follows.From Eq.(7)when the disturbance item J2is neglected,the calculation formula of gravitational acceleration can be simplified as follows:

    So the differential equation of gravitational acceleration is changed to

    where

    Therefore,assume the position error and velocity error in the inertial coordinate system as the state variable,the nAccording to Eqs.(22),(23)and(25),X can be explained as

    where the second item at the right side of the formula is irrelevant to position error and velocity error.The system characteristic equation is

    Thus, the characteristic roots are

    It is clear that the system contains Schuler oscillation period components,and because the characteristic roots contain positive real part,the position errors of three channels in INS are exponentially divergent in the inertial navigation system.

    3.2.Damping loop design and damping coefficient determination

    The navigation error of INS is divergent during the process of reentry section without other information supporting.Therefore,it is necessary to introduce auxiliary information,such as altitude calculated from DDA.Because the damping is realized in the inertial coordinate system at three channels respectively,three dimensional position vectors used in the damping loop is calculated by using the referential height which can be obtained by DDA,and the longitude and latitude obtained by INS.

    The damping loop is designed based on the control the ory.In the original inertial navigation loop,the damping loop is added to change the stability of the system to achieve the purpose of restraining error accumulation and improving navigation accuracy.In this paper,the second-order damping circuit with high accuracy and less damping parameters is designed.The structure is shown in Fig.3.

    The two integral parts of the forward channel are the pure INS and the damping coefficient of the three channels respec-

    within which K1x,K2xare the damping coefficients of x axis,K1y,K2yare the damping coefficients of y axis and K1z,K2zare the damping coefficients of z axis.hrefis auxiliary referential height,and the part in the dotted borders converts the onedimensional height information(^h),longitude(^λ)and the latitude(^L)into a three-dimensional position vector rn,ref.The auxiliary height hrefin the blackout zone is provided by DDA.AccordingtoFig.3,the system differential equations of position and velocity with the second-order damping are given as follows:

    According to Eqs.(22),(28)and(29),after employing the second-order damping loop,the position and velocity error propagation equations of the INS are as follows

    where δrn,refis the position error vector derived by the auxiliary height information,and δrn,refis related to the INS error and the auxiliary height.The propagation equation of the attitude misalignment angle will not change when the damping is added.

    According to Eq.(24),assume that the aircraft altitude h?Re,the n r≈Re,Eq.(24)can be further simplified as

    Then the differential calculation formula of gravitational acceleration is

    bringing Eq.(33)to Eq.(31),the n

    The system characteristic equation is

    The damping coefficient is solved by pole placement method in this paper,where six identical negative poles are assigned,and can ensure the stability of the system.The poles of the damping system are assigned as:

    where τ is time constant reflecting the regulation time and overshoot of the system.From Eqs.(36)and(37),the damping coefficient K1and K2can be calculated.In order to simplify the design of the parameters,assume the damping coefficients are the same,that is K1x=K1y=K1z=K1,K2x=K2y=K2z=K2.K1and K2can be calculated according to control principle,and

    Fig.3 Structure diagram of INS with a second-order damping circuit.

    4.Analysis of DDA simulation

    4.1.Analysis and processing of original data

    The atmospheric model and air drag coefficient model are given in Tables A1 and A2(see Appendix A).The corresponding parameters are obtained by looking up tables.The table of atmospheric model is given by the table of h- ρ(h),while h-H(h)is required by Eq.(7)to obtain h.According to Eq.(6),one can get H(h)=-ln(ρ/ρ0)/h,and the table h- ρ(h)can be transferred to the form of h-H(h).In simulation,according to the altitude calculated by INS,the corresponding atmospheric coefficient is obtained by checking h-H(h)table and carrying on a linear interpolation.

    The drag coefficient is associated with the Mach number and the angle of attack. A two-dimensional table(Ma,α)-CDis obtained by experiments,the n CDcan be obtained by the two-dimensional table looking up and the bilinear interpolation according to the speed and attack angles.If the velocity or the angle of attack is beyond the scope of the table,the boundary value is used to replace the velocity or the angle of attack.

    During the reentry process,the aircraft is only affected by the Earth gravity and the air force,so Γ is the projection of the acceleration vector in the velocity reverse direction and Vris the speed of the aircraft relative to the air,that is,the value of the vehicle speed.The other variables of DDA are the parameters of the aircraft itself.

    4.2.Simulation of DDA under ideal input

    Using the ideal height and acceleration as input,the altitude error of DDA during the whole reentry process is calculated.Then the altitude error is shown in Fig.4(as the altitude of the reentry vehicle decreases monotonically during the whole reentry process,the aircraft height is selected as the abscissa.In this paper,the analysis of DDA is simulated in the blackout zone,and the altitude is 35–90 km).

    Fig.4 Altitude error of DDA.

    It is shown that the accuracy of DDA is the best when the real height is between 20 and 90 km,because at lower altitudes the wind speed has a relatively large effect on the reentry aircraft speed,and at higher altitudes,no significant acceleration produced by residence can be established,so the accuracy of DDA is relatively low.

    In addition,it is clear in Fig.4 that the altitude calculated by DDA in 35–90 km is lower than the actual height,because the re is a constant deviation when using interpolation algorithm to calculate.To furthe r improve the accuracy,the more precise atmospheric model and aerodynamic coefficient model is needed.

    The change regularity of the air density is very complex,and several main factors affecting the density change are:height,the shape of the Earth,the seasons and so on.19,20The reentry environment is also very complex,the drag coefficient of reentry aircraft is affected by angle of attack,Mach number,atmospheric density and rudder angle and so on.Thus,for more accurate atmospheric density models,it is necessary to fully considerate the factors that affect the atmospheric density and aerodynamic coefficient to establish mathematical models.

    4.3.Error analysis of DDA

    The basic principleformula of DDA is given as Eq.(7),the accuracy of the aircraft velocity Vr,the atmospheric coefficient H,the acceleration Γ,the drag coefficient CDare the main factors that affect the DDA accuracy.Assume that ΔH,ΔVr,ΔCDand ΔΓ are the maximum errors on H,Vr,CDand Γ,and the se errors may be caused by the reasons shown in Table 1.

    Among the above reasons,the error analysis of the inertial navigation errors,the INS accuracy can be unified as the error analysis ofinertial navigation.The atmospheric model and the drag coefficient model are based on the method of looking up tables which has certain uncertainty.The wind and the Earth model mainly affect the drag coefficient and the air coefficient,the y are attributed to the model error in this paper.

    4.3.1.Error analysis of INS

    The errors considered in simulation are as follows:

    (1)Installation deviation of INS is 90′.

    (2)Initial position error is 200 m,initial velocity error is 0.5 m/s,and initial pose error is 720′.

    (3)The nonorthogonal installation deviation of gyro is 90′,and the constant drift is 0.1°/h.

    (4)The nonorthogonal installation deviation of accelerometer is 90′,and the constant drift is 10-4g,where g is the gravity acceleration.

    Table 1 Errors causes.

    (5)Standard error of gyro scalefactor is 1 × 10-4(3σ),standard deviation of random noise is 1.4142°/h,where σ is the standard deviation.

    (6)Standard error of accelerometer scale factor is 2 × 10-4(3σ),standard deviation of random noise is 10-4g(3σ).

    In this paper,for the convenience of expression,we use INS installation error,gyro error,and accelerometer error to express the errors considered in simulation.The specific expression is shown in Table 2.

    The error of DDA considering different conditions is shown in Fig.5.Fig.5(a)is the DDA error when the INS installation error and initial altitude error are considered;Fig.5(b)is the DDA error when the INS installation error,initial attitude error and gyro error are considered;Fig.5(c)is the DDA error when the INS installation error,initial attitude error and accelerometer error are considered,and Fig.5(d)is the DDA error when the INS installation error,initial attitude error,gyro error and accelerometer error are considered.

    Fig.6 is the comparison chart of the DDA height errors in the following cases:no error,the INS error only,INS error and gyro error,INS error and accelerometer error,and including the se three kinds of errors.Fig.7 is the comparison chart of the DDA height errors adding different size errors to the gyro.And Fig.8 is the comparison chart of the DDA height errors adding different size errors to the accelerometers.In Figs.7 and 8,the errors added are shown in Table 3.

    It can be shown from Fig.5 that the influence of the INS installation error and initial attitude error is not very large.Fig.6 shows the effect of gyros error and accelerometers error on the DDA error;Fig.7 shows that the effect of gyro error on the DDA error is about 200 m;Fig.8 shows that the effect of the accelerometer error on the DDA error is about 200 m.

    4.3.2.Error analysis of atmospheric coefficient error ΔH

    The atmospheric coefficient will directly affect the accuracy of DDA.The atmospheric model is influenced by altitude,Earth rotation,wind,shape of the Earth and seasons.It is difficult to establish a precise atmospheric model to determine the atmospheric density.This paper uses the method of adding devia-tion to the true value of atmospheric coefficient,and the deviation range is 2%.11–13

    Table 2 Error description.

    Fig.5 Errors of DDA considering different conditions.

    Fig.6 Effect comparison of different type errors on DDA error.

    Fig.7 Effect comparison of different size gyro errors on DDA error.

    Fig.8 Effect comparison of different size accelerometer errors on DDA error.

    Table 3 Errors added in Figs.7 and 8.

    In this paper,the situations that consider random atmospheric coefficient error and the maximum atmospheric coefficient error are analyzed.And the DDA errors are shown in Figs.9 and 10.Fig.9 is the DDA error comparison when different atmospheric coefficient errors are considered,where the big error is 1.02H,the medium error is 1.002H,and the small error is 0.98H.In simulation,the height calculated from INS is replaced by the true height.When the random atmospheric coefficient error is considered to the simulation,the mean value of DDA error is-271.4403 m,and the standard deviation is 702.5892 m during the 35–90 km reentry process.Figs.9 and 10 show that the effect of atmospheric coefficient error on the DDA error is relatively large,and for higher accuracy,it is needed to improve the accuracy of the atmospheric parameters model.

    Fig.9 Effect comparison of different atmospheric coefficient errors on DDA error.

    Fig.10 Effect of random atmospheric coefficient error on DDA error.

    4.3.3.Error analysis of drag coefficient error ΔCD

    Drag coefficient is obtained by experiments under different parameters,it is related to the Mach number,angle of attack and the sideslip angle.The angle of attack and sideslip angle are dependent on INS,and the INS error analysis has been completed in the previous section.Mach number is related to sound velocity,so it is dependent on temperature,atmospheric density and so on.As well as atmospheric coefficient error,the drag coefficient error range in this paper is 10%.

    In this paper,the situations that consider random error and the maximum error of drag coefficient are analyzed.And the DDA errors are shown in Figs.11 and 12.Fig.11 is the DDA error comparison when different drag coefficient errors are considered,where the big error,medium error and small errorare 1.1ΔCD,1.02ΔCDand0.9ΔCDrespectively.The height calculated from inertial navigation system is replaced by the true height.When the random drag coefficient error is considered,the mean value of DDA error is-407.0702 m,and the standard deviation is 445.3268 m during the 35–90 km reentry process.From Figs.11 and 12,it is clear that the effect of drag factor error on the DDA error is relatively large,and for higher accuracy,it is needed to improve the accuracy of the drag parameters model.

    Fig.11 Effect comparison of different drag coefficient errors on DDA error.

    Fig.12 Effect of random drag coefficient error on DDA error.

    4.3.4.Error analysis of comprehensive error

    (1)Consider the drag coefficient error and atmospheric coefficient error,the simulations are shown in Figs.13 and 14.

    (2)Consider INS error,drag coefficient error and atmospheric coefficient error,the simulations are shown in Figs.15 and 16.

    Fig.13 shows the great influence of drag coefficient error and atmospheric coefficient error on the accuracy of DDA.Fig.14 reflects the influence of random drag coefficient error and atmospheric coefficient error on the accuracy of DDA.And the DDA error range is-3000 to 2000 m.In Figs.13 and 15,the big error,medium error and small error are 1.1ΔCD,1.02ΔCDand 0.9ΔCDrespectively.Comparing Figs.13 and 15,it is clear that the INS error has little influence on the accuracy of DDA.Therefore,the main factors affecting the accuracy of DDA are drag coefficient error and atmospheric coefficient error.Fig.16 reflects the effect of INS error,random drag coefficient error and atmospheric coefficient error on the accuracy of DDA.

    Fig.13 Effect comparison of different drag coefficient error and atmospheric coefficient error on DDA error.

    Fig.14 Effect of random drag coefficient error and atmospheric coefficient error on DDA error.

    Fig.15 Effect comparison of INS error,different drag coefficient error and atmospheric coefficient error on DDA error.

    Fig.16 Effect of INS error,random drag coefficient error and atmospheric coefficient error on DDA error.

    4.4.Influence ofiterative algorithm

    The DDA accuracy can be improved by carrying out several DDA measurements as described in Fig.17.As is shown in Fig.17,the altitudes h0,h1,...are used to calculate the atmospheric coefficient,and during the first iteration,the altitude h0is obtained by INS,and during the rest iterations,the altitude h1,h2are obtained by the last DDA.At the last iteration,the output of DDA is the aircraft altitude.

    Fig.17 Method of successive iteration.

    Fig.18 Error comparison of different iterations.

    When considering the maximum error,iterative twice,three times and four times respectively,the simulation is shown in Fig.18(a).When using the atmospheric coefficient table with different precision,iterative twice,three times and four times respectively,the simulation is shown in Fig.18(b)and(c).The accuracy of INS can also affect the convergence of the iteration,and the iteration result after reducing the precision of INS is shown in Fig.18(d)and(e).

    Fig.18(a)shows that the accuracy of three times is better than two,but the accuracy offour times is almost the same as three times.Thus,it is not the more iterations,the higher accuracy.Considering both the computation burden and the accuracy,the optimal number is three times.Fig.18(b)and(c)are the error comparisons of different iterationswhen the data ofatmosphericcoefficienttable reduced by 1/20 and 1/3 respectively.Fig.18(b)and(c)show that if the accuracy of the atmospheric coefficient is reduced to a certain extent,the iterative algorithm will not converge.Fig.18(d)and(e)show that when the precision of INS is reduced,the convergenceofthe iteration willalso be reduced.So in order to ensure the convergence of the iterative algorithm,we must ensure the accuracy of atmospheric coefficient table we used and the accuracy of height calculated by INS.

    5.Simulation analysis of whole reentry trajectory using auxiliary navigation

    The whole reentry process of the aircraft is simulated,the errors of the devices and INS are the same as in Section 4.3.1.Two simulations are carried out,one uses pure INS and the other uses DDA.The time constant τ=80,and the errorcurves of the two casesare shown in Fig.19,where Δsx,Δsy,Δszare the position error of aircraft in x,y,z axis respectively,and Δvx,Δvy,Δvzare the velocity error of aircraft in x,y,z axis respectively.The black curve is the precision of the pure inertial navigation,and the blue curve is the precision of INS+DDA.It is notable that some simulations have been done to find the best τ value.When τ ranges from 70 to 90,the DDA has good navigation result.In Fig.19,the value of τ is set to 80.

    Fig.19(a)and(b)show that the position error curve and velocity error curve of the two cases are very similar;it appears that the re is little influence of DDA on the position error and velocity error.As can be seen in Fig.19(c),when using pure inertial navigation,the average of the altitude error is 912.7 m.While the altitude error is-310.9 m when using DDA,which reveals that the altitude error of the DDA is smaller than that of INS in blackout zone,and it has a certain delay and the trend is more gentle,which is to say that the altitude calculated from DDA has higher accuracy than INS.

    During the reentry process,the blackout zone lasts about 5 min.As can be seen in Fig.19(c),when using pure inertial navigation,the altitude error is about 1 km in 5 min.So in order to ensure the navigation accuracy when the aircraft is out of the blackout zone,it is necessary to use auxiliary method like DDA to restrain the divergence of INS.

    Fig.19 Influence of DDA on navigation accuracy.

    What’s more,after the blackout zone,the altitude error of non-damping is gradually divergent and close to the INS,and in this area,the GPS could be used for the aircraft navigation.

    6.Conclusion

    In this paper,the INS+DDA integrated navigation method has been used in the navigation of the reentry aircraft during the blackout zone.Firstly,the principle of DDA is introduced,and the detailed process is given.Then through the analysis of the wind velocity and the drag factors of air,the applicable condition of DDA is determined as 20–90 km.What’s more,based on the simulation analysis of the influence ofinertial navigation error,atmospheric coefficient and drag coefficient error on DDA,it is concluded that the atmospheric coefficient error and drag coefficient error are the main factors that affect the accuracy of DDA.Iterative times have great influence on the precision of DDA.Simulation analysis shows that the optimal number of the iterations is 3 with respect to both the computation burden and accuracy.Ultimately,the second order damping circuit of altitude is designed,and the simulation analysis is carried out on the whole reentry trajectory.The simulation results show that the altitude accuracy of the solution can be reached by-310 m,which is higher than that of the pure INS.

    Acknowledgements

    This study was supported by the National Natural Science Foundation of China(No.61573059).

    Appendix A

    See Tables A1 and A2.

    Table A1 Atmospheric coefficient.

    Table A2 Drag coefficient.

    1.Zheng FQ.Study on reentry trajectory design and guidance approach for lifting spacecraft[dissertation].Mianyang:China Aerodynamics Research and Development Center;2013[Chinese].

    2.Zhao J,Zhou R.Reentry trajectory optimization for hypersonic vehicle satisfying complex constraints.Chin J Aeronaut 2013;26(6):1544–53.

    3.Bobylev AV,Dyad’kin AA,Kobzev VI,Poedinok VM,Reshetin SN,Suprunenko SN,et al.Problems of control of a re-entry vehicle with a moderate lift-to-drag ratio during its entry into the atmosphere.Cosmic Res 2008;46(1):74–89.

    4.Geng J,Sheng Y,Liu X.Finite-time sliding mode attitude control for a reentry vehicle with blended aerodynamic surfaces and a reaction control system.Chin J Aeronaut 2014;27(4):964–76.

    5.Yang F,Cheng C,Zhang GY,Cheng YM.Integrated navigation method for a sub-orbital vehicle.J Astronaut 2010;31(3):729–33[Chinese].

    6.Li J,Yang B.Study on the application of integrated navigation technology in the reusable launch vehicle.J Shenyang Inst Aeronaut Eng 2007;24(4):9–12[Chinese].

    7.Yang J,Kong LT,Zhang Z.Research on SINS in application of reentry navigation for lifting reentry vehicles.Guidance,navigation and control conference;2014 Aug 9–10;Yantai,China.Piscataway:IEEE Press;2014.

    8.Li J,Yang B.The pivotal technology of re-entry navigation for a reusable launch vehicle[dissertation].Beijing:Beihang University;2007[Chinese].

    9.Wang X,Xia Y.Navigation strategy with the spacecraft communications blackout for Mars entry.Adv Space Res 2015;55(4):1264–77.

    10.Hu HJ,Chen Y,Chen J.The discussion on TT&C technology for spacecraft in blackout area.J Proj,Rockets,Missiles Guid 2012;32(2):197–9,206[Chinese].

    11.Qu X,Fang GP.Introduction and analysis of ‘blackout zone”.Silicon Valley 2010;1(10):149–73[Chinese].

    12.Hu J,Zhang Z.A study on the reentry guidance for a manned lunar return vehicle.Control Theory Appl 2014;31(12):1674–85[Chinese].

    13.Li JP,Lv N,Zhang C.Study on approaches for mitigating hypersonic vehicle communications blackout.Fire Command Control 2012;37(2):155–8[Chinese].

    14.Tian B,Fan W,Zong Q.Integrated guidance and control for reusable launch vehicle in reentry phase.Nonlinear Dynam 2015;80(1–2):397–412.

    15.Viviani A,Pezzella G.Aerodynamic and aerothermodynamic analysis of space mission vehicles.1st ed.Berlin:Springer;2015.p.457–570.

    16.Delaux P,Bouaziz L.Navigation algorithm of the atmospheric re-entry demonstrator. Guidance, navigation, and control conference;1996 Jul 29–31;San Diego,USA.Reston:AIAA;1996.

    17.Pignie G,Delaux P,Bouaziz L,Kebci H,Trinquard F.Parachutes opening triggering algorithm ofthe atmospheric re-entry demonstrator.Reston:AIAA;1996.Report No.:AIAA-1996-3755.

    18.Xu AL,Zhong AH,Sun JH,Yang GR.Vertical structure characteristic from troposphere to low stratosphere in Dali region over southe astern margin of Qinghai-Xizang Plateau.Plateau Meteorol 2016;35(1):77–85[Chinese].

    19.Wang WS.Atmospheric density model in the high-precision autonomous orbit determination of satellites.J XingTai Univ 2008;23(2):89–91[Chinese].

    20.Bai ZX,Bian JC,Chen HB,Chen L.Inertial gravity wave parameters for the lower stratosphere from radiosonde data over China.Sci China Earth Sci 2017;60(2):328–40.

    1 July 2016;revised 15 January 2017;accepted 18 April 2017

    Available online 24 August 2017

    Aided navigation system;

    Blackout zone;

    Drag derived altitude;

    Inertial navigation;

    Integrated navigation;

    Reentry vehicles

    *Corresponding author.

    E-mail address:jing.yang@buaa.edu.cn(J.YANG).

    Peer review under responsibility of Editorial Committee of CJA.

    Production and hosting by Elsevier

    http://dx.doi.org/10.1016/j.cja.2017.08.003

    1000-9361?2017 Production and hosting by Elsevier Ltd.on behalf of Chinese Society of Aeronautics and Astronautics.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

    ?2017 Production and hosting by Elsevier Ltd.on behalf of Chinese Society of Aeronautics and Astronautics.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

    精品久久久久久久末码| 国产欧美日韩精品一区二区| 亚洲第一区二区三区不卡| 午夜日本视频在线| 最近视频中文字幕2019在线8| 日韩欧美精品免费久久| 国产老妇伦熟女老妇高清| 插阴视频在线观看视频| 一级毛片我不卡| 免费搜索国产男女视频| 亚洲婷婷狠狠爱综合网| 午夜免费激情av| 亚洲精品国产av成人精品| 一区二区三区乱码不卡18| 六月丁香七月| 日韩亚洲欧美综合| 三级毛片av免费| av女优亚洲男人天堂| 国产老妇伦熟女老妇高清| 能在线免费观看的黄片| 麻豆av噜噜一区二区三区| 麻豆av噜噜一区二区三区| 国产一级毛片在线| 男女那种视频在线观看| 欧美激情久久久久久爽电影| 亚洲av成人精品一二三区| 一区二区三区四区激情视频| 午夜精品国产一区二区电影 | 91在线精品国自产拍蜜月| 非洲黑人性xxxx精品又粗又长| 啦啦啦观看免费观看视频高清| 九九在线视频观看精品| 亚洲精品色激情综合| 最后的刺客免费高清国语| 99热这里只有精品一区| 在线观看66精品国产| 精品国产三级普通话版| 能在线免费看毛片的网站| 国产精品久久电影中文字幕| 夫妻性生交免费视频一级片| 国产精品国产三级国产av玫瑰| 性插视频无遮挡在线免费观看| 欧美成人免费av一区二区三区| 国产伦在线观看视频一区| 国产亚洲午夜精品一区二区久久 | 成人漫画全彩无遮挡| 非洲黑人性xxxx精品又粗又长| a级一级毛片免费在线观看| 国产熟女欧美一区二区| 国产精品一区二区三区四区久久| 亚洲欧美一区二区三区国产| kizo精华| 日本午夜av视频| 一个人免费在线观看电影| 建设人人有责人人尽责人人享有的 | 精品不卡国产一区二区三区| 最近最新中文字幕免费大全7| 成人无遮挡网站| 国产精品久久电影中文字幕| 欧美日韩一区二区视频在线观看视频在线 | 成人美女网站在线观看视频| 日韩一本色道免费dvd| 变态另类丝袜制服| 久久精品国产亚洲网站| www.av在线官网国产| 亚洲精品成人久久久久久| 在线观看美女被高潮喷水网站| 久久精品久久久久久久性| 国产亚洲精品av在线| 人人妻人人澡欧美一区二区| 日韩av在线大香蕉| 老师上课跳d突然被开到最大视频| 成人亚洲欧美一区二区av| 久久婷婷人人爽人人干人人爱| 久久韩国三级中文字幕| 又黄又爽又刺激的免费视频.| 你懂的网址亚洲精品在线观看 | 久久久久久久久中文| 婷婷色麻豆天堂久久 | 岛国毛片在线播放| 高清日韩中文字幕在线| 久99久视频精品免费| av专区在线播放| 久久热精品热| 一个人免费在线观看电影| 91精品国产九色| 欧美人与善性xxx| 精华霜和精华液先用哪个| 两个人视频免费观看高清| 欧美成人免费av一区二区三区| 国产精品国产高清国产av| 亚洲欧美成人综合另类久久久 | 91午夜精品亚洲一区二区三区| 欧美性猛交黑人性爽| 免费黄色在线免费观看| 少妇熟女aⅴ在线视频| 婷婷色综合大香蕉| 欧美zozozo另类| 亚洲国产精品久久男人天堂| 永久网站在线| 亚洲最大成人av| 99热6这里只有精品| 日韩在线高清观看一区二区三区| 久久精品久久久久久噜噜老黄 | 天天躁夜夜躁狠狠久久av| 久久久久久国产a免费观看| 久久久久性生活片| 看十八女毛片水多多多| 亚洲国产欧洲综合997久久,| 韩国高清视频一区二区三区| 老司机影院毛片| 亚洲自偷自拍三级| 国产精品国产三级国产专区5o | 久久久久九九精品影院| 赤兔流量卡办理| 国产精品福利在线免费观看| 能在线免费看毛片的网站| 久久久精品大字幕| 国产午夜精品久久久久久一区二区三区| 综合色丁香网| 国产高清三级在线| 精品人妻偷拍中文字幕| 免费观看性生交大片5| 国产91av在线免费观看| 大话2 男鬼变身卡| or卡值多少钱| 99久久成人亚洲精品观看| 天堂中文最新版在线下载 | 又粗又爽又猛毛片免费看| 美女xxoo啪啪120秒动态图| 欧美极品一区二区三区四区| 三级国产精品欧美在线观看| 午夜亚洲福利在线播放| 国产大屁股一区二区在线视频| 亚洲中文字幕一区二区三区有码在线看| 一边亲一边摸免费视频| 亚洲欧美精品综合久久99| 国产精品一区二区在线观看99 | 一级毛片久久久久久久久女| videos熟女内射| 国产免费一级a男人的天堂| 国产亚洲最大av| 晚上一个人看的免费电影| 色网站视频免费| 亚洲国产欧洲综合997久久,| 啦啦啦啦在线视频资源| 九九在线视频观看精品| 久久精品综合一区二区三区| 亚洲精品乱码久久久v下载方式| 高清午夜精品一区二区三区| 亚洲精品国产成人久久av| 国产精品蜜桃在线观看| 国产精品久久久久久久电影| 成年女人永久免费观看视频| 亚洲av不卡在线观看| 精品国产一区二区三区久久久樱花 | av女优亚洲男人天堂| 亚洲无线观看免费| 亚洲精品乱久久久久久| 婷婷色av中文字幕| 国产一级毛片七仙女欲春2| 综合色av麻豆| 一卡2卡三卡四卡精品乱码亚洲| 九九在线视频观看精品| 一夜夜www| 尤物成人国产欧美一区二区三区| 精品人妻偷拍中文字幕| av视频在线观看入口| 男女那种视频在线观看| 汤姆久久久久久久影院中文字幕 | 国产人妻一区二区三区在| 欧美+日韩+精品| 亚洲综合色惰| 极品教师在线视频| 日本免费一区二区三区高清不卡| 久久综合国产亚洲精品| 精品国内亚洲2022精品成人| 国产不卡一卡二| 久久国内精品自在自线图片| 久久精品国产亚洲av涩爱| 久久久久久久久大av| 久久国内精品自在自线图片| 大香蕉久久网| 国产高清国产精品国产三级 | 日韩欧美三级三区| 国产v大片淫在线免费观看| 高清午夜精品一区二区三区| 岛国在线免费视频观看| 波多野结衣巨乳人妻| 国产一区二区在线观看日韩| 最近2019中文字幕mv第一页| 在线观看66精品国产| 欧美性猛交黑人性爽| 亚洲精品乱码久久久v下载方式| 男的添女的下面高潮视频| 久久6这里有精品| 日韩视频在线欧美| 男人舔女人下体高潮全视频| 中文字幕精品亚洲无线码一区| 精品酒店卫生间| 五月玫瑰六月丁香| 欧美成人午夜免费资源| av在线观看视频网站免费| 成人无遮挡网站| 亚洲成av人片在线播放无| 亚洲欧美成人综合另类久久久 | 26uuu在线亚洲综合色| 男人和女人高潮做爰伦理| 国产成人一区二区在线| 国产成人精品婷婷| 日本三级黄在线观看| 亚洲精品乱码久久久久久按摩| 丝袜美腿在线中文| 色5月婷婷丁香| 国产精品一及| 一卡2卡三卡四卡精品乱码亚洲| 国产精品爽爽va在线观看网站| 日本一二三区视频观看| 伊人久久精品亚洲午夜| 久久精品综合一区二区三区| 国产综合懂色| 少妇丰满av| 麻豆成人午夜福利视频| 亚洲经典国产精华液单| 听说在线观看完整版免费高清| av线在线观看网站| 国产av在哪里看| 亚洲精品国产av成人精品| av在线亚洲专区| 日韩三级伦理在线观看| 免费黄网站久久成人精品| 久久久久久九九精品二区国产| 亚洲av熟女| 亚洲丝袜综合中文字幕| 久久久久久伊人网av| 精品无人区乱码1区二区| 色网站视频免费| 久久综合国产亚洲精品| 天堂av国产一区二区熟女人妻| 亚洲欧美成人综合另类久久久 | 国产一区二区三区av在线| 一本久久精品| 久久久久久久久久久免费av| 欧美性猛交╳xxx乱大交人| 国产又黄又爽又无遮挡在线| 国产一区二区在线观看日韩| 成年免费大片在线观看| 午夜福利高清视频| 亚洲欧美清纯卡通| 精品酒店卫生间| 久久久久久大精品| 日韩视频在线欧美| 色吧在线观看| 成人亚洲精品av一区二区| 男的添女的下面高潮视频| 卡戴珊不雅视频在线播放| 免费看日本二区| 成人亚洲欧美一区二区av| 国产成人福利小说| 桃色一区二区三区在线观看| 亚洲第一区二区三区不卡| av视频在线观看入口| 麻豆一二三区av精品| 久久婷婷人人爽人人干人人爱| 韩国av在线不卡| 欧美日韩一区二区视频在线观看视频在线 | 在线观看66精品国产| 欧美性猛交╳xxx乱大交人| 日本一本二区三区精品| 久久这里只有精品中国| 晚上一个人看的免费电影| 精品人妻偷拍中文字幕| 1024手机看黄色片| 熟女电影av网| 毛片一级片免费看久久久久| 亚洲欧美日韩高清专用| www.av在线官网国产| 国产高清国产精品国产三级 | 少妇裸体淫交视频免费看高清| 啦啦啦韩国在线观看视频| 国产精品久久久久久精品电影| 男人狂女人下面高潮的视频| 久久99精品国语久久久| 国产中年淑女户外野战色| 亚洲电影在线观看av| 毛片一级片免费看久久久久| 女人久久www免费人成看片 | 97超碰精品成人国产| 五月伊人婷婷丁香| 亚洲国产欧洲综合997久久,| 欧美日韩一区二区视频在线观看视频在线 | 99热网站在线观看| 精品一区二区三区视频在线| kizo精华| 久久精品国产自在天天线| 春色校园在线视频观看| 日韩在线高清观看一区二区三区| 亚洲av免费在线观看| 国产精品伦人一区二区| 免费无遮挡裸体视频| 高清av免费在线| 免费黄网站久久成人精品| 日韩一本色道免费dvd| 亚洲高清免费不卡视频| 国产激情偷乱视频一区二区| 1000部很黄的大片| 国产精品精品国产色婷婷| 欧美3d第一页| 成人高潮视频无遮挡免费网站| 99热网站在线观看| 国产精品美女特级片免费视频播放器| 伦精品一区二区三区| 国产精品一及| 51国产日韩欧美| 免费看光身美女| 精品一区二区免费观看| 水蜜桃什么品种好| 亚洲最大成人中文| 精品午夜福利在线看| 国产乱来视频区| av福利片在线观看| 亚洲国产欧美人成| 成年av动漫网址| 哪个播放器可以免费观看大片| 青春草国产在线视频| 男女啪啪激烈高潮av片| 三级毛片av免费| 欧美一区二区精品小视频在线| 春色校园在线视频观看| 午夜福利在线观看免费完整高清在| 能在线免费看毛片的网站| 国语自产精品视频在线第100页| 一夜夜www| 久久久久久久久久成人| 最近中文字幕高清免费大全6| 中文字幕久久专区| 免费观看a级毛片全部| 亚洲四区av| 亚洲无线观看免费| 久久精品久久久久久噜噜老黄 | 国产免费一级a男人的天堂| 狂野欧美激情性xxxx在线观看| 亚洲久久久久久中文字幕| 久久精品91蜜桃| 久久久久久大精品| 亚洲美女搞黄在线观看| 纵有疾风起免费观看全集完整版 | 亚州av有码| 欧美丝袜亚洲另类| 一级黄色大片毛片| 淫秽高清视频在线观看| 日本三级黄在线观看| 久久久亚洲精品成人影院| 卡戴珊不雅视频在线播放| 国产精品野战在线观看| 国产精品av视频在线免费观看| 亚洲欧美一区二区三区国产| 国产乱人视频| 乱系列少妇在线播放| 欧美精品一区二区大全| 亚洲欧美一区二区三区国产| 国产av在哪里看| 国产精品国产三级国产专区5o | 国产在线一区二区三区精 | 久久精品久久久久久久性| 人人妻人人澡欧美一区二区| 久久草成人影院| 国产一区二区三区av在线| 婷婷色av中文字幕| 最新中文字幕久久久久| 亚洲精品一区蜜桃| 最近最新中文字幕免费大全7| 老女人水多毛片| 久久久久久久久大av| videos熟女内射| 国产精品,欧美在线| 亚洲精品一区蜜桃| 亚洲自偷自拍三级| 看片在线看免费视频| 丰满乱子伦码专区| 麻豆成人av视频| 日韩成人av中文字幕在线观看| 国产探花极品一区二区| 成人午夜高清在线视频| av福利片在线观看| 99久国产av精品国产电影| 国产免费福利视频在线观看| 亚洲欧美一区二区三区国产| 国产亚洲5aaaaa淫片| 一级黄片播放器| 国产av在哪里看| 超碰97精品在线观看| av国产免费在线观看| 国产精品久久久久久久久免| 国产精品av视频在线免费观看| 国产av一区在线观看免费| 高清日韩中文字幕在线| 99热这里只有是精品50| 国产一区二区在线av高清观看| a级一级毛片免费在线观看| 亚洲最大成人av| 亚洲性久久影院| 国产黄片视频在线免费观看| 欧美人与善性xxx| 两个人视频免费观看高清| 久久久久久久久久久免费av| 免费观看的影片在线观看| 成人毛片60女人毛片免费| 欧美xxxx性猛交bbbb| 国产精品国产三级专区第一集| 美女高潮的动态| 日韩中字成人| 中文在线观看免费www的网站| 国产一区亚洲一区在线观看| 免费一级毛片在线播放高清视频| 国产精品女同一区二区软件| 九草在线视频观看| 亚洲欧美清纯卡通| 久久精品人妻少妇| 国产精品久久久久久久久免| 美女国产视频在线观看| 国产精品麻豆人妻色哟哟久久 | 精品人妻熟女av久视频| 亚洲欧美精品专区久久| 久久精品国产亚洲av天美| 亚洲综合色惰| 欧美潮喷喷水| 青青草视频在线视频观看| 国产毛片a区久久久久| 好男人视频免费观看在线| 国产午夜福利久久久久久| 欧美日韩精品成人综合77777| 欧美成人精品欧美一级黄| 婷婷六月久久综合丁香| 成人综合一区亚洲| 舔av片在线| 一个人看视频在线观看www免费| 亚洲欧美日韩卡通动漫| 在线天堂最新版资源| 久久久色成人| 欧美日本亚洲视频在线播放| 国产成人aa在线观看| 亚洲精品aⅴ在线观看| 一区二区三区乱码不卡18| 亚洲成色77777| 99久久无色码亚洲精品果冻| 伦理电影大哥的女人| 国产在视频线精品| 日韩av在线免费看完整版不卡| 超碰97精品在线观看| 婷婷六月久久综合丁香| 99热精品在线国产| 天美传媒精品一区二区| 欧美最新免费一区二区三区| 成人av在线播放网站| 亚洲真实伦在线观看| 一卡2卡三卡四卡精品乱码亚洲| 亚洲av.av天堂| 日韩一区二区视频免费看| 一卡2卡三卡四卡精品乱码亚洲| 国产精品久久久久久精品电影| 欧美97在线视频| 最近视频中文字幕2019在线8| 黄色配什么色好看| 男人和女人高潮做爰伦理| 午夜福利在线在线| 久久精品国产亚洲av天美| 欧美精品国产亚洲| 国产精品不卡视频一区二区| 日日摸夜夜添夜夜爱| 中文亚洲av片在线观看爽| 在线免费十八禁| av在线观看视频网站免费| 欧美高清性xxxxhd video| 国产欧美另类精品又又久久亚洲欧美| 久久亚洲国产成人精品v| 一二三四中文在线观看免费高清| 精品一区二区三区人妻视频| 欧美3d第一页| 最近最新中文字幕免费大全7| 午夜视频国产福利| 99国产精品一区二区蜜桃av| 日本黄大片高清| 韩国av在线不卡| 亚洲av成人精品一二三区| 夜夜爽夜夜爽视频| 国产精品女同一区二区软件| 最近中文字幕2019免费版| 能在线免费观看的黄片| 色噜噜av男人的天堂激情| 天堂网av新在线| 一级二级三级毛片免费看| 国语自产精品视频在线第100页| 精品一区二区三区视频在线| 免费观看a级毛片全部| 国产高清三级在线| 国产精品女同一区二区软件| 全区人妻精品视频| 建设人人有责人人尽责人人享有的 | 亚洲av熟女| 永久网站在线| 国产精品伦人一区二区| 97超视频在线观看视频| 在线播放无遮挡| 久久午夜福利片| 日韩人妻高清精品专区| 九九爱精品视频在线观看| 天堂√8在线中文| 特级一级黄色大片| 99久久精品国产国产毛片| 国产av在哪里看| av黄色大香蕉| 中国国产av一级| 色网站视频免费| 久久人人爽人人片av| 久久婷婷人人爽人人干人人爱| 国产激情偷乱视频一区二区| 91精品伊人久久大香线蕉| 日韩欧美 国产精品| 一级av片app| 亚洲自偷自拍三级| 亚洲最大成人中文| 亚洲国产色片| 免费播放大片免费观看视频在线观看 | 麻豆久久精品国产亚洲av| 国产精品伦人一区二区| 欧美精品国产亚洲| www.av在线官网国产| 精品一区二区三区视频在线| 国产一区二区在线观看日韩| 亚洲丝袜综合中文字幕| 国产成人a∨麻豆精品| 久久热精品热| 国产在线一区二区三区精 | 啦啦啦观看免费观看视频高清| 久久99热这里只频精品6学生 | 亚洲高清免费不卡视频| 波多野结衣高清无吗| 直男gayav资源| 成人特级av手机在线观看| 老司机影院毛片| 秋霞在线观看毛片| 最近最新中文字幕免费大全7| 欧美日韩国产亚洲二区| 毛片女人毛片| 3wmmmm亚洲av在线观看| 国产激情偷乱视频一区二区| 久久午夜福利片| 国产av在哪里看| 国产色爽女视频免费观看| 天美传媒精品一区二区| 边亲边吃奶的免费视频| 亚洲精品乱码久久久v下载方式| 男女国产视频网站| 在线观看66精品国产| 日韩一区二区三区影片| 我的女老师完整版在线观看| 熟妇人妻久久中文字幕3abv| 男人和女人高潮做爰伦理| eeuss影院久久| 国产伦精品一区二区三区四那| 六月丁香七月| 免费看a级黄色片| 久久人妻av系列| 日韩成人av中文字幕在线观看| 人妻夜夜爽99麻豆av| 成年av动漫网址| 狂野欧美激情性xxxx在线观看| 精品国产露脸久久av麻豆 | 乱码一卡2卡4卡精品| 欧美一区二区亚洲| 永久网站在线| 99热这里只有精品一区| 国产高清国产精品国产三级 | 蜜臀久久99精品久久宅男| 国产av一区在线观看免费| 乱人视频在线观看| 九九在线视频观看精品| 超碰97精品在线观看| 久久人人爽人人爽人人片va| 大香蕉久久网| 99热6这里只有精品| 高清日韩中文字幕在线| 国产免费男女视频| 亚洲欧美中文字幕日韩二区| 亚洲怡红院男人天堂| 亚洲国产日韩欧美精品在线观看| 午夜免费激情av| 97人妻精品一区二区三区麻豆| 久久久精品大字幕| 91久久精品国产一区二区成人| 久久精品人妻少妇| 老女人水多毛片| 69人妻影院| 中文天堂在线官网| 尤物成人国产欧美一区二区三区| 超碰av人人做人人爽久久| 我的老师免费观看完整版| 欧美另类亚洲清纯唯美| or卡值多少钱| 亚洲欧美精品自产自拍| 精品99又大又爽又粗少妇毛片| 久久久久久国产a免费观看| 小说图片视频综合网站| 久久精品国产99精品国产亚洲性色| 欧美精品国产亚洲| 一个人看的www免费观看视频| 99热这里只有是精品在线观看| 亚洲18禁久久av| 白带黄色成豆腐渣| 建设人人有责人人尽责人人享有的 | 1024手机看黄色片| 国产伦精品一区二区三区视频9| 欧美高清成人免费视频www| 18禁在线播放成人免费|