Normal gravity model for inertial navigation of a hypersonic boost-glide vehicle

Kai CHEN,Cheng-zhi ZENG,Sen-sen PEI,Wen-chao LIANG

School of Astronautics,Northwestern Polytechnical University,Xi’an 710072,China

Abstract:The normal gravity model of a hypersonic boost-glide vehicle in near space is studied in this paper with the aim of alleviating the influence of the gravity model error on the precision of the inertial navigation system(INS)during flight.First,a spherical harmonic model of the Earth’s gravitational field is introduced and the normal gravity of the Earth is derived from it.Then,the coordinate transformation needed for the application of the gravity model to the near-space navigation algorithm is formulated.Subsequently,the gravity disturbance in near space and the impact of J2 and J4 gravity truncation errors are analyzed.Finally,different normal gravity models and different precisions of inertial measurement unit (IMU) are exploited to simulate the near-space navigation algorithm.Based on this,the influence of the independent and combined effects caused by the interference factors is analyzed,and the applicable conditions of the normal gravity model are discussed.

Key words:Hypersonic boost-glide vehicle;Inertial navigation system(INS);Normal gravity;Gravity disturbance

1 Introduction

The speed of a hypersonic vehicle is usually above Mach 5,which gives such vehicles great value in commercial applications and in scientific research.Traditional subsonic aircraft need 10 to 20 h for a long flight for which,in theory,hypersonic vehicles need only 2 to 3 h (Bykerk et al.,2020).An inertial navigation system can provide the attitude,velocity,and position of the aircraft at a high sampling frequency.Thus,it is an essential part of the navigation of nearspace aircraft(Chen et al.,2020a).The inertial naviga‐tion system and global navigation satellite system(INS/GNSS) integrated navigation system is commonly used on near-space vehicles (Chen et al.,2020b).For example,NASA’s X-43A adopts Honeywell’s H-764 series strap-down inertial navigation system and global positioning system (SINS/GPS) integrated navigation system(Bahm et al.,2005),SHEFEX-2 uses the SINS/GPS integrated navigation system(Steffes,2013),and the HTV-2 aircraft uses the SINS/GPS tightly-coupled integrated navigation system(Walker et al.,2008).

Improving the precision of the INS is of great significance for the performance of the integrated nav‐igation system.In the INS,gravity is an inertial force which cannot be measured by inertial measurement units (IMUs),and should be calculated by gravity models (Chatfield et al.,1975).The normal gravity model based on the spheroid assumption is conve‐nient to use and can meet the requirements of most engineering applications (Hsu,1996).Hence,it is cur‐rently the most widely used gravity model.However,with the improvement of the measurement precision of IMUs,the model error of the normal gravity model gradually becomes one of the main sources of error in inertial navigations.

With the development of geodesy technologies,high-precision gravity model databases,such as EGM2008 and EIGEN-6C4,were established after the launch of satellite gravimetry missions,including gravity recovery and climate experiment (GRACE) and gravity field and steady-state ocean circulation explorer(GOCE)(Tsoulis and Patlakis,2014;Chen et al.,2018;Slobbe et al.,2019).These models can model the gravi‐tational field of the Earth with quite high precision(Pavlis et al.,2012;Foerste et al.,2014).However,although the high-order spherical harmonic models offer global high precision,they suffer from a large computation load (Guy and Keith,1998;Claessens and Hirt,2015).Younis and Becker(2013)proposed a spherical cap harmonic expansion method and estab‐lished a local model to simplify the computation.How‐ever,this model fails to meet the long-distance naviga‐tion requirements of some hypersonic vehicles.

Since the precision improvement of the gravity models is always accompanied by higher complexity and more stringent performance requirements of the navigation computer,it is essential to investigate the applicable conditions of the gravity models with dif‐ferent precisions.Zhu et al.(2020)researched a method of gravity disturbance compensation to improve the attitude measurement accuracy,but the navigation speeds and position errors were not discussed.Chang et al.(2019) studied a compensation algorithm based on the INS error model,which selectively compensated for the gravity error according to whether the acceler‐ometer’s constant bias and gravity offset each other.These studies conducted only ground verification without full consideration of the altitude channels.Besides,inertial navigation has the characteristic of altitude channel divergence.In addition,the gravity calculation error caused by the altitude error makes the altitude error more serious,forming a positive feed‐back.Near-space hypersonic vehicles use the launchcentered inertial frame navigation algorithm,where the divergence in the altitude direction is projected to the three channels of the launch-centered inertial frame.Currently,the influence of gravity model errors on the launch-centered inertial frame navigation has not been explored in depth.

This paper describes research on a hypersonic boost-glide vehicle in near space.The EIGEN-6C4 highorder spherical harmonic model is directly applied to generate a flight trajectory with the dynamic model.Meanwhile,different gravity models are used in navi‐gation.The measurement precision of different IMUs is exploited in navigation simulation.The navigation error changes under the coupling effect of the gravity model error and the IMU error.These are analyzed to formulate the usage conditions of the normal gravity model.

The rest of this paper is organized as follows.Section 2 introduces the mathematical foundation of the Earth’s spherical harmonic gravity model.Section 3 describes the high-precision strap-down inertial navigation algorithm and the conversion of the gravity model from the respective reference system to the navigation reference system.Section 4 simulates the global distribution of gravitational disturbance in near space,and the simulation results are analyzed.Section 5 presents the simulated flight trajectory,the gravita‐tional disturbance in the flight area,and the naviga‐tion simulation results under different gravity models and IMU precisions.

2 Earth’s gravitational field model

2.1 Spherical harmonic model of the Earth’s gravitational field

The spherical harmonic model of the Earth’s gravitational field is the gravity model used by the EGM2008 and EIGEN-6C4 databases.These databases adopt a variety of measurement methods,such as altimetry satellites and gravity satellites.The spherical harmonic model built by long-term measurement can achieve a better precision.The spherical harmonic model adopts the spherical coordinate system,and the expression of the gravitational potential at any point(ρ,φ,λ) outside the Earth is expressed as follows(Pavlis et al.,2012):

whereφis the geocentric latitude;λis the longitude;ρis the radius of the geocentric vector;Gis the gravita‐tional constant;Mis the mass of the Earth;Reis the long radius of the Earth;the coefficientsCn,kandSn,kare provided by databases,such as EIGEN-6C4;Pn,k(sinφ)is the associated Legendre polynomial.

The solution of the gravitational vector relies on the gravitational potential function to determine the gradient.The gradient expression in the spherical coordinate system is

wherei,j,andkare the unit vectors along the tangent directions of the three coordinate lines of the spherical coordinate system.

After the partial derivative of the gravitational potentialVis obtained,the expression of the gravita‐tional vector is represented as follows:

The three components of the vector,i.e.gE,gU,andgN,are projected onto the three directions of the tangent of the meridian,the tangent of the latitude,and the normal on the sphere,respectively.

2.2 Normal gravity model in near space

Normal gravity is the gravitational field excited by simplifying the real Earth to a reference ellipsoid.Due to its similarity to the real gravity field,this model is commonly used in navigation,and the basic pa‐rameters of the reference ellipsoid can be provided by Earth models such as WGS84.The reference ellipsoid can be described by a clear and concise mathematical formula,which makes the normal gravity model suit‐able for engineering applications.

Eq.(1) can be exploited to describe the external gravitational potential of any mass body.It can be simplified by the application of the properties of the reference ellipsoid.The reference ellipsoid has a uni‐form mass distribution,and its normal gravity does not change with the longitudeλ.The coefficient of Eq.(1)varying with longitude should be equal to zero.That is,whenk≠0,the coefficientsCn,kandSn,kare both zero;otherwise,Cn,kcos(kλ)+Sn,ksin(kλ) becomes a constant,Jn.Pn,0(sinφ) is a Legendre polynomial and is denoted asPn(sinφ).Based on the above analysis,the gravita‐tional potential function of Eq.(1)can be expressed as

Based on the research of Chachan and Stevenson(2019),

In addition,the reference ellipsoid has the characteristic of north-south symmetry.Ifnis an odd number,sinn-2mφis an odd function with respect toφ.Hence,Pn(sinφ) is an odd function with respect toφ.These terms do not conform to the symmetry of normal gravity about the equator,so the coefficients of these terms in the gravitational field are all equal to zero.

The obtainedJ2ninfinite series expression of normal gravity is as follows:

The coefficientJ2nis derived from the basic parameters of the reference ellipsoid.When theJ2nmodel is exploited to solve the gravitational vector,it is generally converted to an Earth-centered Earthfixed frame.The conversion is expressed as follows:

whereP2n(z/ρ) is the Legendre function.According to Eq.(5),values ofP2n(z/ρ) corresponding ton=1,2,3 are presented as follows:

By substituting the above equations and truncating the higher-order terms,theJ6gravitational potential function can be obtained,and is represented in the following form:

The gradient of the gravitational potential function is the gravitational vector.The gradient in the rectan‐gular coordinate system is:

The partial derivative ofVis actually the partial derivative ofzm/ρnwith different powers in Eq.(9).It can be seen that the results of the partial derivative ofzandx,yfor the term pair are different.The situation wherezm/ρnobtains partial derivatives ofxandyrespectively is as follows:

SubstitutingP2,P4,andP6,the gradient ofV(ρ,φ)can be determined.Then,the expression of theJ6gravitational potential is

If all terms withJ6coefficients in the above equation are ignored,theJ4model is obtained.If all the terms withJ4coefficients are ignored,theJ2model is obtained.TheJ2model is a first-order truncation of theJ2nmodel,and it is a normal gravity model commonly used in the traditional aerospace field for an altitude greater than 200 km.However,in the range of 20–100 km,the influence of higher-order terms in Eq.(13),such asJ4,can increase,and the influence of gravitational disturbance can also increase.Therefore,the use of theJ2model in near space may cause a larger gravity model error.The truncation error of theJ2model and whether the error of the normal gravity model is accepted by the navigation system are issues that need to be further discussed.

3 Results and discussion

This section introduces the algorithm of inertial navigation and the coordinate transformation re‐quired by the use of the gravity model in a navigation system.

3.1 High-precision strap-down inertial navigation algorithm

IMUs are exploited by inertial navigation to measure the angular velocity and specific force of the vehicle,and are integrated to obtain the attitude,position,and velocity of the vehicle.The velocity differential equation for navigation in the launchcentered inertial frame is(Chen et al.,2020a)

wherefais the specific force of the carrier in the launch-centered inertial frame,gais the gravity received by the carrier in the launch-center inertial frame,Vais the velocity of the vehicle,Pais the position of the vehicle,is the attitude matrix,andis the antisymmetric matrix of vehicle’s angular velocitycan be obtained through calculation on the navigation computer.Therefore,the difference betweengacalculated by the gravity model used in navigation and the real gravity results in an error in.Then,errors are introduced to the velocity and position of the carrier.

The velocity increment of gravityis then used to update the velocity and position.The gravity model error is converted into velocity and position errors.

3.2 Conversion of the gravity vector coordinate system

Gravityis solved by the gravity model in‐troduced in Section 2.Meanwhile,the spherical har‐monic function introduced in Section 2 is described in the geocentric spherical coordinate frame,and normal gravity is described in the Earth-centered Earth-fixed frame.Coordinate conversion is required when the gravity models are used in a launch-centered inertial frame.Gravity is transformed to the launch-centered inertial frame (aframe) for description (Chen et al.,2020a).

The three components of gravity obtained by cal‐culating the gradient through Eq.(2),i.e.gE,gN,andgUare respectively along the tangent directions of the three coordinate lines in the spherical coordinate frame.Since the spherical coordinate frame is an orthogonal surface coordinate frame,the three tangents constitute a rectangular coordinate frame,which is denoted as thelframe.Thelframe and theaframe respectively use the geocentric latitude and the geographic lati‐tude,and these two latitudes are fixed to the Earth and to the inertial space,respectively.Furthermore,thex-axis of theaframe points to the target instead of east,which makes the rotation between the two refer‐ence frames complicated.To address this problem,the launch Earth-centered inertial frame (tframe) is used as the transitional reference frame to perform two con‐versions:first from thelframe to thetframe,and then from thetframe to theaframe.

The transformation matrix for conversion from thelframe to thetframe is:

whereωieis the angular velocity of the Earth’s rotation,andtis the time calculated from the moment of launch.

The transformation matrix for conversion from thetframe to theaframe is:

whereλ0is the longitude of the launch point,B0is the geographic latitude of the launch point,andA0is the launch direction.

The rotation matrix for conversion from thelframe to theaframe can be obtained by multiplying the two matrices in Eqs.(16)and(17):

The normal gravity calculated by Eq.(13) is described under the Earth-centered Earth-fixed frame(eframe).Theeframe and thetframe only differ in respect of the Earth’s rotation.The transformation matrix is:

The transformation matrix for conversion from the Earth-centered Earth-fixed frame to the launchcentered inertial frame can be obtained by multiplying the two matrices in Eq.(19):

Following the above computation,the coordinate conversion of the gravity model is completed.The gravity model can be applied in the navigation algorithm.

4 Simulation analysis of the gravity model in near space

4.1 Analysis of Earth’s real gravitational disturbance

To analyze the impact of gravitational distur‐bance,the 2190-order EIGEN-6C4 spherical harmonic model and the WGS-84 ellipsoid normal gravity model are adopted in the experiment.Specifically,these two models output the latitudes of 80°S–80°N and 180°W–179°E.The surface plot of the three components of the Earth’s gravity is shown on the grid of 1°×1°,and the three-component curve of normal gravity is at lati‐tudes of 80°S–80°N.

Through comparison,it is found that the Earth’s gravity and normal gravity at the poles are slightly larger than those at the equator.Also,these two gravity models have the same trends of changes with latitude.This change is caused by the Earth’s oblateness,and indicates that the normal gravity model can reflect the impact of that oblateness.

Some small gravitational disturbances can also be observed on the Earth’s gravity surface.Gravita‐tional disturbance refers to the difference between the Earth’s gravity and the normal gravity.These distur‐bances are caused by uneven distribution of the Earth’s mass and topographical fluctuations.The modeling of these disturbances is very complicated.

An altitude of 20 km is taken in calculating the global gravitational disturbance,and the EIGEN-6C4 gravity are drawn in Fig.2 and the ellipsoidal gravity are drawn in Fig.3.

Fig.1 Numerical updating algorithm

Fig.2 Gravity data at 20-km altitude provided by EIGEN-6C4 model:(a) east component;(b) north component;(c)upward component

Fig.3 Ellipsoidal gravity at 20-km altitude:(a) east component;(b)north component;(c)upward component

The isoline map of the global gravitational disturbance modulus values is drawn in Fig.4,where mgis the unit of measurement,andg=9.806 65 m/s2(mg=10−3g).The three isolines in Fig.4 are 0.02 mg,0.06 mg,and 0.10 mg,indicating that gravitational disturbances are widespread everywhere in the world.The isoline corresponding to 0.06 mgcovers most areas;in few areas is it greater than 0.10 mg.Besides,places with large gravitational disturbances tend to have large terrain undulations,such as the Qinghai-Tibet Plateau near 40°N and 100°E.Additionally,the shape of the isolines in Fig.4 is tortuous,indicating the obvious effect caused by the high-frequency components of the low-altitude gravitational disturbance.

Fig.4 Isoline map of the global gravitational disturbance modulus values at 20-km altitude

Only 16.67%of total area has a gravitational dis‐turbance greater than 0.05 mg,only 9.45%of the area has a disturbance greater than 0.06 mg(Fig.5,area in black),and only 1.68%has a disturbance greater than 0.10 mg(Fig.6).

Fig.5 Area with gravitational disturbance larger than 0.06 mg

Fig.6 Area with gravitational disturbance larger than 0.10 mg

Then,the altitudes of 50 km and 100 km are taken to draw isoline maps shown in Fig.7.

Fig.7 Global gravitational disturbance isoline maps at altitudes of 50 km(a)and 100 km(b)

It can be seen from Fig.5 that the magnitude of gravitational disturbance attenuates with the increase of altitude.Meanwhile,at altitudes of 50 km and 100 km,the areas with gravitational disturbance over 0.10 mgand the areas surrounded by the 0.06-mgisoline are significantly less than those at 20-km altitude.In addi‐tion,the gravitational disturbance isolines are rela‐tively smooth,indicating that the high-frequency com‐ponents are significantly reduced.However,even in the top layer of near space at 100 km,areas with grav‐itational disturbance larger than 0.06 mgdo still exist.

These results indicate that,for an accelerometer with a constant bias of 0.10 mg,the model error of the normal gravity model cannot be ignored in near space.However,for accelerometers with precisions equal or greater than 0.01 mg,the gravity model error may become one of the main sources of error.In this case,the improvement of the IMU precision fails to effec‐tively improve the precision of the inertial navigation.This conclusion will be verified through an inertial navigation algorithm simulation in the next section.

4.2 Comparison of different normal gravity models

Next,the selection basis of the normal gravity model is discussed,and whether the simplification of normal gravity byJ2andJ4models significantly reduces the precision of the normal gravity model is determined.The truncation errors of the gravity model are output,and are illustrated in Fig.8.

Fig.8 J2(a)and J4(b)truncation errors at 20-km altitude.δg:gravitation error

At the altitude of 20 km,theJ2andJ4models are directly compared with the precision normal gravity model.The truncation errors caused by the simplification of theJ2model andJ4model are of the magnitude of 0.01 mgand 10−5mg,respectively.The truncation error of theJ2model has a similar magnitude to that of the gravity disturbance.Therefore,the application of theJ2model may reduce the precision of the normal gravity model to a certain extent.The truncation error of theJ4model is much smaller than the gravity disturbance in near space,and so can be ignored.

5 Navigation simulation verification in near space

The gravity model is verified in this section with the strap-down inertial navigation algorithm of a hy‐personic vehicle in near space.The navigation coordinate system exploits the launch-centered inertial frame.

5.1 Flight trajectory

A boost-glide flight trajectory,which is one of the classical flight paths of hypersonic vehicles,is adopted.Specifically,it firstly boosts up to 67 km altitude,then dives and glides at a high speed around an altitude of 30 km(Chen et al.,2020c).The EIGEN-6C4 model of 2190 orders uses a spherical harmonic model to generate trajectories.The 3D trajectory of the flight trajectory is illustrated in Fig.9.

Fig.9 Simulated trajectory

The area covered by the flight trajectory is roughly located at 103–109°E and 20–35°N.As shown in Fig.10,this area is representative of an area with large gravitational disturbance.The isoline map of gravitational disturbance at 20-km altitude is output,and four isolines corresponding to gravitational distur‐bances of 0.025 mg,0.05 mg,0.075 mg,and 0.10 mgare drawn.The largest gravitational disturbance of this area is between 0.05 and 0.10 mg,which is in the 16.7%of areas that show gravitational disturbances larger than 0.05 mg.

Fig.10 Gravitational disturbance isoline map at 20-km altitude near the simulated trajectory

5.2 Impacts of gravity model errors

Based on the above flight trajectory,the gravity models includingJ2,J4,andJ6,and the EIGEN-6C4 models of maximum ordersN=18 andN=200 are used as the navigation gravity models in the navigation solution.The inertial navigation simulation is performed to obtain the position and velocity error images,which are shown in Fig.11.

Fig.11 Navigation error curves with only gravity model errors:(a)position error;(b)velocity error

The error curves ofJ2andJ4are close,though some difference exists.Meanwhile,the error curves ofJ4andJ6almost coincide.These results indicate that the truncation error of theJ2model used in near space is smaller than the model error,and the addition of theJ4term can improve navigation precision.The influence of theJ6and terms of higher order can be neglected completely,which is consistent with the result of direct comparison.The navigation error of theJ6model is caused by gravity disturbance,and it is critical to model that disturbance to obtain a gravity model with higher precision.Through the two simulations with the spherical harmonic gravity of order 18 and order 200,the navigation errors are reduced while the central processing unit (CPU) time increases.The final velocity and position errors are shown in Table 1.

As the navigation error of INS increases with time,INS is often integrated with GNSS or another navigation system rather than working alone.When the GNSS signal is lost for an INS/GNSS integrated navigation system,the security of the vehicle depends on the INS.For example,if the vehicle is landing when the GNSS signal is lost,the navigation error caused by gravity error will increase to 20 m in about 200 s when using theJ4model and in about 800 s when using the 200-order EIGEN model (only the gravity model error is considered).Considering the width of runways is～20 to 50 m,the vehicle using higher accuracy gravity may get about 800-s landing time in emergency (considering gravity error only)with INS working alone,while theJ4model has about 200 s only.Hypersonic vehicles fly at higher speeds and altitudes than traditional aircraft.More time is to be spent in landing,which justifies the use of the highaccuracy gravity model.

Since theJ2model and theJ4model are still commonly used,it is necessary to investigate their applicable conditions.TheJ4model is used for navigation on the aforementioned flight trajectory,and only accelerometer errors of different levels in Table 2 are added to the model.

The results are shown in Fig.12.Navigation position errors and velocity errors are summarized and listed in Table 3.

Fig.12 Navigation errors caused by the combined effect of normal gravity model errors and accelerometer errors with different levels:(a)position error;(b)velocity error

Table 1 Navigation errors with only the gravity model error

Table 2 Levels of accelerometer errors

Table 3 Navigation position errors and velocity errors

The results indicate that,after the accelerometer constant bias of 100 μgis added,the velocity errorand the position error caused by the combined effect of the accelerometer error and the gravity model are 2.033 m/s and 1019.98 m,respectively.Compared with the independent effect of the gravity model error,the values of these two errors are significantly increased.In this case,the constant accelerometer bias is the main cause of the errors.As the measurement error of the accelerometer decreases,the navigation position error and the velocity error gradually con‐verge to the navigation error curves without IMU errors.Since the accelerometer constant bias and the gravity error have equivalent impacts on the navigation error,adding these two kinds of errors simultaneously may cause the errors to offset each other in a single simulation,thus reducing or increasing the total error.In this case,it can be seen from the simulation result that when the IMU error is less than 10 μg,the gravity model error can be considered as the main cause of navigation error.This result is consistent with the result obtained by direct analysis of gravity disturbance.

Finally,the measurement errors of the gyroscope and accelerometer are added simultaneously to the model,to constitute a more complete IMU error model.In recent years,some high-accuracy IMU products have been used in INS.For example,QA-3000 manu‐factured by Honeywell(Morristown,USA)is a quartz flexible accelerometer with bias repeatability better than 40 μg/h(Zhang et al.,2021).Becka et al.(2008)and Todorokihara et al.(2018) introduced the prototype of an accelerometer with its bias repeatability better than 1 μg/h.The IMU error parameters used for simu‐lation verification are listed in Table 4,and the results are shown in Fig.13.

Fig.13 Navigation errors caused by the combined effect of IMU errors and different gravity models:(a) position error;(b)velocity error

Table 4 IMU error model used in simulation

The results indicate that,for the INS using highprecision IMUs,the error introduced by the normal gravity model cannot be ignored.Under the simula‐tion conditions in this study and considering both the gravity model errors and the IMU errors,the naviga‐tion error is nearly doubled compared to that obtained by considering the IMU errors alone.Hence,under the simulation conditions,the errors of the normal gravity model have a greater impact on the navigation errors.The establishment of a more accurate gravitymodel is essential for improving navigation precision.In addition,this study uses the data from the maximum orders 18 and 200 of EIGEN-6C4 in the simulations.In theory,the 200-order model can reduce the gravity model error to below 1 μg.However,under the given simulation conditions,the navigation error added by the 200-order model to the IMU error model’s error is negligible.

6 Conclusions

Gravitational disturbances ranging from 0.01 mgto 0.10 mgexist in most areas of near space,and the module value of disturbance is large.Meanwhile,the impact of a high-frequency component of the lowaltitude disturbance is obvious.Gravitational disturbance at high altitudes is smaller than that at low altitudes,and its changes are less marked.

As a normal gravity model,theJ4model is of high precision in near space,and the effect of theJ6term for correction can be completely ignored.The correction of theJ4term itself in near space is of the magnitude of 0.01 mg,and this magnitude is the same as that of the gravity disturbance in most regions of the world.The use of theJ4term for correction has a positive effect on reducing navigation errors,but the gravity disturbance in near space has a similar magni‐tude to theJ4term.Therefore,only adding theJ4term for correction usually fails to meet the navigation system’s requirement for a high-precision gravity model.To further improve the precision of the gravity model,it is critical to model the gravity disturbance and adopt a model that describes the Earth’s gravity with more accuracy.

For INS,theJ2orJ4normal gravity model is ap‐plicable if the used accelerometer constant bias is above 0.10 mg.When the gyroscope has sufficiently high precision and the precision of the accelerometer reaches 10 μgor higher,the model error of the normal gravity model becomes the main cause of error,and theJ2orJ4model is no longer applicable.This conclu‐sion is verified by simulations without gyroscope error and with a gyroscope constant bias of 0.003°/h.

Author contributions

Kai CHEN designed the research.Cheng-zhi ZENG pro‐cessed the corresponding data and wrote the first draft of the manuscript.Sen-sen PEI helped to organize the manuscript.Wen-chao LIANG revised and edited the final version.

Conflict of interest

Kai CHEN,Cheng-zhi ZENG,Sen-sen PEI,and Wenchao LIANG declare that they have no conflict of interest.