CN108225310A - A kind of Gravity-aided navigation path planning method - Google Patents

A kind of Gravity-aided navigation path planning method Download PDF

Info

Publication number
CN108225310A
CN108225310A CN201711414134.8A CN201711414134A CN108225310A CN 108225310 A CN108225310 A CN 108225310A CN 201711414134 A CN201711414134 A CN 201711414134A CN 108225310 A CN108225310 A CN 108225310A
Authority
CN
China
Prior art keywords
flight path
inertial navigation
random
gravity
error
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
CN201711414134.8A
Other languages
Chinese (zh)
Other versions
CN108225310B (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.)
707th Research Institute of CSIC
Original Assignee
707th Research Institute of CSIC
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 707th Research Institute of CSIC filed Critical 707th Research Institute of CSIC
Priority to CN201711414134.8A priority Critical patent/CN108225310B/en
Publication of CN108225310A publication Critical patent/CN108225310A/en
Application granted granted Critical
Publication of CN108225310B publication Critical patent/CN108225310B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20Instruments for performing navigational calculations

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

The present invention relates to a kind of Gravity-aided navigation path planning method, technical characterstic is:Include the following steps:Step 1 obtains the adaptation gravity anomaly Background in area and its error to standard deviation g by investigatingm;Investigate the performance parameter of each component devices of inertia/gravity integrated navigation system;Step 2 utilizes search --- and traversal flight path algorithm completes the construction and traversal of various possible random flight paths;Gravity measurement error to standard deviation g in step 3, the whole alternative flight paths of calculatingcWith inertial navigation position error diverging value δ during matching;Step 4, calculate whole alternative flight paths estimate matching positioning accuracy, the matching highest alternative flight path of positioning accuracy is as the planning flight path in the adaptation area.The present invention simulates submarine gravity survey process using gravity measurement error model, and the operating mode navigated by water under water using inertial navigation position error calculation method analog carrier can reappear practical underwater matching operating mode, so as to ensure that the reliability for filtering out flight path.

Description

A kind of Gravity-aided navigation path planning method
Technical field
The invention belongs to inertia/gravity integrated navigation technology fields, are related to path planning method, and especially a kind of gravity is auxiliary Assistant director of a film or play's boat path planning method.
Background technology
Gravity-aided navigation trajectory planning refers to filter out to carry out the boat of inertia/gravity integrated navigation in ocean Mark, integrated navigation system carry out matching positioning calculation on the flight path, high-precision location information can be provided, so as to fulfill right The biharmonic correction of inertial navigation.
At present, common path planning method such as A* algorithms, mostly with bathymetric data, gravity field feature statistic, as flight path The evaluation criterion of screening.These algorithms are mostly taking human as threshold value of the experience index of setting as evaluation criterion, and randomness is stronger, nothing The matching positioning accuracy of flight path that method guarantee filters out, poor reliability.And matching positioning accuracy can not be provided and united with gravity field feature Quantitative relationship between metering, intuitive are poor.
Invention content
It is an object of the invention to overcome the deficiencies of the prior art and provide a kind of intuitive and reliable and can consider carrier Underwater navigation operating mode, and to match positioning accuracy as the adaptation path planning method of screening and assessment standard.
The present invention solves its technical problem and following technical scheme is taken to realize:
A kind of Gravity-aided navigation path planning method, includes the following steps:
Step 1 obtains the adaptation gravity anomaly Background in area and its error to standard deviation g by investigatingm;Inertial navigation, again is set The operating mode of power instrument and the sail mode of carrier;The performance parameter of each component devices of inertia/gravity integrated navigation system is investigated, And inertial navigation position error maximum value Y during matching;
Under the conditions of step 2, the parameter set in step 1 and operating mode, search is utilized --- traversal flight path algorithm is completed The construction and traversal of various possible random flight paths:Degree calculating side is enriched with reference to inertial navigation position error calculation method and gravity field feature Method, under the conditions of filtering out random length, the arbitrary speed of a ship or plane in adaptation area, random flight path the abundantest Features of Gravity Anomalies GL, And calculate its match time length Tp, alternative flight path collection is formed, so as to complete the construction and traversal of various possible random flight paths;
Step 3 utilizes gravity measurement error model, under the conditions of calculating the parameter set in step 1 and operating mode, whole Gravity measurement error to standard deviation g on alternative flight pathc;Using inertial navigation position error calculation method, the parameter that step 1 is set is calculated Under the conditions of operating mode, inertial navigation position error diverging value δ during the matching on whole alternative flight paths;
Step 4:By the T of all alternative flight pathsp、gm、gc, GL, δ substitute into matching positioning accuracy predictor formula, calculate whole Alternative flight path estimates matching positioning accuracy, and the matching highest alternative flight path of positioning accuracy is as the planning boat in the adaptation area Mark.
Moreover, setting inertial navigation, the operating mode of gravimeter and the sail mode of carrier of the step 1, including following interior Hold:
(1) using the inertial navigation under horizontal damping operating mode, guiding carrier is at the uniform velocity sailed through to;
(2) using the gravimeter being operated under horizontal damping operating mode, operation is measured;
(3) information provided using the inertial navigation under horizontal damping operating mode carries out gravity anomaly correction.
Moreover, the performance parameter of each component devices of investigation inertia/gravity integrated navigation system of the step 1, including with Lower content:
(1) the components accuracy index of target inertial navigation:ΔAx、ΔAxr、ΔAy、ΔAyrInertial navigation east orientation, north orientation water are represented respectively The equivalent constant value zero bias of Ping Jiabiao and equivalent random error;εx、εxr、εy、εyrDo not represent the horizontal gyro of inertial navigation east orientation, north orientation etc. Imitate constant value drift and equivalent random error;
(2) the horizontal damping parameter of target inertial navigation:w1~w10For horizontal damping net coefficients, above-mentioned parameter can be by right Inertial navigation performance is investigated to obtain;
(3) the components accuracy index of target gravimeter:ΔAxg、ΔAxgr、ΔAyg、ΔAygrRepresent respectively gravimeter east orientation, North orientation level adds the equivalent constant value zero bias of table and equivalent random error;εxg、εxgr、εyg、εygrGravimeter east orientation, north are represented respectively Equivalent constant value drift and equivalent random error to horizontal gyro;
(4) the horizontal damping parameter of target gravimeter:wg1~wg10For horizontal damping net coefficients, above-mentioned parameter can lead to It crosses and gravimeter is investigated to obtain.
(5) the total accuracy of sounding index of latent instrument is surveyed in target depth measurement:δ h can use hand to survey latent instrument error, value from latent instrument is surveyed It is found in volume;
(6) target Doppler tachometer rate accuracy index:δVry、δVrxRepresent the north orientation speed that Doppler log provides Degree, east orientation speed error, value can be found from Doppler log service manual.
Moreover, the step 2 is using search --- traversal flight path algorithm complete it is various may random flight paths constructions and time The specific steps gone through include:Include the following steps:
(1) shape for being adapted to area is standardized:
Whole marginal points in some adaptation area are traversed, the maximum value for calculating space D is, Dis between any two points is It is adapted to the external diameter of a circle in area, corresponding two marginal points of maximum space D is are the point of contact for being adapted to area and circumscribed circle, two point of contacts The average value of coordinate is the central coordinate of circle of circumscribed circle, and the region in the range of circumscribed circle is regarded as the adaptation area after standardization;
(2) whether the size in judgement adaptation area is enough:
When the diameter Dis in the adaptation area after standardization is more than 2*Y nautical miles, the Features of Gravity Anomalies on actual position is enough It is abundant, therefore the search spread of flight path can be carried out in the adaptation area.The diameter Dis in adaptation area after standardization is less than 2* At Y nautical miles, the true flight path of carrier may be driven out to adaptation area, can not ensure to match positioning accuracy, therefore this kind of adaptation area is not It can use, directly give up;
(3) displacement error band traversal adaptation area:
In order to make full use of the Features of Gravity Anomalies in adaptation area, random flight path should be located on circumscribed circle diameter, at this time root According to the location diagram between circumscribed circle and random flight path, as X < Dis-2*Y, random flight path can be along circumscribed circle diameter Mobile, X represents random flight path length, and B represents translation position of the random flight path on circumscribed circle diameter, and K represents circumscribed circle diameter The course of random flight path;As the possibility value to B traverses, error band can be moved along the circumscribed circle diameter on a certain course It is dynamic, the adaptation area near Landfill covering this diameter;Then the possibility value of K is traversed, with the rotation of circumscribed circle diameter Turn, error band can completely cover circular adaptation area, so as to fulfill the traversal of flight path random to error band;
(4) search for --- the logical order of traversal:
Search --- the value of random flight path length X and speed of a ship or plane V are first solidifying in ergodic process, is X by length0The speed of a ship or plane is V0Random flight path corresponding to error band, adaptation area diametrically translate and as diameter rotates, while positioned according to inertial navigation Error calculation method solves inertial navigation position error;And level calculating method is enriched using gravity field feature, it is X to calculate whole length0、 The speed of a ship or plane is V0Random flight path GL values, it is X to take a random flight path of GL maximums as length0, speed of a ship or plane V0Under the conditions of it is standby Flight path is selected, and calculates the T of the flight pathp=X0/V0
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, filter out arbitrary flight path length, the arbitrary speed of a ship or plane Under the conditions of GL maximums random flight path, form alternative flight path collection, and calculate the T of all alternative flight pathsp=X/V, so as to complete The construction and traversal of various possible random flight paths.
Moreover, (4) step of the step 2 solves the specific of inertial navigation position error according to inertial navigation position error calculation method Step includes:
1. random flight path length is set first as X0In the sea, speed of a ship or plane V0Section;
2. solve the latitude and longitude coordinates of random flight path;
Set the initial heading K of search spread0Degree, initial translation position B0In the sea;It is K by course0Adaptation area diameter with The intersection point of circumscribed circle translates B along circumscribed circle diameter0In the sea, you can obtain the starting point of random flight path.It is pointed out from random flight path Hair is K along course0Diameter, intercept out length as X0Line segment in the sea, as required random flight path;It is positioned in view of matching The resolving frequency of system is 1HZ, therefore using linear interpolation algorithm, random flight path interpolation is become Tp=X0/V0A sampled point Position sequence, you can obtain carrier along random flight path navigate by water when, the latitude and longitude coordinates of any time;
3. solve carrier along random flight path navigate by water when, rotation and acceleration that gyro and Jia Biao sensitivities arrive;
Carrier is calculated first when being run on above-mentioned random flight path, target inertial navigation adds table and gyro sensitivity arrives acceleration and Rotation information is shown below:
ωx=-vy/R
In above formula, Ax、Ay、AzRepresent Jia Biao sensitive acceleration arrived under geographic coordinate system;ωx、ωy、ωzRepresent gyro The sensitive angular velocity of rotation arrived under geographic coordinate system, to treat resolving amount;
vx=V0*cos(K0)、vy=V0*sin(K0)、East orientation speed, the north orientation speed of above-mentioned random flight path are represented respectively Degree, latitude;
Ω represents rotational-angular velocity of the earth, and R represents earth radius, is constant.
G represents local absolute gravity value, can be calculated with the latitude combination WGS-84 ellipsoidal models of above-mentioned random flight path;
vzThe vertical velocity of carrier is represented, can be obtained from the data for the experimental record that goes to sea.
4. solve inertial navigation position error;
After the inertance element rotation that sensitivity arrives on random flight path and acceleration is got, it can be used to using horizontal damping The fundamental equation of guiding systems on computers simulates inertial navigation output, and longitude, the latitude of inertial navigation output are that carrier is true Flight path, the longitude of inertial navigation output, longitude, the latitude of latitude and random flight path subtract each other, you can obtain inertial navigation position error.Level resistance The fundamental equation of Buddhist nun's inertial navigation system is as follows:
Apx=Ax-β*g+ΔAx+ΔAxr
Apy=Ay+α*g+ΔAy+ΔAyr
Apx、ApyAdd east orientation acceleration, the north orientation acceleration of table output for inertial navigation level under geographic coordinate system;α, β generation respectively Table inertial navigation inertial platform deviates geographical horizontal attitude error angle;γ represents the course error that inertial reference calculation goes out;Vcx、Vcy λcThe carrier east orientation speed, north orientation speed, latitude, longitude that inertial reference calculation goes out are represented respectively;wcx、wcy、wczRepresentative is applied to gyro Torque, above-mentioned parameter is treats resolving amount;
Vry、VrxNorth orientation speed, the east orientation speed that Doppler log provides are represented, it can be by random flight path velocity vx、vyUpper superposition Doppler log range rate error δ Vry、δVrxIt obtains;
In the case that inertial navigation is operated in interior horizontal damping operating mode, K1=0, A=B=1, C=D=E=0, inertial navigation are operated in In the case of outer level damp operating mode, K1=1, A=B=1, C=D=E=0;
ua~ufFor the intermediate variable in solution process;
Above-mentioned equation group is only made of first-order linear homogeneous differential equation and linear function, by Runge-Kutta method It solves;Wherein λc、φcThe as true longitude and latitude of carrier subtracts each other its longitude with random flight path, latitude, you can obtain inertial navigation Position error.
Moreover, the gravity field feature that (4) step of the step 2 uses enriches level calculating method, calculating whole length is X0, speed of a ship or plane V0The specific steps of GL values of random flight path include:
1. true flight path is translated;
By by any one true flight path of certain random flight path generation along weft respectively east, west translation 0.5 it is extra large In, using bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
By by any one true flight path of certain random flight path generation along warp respectively south, north translation 0.5 it is extra large In, using bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
2. calculate the gravity anomaly gradient on true flight path:
On the true flight path, the gradient sequence Tx in longitudinal, the gradient sequence Ty on latitude direction, such as following formula institute Show:
3. calculate the GL of random flight path:
To any one random flight path, the true flight path quantity by the emulation generation of inertial navigation position error calculation method is 100, then all the average value of gradient variance can be defined as follows in true flight path longitudinal:
All the average value of gradient variance can be defined as follows on true flight path latitude direction:
Wherein, Txi、TyiLongitude and latitude gradient sequence on respectively i-th true flight path.
Then the Features of Gravity Anomalies of corresponding random flight path enriches degree GL, is represented by:
Level calculating method and inertial navigation position error calculation method are enriched with reference to above-mentioned gravity field feature, you can are calculated complete Minister's degree is X0The speed of a ship or plane is V0Random flight path GL values, a random flight path of GL maximums is that length is X0, speed of a ship or plane V0Item Alternative flight path under part.
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, you can filter out arbitrary flight path length, arbitrary The random flight path of GL maximums, forms alternative flight path collection under the conditions of the speed of a ship or plane.
Moreover, the gravity measurement error model that the step 3 utilizes, is shown below:
In above formula, g represents local absolute gravity value, can be calculated with the latitude combination WGS-84 ellipsoidal models of alternative flight path It arrives;Represent the latitude of alternative flight path;veRepresent the east orientation speed of alternative flight path;vnRepresent the north orientation speed of alternative flight path.Ω generations Table rotational-angular velocity of the earth, R represent earth radius, A1=978032.67714, B1=0.00193185138639, C1=- 0.00669437999013, D1=-0.3086 is constant;
In above formula, δ gcFor gravity measurement error to be solved, error source is divided into three classes:δ h take to survey latent instrument error Value can be found from surveying in latent instrument service manual;δaE、δaNδve、δvnThe error of information is provided for inertial navigation;α, β are gravimeter The misaligned angle of the platform;For rear two classes error source, horizontal damping inertial navigation fundamental equation combination Runge-Kutta method can be utilized to solve meter It obtains.Whole error sources are substituted into formula, you can acquire δ gcNumerical solution.Gravity error standard deviation g is surveyed during then matchingc =std (δ gc)。
Moreover, the step 3 utilizes inertial navigation position error calculation method, parameter and Working mould that step 1 is set are calculated Under the conditions of formula, all inertial navigation position error diverging value δ during the matching on alternative flight path, is shown below:
In above formula, subscript t ∈ [t1, t2], t1And t2Represent the matched beginning and ending time.λctIt represents using described in step 2 The longitude and latitude on the true flight path of carrier that inertial navigation position error calculation method is calculated;λtRepresent alternative flight path longitude and latitude. Inertial navigation position error diverging value δ=max (δ during then matchingt)。
Moreover, the matching positioning accuracy predictor formula of the step 4, is shown below:
Inertial navigation position error diverging value calculating method, meter during the gravity measurement error model that is provided using step 3, matching Calculate the g of all alternative flight pathsc, δ, and the whole being calculated with reference to step 2 is chosen GL, T of flight pathpAnd investigate what is obtained gm, substitute into matching positioning accuracy predictor formula, you can obtain whole alternative flight paths estimates matching positioning accuracy, all alternative boat The matching highest flight path of positioning accuracy, the planning flight path as finally obtained are estimated in mark.
The advantages and positive effects of the present invention are:
1st, the present invention is that one kind considers carrier underwater navigation operating mode, and commented as screening flight path to match positioning accuracy The adaptation path planning method of valency index, can directly give the matching highest flight path of positioning accuracy, and intuitive is strong.
2nd, screening process of the invention, has fully considered the characteristics of carrier navigates by water under water, has utilized gravity measurement error mould Pattern intends submarine gravity survey process, and the operating mode navigated by water under water using inertial navigation position error calculation method analog carrier can Reappear actual match operating mode, so as to ensure that the reliability for filtering out flight path.
3rd, the present invention is using search --- and traversal flight path algorithm enables screening process completely to cover adaptation area, so as to protect The optimality of the selection result is demonstrate,proved.
Description of the drawings
Fig. 1 is the logarithm amplitude-frequency response characteristic curve graph of the horizontal damping network of the present invention;
Fig. 2 is the ins error circle schematic diagram of the present invention;
Fig. 3 is the ins error band schematic diagram of the present invention;
Fig. 4 is the adaptation area standardization schematic diagram of the present invention;
Fig. 5 is the circumscribed circle of the present invention and the spatial position figure of error band;
Fig. 6 is the location diagram between the circumscribed circle of the present invention and planning flight path;
Fig. 7 is the process chart of the Path Planning of the present invention;
Inertial navigation position error diverging value figure during Fig. 8 is matching of the invention.
Specific embodiment
The embodiment of the present invention is described in further detail below in conjunction with attached drawing:
The purpose of trajectory planning is to filter out the highest track of matching precision, and background technology is gravity field data library adaptation Property assay method, the Gravity Matching precision estimation formula obtained using the evaluation method,
It is shown below:
In order to establish Path Planning, it is necessary first to matching precision predictor formula be analyzed, find out influence matching The parameter of precision, in the parameter occurred in above formula:
Match time length Tp, length X, the speed of a ship or plane V to planning flight path are related, Tp=X/V;
Inertial navigation position error diverging value δ during matching, length X, the speed of a ship or plane V to planning flight path are related, determine using inertial navigation Position error calculation method is calculated;
Gravity gradient standard deviation GL during matching on the true flight path of carrier, characterizes gravity anomaly on the true flight path of carrier The abundant degree of feature, value is related to the length X, course K, location B, the speed of a ship or plane V that plan flight path, available to search Rope --- traversal flight path algorithm combination inertial navigation position error calculation method, gravity field feature are enriched level calculating method and are calculated;
Gravity error standard deviation g is surveyed during matchingc, with plan the length X of flight path, speed of a ship or plane V, flight path shape whether be Whether straight line, headway stablize correlation, are calculated using gravity measurement error model.
It can be seen that the factors such as the location of flight path, shape, course, length, speed of a ship or plane, it can be to finally matching positioning Precision has an impact.The present invention is based on matching precision predictor formulas, to obtain high-precision matching position location as target, carry out boat Mark programme designs, and gives the method for determining the factors such as planning flight path position, shape, course, length, the speed of a ship or plane.
A kind of Gravity-aided navigation path planning method, as shown in fig. 7, comprises following steps:
Step 1 obtains the adaptation gravity anomaly Background in area and its error to standard deviation g by investigatingm;Inertial navigation, again is set The operating mode of power instrument and the sail mode of carrier;The performance parameter of each component devices of inertia/gravity integrated navigation system is investigated, And inertial navigation position error maximum value Y during matching;
(1) inertial navigation, the operating mode of gravimeter and the sail mode of carrier are set:
The realization of Gravity Matching depends on high-precision real gravity anomaly information, and gravity anomaly information is by relative gravity What instrument measured.In the signal of relative gravity instrument output, other than gravity anomaly information, also generated comprising carrier movement additional The distracters such as centrifugal acceleration, Corioli's acceleration, normal gravity, correspondingly, needing the output signal work to relative gravity instrument Normal field corrects and Etvs corrections, to complete the separation to distracter.These corrections are needed with carrier latitude information, level Based on velocity information.When carrier navigates by water under water, continuous, stable, prolonged speed, reference position information, generally by Inertial navigation provides, it can be seen that, the computational accuracy that the precision of inertial navigation output information will directly affect distracter, and then influence actual measurement weight The precision of power exception.
In order to promote the precision that inertial navigation provides information, inertial navigation need to be arranged under horizontal damping operating mode, horizontal damping network Logarithm amplitude-frequency response characteristic, as shown in Figure 1, horizontal damping operating mode can apply noise strong attenuation, promote what inertial navigation provide The precision of information, the horizontal attitude precision of gravimeter platform, and then the precision of real gravity anomaly is promoted, therefore, advised in flight path During drawing, inertial navigation and gravimeter need to be operated under horizontal damping operating mode.The function of the damping network mainly has following two side Face:On the one hand, which resolves in inertial navigation level and increases pole on channel, can eliminate the Schuler in the information that inertial navigation provides Periodic vibration;On the other hand, which can be to random walk part, the sea of gyroscopic drift random walk part plus table zero bias The high frequency impulse noise that the factors such as wave introduce inertial navigation system is damped, and accumulation of the inertial navigation to noise response is eliminated, so as to be promoted Inertial navigation provides the precision of information.But the introducing of horizontal damping network destroys the Schuler regularization condition of inertial navigation, it is violent in carrier Under maneuvering condition, overshoot can occur for the information that inertial navigation provides, and can reduce its precision instead, therefore in order to play damping network to used The inhibition of error is led, carrier is needed under inertial navigation guiding, the state for the direct route that remains a constant speed.
In addition, relative gravity instrument is substantially also an inertial navigation system.Maneuvering condition, the horizontal damping work at the uniform velocity sailed through to Condition creates one for promoting the horizontal attitude precision of relative gravity instrument electromechanics platform or mathematical platform, and then for gravimeter A stable working environment, is also advantageous.
Summary theory analysis, for guarantee matching precision, the equipment of carrier and its carrying need to be operated in following state:
Using the inertial navigation under horizontal damping operating mode, guiding carrier is at the uniform velocity sailed through to;
Using the gravimeter being operated under horizontal damping operating mode, operation is measured;
The information provided using horizontal damping inertial navigation carries out gravity anomaly correction.
(2) performance parameter of each component devices of inertia/gravity integrated navigation system, including the following contents:
The components accuracy index of target inertial navigation:ΔAx、ΔAxr、ΔAy、ΔAyrInertial navigation east orientation is represented respectively, north orientation level adds The equivalent constant value zero bias of table and equivalent random error;εx、εxr、εy、εyrThe equivalent normal of the horizontal gyro of inertial navigation east orientation, north orientation is not represented Value drift and equivalent random error;
The horizontal damping parameter of target inertial navigation:w1~w10For horizontal damping net coefficients, above-mentioned parameter can be by used Performance is led to be investigated to obtain.
The components accuracy index of target gravimeter:ΔAxg、ΔAxgr、ΔAyg、ΔAygrGravimeter east orientation, north are represented respectively To horizontal plus table equivalent constant value zero bias and equivalent random error;εxg、εxgr、εyg、εygrGravimeter east orientation, north orientation are represented respectively The equivalent constant value drift of horizontal gyro and equivalent random error;
The horizontal damping parameter of target gravimeter:Wg1~Wg10For horizontal damping net coefficients, above-mentioned parameter can be by right Gravimeter is investigated to obtain.
The total accuracy of sounding index of latent instrument is surveyed in target depth measurement:δ h can be from the latent instrument service manual of survey to survey latent instrument error, value In find;
Target Doppler tachometer rate accuracy index:δVry、δVrxRepresent north orientation speed that Doppler log provides, East orientation speed error, value can be found from Doppler log service manual.
Under the conditions of step 2, the parameter set in step 1 and operating mode, search is utilized --- traversal flight path algorithm is completed The construction and traversal of various possible random flight paths:Degree calculating side is enriched with reference to inertial navigation position error calculation method and gravity field feature Method, under the conditions of filtering out random length, the arbitrary speed of a ship or plane in adaptation area, random flight path the abundantest Features of Gravity Anomalies GL, And calculate its match time length Tp, alternative flight path collection is formed, so as to complete the construction and traversal of various possible random flight paths;
(1) step 2 described search --- traversal flight path algorithm specific implementation step is as follows:
It is guided when carrier matches under water by inertial navigation, along planning flight path navigation, the presence of inertial navigation position error can cause Carrier real trace deviates planning flight path, it is contemplated that it is fixed that the gravity gradient standard deviation GL in carrier real trace directly affects matching Position precision, therefore during trajectory planning, need to calculate the GL values for the true track being located in flight path near zone, it is suitable as it Evaluation criterion with property;Again by the method for search spread, planning flight path is filtered out from arbitrary flight path.It can be seen that:Arbitrarily The suitability of flight path is determined by the GL values for being located at the true track in flight path near zone:To the search spread of flight path, essence On be search spread to flight path near zone.
Provide search step by step below --- the implementation method of traversal flight path algorithm.
1) shape for being adapted to area is standardized:
In general, before trajectory planning is carried out, alternative marine site can be artificially divided into multiple portions, calculated respectively each GL in part delimit the part that GL is larger to be adapted to area.For the irregular adaptation area of shape, for the ease of search through It goes through, needs to carry out standardization processing to the shape for being adapted to area:Traverse some adaptation area whole marginal points, calculate any two points it Between the maximum value of space D is, Dis be to be adapted to the external diameter of a circle in area, corresponding two marginal points of maximum space D is are Area and the point of contact of circumscribed circle are adapted to, the average value of two point of contact coordinates is the central coordinate of circle of circumscribed circle.In the range of circumscribed circle Region is regarded as the adaptation area after standardization, as shown in figure 4, the method for original adaptation area's circumscribed circle is solved, after ensureing to standardize Adaptation area while be fully contemplated by arbitrary shape original adaptation area, it is ensured that can make full use of gravitational field during trajectory planning Feature.
2) whether the size in judgement adaptation area is enough:
Using inertial navigation indicating positions as the center of circle, using the maximum value Y of the position error of inertial navigation as radius, a circle of uncertainty is made, Then any time in the matching process, the actual position of carrier is certain to include in the circle, as shown in Fig. 2, being positioned with inertial navigation Max value of error Y is radius, using inertial navigation indicating positions as the center of circle, makes inertial navigation position error circle, carrier actual position one positions Inside the circle of uncertainty.During matching, carrier needs to sail through under inertial navigation guiding, with matched progress, inertial navigation indicating positions By by elongating as straight line, correspondingly, the error fenestra near inertial navigation indicating positions can be elongated to as a rectangular error band, The width of the error band is 2*Y nautical miles, at this time the real trace of carrier, the error band being still located near inertial navigation instruction flight path Interior, as shown in figure 3, during matching, carrier needs to sail through under inertial navigation guiding, with matched progress, inertial navigation indicating positions quilt It is elongated by as straight line, correspondingly, the error fenestra near inertial navigation indicating positions can be elongated to as a rectangular error band, it should The width of error band is 2*Y, and carrier actual position one is positioned inside the error band.
In order to ensure that real trace can be located in adaptation area to ensure that error band can be fallen completely in adaptation area Portion.When the diameter Dis in the adaptation area after standardization is more than 2*Y nautical miles, the adaptation area after standardization being capable of complete overlay errors Band it ensure that the Features of Gravity Anomalies on actual position is abundant enough, therefore can carry out searching for flight path in the adaptation area Rope traverses.When the diameter Dis in the adaptation area after standardization is less than 2*Y nautical miles, the true flight path of carrier may be driven out to adaptation Area can not ensure to match positioning accuracy, therefore this kind of adaptation area can not use, and directly give up.
3) displacement error band traversal adaptation area:
By the sail mode of carrier during matching it is found that matching during inertial navigation instruction flight path should with plan merged, because This, the real trace of carrier, one is positioned in the error band near planning flight path.It can be seen that error band can characterize flight path The region of true track is nearby included, is substantially to inertial navigation position error band to the search spread of flight path and its near zone Search spread.
When external diameter of a circle Dis is more than 2*Y nautical miles, shown in spatial relation Fig. 5 of circumscribed circle and error band, figure Middle inertial navigation instruction track (planning flight path) and the circle center distance of circumscribed circle are d.Due to can not be determined in advance most before trajectory planning The length of excellent flight path, in order to make full use of the Features of Gravity Anomalies in adaptation area, to promote the success of trajectory planning to greatest extent Rate is maximum it is necessary to try to ensure that the area that error band is fallen into adaptation area, and error band falls into the area of circumscribed circle, such as following formula institute Show:
By S to d derivations, as d=0, area S obtains maximum value, this shows that inertial navigation instruction track (planning flight path) is located at When on circumscribed circle diameter, the area that error band falls into adaptation area is maximum.Therefore, planning flight path will be searched on circumscribed circle diameter It arrives, correspondingly, the search spread of error band will be carried out along adaptation area's circumscribed circle diameter.
When plan flight path be located on circumscribed circle diameter when, circumscribed circle and plan flight path between position relationship as shown in fig. 6, In figure, X represents random flight path length, and as X < Dis-2*Y, random flight path can be moved along circumscribed circle diameter, be represented with B Random translation position of the flight path on circumscribed circle diameter, K represent the course of circumscribed circle diameter (random flight path).With to B can Energy value is traversed, and error band can be moved along the circumscribed circle diameter on a certain course, near Landfill covering this diameter It is adapted to area.Then the possibility value of K is traversed, with the rotation of circumscribed circle diameter, error band can completely cover circle The adaptation area of shape, so as to fulfill the traversal to error band (random flight path).Under normal circumstances, the traversal step-length of course K is 1 degree, The traversal step-length of position B is 1 nautical mile.
4) search for --- the logical order of traversal:
It needs to be determined completely by length X, speed of a ship or plane V, course tetra- parameters of K peace pan positions B, and each in view of random flight path A parameter is not quite similar, therefore searching for the Effect Mode for matching positioning accuracy --- in ergodic process, need reasonable arrangement The screening sequence of parameter determines the optimal value of flight path parameter one by one.
The discussion that is saved by background technology one has an effect on δ it is found that flight path length X and speed of a ship or plane V not only influence GL;And flight path Translation position B and course K only influences GL, does not influence δ.In line with the principle of Stepwise Screening, it is fixed to matching that this algorithm excludes δ first The influence of position precision namely curing flight path length X and speed of a ship or plane V, are X by length0The speed of a ship or plane is V0Random flight path corresponding to error Band diametrically translates and as diameter rotates in adaptation area, to the random flight path of arbitrary course K, arbitrary intercept B, utilizes third The inertial navigation position error calculation method provided is saved, 100 carrier real traces are constructed in the random error band corresponding to it, The GL in the whole 100 articles of carrier real traces of level calculating method calculating is enriched, and will using Section of five gravity field feature provided The mean value of 100 real trace GL, the evaluation index as the random flight path suitability.All in random flight path, GL maximums One random flight path, as length are X0In the sea and the speed of a ship or plane is V0Under the premise of, the alternative flight path in the adaptation area, T at this timep=X0/ V0
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, filter out arbitrary flight path length, the arbitrary speed of a ship or plane Under the conditions of GL maximums random flight path, form alternative flight path collection, and calculate the T of all alternative flight pathsp=X/V.So as to complete The search of various possible random flight paths and traversal.The value of random flight path length X can be 20 nautical miles or 30 nautical miles, speed of a ship or plane V Value can be 7 section or 10 section.
(2) inertial navigation position error calculation method specific implementation step is as follows described in step 2:
Degree direct influence matching positioning accuracy is enriched in view of the gravity field feature on the true flight path of carrier, in order to right The suitability of certain random flight path is analyzed, and needs to construct true flight path in its corresponding inertial navigation position error band, very Real flight path can be obtained by being superimposed inertial navigation position error on random flight path.The calculation method of inertial navigation position error is as follows:
Random flight path length is set first as X0In the sea, speed of a ship or plane V0Section.
1) latitude and longitude coordinates of random flight path are solved:
Set the initial heading K of search spread0Degree, initial translation position B0In the sea.It is K by course0Adaptation area diameter with The intersection point of circumscribed circle translates B along circumscribed circle diameter0In the sea, you can obtain the starting point of random flight path.It is pointed out from random flight path Hair is K along course0Diameter, intercept out length as X0Line segment in the sea, as required random flight path.It is positioned in view of matching The resolving frequency of system is 1HZ, therefore using linear interpolation algorithm, random flight path interpolation is become Tp=X0/V0A sampled point Position sequence, you can obtain carrier along random flight path navigate by water when, the latitude and longitude coordinates of any time.
2) solve carrier along random flight path navigate by water when, rotation and acceleration that gyro and Jia Biao sensitivities arrive:
Due to inertial navigation system be non-linear, time-varying system, inertial navigation position error diverging value δ during can not being matched Analytic solutions.In order to determine the value of δ, need to be first depending on the latitude of alternative flight path, east orientation speed, north orientation speed, resolve used Property the element sensitive rotation arrived and acceleration under geographic coordinate system.It is shown below:
ωx=-vy/R
In above formula, Ax、Ay、AzRepresent Jia Biao sensitive acceleration arrived under geographic coordinate system;ωx、ωy、ωzRepresent gyro The sensitive angular velocity of rotation arrived under geographic coordinate system, to treat resolving amount.
vx、vyThe east orientation speed, north orientation speed, latitude of alternative flight path are represented respectively.
Ω represents rotational-angular velocity of the earth, and R represents earth radius, is constant.
G represents local absolute gravity value, can be calculated with the latitude combination WGS-84 ellipsoidal models of alternative flight path.
vzThe vertical velocity of carrier is represented, can be obtained from the data for the experimental record that goes to sea.
3) inertial navigation position error is solved:
After the inertance element rotation that sensitivity arrives on random flight path and acceleration is got, it can be used to using horizontal damping The fundamental equation of guiding systems on computers simulates inertial navigation output.Longitude, the latitude of inertial navigation output are that carrier is true Flight path can be used to calculate the GL of random flight path;The longitude of inertial navigation output, longitude, the latitude of latitude and random flight path subtract each other, you can Obtain inertial navigation position error δ.The fundamental equation of horizontal damping inertial navigation system is as follows:
Apx=Ax-β*g+ΔAx+ΔAxr
Apy=Ay+α*g+ΔAy+ΔAyr
A in equation group (2)px、ApyEast orientation acceleration, the north orientation for table being added to export for inertial navigation level under geographic coordinate system accelerate Degree;α, β represent inertial navigation inertial platform and deviate geographical horizontal attitude error angle respectively;γ represents the course that inertial reference calculation goes out and misses Difference;Vcx、VcyλcThe carrier east orientation speed, north orientation speed, latitude, longitude that inertial reference calculation goes out are represented respectively;wcx
wcy、wczThe torque that inertial navigation system is applied to gyro is represented, above-mentioned parameter is treats resolving amount.
ΔAx、ΔAxr、ΔAy、ΔAyrInertial navigation east orientation is represented respectively, north orientation level adds the constant value zero bias and random error of table; εx、εxr、εy、εyrThe constant value drift and random error of the horizontal gyro of inertial navigation east orientation, north orientation are not represented;w1~w10For horizontal damping net Network coefficient, above-mentioned parameter can be by being investigated to obtain to inertial navigation.
Vry、VrxNorth orientation speed, the east orientation speed that Doppler log provides are represented, the speed in alternative flight path can be passed through Doppler log range rate error is superimposed on degree to obtain.
In the case that inertial navigation is operated in interior horizontal damping operating mode, K1=0, A=B=1, C=D=E=0, inertial navigation are operated in In the case of outer level damp operating mode, K1=1, A=B=1, C=D=E=0.
ua~ufFor the intermediate variable in solution process.
Equation group (2) is only made of first-order linear homogeneous differential equation and linear function, can be asked by Runge-Kutta method Solution.Wherein λc、φcThe as true longitude and latitude of carrier subtracts each other its longitude with random flight path, latitude, you can obtain inertial navigation and determine Position error.
It is random in view of inertial navigation components error, therefore carrier actual position is likely located in error band during matching Any position.Operating mode is really matched in order to simulate, during trajectory planning, to any one random flight path, by random The mode of producing element error, you can construct 100 true flight paths in the range of the corresponding error band of the random flight path, be used for The Features of Gravity Anomalies for calculating the random flight path enriches degree.
(3) gravity field feature that the step 2 uses enriches level calculating method, includes the following steps:
In view of Gravity Matching location algorithm using the tendency of gravity anomaly as observed quantity, therefore Features of Gravity Anomalies rises Volt directly reflects the abundant degree of Features of Gravity Anomalies, and the only violent region of gravity anomaly fluctuating could be that matching positioning is calculated Method provides effective observed quantity.The fluctuating of gravity anomaly can be characterized with the standard deviation GL of gravity anomaly gradient, and GL is bigger, Prove that the gravity anomaly gradient each opposite sex of the gravitational field on each position, direction is stronger, correspondingly, the fluctuating of gravity anomaly is also got over Acutely.
On gravity map, by searching for --- traversal flight path algorithm generates random flight path, and the location sets on flight path areIn view of latent device under water by inertial navigation guiding navigation, therefore the real trace of carrier will necessarily deviate at random Flight path, using horizontal damping inertial navigation fundamental equation, based on random flight path, you can construct carrier real trace.
By by any one true flight path of certain random flight path generation along weft respectively east, west translation 0.5 it is extra large In, using bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
By by any one true flight path of certain random flight path generation along warp respectively south, north translation 0.5 it is extra large In, using bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
Then on the true flight path, the gradient sequence Tx in longitudinal, the gradient sequence Ty on latitude direction, such as following formula institute Show:
To any one random flight path, the true flight path quantity by the emulation generation of inertial navigation position error calculation method is 100, then all the average value of gradient variance can be defined as follows in true flight path longitudinal:
All the average value of gradient variance can be defined as follows on true flight path latitude direction:
Wherein, Txi、TyiLongitude and latitude gradient sequence on respectively i-th true flight path.
Then the Features of Gravity Anomalies of corresponding random flight path enriches degree GL, is represented by:
Level calculating method and inertial navigation position error calculation method are enriched with reference to above-mentioned gravity field feature, you can are calculated complete Minister's degree is X0The speed of a ship or plane is V0Random flight path GL values, a random flight path of GL maximums is that length is X0, speed of a ship or plane V0Item Alternative flight path under part.
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, you can filter out arbitrary flight path length, arbitrary The random flight path of GL maximums, forms alternative flight path collection under the conditions of the speed of a ship or plane.
Step 3:The gravity measurement error model designed using the present invention calculates parameter and operating mode that step 1 is set Under the conditions of, the gravity measurement error to standard deviation g on whole alternative flight pathsc;It is resolved using the inertial navigation position error that the present invention designs Method, under the conditions of calculating the parameter and operating mode that step 1 sets, whole inertial navigation position errors during the matching on alternative flight paths Diverging value δ.
(1) the gravity measurement error model that the step 3 uses, including the following contents:
Since Gravity Matching location algorithm is using the fluctuating of gravity anomaly as observed quantity, matching starts preceding gravity sensitive The accumulated error that device drift generates, will not have an impact matching precision;The matching duration is usually no more than 3 hours, herein Period, the drift error accumulated in gravimeter output signal do not exceed ± 0.05mgal;In addition, environment temperature in latent device Stablize relatively, the temperature control system of gravimeter itself equipment, can be by the temperature of gravimeter near zone under latent device environment Control is within ± 0.1 DEG C, therefore the gravity measurement error as caused by temperature change, can be neglected.Rejecting above-mentioned error Afterwards, actual measurement gravity error can be divided into following three types according to its source:
Under horizontal damping working condition, gravimeter platform coordinate system can deviate local geographic coordinate system, cause instrument sensitive Axis deviates gravimetric plumb line direction, generates absolute gravity measurement error.
Carrier is using inertial navigation guiding navigation, and therefore, horizontal acceleration inevitably has manipulation error, the error meeting By the attitude error of gravimeter platform, it is coupled into gravity sensitive axis direction.
Horizontal damping inertial navigation velocity error, horizontal damping inertial navigation latitude position error survey the latent deep error that latent instrument obtains, and lead The gravity anomaly correction error of cause.
In this section, modeling description is carried out to this three classes error first, establishes gravity measurement error model;Then to model packet The various error sources contained are analyzed, and provide the computational methods for solving error source numerical solution.
For the gravimeter being operated under horizontal damping operating mode, it is assumed that the unit vector of its sensitive axis direction is usedIt represents,It is assumed that its east orientation, north orientation platform error angle are respectively:WithThe equal very little in the two angles, therefore vectorWith local ground Reason coordinate system day is to the amplitude of angle:
Actual measurement gravity error, is represented by caused by being deflected by gravimeter platform stance: Its In,For the absolute gravity vector at alternative flight path, can be calculated with the latitude combination WGS-84 ellipsoidal models of alternative flight path It arrives.The value of error is as shown in right formula:Due toTherefore very little, can be approximately consideredWithSide To consistent, it can thus be concluded that:
Above formula is measurement error, the scalar shape being converted into as unit of mgal caused by the deflection of gravimeter platform stance Formula is as follows:
For the carrier using inertial navigation guiding navigation under horizontal damping operating mode, it is assumed that its east orientation, north orientation acceleration are respectively:WithThen carrier east orientation acceleration is projected as on gravity sensitive axis:Direction is by the right-hand rule It obtains, value is: In formulaForFrom toAngle, it is suitable counterclockwise for just Hour hands are negative.Likewise, north orientation acceleration is projected as on gravity sensitive axis:Determined by the right hand in direction It then obtains, value is:In formulaForFrom toAngle, the inverse time Needle is is clockwise negative just.
WithVector sum be:Its value is shown below:
Above formula is coupled into the measurement error of gravity sensitive axis for carrier levels acceleration, be converted into using mgal as The scalar form of unit is as follows:
δg4=(δ aE*β-δaN*α)*105 (4)
Etvs corrections formula, normal field correction formula and vertical acceleration correction formula as unit of mgal are such as Under:
Wherein A1=978032.67714, B1=0.00193185138639,
C1=-0.00669437999013, D1=-0.3086, T is the carrier height sampling period, in general T=1s, Ω For rotational-angular velocity of the earth, R is earth radius, veFor the east orientation speed of alternative flight path, vnFor the north orientation speed of alternative flight path, For the latitude of alternative flight path, h is to survey the carrier that latent instrument provides to dive depth.Will three formula above, respectively to latitude, east speed, north speed, Height derivation can obtain:
δgh=(δ h (t-1)+δ h (t+1) -2* δ h (t)) * 105 (7)
Formula (3)~(7) are combined, you can obtain gravity measurement error model, be shown below:
Error source in above formula is divided into three classes:δ h can be found to survey latent instrument error, value from surveying in latent instrument service manual; δaE、δaNδve、δvnThe error of information is provided for inertial navigation;α, β are gravimeter the misaligned angle of the platform;For rear two classes error Source can utilize the solution of horizontal damping inertial navigation fundamental equation combination Runge-Kutta method to be calculated.Whole error sources are substituted into Formula, you can acquire δ gcNumerical solution.Gravity error standard deviation g is surveyed during then matchingc=std (δ gc)。
(2) inertial navigation position error diverging value δ, computational methods are as follows during the matching described in the step 3:
Inertial navigation position error transmitting case is as shown in figure 8, by inertial navigation fundamental equation using Long Ze --- Ku Tafa is resolved It arrives.
In above formula, subscript t ∈ [t1, t2], t1And t2Represent the matched beginning and ending time.λctIt represents using described in step 2 The longitude and latitude on the true flight path of carrier that inertial navigation position error calculation method is calculated;λtRepresent alternative flight path longitude and latitude. Inertial navigation position error diverging value δ=max (δ during then matchingt)。
Step 4:By the T of all alternative flight pathsp、gm、gc, GL, δ substitute into matching positioning accuracy predictor formula, can calculate complete The alternative flight path in portion estimates matching positioning accuracy, matches the highest alternative flight path of positioning accuracy, the as planning in the adaptation area Flight path.
The matching positioning accuracy predictor formula of the step 4, is shown below:
Inertial navigation position error diverging value calculating method, meter during the gravity measurement error model that is provided using step 3, matching Calculate the g of all alternative flight pathsc、δ.The whole being calculated with reference to step 2 is chosen GL, T of flight pathpAnd investigate what is obtained gm, substitute into matching positioning accuracy predictor formula, you can obtain whole alternative flight paths estimates matching positioning accuracy, all alternative boat The matching highest flight path of positioning accuracy, the planning flight path as finally obtained are estimated in mark.
It is emphasized that embodiment of the present invention is illustrative rather than limited, therefore present invention packet Include the embodiment being not limited to described in specific embodiment, it is every by those skilled in the art according to the technique and scheme of the present invention The other embodiment obtained, also belongs to the scope of protection of the invention.

Claims (9)

1. a kind of Gravity-aided navigation path planning method, it is characterised in that:Include the following steps:
Step 1 obtains the adaptation gravity anomaly Background in area and its error to standard deviation g by investigatingm;Inertial navigation, gravimeter are set The sail mode of operating mode and carrier;Investigate each component devices of inertia/gravity integrated navigation system performance parameter and With period inertial navigation position error maximum value Y;
Under the conditions of step 2, the parameter set in step 1 and operating mode, search is utilized --- traversal flight path algorithm is completed various The construction and traversal of possible random flight path:Level calculating method is enriched with reference to inertial navigation position error calculation method and gravity field feature, Under the conditions of random length, the arbitrary speed of a ship or plane being filtered out in adaptation area, random flight path the abundantest Features of Gravity Anomalies GL, and count Calculate its match time length Tp, alternative flight path collection is formed, so as to complete the construction and traversal of various possible random flight paths;
Step 3, using gravity measurement error model, it is all alternative under the conditions of calculating the parameter set in step 1 and operating mode Gravity measurement error to standard deviation g on flight pathc;Using inertial navigation position error calculation method, parameter and work that step 1 is set are calculated Under the conditions of operation mode, inertial navigation position error diverging value δ during the matching on whole alternative flight paths;
Step 4, the T by all alternative flight pathsp、gm、gc, GL, δ substitute into matching positioning accuracy predictor formula, calculate all alternative Flight path estimates matching positioning accuracy, and the matching highest alternative flight path of positioning accuracy is as the planning flight path in the adaptation area.
2. a kind of Gravity-aided navigation path planning method according to claim 1, it is characterised in that:The step 1 Inertial navigation, the operating mode of gravimeter and the sail mode of carrier are set, including the following contents:
(1) using the inertial navigation under horizontal damping operating mode, guiding carrier is at the uniform velocity sailed through to;
(2) using the gravimeter being operated under horizontal damping operating mode, operation is measured;
(3) information provided using the inertial navigation under horizontal damping operating mode carries out gravity anomaly correction.
3. a kind of Gravity-aided navigation path planning method according to claim 1 or 2, it is characterised in that:The step 1 Each component devices of investigation inertia/gravity integrated navigation system performance parameter, including the following contents:
(1) the components accuracy index of target inertial navigation:ΔAx、ΔAxr、ΔAy、ΔAyrInertial navigation east orientation is represented respectively, north orientation level adds The equivalent constant value zero bias of table and equivalent random error;εx、εxr、εy、εyrThe equivalent normal of the horizontal gyro of inertial navigation east orientation, north orientation is not represented Value drift and equivalent random error;
(2) the horizontal damping parameter of target inertial navigation:w1~w10For horizontal damping net coefficients, above-mentioned parameter can be by inertial navigation Performance is investigated to obtain;
(3) the components accuracy index of target gravimeter:ΔAxg、ΔAxgr、ΔAyg、ΔAygrGravimeter east orientation, north orientation are represented respectively The equivalent constant value zero bias and equivalent random error of level plus table;εxg、εxgr、εyg、εygrGravimeter east orientation, north orientation water are represented respectively The equivalent constant value drift of flat gyro and equivalent random error;
(4) the horizontal damping parameter of target gravimeter:wg1~wg10For horizontal damping net coefficients, above-mentioned parameter can be by right Gravimeter is investigated to obtain;
(5) the total accuracy of sounding index of latent instrument is surveyed in target depth measurement:δ h can be from the latent instrument service manual of survey to survey latent instrument error, value It finds;
(6) target Doppler tachometer rate accuracy index:δVry、δVrxRepresent north orientation speed, the east that Doppler log provides To velocity error, value can be found from Doppler log service manual.
4. a kind of Gravity-aided navigation path planning method according to claims 1 or 2 or 3, it is characterised in that:The step Rapid 2 using search --- and traversal flight path algorithm completes the various possible random constructions of flight path and the specific steps of traversal include:
(1) shape for being adapted to area is standardized:
Whole marginal points in some adaptation area are traversed, the maximum value for calculating space D is, Dis between any two points is to be adapted to The external diameter of a circle in area, corresponding two marginal points of maximum space D is are the point of contact for being adapted to area and circumscribed circle, two point of contact coordinates Average value be circumscribed circle central coordinate of circle, by the region in the range of circumscribed circle be regarded as standardization after adaptation area;
(2) whether the size in judgement adaptation area is enough:
When the diameter Dis in the adaptation area after standardization is more than 2*Y nautical miles, the Features of Gravity Anomalies on actual position is rich enough Richness, therefore the search spread of flight path can be carried out in the adaptation area;The diameter Dis in adaptation area after standardization is less than 2*Y Nautical mile when, the true flight path of carrier may be driven out to adaptation area, can not ensure to match positioning accuracy, therefore this kind of adaptation area cannot It is enough to use, directly give up;
(3) displacement error band traversal adaptation area:
In order to make full use of the Features of Gravity Anomalies in adaptation area, random flight path should be located on circumscribed circle diameter, at this time according to outer The location diagram between circle and random flight path is connect, as X < Dis-2*Y, random flight path can be moved along circumscribed circle diameter, X represents random flight path length, and B represents translation position of the random flight path on circumscribed circle diameter, and K represents circumscribed circle diameter and navigates at random The course of mark;As the possibility value to B traverses, error band can be moved along the circumscribed circle diameter on a certain course, most The adaptation area near this diameter is covered eventually;Then the possibility value of K is traversed, with the rotation of circumscribed circle diameter, missed Difference band can completely cover circular adaptation area, so as to fulfill the traversal of flight path random to error band;
(4) search for --- the logical order of traversal:
Search --- the value of random flight path length X and speed of a ship or plane V are first solidifying in ergodic process, is X by length0The speed of a ship or plane is V0's Error band corresponding to random flight path diametrically translates and as diameter rotates, while in adaptation area according to inertial navigation position error Calculation method solves inertial navigation position error;And level calculating method is enriched using gravity field feature, it is X to calculate whole length0, the speed of a ship or plane For V0Random flight path GL values, it is X to take a random flight path of GL maximums as length0, speed of a ship or plane V0Under the conditions of alternative boat Mark, and calculate the T of the flight pathp=X0/V0
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, filter out arbitrary flight path length, arbitrary speed of a ship or plane condition The random flight path of lower GL maximums, forms alternative flight path collection, and calculates the T of all alternative flight pathsp=X/V, it is various so as to complete The construction and traversal of possible random flight path.
5. a kind of Gravity-aided navigation path planning method according to claim 4, it is characterised in that:The step 2 The specific steps that (4) step solves inertial navigation position error according to inertial navigation position error calculation method include:
1. random flight path length is set first as X0In the sea, speed of a ship or plane V0Section;
2. solve the latitude and longitude coordinates of random flight path;
Set the initial heading K of search spread0Degree, initial translation position B0In the sea;It is K by course0Adaptation area diameter with it is external Round intersection point translates B along circumscribed circle diameter0In the sea, you can obtain the starting point of random flight path;From random flight path starting point, edge Course is K0Diameter, intercept out length as X0Line segment in the sea, as required random flight path;In view of matching alignment system Resolving frequency for 1HZ, therefore using linear interpolation algorithm, random flight path interpolation is become into Tp=X0/V0The position of a sampled point Sequence, you can obtain carrier along random flight path navigate by water when, the latitude and longitude coordinates of any time;
3. solve carrier along random flight path navigate by water when, rotation and acceleration that gyro and Jia Biao sensitivities arrive;
Carrier is calculated first when being run on above-mentioned random flight path, target inertial navigation adds table and gyro sensitivity arrives acceleration and rotation Information is shown below:
ωx=-vy/R
In above formula, Ax、Ay、AzRepresent Jia Biao sensitive acceleration arrived under geographic coordinate system;ωx、ωy、ωzGyro is represented on ground The sensitive angular velocity of rotation arrived under reason coordinate system, to treat resolving amount;
vx=V0*cos(K0)、vy=V0*sin(K0)、The east orientation speed, north orientation speed, latitude of above-mentioned random flight path are represented respectively Degree;
Ω represents rotational-angular velocity of the earth, and R represents earth radius, is constant;
G represents local absolute gravity value, can be calculated with the latitude combination WGS-84 ellipsoidal models of above-mentioned random flight path;
vzThe vertical velocity of carrier is represented, can be obtained from the data for the experimental record that goes to sea;
4. solve inertial navigation position error;
After the inertance element rotation that sensitivity arrives on random flight path and acceleration is got, horizontal damping inertial navigation system can be utilized The fundamental equation of system on computers simulates inertial navigation output, and longitude, the latitude of inertial navigation output are that carrier really navigates Mark, the longitude of inertial navigation output, longitude, the latitude of latitude and random flight path subtract each other, you can obtain inertial navigation position error;Horizontal damping The fundamental equation of inertial navigation system is as follows:
Apx=Ax-β*g+ΔAx+ΔAxr
Apy=Ay+α*g+ΔAy+ΔAyr
Apx、ApyAdd east orientation acceleration, the north orientation acceleration of table output for inertial navigation level under geographic coordinate system;α, β represent used respectively It leads inertial platform and deviates geographical horizontal attitude error angle;γ represents the course error that inertial reference calculation goes out;Vcx、VcyλcPoint The carrier east orientation speed, north orientation speed, latitude, longitude that inertial reference calculation goes out are not represented;wcx、wcy、wczRepresent the torsion for being applied to gyro Square, above-mentioned parameter is treats resolving amount;
Vry、VrxNorth orientation speed, the east orientation speed that Doppler log provides are represented, it can be by random flight path velocity vx、vy Upper superposition Doppler log range rate error δ Vry、δVrxIt obtains;
In the case that inertial navigation is operated in interior horizontal damping operating mode, K1=0, A=B=1, C=D=E=0, inertial navigation are operated in outer level In the case of damping operating mode, K1=1, A=B=1, C=D=E=0;
ua~ufFor the intermediate variable in solution process;
Above-mentioned equation group is only made of first-order linear homogeneous differential equation and linear function, can be solved by Runge-Kutta method; Wherein λc、φcThe as true longitude and latitude of carrier subtracts each other its longitude with random flight path, latitude, you can obtains inertial navigation positioning Error.
6. a kind of Gravity-aided navigation path planning method according to claim 4, it is characterised in that:The step 2 (4) step enriches level calculating method using gravity field feature, and it is X to calculate whole length0, speed of a ship or plane V0Random flight path GL values Specific steps include:
1. true flight path is translated;
By by any one true flight path of certain random flight path generation along weft respectively east, west translate 0.5 nautical mile, it is sharp With bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
By by any one true flight path of certain random flight path generation along warp respectively south, 0.5 nautical mile of north translation, it is sharp With bilinear interpolation, gravity anomaly value sequence on this two flight paths can be obtained, respectively:
2. calculate the gravity anomaly gradient on true flight path:
On the true flight path, the gradient sequence Tx in longitudinal, the gradient sequence Ty on latitude direction are shown below:
3. calculate the GL of random flight path:
To any one random flight path, the true flight path quantity by the emulation generation of inertial navigation position error calculation method is 100, Then all the average value of gradient variance can be defined as follows in true flight path longitudinal:
All the average value of gradient variance can be defined as follows on true flight path latitude direction:
Wherein, Txi、TyiLongitude and latitude gradient sequence on respectively i-th true flight path;
Then the Features of Gravity Anomalies of corresponding random flight path enriches degree GL, is represented by:
Level calculating method and inertial navigation position error calculation method are enriched with reference to above-mentioned gravity field feature, you can whole length are calculated It spends for X0The speed of a ship or plane is V0Random flight path GL values, a random flight path of GL maximums is that length is X0, speed of a ship or plane V0Under the conditions of Alternative flight path;
And then random flight path length X and speed of a ship or plane V is changed, according to the method described above, you can filter out arbitrary flight path length, the arbitrary speed of a ship or plane Under the conditions of GL maximums random flight path, form alternative flight path collection.
7. a kind of Gravity-aided navigation path planning method according to claims 1 or 2 or 3, it is characterised in that:The step The rapid 3 gravity measurement error models utilized, are shown below:
In above formula, g represents local absolute gravity value, can be calculated with the latitude combination WGS-84 ellipsoidal models of alternative flight path; Represent the latitude of alternative flight path;veRepresent the east orientation speed of alternative flight path;vnRepresent the north orientation speed of alternative flight path;Ω represents ground Revolutions angular speed, R represent earth radius, A1=978032.67714, B1=0.00193185138639, C1=- 0.00669437999013, D1=-0.3086 is constant;
In above formula, δ gcFor gravity measurement error to be solved, error source is divided into three classes:For δ h to survey latent instrument error, value can It is found from surveying in latent instrument service manual;δaE、δaNδve、δvnThe error of information is provided for inertial navigation;α, β are gravimeter platform Misalignment;For rear two classes error source, the solution of horizontal damping inertial navigation fundamental equation combination Runge-Kutta method can be utilized to calculate It arrives;Whole error sources are substituted into formula, you can acquire δ gcNumerical solution;Gravity error standard deviation g is surveyed during then matchingc= std(δgc)。
8. a kind of Gravity-aided navigation path planning method according to claims 1 or 2 or 3, it is characterised in that:The step It is all alternative to navigate under the conditions of rapid 3 parameter set using inertial navigation position error calculation method, calculating step 1 and operating mode Inertial navigation position error diverging value δ, is shown below during matching on mark:
In above formula, subscript t ∈ [t1, t2], t1And t2Represent the matched beginning and ending time;λctIt represents and utilizes inertial navigation described in step 2 The longitude and latitude on the true flight path of carrier that position error calculation method is calculated;λtRepresent alternative flight path longitude and latitude;Then With period inertial navigation position error diverging value δ=max (δt)。
9. a kind of Gravity-aided navigation path planning method according to claims 1 or 2 or 3, it is characterised in that:The step Matching positioning accuracy predictor formula in rapid 4, is shown below:
Inertial navigation position error diverging value calculating method, calculates during the gravity measurement error model that is provided using step 3, matching All g of alternative flight pathc, δ, and the whole being calculated with reference to step 2 is chosen GL, T of flight pathpAnd the g that investigation obtainsm, Substitute into matching positioning accuracy predictor formula, you can obtain whole alternative flight paths estimates matching positioning accuracy, all alternative flight paths In estimate matching the highest flight path of positioning accuracy, the planning flight path as finally obtained.
CN201711414134.8A 2017-12-22 2017-12-22 Gravity-assisted navigation track planning method Active CN108225310B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711414134.8A CN108225310B (en) 2017-12-22 2017-12-22 Gravity-assisted navigation track planning method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711414134.8A CN108225310B (en) 2017-12-22 2017-12-22 Gravity-assisted navigation track planning method

Publications (2)

Publication Number Publication Date
CN108225310A true CN108225310A (en) 2018-06-29
CN108225310B CN108225310B (en) 2020-08-25

Family

ID=62648578

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711414134.8A Active CN108225310B (en) 2017-12-22 2017-12-22 Gravity-assisted navigation track planning method

Country Status (1)

Country Link
CN (1) CN108225310B (en)

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109059964A (en) * 2018-09-19 2018-12-21 中国船舶重工集团公司第七0七研究所 A kind of inertial navigation based on gravity peak and the double calibration methods of gravity measurement
CN109085655A (en) * 2018-09-19 2018-12-25 中国船舶重工集团公司第七0七研究所 A kind of underwater platform gravity measurement scheme and verification method
CN109141426A (en) * 2018-08-10 2019-01-04 中国空间技术研究院 A kind of preferred method in subaqueous gravity matching navigation adaptation area
CN109269526A (en) * 2018-07-16 2019-01-25 哈尔滨工程大学 Rotary grid inertial navigation horizontal damping method based on damping network
CN110031000A (en) * 2019-05-21 2019-07-19 北京理工大学 A kind of evaluation method of Method in Gravity Aided INS region suitability
CN110196050A (en) * 2019-05-29 2019-09-03 哈尔滨工程大学 A kind of vertical height of Strapdown Inertial Navigation System and speed measurement method
CN111006694A (en) * 2019-12-29 2020-04-14 北京理工大学 Design method of long-endurance inertial navigation system track generator based on track planning
CN111426313A (en) * 2020-04-26 2020-07-17 中国人民解放军61540部队 Submarine navigation method and system based on gravity beacon
CN111473790A (en) * 2020-04-26 2020-07-31 中国人民解放军61540部队 Submarine navigation method and system of gravity beacon along track
CN111721300A (en) * 2020-06-30 2020-09-29 清华大学 Gravity beacon navigation method and system based on artificial intelligence algorithm
CN112541046A (en) * 2020-11-30 2021-03-23 中国电子科技集团公司第二十八研究所 Co-occurrence target monitoring method based on time and space
CN112729275A (en) * 2021-01-08 2021-04-30 中国船舶重工集团公司第七0七研究所 Satellite inversion chart gravity adaptation area selection method utilizing factor analysis
CN113503871A (en) * 2021-05-19 2021-10-15 北京理工大学 Gravity matching method based on correlation filtering
CN113514059A (en) * 2021-05-19 2021-10-19 北京理工大学 Gravity-assisted inertial navigation system simulation platform
CN113587923A (en) * 2021-05-31 2021-11-02 中国人民解放军61540部队 Submersible vehicle positioning method and system for screening multi-dimensional gravity gradient lighthouse matching area
CN114136320A (en) * 2021-11-19 2022-03-04 中国船舶重工集团公司第七0七研究所 Ocean multi-physical-field parameter feature positioning fusion method based on feature complementary characteristics
US11268813B2 (en) 2020-01-13 2022-03-08 Honeywell International Inc. Integrated inertial gravitational anomaly navigation system

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102128625A (en) * 2010-12-08 2011-07-20 北京航空航天大学 Initial matching method for use in gravimetric map matching in gravity-aided inertial navigation system
CN103604430A (en) * 2013-11-06 2014-02-26 哈尔滨工程大学 Marginalized cubature Kalman filter (CKF)-based gravity aided navigation method
CN105043388A (en) * 2015-06-29 2015-11-11 中国船舶重工集团公司第七0七研究所 Vector search iterative matching method based on inertia/gravity matching integrated navigation
CN105157703A (en) * 2015-06-03 2015-12-16 北京理工大学 Evaluation method of navigability of gravity-assisted inertial navigation adaptation zone
CN107289973A (en) * 2017-06-22 2017-10-24 北京理工大学 A kind of gravitational field suitability determination methods in Gravity Matching navigation

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102128625A (en) * 2010-12-08 2011-07-20 北京航空航天大学 Initial matching method for use in gravimetric map matching in gravity-aided inertial navigation system
CN103604430A (en) * 2013-11-06 2014-02-26 哈尔滨工程大学 Marginalized cubature Kalman filter (CKF)-based gravity aided navigation method
CN105157703A (en) * 2015-06-03 2015-12-16 北京理工大学 Evaluation method of navigability of gravity-assisted inertial navigation adaptation zone
CN105043388A (en) * 2015-06-29 2015-11-11 中国船舶重工集团公司第七0七研究所 Vector search iterative matching method based on inertia/gravity matching integrated navigation
CN107289973A (en) * 2017-06-22 2017-10-24 北京理工大学 A kind of gravitational field suitability determination methods in Gravity Matching navigation

Cited By (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109269526A (en) * 2018-07-16 2019-01-25 哈尔滨工程大学 Rotary grid inertial navigation horizontal damping method based on damping network
CN109269526B (en) * 2018-07-16 2022-06-17 哈尔滨工程大学 Rotary grid inertial navigation horizontal damping method based on damping network
CN109141426B (en) * 2018-08-10 2020-11-10 中国空间技术研究院 Method for matching navigation adaptation area by underwater gravity
CN109141426A (en) * 2018-08-10 2019-01-04 中国空间技术研究院 A kind of preferred method in subaqueous gravity matching navigation adaptation area
CN109085655A (en) * 2018-09-19 2018-12-25 中国船舶重工集团公司第七0七研究所 A kind of underwater platform gravity measurement scheme and verification method
CN109085655B (en) * 2018-09-19 2020-06-19 中国船舶重工集团公司第七0七研究所 Underwater platform gravity measurement scheme and verification method
CN109059964A (en) * 2018-09-19 2018-12-21 中国船舶重工集团公司第七0七研究所 A kind of inertial navigation based on gravity peak and the double calibration methods of gravity measurement
CN110031000A (en) * 2019-05-21 2019-07-19 北京理工大学 A kind of evaluation method of Method in Gravity Aided INS region suitability
CN110196050A (en) * 2019-05-29 2019-09-03 哈尔滨工程大学 A kind of vertical height of Strapdown Inertial Navigation System and speed measurement method
CN110196050B (en) * 2019-05-29 2022-11-18 哈尔滨工程大学 Vertical height and speed measuring method of strapdown inertial navigation system
CN111006694A (en) * 2019-12-29 2020-04-14 北京理工大学 Design method of long-endurance inertial navigation system track generator based on track planning
US11268813B2 (en) 2020-01-13 2022-03-08 Honeywell International Inc. Integrated inertial gravitational anomaly navigation system
CN111473790A (en) * 2020-04-26 2020-07-31 中国人民解放军61540部队 Submarine navigation method and system of gravity beacon along track
CN111426313A (en) * 2020-04-26 2020-07-17 中国人民解放军61540部队 Submarine navigation method and system based on gravity beacon
CN111721300A (en) * 2020-06-30 2020-09-29 清华大学 Gravity beacon navigation method and system based on artificial intelligence algorithm
CN112541046A (en) * 2020-11-30 2021-03-23 中国电子科技集团公司第二十八研究所 Co-occurrence target monitoring method based on time and space
CN112729275A (en) * 2021-01-08 2021-04-30 中国船舶重工集团公司第七0七研究所 Satellite inversion chart gravity adaptation area selection method utilizing factor analysis
CN113503871A (en) * 2021-05-19 2021-10-15 北京理工大学 Gravity matching method based on correlation filtering
CN113514059A (en) * 2021-05-19 2021-10-19 北京理工大学 Gravity-assisted inertial navigation system simulation platform
CN113514059B (en) * 2021-05-19 2024-02-13 北京理工大学 Gravity-assisted inertial navigation system simulation platform
CN113503871B (en) * 2021-05-19 2024-04-16 北京理工大学 Gravity matching method based on correlation filtering
CN113587923A (en) * 2021-05-31 2021-11-02 中国人民解放军61540部队 Submersible vehicle positioning method and system for screening multi-dimensional gravity gradient lighthouse matching area
CN113587923B (en) * 2021-05-31 2024-04-26 中国人民解放军61540部队 Submersible positioning method and system for screening matching areas of multidimensional gravity gradient lighthouse
CN114136320A (en) * 2021-11-19 2022-03-04 中国船舶重工集团公司第七0七研究所 Ocean multi-physical-field parameter feature positioning fusion method based on feature complementary characteristics

Also Published As

Publication number Publication date
CN108225310B (en) 2020-08-25

Similar Documents

Publication Publication Date Title
CN108225310A (en) A kind of Gravity-aided navigation path planning method
CN104197927B (en) Submerged structure detects robot real-time navigation system and method
US9551561B2 (en) Determining location using magnetic fields from AC power lines
Chutia et al. A review of underwater robotics, navigation, sensing techniques and applications
Gade NavLab, a generic simulation and post-processing tool for navigation
CN109556632A (en) INS/GNSS/polarization/geomagnetic integrated navigation alignment method based on Kalman filtering
CN109443379A (en) A kind of underwater anti-shake dynamic alignment methods of the SINS/DVL of deep-sea submariner device
CN103743395B (en) The compensation method of time delay in a kind of inertia-gravity coupling integrated navigation system
CN104390646B (en) The location matching method of underwater hiding-machine terrain aided inertial navigation system
CN108225308A (en) A kind of attitude algorithm method of the expanded Kalman filtration algorithm based on quaternary number
CN106979780B (en) A kind of unmanned vehicle real-time attitude measurement method
Song et al. Neural-network-based AUV navigation for fast-changing environments
CN105242682B (en) Target drone target signature measurement system
CN107990910A (en) A kind of naval vessel Large azimuth angle Transfer Alignment based on volume Kalman filtering
CN105043415A (en) Inertial system self-aligning method based on quaternion model
CN106979778A (en) A kind of localization method, device and mobile terminal
CN106017460B (en) A kind of underwater hiding-machine navigation locating method of terrain aided inertial navigation tight integration
CN106840211A (en) A kind of SINS Initial Alignment of Large Azimuth Misalignment On methods based on KF and STUPF combined filters
CN109596144A (en) Initial Alignment Method between GNSS location assists SINS to advance
CN102116634A (en) Autonomous dimensionality reduction navigation method for deep sky object (DSO) landing detector
CN107576327A (en) Varistructure integrated navigation system design method based on Observable degree analysis of Beidou double
CN109540135A (en) The method and device that the detection of paddy field tractor pose and yaw angle are extracted
CN107525502A (en) A kind of method for improving submarine navigation device inertia terrain match navigation mean accuracy
CN105988129A (en) Scalar-estimation-algorithm-based INS/GNSS combined navigation method
US20140249750A1 (en) Navigational and location determination system

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