WO2009028929A1  Device and method for calculating position of mobile station  Google Patents
Device and method for calculating position of mobile station Download PDFInfo
 Publication number
 WO2009028929A1 WO2009028929A1 PCT/NL2007/050427 NL2007050427W WO2009028929A1 WO 2009028929 A1 WO2009028929 A1 WO 2009028929A1 NL 2007050427 W NL2007050427 W NL 2007050427W WO 2009028929 A1 WO2009028929 A1 WO 2009028929A1
 Authority
 WO
 WIPO (PCT)
 Prior art keywords
 observables
 integer
 solution
 mobile station
 test
 Prior art date
Links
 230000001133 acceleration Effects 0 description 1
 238000009825 accumulation Methods 0 description 27
 230000006399 behavior Effects 0 description 2
 238000005452 bending Methods 0 description 1
 238000004422 calculation algorithm Methods 0 description 1
 238000004364 calculation methods Methods 0 description 19
 239000000969 carrier Substances 0 abstract claims description 49
 238000004891 communication Methods 0 description 10
 238000004590 computer program Methods 0 claims description 8
 230000000875 corresponding Effects 0 description 4
 230000000593 degrading Effects 0 description 1
 230000001419 dependent Effects 0 description 1
 238000009826 distribution Methods 0 description 3
 230000000694 effects Effects 0 description 3
 239000000284 extracts Substances 0 description 1
 239000010410 layers Substances 0 description 2
 239000011133 lead Substances 0 description 3
 239000011159 matrix materials Substances 0 description 26
 230000015654 memory Effects 0 description 3
 238000000034 methods Methods 0 description 47
 239000000203 mixtures Substances 0 description 2
 230000004048 modification Effects 0 description 3
 238000006011 modification Methods 0 description 3
 230000000051 modifying Effects 0 claims description 2
 239000011295 pitch Substances 0 description 2
 238000003908 quality control methods Methods 0 description 2
 230000003068 static Effects 0 description 1
 238000000528 statistical tests Methods 0 claims description 13
 230000001131 transforming Effects 0 description 5
 239000011701 zinc Substances 0 description 1
Classifications

 G—PHYSICS
 G01—MEASURING; TESTING
 G01S—RADIO DIRECTIONFINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCEDETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
 G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
 G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
 G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting timestamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
 G01S19/42—Determining position
 G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
 G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Leastsquares AMBiguity Decorrelation Adjustment] method
Abstract
Description
DESCRIPTION
DEVICE AND METHOD FOR CALCULATING POSITION OF MOBILE STATION
TECHNICAL FIELD
The present invention relates to a device and a method for calculating the position of the mobile station.
BACKGROUND ART
Carrier phase ambiguity resolution is the key to fast and highprecision GNSS (Global Navigation Satellite System) positioning and navigation. It applies to a great variety of current and future models of GPS, modernized GPS and Galileo.
Ambiguity resolution is the process of resolving the unknown cycle ambiguities of double difference carrier phase data (or single difference carrier phase data) as integers . Various methods using various models are known or proposed to estimate the ambiguity. Integer ambiguity estimation using the leastsquares method, the Kalman filter, or the like is wellknown.
Despite the differences in application of the various methods, their ambiguity resolution problems are intrinsically the same. One of the problems is that they cannot provide the user or analyst with tools to evaluate the quality of the integer solution. It is of importance to be able to evaluate the quality of the integer solution, since unsuccessful ambiguity resolution, when remaining unnoticed, may lead to unacceptable errors in the positioning result.
Unsuccessful ambiguity resolution may be caused by various error factors such as outliers in code data, slips in carrier phase data (referred to as cycleslip) , or the like. Some of these error factors cannot be removed even when double difference carrier phase data are used. The error factors may also be hidden in the integer ambiguity estimating process. Therefore, it is of importance to be able to identify or narrow the error factors if it is likely that there are errors in the positioning result or lack of reliability in the observation data.
Papers titled "Quality Control in Integrated Navigation Systems" and ^{W}AN INTEGRITY AND QUALITY CONTROL PROCEDURE FOR USE IN MϋLTI SENSOR INTEGRATION" published in 1990s by P. J.G. Teunissen discloses methods of evaluating the quality of the ambiguity resolution.
The disclosed methods are directed at evaluating the quality of the ambiguity resolution conducted by using the recursive Kalman filter algorithm. Of course, it is of importance to be able to evaluate the quality of the ambiguity resolution in such a multiepoch application; however, the multiepoch application has a disadvantage in that the ambiguity resolution at each epoch is affected by observables obtained at the previous epochs. Thus, according to these methods, it is difficult to evaluate the quality of the observables (and thus the quality of the ambiguity resolution conducted by using the observables) on epoch basis.
In particular, in a single epoch application where the ambiguity resolution is conducted by using the observables obtained at a single epoch, it is of importance to be able to evaluate the quality of the integer solution on epoch basis, because in the single epoch application the quality of the integer solution may be significantly affected by observables obtained at each epoch.
DISCLOSURE OF INVENTION It is an object of the present invention to provide a device and a method for calculating the position of a mobile station which can appropriately evaluate reliability of observables on epoch basis.
In order to achieve the abovementioned objects, according to one aspect of the present invention a method for calculating the position of a mobile station using an integer ambiguity is provided, comprising; a testing unit configured to evaluate reliability of observables on epoch basis, said observables including at least code data and carrier phase obtained at a single epoch, wherein said testing unit includes a first testing unit configured to test residuals as a whole, and a second testing unit configured to test one or more items of the residuals separately.
In this aspect, the first testing unit may be configured to determine whether the sum of squares of the residuals is larger than a reference value, and the second testing unit is configured to perform the^{*} test . if the sum of squares of the residuals is larger than the reference value.
Preferably, the residuals are derived using the integer ambiguity which is estimated based on the observables . Preferably, the device further includes a third testing unit configured to evaluate reliability of the integer ambiguity estimated based on the observables.
Preferably, the device further a reliability outputting unit configured to output the reliability of the calculated position of the mobile station, said reliability being derived based on the test result.
Preferably, items of the observables to be used for estimating the integer ambiguity are selected according to the test result.
Preferably, the integer ambiguity is estimated using the observables obtained at a single epoch.
Preferably, the integer ambiguity is estimated after modifying a float solution according to the test result.
Preferably, if the test result at a certain epoch does not meet predetermined criteria, the integer ambiguity at the epoch is estimated using an item of observables obtained at another epoch.
According to another aspect of the present invention a method for calculating a position of a mobile station using an integer ambiguity is provided, comprising; a float solution testing unit configured to evaluate reliability of a float solution of the integer ambiguity derived based on observables which are obtained at a single epoch; and an integer solution testing unit configured to evaluate reliability of an integer solution of the integer ambiguity estimated based the float solution.
In this aspect, preferably, if the test result on the float solution doesn't meet predetermined criteria, the estimation of the integer solution based on said float solution is ceased. Preferably, the tests are performed separately on the float solution derived based on double difference observables and on the float solution derived based on single difference observables. Preferably, float solutions to be used for estimating the integer ambiguity are selected according to the test result.
Preferably, the tests are performed separately on the integer solution derived based on double difference observables and on the integer solution derived based on single difference observables.
Preferably, integer solutions to be used for calculating the position of a mobile station are selected according to the test result.
According to another aspect of the present invention, a method for calculating a position of a mobile station using an integer ambiguity is provided, comprising; a float solution testing step of evaluating reliability of a float solution of the integer ambiguity derived based on observables which are obtained at a single epoch; and an integer solution testing step of evaluating reliability of an integer solution of the integer ambiguity estimated based the float solution.
According to another aspect of the present invention, a computer program is provided for calculating a position of a mobile station using an integer ambiguity, said computer program causing a computer to execute the following step; a testing step for evaluating reliability of observables including at least code data and carrier phase data, said observables being obtained at a single epoch, wherein said testing step includes a first testing step in which residuals of variables with respect to the observables as a whole are evaluated, and a second testing step in which a statistical test for detecting at least one of error factors is performed.
According to another aspect of the present invention a computer program is provided for calculating the position of a mobile station using an integer ambiguity, said computer program causing a computer to execute the following step: a float solution testing step of evaluating reliability of a float solution of the integer ambiguity derived based on observables obtained at a single epoch; and an integer solution testing step of evaluating reliability of an integer solution of the integer ambiguity estimated based the float solution.
BRIEF DESCRIPTION OF THE DRAWINGS
These and other objects, features, and advantages of the present invention will become more apparent from the following detailed description of preferred embodiments given with reference to the accompanying drawings, in which:
Fig. 1 is a schematic diagram of a carrier phase GPS positioning device according to an embodiment of the present invention; Fig. 2 is a diagram showing a configuration of the carrier phase GPS positioning device in FIG. 1;
Fig. 3 is a block diagram showing an embodiment of a carrier phase GPS positioning device 34 installed in the mobile station 30 according to the present invention; Fig. 4 is a diagram illustrating the definitions of coordinate systems used in descriptions;
Fig. 5 is a flowchart illustrating the input test according to an embodiment of the present embodiment; Fig. 6 is a flowchart illustrating a method of giving reliability to the calculated position of the mobile station 30;
Fig. 7 is a flowchart illustrating a method of calculating the position of the mobile station 30 for the epoch when the flag 2 and/or the flag 1 is set to "1";
Fig. 8A shows the curve of the positions calculated according to the conventional single epoch application; Fig. 8B shows the curve of the positions calculated according to an embodiment of the present invention;
Fig. 9 is a flowchart illustrating a method of evaluating the quality of the calculated position of the mobile station 30 according to another embodiment of the present invention; and
Fig. 10 is a flowchart illustrating a method of improving the accuracy in the calculated position of the mobile station 30 based on the test results according to an embodiment of the present invention.
BEST MODE FOR CARRYING OUT THE INVENITON
Hereafter, the preferred embodiments according to the present invention are explained with reference to the drawings.
FIG. 1 is a schematic diagram of a carrier phase GPS positioning system according to an embodiment of the present invention.
As illustrated in FIG. 1, the carrier phase GPS positioning system includes GPS satellites 10 orbiting around the earth, a reference station 20 located at a fixed position (known position) , and a mobile station 30 that is on the earth and is able to move on the earth. Each of the GPS satellites 10 broadcasts navigation messages toward the earth continuously. The navigation messages include orbital information of the corresponding GPS satellite 10, clock correction values, and correction coefficients of the ionospheric layer. The navigation messages are spread using a PRN code, such as C/A code or P code. At present, the C/A code is carried on a Ll carrier (frequency: 1575.42 MHz) and the P codes are carried on the Ll carrier and a L2 carrier (frequency: 1227,6 MHz), and are broadcast toward the earth.
Presently, there are 24 GPS satellites orbiting the earth at an altitude of 20,000 km. Every four GPS satellites are equally arranged on one of six orbital planes of the earth, which are inclined 55 degrees relative to each other. Therefore, at least five satellites are always observable from a position as long as the position is open to the sky, no matter where the position is on the earth.
In FIG. 2, the mobile station 30 has a GPS receiver 32. In the GPS receiver 32, there is an oscillator (not illustrated) having an oscillating frequency equal to the carrier frequency of the GPS satellite 10. The GPS receiver 32 converts an electromagnetic wave, which is emitted from the GPS satellite 10 and is received by the GPS receiver 32 via a GPS antenna 32a, performs C/A code synchronization using the C/A codes generated in the GPS receiver 32, and extracts the navigation messages .
The GPS receiver 32 calculates a carrier phase accumulation value Φi_{U} of the carrier waves from the GPS satellites 1Oi. Here, in the phase accumulation value Φi_{u}, the subscript i(=l, 2, ...) represents the numbers assigned to the GPS satellite 1Oi, and the subscript u represents _Q—
that the accumulation value is calculated on the side of the mobile station 30.
The phase accumulation value Φχ_{u} can be described as the difference between a phase Θi_{U}(t) of the oscillator at the time t of receiving the carrier wave and a phase Θj_{.u}(tτ) of the carrier wave when the satellite signal from the GPS satellite 1Oi is generated, as shown by the following formula (1) .
Here, τ_{u} represents travel time from the GPS satellite 10 to the GPS receiver 32, and ε i_{U} represents noise (uncertainty) . Further, at the time when starting observing the phase difference, the GPS receiver 32 can accurately determine the carrier phase within one wavelength of the carrier wave, but cannot determine what number of the wavelengths the present wavelength is. For this reason, in the phase accumulation value Φj_{.u}(t), as shown in the formula (1), there is an uncertainty factor Ni_{U}, known as "integer ambiguity".
The phase accumulation value Φi_{U} can be calculated for both of .the Ll wave and the L2 wave. In this case, two phase accumulation values ~Φ_{LU} can be obtained at each epoch.
The GPS receiver 32 also calculates a pseudo range Pi_{U}(t) based on the C/A codes carried by the carrier waves from the GPS satellites 1Oi. The pseudo range p i_{u}(t) calculated here includes errors such as a range error, as shown by the following formula (11) .
Here, c denotes speed of light, and b, which is also referred to as a clock bias, corresponds to a range error due to a clock error in GPS receiver 32.
Similarly, the pseudo ranges Pi_{U}(t) can be measured using the P codes carried by the Ll wave and the L2 wave. In this case, two pseudo ranges Pi_{U}(t) based on the P codes can be obtained at each epoch.
The mobile station 30 also includes a communication device 33, such as a mobile phone. As described below, the communication device 33 is capable of communicating with a communication facility 23 installed on the reference station 20 side, such as a base station for mobile phones, by bidirectional communication.
A GPS receiver 22 having a GPS antenna 22a is installed in the reference station 20. The GPS receiver 22, the same as the GPS receiver 32 in the mobile station 30, calculates a carrier phase accumulation value Φi_{b} at time t based on the carrier waves from the GPS satellites 10i, as shown by the following formula (2) .
Here, Nu_{3} is an integer ambiguity, and e _{H3} represents noise (uncertainty) . In the ptiase accumulation value Φib, the subscript b represents that the accumulation value is calculated on the side of the reference station 20.
Similarly, the phase accumulation value Φ_{ϋ} can be measured using both of the Ll wave and the L2 wave. In this case, two phase accumulation values Φjjj can be obtained at each epoch.
The GPS receiver 22 also calculates a pseudo range pa_{>}{t) based on the C/A codes carried by the carrier waves from the GPS satellites 10i, as shown by the following formula (21) .
Similarly, the pseudo range p^ (t) can be measured using the P codes carried by the Ll wave and the
L2 wave. In this case, two the pseudo ranges Pϋ_{>} (t) based on the P codes can be obtained at each epoch. The reference station 20 transmits observation data including the obtained carrier phase accumulation value Φi_{b} and the pseudo range p a, to the mobile station 30 via the communication facility 23. More than one reference station 20 may be installed in a specified region. As illustrated in FIG. 2, each of the reference stations 20 may be connected to one or more communication facilities 23 through the Internet or other networks, or a communication facility 23 may be installed in each of the reference stations 20. In the former case, as long as the mobile station 30 is able to communicate with the communication facility 23, the mobile station 30 can obtain the information received by each of the reference stations 20.
FIG. 3 is a block diagram showing an embodiment of a carrier phase GPS positioning device 34 installed in the mobile station 30 according to an embodiment of the present invention.
The carrier phase GPS positioning device 34 of the present embodiment includes a calculation unit 40, which is connected to the GPS receiver 32 and the communication device 33, and further, to various sensors 50 provided in the mobile station 30. The calculation unit 40 may also be installed in the GPS receiver 32. When the mobile station is a vehicle, the GPS receiver 32, the calculation unit 40 and/or the communication device 33 may also be mounted in a navigation device.
The calculation unit 40 consists mainly of a 5 microcomputer. A microcomputer consists of a CPU for processing data, memory in which computer programs and data for implementing the process described below or processed data are stored, interfaces, and the like.
The calculation unit 40, based on the orbital
10 information in the navigation messages received by the GPS receiver 32, calculates positions (X±(t), Yi (t) , Zi (t) ) of all observable GPS satellites 1Oi at time t in a global coordinate system.
FIG. 4 is a diagram illustrating the
15 definitions of coordinate systems used in the following description. FIG. 4 shows relationships between the global coordinate system, a local coordinate system, and a body coordinate system.
The body coordinate system is defined on the
20 body of the vehicle.
Because movement of each of the GPS satellites 10 is confined to an orbital plane passing through the center of gravity of the earth, and the orbit of each of the GPS satellites 10 is an ellipse with the center of
25 gravity of the earth as a focus, positions of each of the GPS satellites 10 in the orbital plane can be calculated by successive numerical solutions of Kepler's equation. Because the orbital planes of each of the GPS satellites and the equatorial plane in the global coordinate system
30. satisfy a rotational transformation relationship, positions (Xi(t), Yi(t), Zi(t)) of the GPS satellites 10 at the time t of receiving the carrier waves can be calculated by three dimensional rotational coordinate transformation of the positions of the GPS satellites 10 on the orbital planes.
The calculation unit 40, based on the output signals of the various sensors 50 input periodically, calculates quantities related to movement of the mobile station 30.
For example, if the mobile station 30 is a vehicle, the calculation unit 40 calculates the speed Vx (t) (speed in the forward and backward directions) and Vy (t) (speed in the right and left directions) at the time t of receiving the carrier wave based on output signals from various sensors 50, for example, two wheel speed sensors mounted on the driven wheels of the vehicle, a steering sensor, a yaw rate sensor, a left and right G acceleration sensor, and an azimuth meter.
Because the speed vector (Vx (t) , Vy (t) ) of the vehicle is defined in the body coordinate system whose origin is on the body of the vehicle, it is necessary for the calculation unit 40 to transform the speed vector (Vx (t), Vy (t) ) from the body coordinate system to the global coordinate system via the local coordinate system. Usually, the rotational transformation of coordinates can be performed by using Euler angles. In the present embodiment, the transformation from the body coordinate system to the local coordinate system is performed using only a yaw angle φ (t) since the roll angle and pitch angle are small. Depending on the situation, the roll angle and pitch angle may also be considered, or the yaw angle may also be ignored. The transformation from the local coordinate system to the global coordinate system is performed by using the longitude φ (t) and latitude λ (t) of the position of the vehicle.
Such dynamic information indicating dynamic behavior of the mobile station 30 may be used to generate or modify models for calculating the position of the mobile station 30, or may be used to modify the calculated position of the mobile station 30. The calculation unit 40 performs an input test on the observation data obtained at the mobile station 30 and the reference station 20, or other sensor data from the various sensors 50. The observation data include the carrier phase accumulation values Φi_{U}, Φϋ_{>}/ and the pseudo ranges p i_{U}, p ϋ.
The input test is performed for checking reliability of the observation data, as described below. The input test is performed in two stages. The first stage is for diagnosing whether an unspecified model error has occurred. If the error indeed present, then the second stage is initiated for detecting a potential source (error factor) of the model error detected in the first stage .
Fig. 5 is a flowchart illustrating the input test according to an embodiment of the present embodiment. The process routine shown in Fig. 5 is performed at each epoch independently. The process routine shown in Fig. 5 is executed by the calculation unit 40.
In step SlOO, a residual vector v consisting of residuals is derived using an appropriate model. The various models known or proposed are applicable. In this embodiment, the simple model referred to as "geometryfree model" and/or "geometrybased model" is used. The geometryfree model is given as follows.
Or the
follows.Here, Φ (t) is an observable related to the carrier phase accumulation value, and p (t) is an observable related to the pseudo range. D_{p}(t), D_{φ} (t) are each an unknown quantity related to a stationsatellite range. N is an unknown quantity related to the integer ambiguity. λ denotes a wavelength of the carrier wave. n_{p} (t) and n_{φ}(t) each denotes an observation error, including the clock error in the case of the former equation.
According to this embodiment, since the model is established by using the pseudo ranges pi_{u}, p_{ϋ},in addition to the carrier phase accumulation values Φi_{U/} Φ i_{b} , it becomes possible to have the redundancy that is required to accurately derive the optimal solutions, even if the number of the GPS satellites 10 that is being tracked is relatively small. In other words, it becomes possible to perform the integer resolution on epoch basis, even if the number of the GPS satellites 10 that is being tracked is relatively small.
In one embodiment, the calculation unit 40 calculates^{"}double differences of the observation data, that is, double differences of the carrier phase accumulation values and the pseudo ranges.
The double difference of the carrier phase accumulation values related to the GPS satellites 10_{j} and 1O_{h} (j is not equal to h) at time t (a certain epoch) can be expressed by the following formula (3) .
The double difference of the pseudo ranges related to the GPS satellites 1Oj and 1Oh (j is not equal to h) at time t (a certain epoch) can be expressed by the following formula (4) .
In this embodiment, the phase accumulation value Φ_{jMw} of the carrier phase accumulation values is used as observables Φ (t) in the model, and the double difference P_{jhbu}Of pseudo ranges is used as observables p (t) in the model .
In this embodiment, the leastsquares method is applied to derive the residual vector v . The least squar method is iven as follows.
Here, y denotes a vector of observables consisting of Φjhbu and Pjhbu_{f} x denotes a vector of unknown quantities consisting of D_{p} (t) , D_{φ} (t) , and N, n denotes an observation noise including errors which cannot be represented even by approximations, and W denotes an design matrix. The design matrix W is given as follows.
Here, W is a combination of an design matrix WX for the observables Φjhbu and an design matrix Wl for the observables Pjhbu The design matrix W may be adapted by the solutions of the variables x in the case of the design matrix W being dependent on the variables x. In the case of double difference the same 5 satellites being observable at the mobile station 30 and the reference station 20 and observation data being generated from the Ll wave and the L2 wave, the number of the observables is 12, and the number of the unknown quantities is 11.
As is well known, the optimum solution x of x is given as follows.
As a result, the residual vector v is given as follows;
that in order to derive the residual vector v the single differences of the observation data or the observation data itself can be used as observables. In the case of the double phase difference, it is possible to eliminate influence of the initial phase of oscillators in the GPS receivers 22 and 32, and clock uncertainties. In the case of the single phase difference, it is possible to eliminate influence of the initial phase of oscillators in the GPS satellite 10, and the GPS clock error.
It is also noted that in order to derive the residual vector v , the combination of the single differences of the observation data and the double differences of the observation data can be used as observables.
In single difference application, the single difference of the carrier phase accumulation values related to the GPS satellites 1Oj at time t can be expressed by the following formula (5) .
The single difference of the pseudo ranges related to the GPS satellites 1Oj at time t can be expressed by the following formula (6) .
It is also noted that in order to derive the residual vector v , another model can be used. In another embodiment, a geometrybased model is used. In this embodiment, the distance between the GPS satellite 1Oi and the GPS receiver 22 or 32 equals the wavelength λ of the carrier wave multiplied by the phase accumulation value, and the double difference Φjh_{bu} of the phase accumulation satisfies the following formula (7) .
In formula (7), [Xj_{3} (t) , Y_{b}(t), Z_{b}(t)] are coordinates (known) of the reference station 20 at time t in the global coordinate system, [X_{u}(t), Y_{u}(t), Z_{u}(t)] are coordinates (unknown) of the mobile station 30 at time t, and [X,(t)_{f} Yj (t), Zj (t)] and [X_{1}, (t), Y_{h}(t), Z_{h}(t)] are coordinates of the GPS satellites 10_{j} and 1O_{h} at time t calculated by the calculation unit 40. Njht_{>u} represents the double difference of the integer ambiguity, that is, Nj_{h}__{>u}
he double difference Pj_{h}buOf the phase accumulation satisfies the following formula (8) .
In formula (8), [X_{b}(t), Y_{b}(t), Z_{b}(t)], [X_{u}(t), Yu(t), Z_{u}(t)], [Xj (t), Yj (t), Zj (t)] and [X_{h}(t), Y_{h}(t), Z_{h} (t) ] are the same as mentioned previously. bj_{hbU} represents the double difference of the clock biases. Similarly, in this embodiment, the least squares method is applied to derive the residual vector v The leastsquares method is given as follows.
Here, the vector of observables y consists of the double difference Φj_{ht>}u of the phase accumulation
(refer to formula (7)) and the double difference p ^_{u} of pseudo ranges (refer to formula (8)). x denotes a vector of unknown quantities consisting of variables X_{11}, Y_{u}, Z_{u}, and the double phase difference Njh_{b}u of the integer ambiguity.
Since the vector of observables y is nonlinear relative to the variables X_{u}, Y_{u}, Z_{n}, the items in the formulas (7) and (8) are partially differentiated (linearized) relative to X_{u}, Y_{u}, Z_{u}. Sirtiilarly, in the embodiment, the optimum solution x of x and the residual vector v are given as follows .
In step SIlO, a test statistic T_{1n} for testing the overall validity of observables is derived and model errors are checked using the test statistic T_{1n} . The test statistic T_{01} , which is the sum of squares of the residuals, is given as follows.
Here, Q_{v} denotes an error matrix (i.e., a variancecovariance matrix of residuals) . m denotes the number of unknown quantities. c_{ip} (i =l,....,<gr, ) is a vector for identifying an expected error factor e_{t} ( /=1,....,#, ) that could contribute to the test statistic T_{m} .
In the case where the observables consist of the pseudo ranges based on the C/A code, the P code in the Ll wave, the P code in the L2 wave, and the carrier phase accumulation values based on the Ll wave and the L2 wave, and the geometryfree model is used, the variable vector x for a minimum unit (i.e., one pair of the GPS satellites 10j and 1Oh ) can be represented as follows.
Then, residual vector v is given as follows.
he residual vector v used for calculating is modified by omitting the items associated with integer ambiguity as follows. e
Here, represents a vector for identifying an error in Φ_{u}, that is, error in the carrier phase data of the Ll wave, such as Φ (t) for the
Ll wave . epresents a vector for
identifying an error in Φ^ , that is, error in the carrier phase data of the L2 wave, such as Φ (t) for the L2 wave.represents a vector for identifying
errors that occur simultaneously in Φ_{Ll} and Φ_{L2}.represents a vector for
identifying an error proportional to the ratio between the frequency Ll and the frequency L2, such as an ionosphericerror. represents a vector for
identifying an error proportional to a ratio 9/7.represents a vector for identifying an
error in code data (pseudo ranges based on the C/A code) .In step S115, a test is performed on the test statistic T_{m} . In one embodiment, the test for testing the overall validity of the hypothesis is performed as follows .
Here, χ^{2} _{a}(b,0) is the upper a probability point of the central % ^{2}distribution with b degrees of freedom.
If T_{m}≥χl(b,0) , it is concluded that unspecified model errors have occurred, and then the test process of Fig. 5 proceeds to step 120. Otherwise it is concluded that no unspecified model errors have occurred and therefore the test process finishes.
In step S120, p indicating the number of repetitions is incremented as follows. p = p +l.
In step S130, a test statistic t_{lp} is derived for each of the expected error factors by applying each corresponding vector C_{1} , and the maximum test statistic
If_{1}J is selected as follows.
In step S140, a test is performed on the maximum
test statistic W_{j} A ■ In one embodiment, the test is
performed using a normal distribution with variance 0 and mean 1 as follows.
est
process of Fig. 5 finishes, concluding that all errors have been found. In this case, the vector c_{Jp} related tothe maximum test statistic f_{Λj}_ selected in step 130 is
considered to indicate the error factor. In other words, it is concluded that the vector c_{Jp} related to the maximum
test statistic k J indicates the error factor (error
type) that has contributed to a deterioration of the test statistic T_{1n} . Otherwise the test process proceeds to step 150.
In step S150, the redundancy is checked, as follows .
If m≤p, the test process of Fig. 5 finishes, concluding that there have been too many errors in the observation data to be identified. Otherwise the test process proceeds to step 160. In step S160, in order to eliminate or reduce
the effect of the maximum test statistic f,J on the test
statistic T_{n} , the test statistic T_{1n} is modified as follows .
In step S170, the test is performed again on the modified test statistic T^{*} _{m} . In one embodiment, the test is performed as follows.
Here, F_{a}{m,∞,ϋ) is the upperα probability point of the central Fdistribution with m, ∞ degrees of freedom. If T'm≥T Thr, it is concluded that there may be error factors other than the error factor specified by the vector C_{j p}, and the test process of Fig. 5 proceeds to step 180. Otherwise it becomes clear that the error factor specified by the vector c_{Jp} is the last error factor to be found, and thus the test process finishes, concluding that all errors have been found.
In step S180, in order to search for the error factors in the different direction, the vector c_{ip} is orthogonalized to be used for the next routine, as follows;
Then, the test process of Fig. 5 returns to step
120 with the orthogonalized vector c_{ip} , and continues recursively until it is concluded that all errors have been found or there have been too many errors.
According to the embodiment, it becomes possible to evaluate the quality of the observation data as a whole as well as each of the error factors separately.
It is noted that if the four GPS satellites can be tracked, the aforementioned tests may be performed on the observation data for three pairs of the GPS satellites (if one of GPS satellites is regarded as a reference satellite) . In this case, it is also possible to perform the aforementioned tests simultaneously for 3 pairs of the GPS satellites by integrating the corresponding vectors and the matrixes. In this embodiment, it is also advantageous to modify the float solution and/or the variance according to the test result. The modification is implemented by using matrixes Δ , Q_{A} , and K . The matrixes Δ , β_{Δ} , and K are given as follows.
Here, c_{J} ^{T} _{p} is the vector that maximizes the test
statistic W_{j} A (see step 130 in Fig. 5) . Q_{x} is a variance
covariance matrix of the estimated x, and is represented as follows.
If the error factor specified by cj_{p} is outliers
in the observation data (typically, in the code data) , the optimum solution (float solution) Jc and the variance covariance matrix of the observables Q_{v} are modified by being derived as follows.
On the other hand, if the error factor specified
by c^_{p} is slip in the carrier phase data (i.e. cycleslip),
x and Q_{v} are modified by being derived as follows .
Here, [Δ] represents roundedoff components of the matrix Δ .Next, referring to Fig. 6, a method of giving reliability representative of the quality of the calculated position of the mobile station 30 to the output result is described. The process routine shown in Fig. 6 is performed at each epoch independently. The process routine shown in Fig. 6 is executed by the calculation unit 40.
In step 200, two flags (flag 1 and flag 2) are provided. Both flags have initial values ^{M}0" indicating that the quality of the calculated position of the mobile station 30 meets predetermined criteria.
In step 210, the data to be used to calculate the position of the mobile station 30 are received (input) The input data may include the observation data, such as carrier phase accumulation values and pseudo ranges, and the dynamic information indicating dynamic behavior of the mobile station 30 derived based on the output signals from the various sensors 50. In step 220, the test is performed on the input data. The test is performed according to the method described in Fig. 5.
In step 220, if the test result doesn't meet predetermined criteria, the flag 1 is set to Λ^{Λ}1" indicating that there may be an error in the input data that could lead to an unsuccessful ambiguity resolution. For example, if in the step 150 (see Fig. 5) it is concluded that there have been too many errors in the observation data, the flag 1 is set to "1". Otherwise, the process proceeds to step 230.
In step 230, the ambiguity resolution is performed. The ambiguity resolution process includes the process to derive the float solution and the process to estimate the integer solution (also referred to as ^{Λ}fixed solution' ) based on the float solution. In the process to derive the .float solution, a leastsquares method or a modified Kalman filter (described later) or the like can be used. In the process to estimate the integer ambiguity, the LAMBDA, for example, can be used, which decorrelates the integer ambiguities, and narrows the searching space of the integer solutions so as to facilitate finding the integer solution.
In step 240, an output test is performed on the derived float solution and/or the estimated integer ambiguity. The output test may be performed using conditional variance, a ratio test, or the like. The probability of estimating the correct integers is given as follows .
Here, σ_{v} denotes diagonal components (variance components) of the error matrix Q_{v} .
The ratio test is given as follows.
Here, Q_{a}, denotes an appropriate weight matrix, a denotes a subset of the float solution of the variables x, a_{s} denotes the integer solution of the variables x with the highest probability and 5_{S2} denotes the integer solution with the second highest probability. The optimum solution Jc or the modified solution Jc^{*} can be used as a. It is noted that other ways of conducting the output test are also applicable. For example, the output test using the innovation or residual vector v and/or the covariance matrix Q_{x} is applicable.
In step 240, if the test result of the output test doesn't meet the predetermined criteria, the flag 2 is set to "1" indicating that there may be an error in the estimated integer ambiguity that could lead to incorrect position determination of the mobile station 30. For
example, the flag 2 is set to
W η //
J .
After the integer ambiguity is resolved, the position (coordinates) of the mobile station 30 can be determined using various differential positioning methods, for example.
In step 250, the calculated position of the mobile station 30 together with the flags are output as a result. The position of the mobile station 30 obtained in this way may be used in various controls or be presented as information, for example, it may be output and displayed on a screen of a navigation device, or displayed in a map shown on the screen of a mobile phone. In this way, it becomes possible for a user or analyst to check the quality of the calculated position by checking the value of the flags. For example, if the flag 1 is set to "0" and the flag 2 is set to "1", user or analyst may understand there is a possibility of errors in the calculated position as well as that such error would result from the ambiguity resolution process. It is noted that more flags can be used to indicate the type of error that has been found by the vector c_{ip} , for example.
Next, referring to Fig. 7, one embodiment of a method of calculating the position of the mobile station 30 at the epoch in which the flag 2 and the flag 1 are set to "1" is described. The process routine shown in Fig. 7 is executed by the calculation unit 40. The process routine shown in Fig. 7 may be performed before the step 250 in Fig. 6 in the case where the flag 2 and the flag 1 are set to "1". It is noted that the process routine shown in Fig. 7 may be performed in the case where either the flag 2 or the flag 1 is set to "1".
In step 300, Ar—Ϊ < 0 is checked. Here, k is the current epoch number and J = O_{5}I,...., (initial value of i is 0) . If k~i>0, the process proceeds to step 310. Otherwise the process finishes, setting the flag 1 to "ki(=0)".
This indicates that the process routine has ended without any modification.
In step 310, the data having been obtained at the epoch ki is read out from the memory. In other words, the previous epoch data are read out from the memory.
In step 320, the ambiguity resolution is performed based on the data obtained from epoch ki to epoch k. In this case, the ambiguity resolution may be performed by the observation data obtained over more than 2 epochs. Therefore, the recursive Kalman filter can be applied to the ambiguity resolution. If applying the Kalman filter, the following equations can be obtained. For updating of the epoch,
For updating of the observables,
Here, Q and^{'} R represent the covariance matrix of the external noise and the covariance matrix of the observation noise, respectively. Ki denotes Kalman gain (corrective gain). P(k)^{(+>} and P(k)^{M} are the covariance matrix of the expected errors, and the covariance matrix of the estimated errors, respectively. The formulae (9) and (12) are filter equations, and the formulae (10) and (13) are covariance equations. Here, the superscript ^{("}' and ^{!+!} indicate time before and after the updating, respectively. In the static model, the item U(kl) derived from the dynamic information in formula (9) does not exist. In the case of the modified Kalman filter for the single epoch application, the P(k) is initialized on epoch basis and thus the covariance is not took^{*} over for^{^}the next epoch. It is noted that the ambiguity resolution may be performed by the leastsquare method using observation data obtained at more than two epochs .
In step 330, the statistical test is performed on the float solution and/or the integer ambiguity which are derived or estimated based on the observation data obtained over more than 2 epochs. The method of the statistical test may be the same as described with reference to step 240' in Fig. 6. The test using the innovation (y(k) W(k) *x(k) ^{l~}}) or residual vector v is also applicable. Further, the input test as described with reference to Fig. 5 may be performed on the observation data obtained over more than 2 epochs. If the result of the statistical test doesn't meet the predetermined criteria, the process returns to step 300 and continues recursively until the result of the statistical test meets the predetermined criteria, or ki≤Q . In step 340, the calculated position of the mobile station 30 together with the flags are output as a result. In this way, it becomes possible for the user or analyst to check the quality of the calculated position by checking the value of the flags. According to the embodiment, while the ambiguity resolution may be performed on epoch basis, the ambiguity resolution on a multiepoch basis is performed in exceptional cases where the flag 2 and/or the flag 1 is set to "1". Therefore, it is possible to output the calculated position of the mobile station 30 with high accuracy without degrading the advantage of the single epoch application.
It is noted that it is also possible to use more than two GPS antennas 32a, which are preferably provided at different places on the mobile station 30, in order to have more redundancy to ensure reliability.
Fig. 8A graphs the positions of the mobile station 30 at the plural epochs determined according to the conventional single epoch application. Fig. 8B graphs the positions at the plural epochs determined according to the embodiment shown in Fig. 7 of the present invention. The same observation data were used for both applications. In the conventional single epoch application, at 2 epochs, there are hops in the determined position with respect to the neighboring epochs, as shown in Fig. 8A. To the contrary, according to the embodiment shown in Fig, 7 of the present invention, there are no such hops in the determined position, as shown in Fig. 8B. The hop resulting from the epoch kl on the left side in Fig. 8A is avoided by using the observation data obtained over two epochs (epochs kl and kl1) . The hop resulting from the epoch k2 on the right side in Fig. 8A is avoided by using the observation data obtained over 6 epochs
(epochs k2 and kl6) . It is noted that in Fig. 8B the positions of the mobile station 30 at other epochs were determined by the observation data obtained at the single corresponding epoch. Next, referring to Fig. 9, a method of evaluating the quality of the calculated position of the mobile station 30 according to another embodiment of the present invention is described. The process routine shown in Fig. 9 is performed at each epoch independently. The process routine shown in Fig. 9 is executed by the calculation unit 40.
In step 400', two flags (flag 3 and flag 4) are provided. The flag 3 indicates reliability of the derived float solutions derived. The flag 4 indicates reliability of the derived fixed solutions. Both flags have initial values "0" indicating that the reliability is relatively high.
In step 410, the float solutions are derived based on the observation data. The float solutions can be derived using the leastsquares method or the like.
In step 420, the innovation or residual vector v and/or the covariance matrix Q_{x} are derived using the derived float solutions. In step 430, the statistical test is performed on the derived float solutions by checking the innovation or residual vector v and/or the covariance matrix Q_{x} . The statistical test according to the Fig. 5 may be performed by using the residual vector v . The statistical test may be performed by checking the conditional variance.
If the test result doesn't meet the predetermined criteria, the flag 3 is set to "1" indicating that there may be an error in the derived float solutions. In this case, the derived float solutions are not reliable enough to be used for the integer estimating process (step 440) . Thus, the process at this epoch finishes by skipping the integer estimating process. Since the integer estimating process requires a relatively heavy calculation work load, it is useful to check the quality before the integer estimating process. It is noted that it is also possible to rederive the float solutions using the previous epoch data, as is described referring to Fig. 7. On the other hand, if the test result meets the predetermined criteria, the process proceeds to step 440.
In step 440, the fixed solutions are estimated from the derived float solutions. The fixed solutions can be estimated using the LAMBDA or the like.
In step 450, the innovation or residual vector v and/or the covariance matrix Q_{x} is derived using the estimated fixed solutions.
In step 460, the statistical test is performed on the estimated fixed solutions by checking innovation or residual vector v and/or the covariance matrix Q_{x} . The statistical test according to the Fig. 5 may be performed by using the residual vector v . The statistical test may be performed by checking the conditional variance and/or by using the ratio test as described previously.
If the test result doesn't meet the predetermined criteria, the flag 4 is set to "1" indicating that there may be an error in the estimated fixed solutions, and the process proceeds to step 470. If the test result meets the predetermined criteria, the process proceeds to step 470 without changing the values of the flags. In step 470, the calculated position of the mobile station 30 together with the flags are output as a result. In this way, it becomes possible for a user or analyst to check the quality of the calculated position by checking the value of the flags. Further, if the flag 3 is "1", the user or analyst can understand that the integer estimation is skipped because of the lack of reliability in the solutions.
According to the embodiment, since the reliability of the estimated fixed solutions is checked in addition to the reliability of the derived float solutions, it becomes possible to detect the errors that could occur in the integer estimating process.
Next, referring to Fig. 10, a method of improving the accuracy in the calculated position of the mobile station 30 based on the test results according to an embodiment of the present invention is described. The process routine shown in Fig. 10 is performed at each epoch independently. The process routine shown in Fig. 10 is executed by the calculation unit 40. In step 500, the ambiguity resolution is performed using the single difference (SD) observation data (see the formulas (5) and (6) ) .
In step 510, the innovation or residual vector v and/or the covariance matrix Q_{x} are derived using the fixed solutions and/or the float solution derived at step 500, as is described referring to Fig. 9.
In step 520, the ambiguity resolution is performed using the double difference (DD) observation data (see the formulas (3) and (4)).
In step 530, the innovation or residual vector v and/or the covariance matrix Q_{x} are derived using the fixed solutions and/or the float solution derived at step 520, as is described referring to Fig. 9.
In step 540, the reliability of the estimated fixed solutions derived at step 500 is evaluated by checking innovation or residual vector v and/or the covariance matrix Q_{1} derived at 510. Similarly, the reliability of the estimated fixed solutions derived at step 520 is evaluated by checking innovation or residual vector v and/or the covariance matrix Q_{x} derived at step 530. Then, the reliabilities for both of the ambiguity resolutions are compared. In step 550, the more reliable estimated fixed solutions are selected to calculate the position of the mobile station 30. Then, the calculated position of the mobile station 30 is output as a result. In this way, it becomes possible to have more redundancy to ensure reliability and to appropriately select a reliable output by using the test result.
In this embodiment, it is also possible to derive the conditional variance or the ratio used in the ratio test for each of the ambiguity resolutions and then compare them so as to output the more reliable estimated fixed solutions.
In this embodiment, it is also possible to compare the reliability of the float solutions based on single difference observation data with the reliability of the float solutions based on double difference observation data before the integer estimating process. In this case, the more reliable float solutions are used to estimate the fixed solutions. Since the integer estimating process requires a relatively heavy calculation work load, it is useful to compare the reliability levels of the float solutions before the integer estimating process. It is noted that it is also possible to calculate the attitude or the azimuth angle of the mobile station 30 based on the observation data, as is the case with the position of the mobile station 30. The above described embodiments can be applicable in such an application. The present invention is disclosed with reference to the preferred embodiments. However, it should be understood that the present invention is not limited to the abovedescribed embodiments, and variations and modifications may be made without departing from the scope of the present invention.
For example, in the above embodiments, the influence of the ionospheric layer refraction effect, tropospheric bending effect, and multipath are not considered^{"} in the model, but it^{*} is also possible to modόl these error factors.
In the above embodiment, a vehicle is illustrated as an example of the mobile station 30. The mobile station 30 may also include a folk lift or a robot with the receiver 32 and the calculation unit 40, for example.
Claims
Priority Applications (1)
Application Number  Priority Date  Filing Date  Title 

PCT/NL2007/050427 WO2009028929A1 (en)  20070829  20070829  Device and method for calculating position of mobile station 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

PCT/NL2007/050427 WO2009028929A1 (en)  20070829  20070829  Device and method for calculating position of mobile station 
Publications (1)
Publication Number  Publication Date 

WO2009028929A1 true WO2009028929A1 (en)  20090305 
Family
ID=39446121
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

PCT/NL2007/050427 WO2009028929A1 (en)  20070829  20070829  Device and method for calculating position of mobile station 
Country Status (1)
Country  Link 

WO (1)  WO2009028929A1 (en) 
Cited By (3)
Publication number  Priority date  Publication date  Assignee  Title 

CN102353969A (en) *  20110902  20120215  东南大学  Method for estimating phase deviation in precise singlepoint positioning technology 
WO2015145718A1 (en) *  20140328  20151001  三菱電機株式会社  Positioning device 
US10393879B2 (en)  20140328  20190827  Mitsubishi Electric Corporation  Global positioning device 
Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

US6127968A (en) *  19980128  20001003  Trimble Navigation Limited  Onthefly RTK positioning system with single frequency receiver 
US6313788B1 (en) *  19980814  20011106  Seagull Technology, Inc.  Method and apparatus for reliable interantenna baseline determination 
US7148843B2 (en) *  20030702  20061212  Thales North America, Inc.  Enhanced real time kinematics determination method and apparatus 
US20070120733A1 (en) *  20051003  20070531  Trimble Navigation Limited  MultipleGNSS and FDMA high precision carrierphase based positioning 

2007
 20070829 WO PCT/NL2007/050427 patent/WO2009028929A1/en active Application Filing
Patent Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

US6127968A (en) *  19980128  20001003  Trimble Navigation Limited  Onthefly RTK positioning system with single frequency receiver 
US6313788B1 (en) *  19980814  20011106  Seagull Technology, Inc.  Method and apparatus for reliable interantenna baseline determination 
US7148843B2 (en) *  20030702  20061212  Thales North America, Inc.  Enhanced real time kinematics determination method and apparatus 
US20070120733A1 (en) *  20051003  20070531  Trimble Navigation Limited  MultipleGNSS and FDMA high precision carrierphase based positioning 
NonPatent Citations (1)
Title 

ODIJK D ET AL: "Recursive Detection, Identification and Adaptation of Model Errors for Reliable HighPrecision GNSS Positioning and Attitude Determination" RECENT ADVANCES IN SPACE TECHNOLOGIES, 2007. RAST '07. 3RD INTERN ATIONAL CONFERENCE ON, IEEE, PI, 1 June 2007 (20070601), pages 624629, XP031123377 ISBN: 9781424410569 * 
Cited By (5)
Publication number  Priority date  Publication date  Assignee  Title 

CN102353969A (en) *  20110902  20120215  东南大学  Method for estimating phase deviation in precise singlepoint positioning technology 
WO2015145718A1 (en) *  20140328  20151001  三菱電機株式会社  Positioning device 
JPWO2015145718A1 (en) *  20140328  20170413  三菱電機株式会社  Positioning device 
US10371820B2 (en)  20140328  20190806  Mitsubishi Electric Corporation  Positioning device 
US10393879B2 (en)  20140328  20190827  Mitsubishi Electric Corporation  Global positioning device 
Similar Documents
Publication  Publication Date  Title 

US6229479B1 (en)  Relative position measuring techniques using both GPS and GLONASS carrier phase measurements  
US7522099B2 (en)  Position determination using carrier phase measurements of satellite signals  
US5495257A (en)  Inverse differential corrections for SATPS mobile stations  
EP1430272B1 (en)  Hybrid inertial navigation system with improved integrity  
US6975266B2 (en)  Method and apparatus for locating position of a satellite signal receiver  
CA2550600C (en)  A method for combined use of a local rtk system and a regional, widearea, or global carrierphase positioning system  
US6424915B1 (en)  System for determining the heading and/or attitude of a body  
US7133772B2 (en)  Method and apparatus for navigation using instantaneous Doppler measurements from satellites  
EP0856747A1 (en)  Method and apparatus for attitude determination utilizing an inertial measurement unit and a plurality of satellite transmitters  
US7937190B2 (en)  Apparatus for an automated aerial refueling boom using multiple types of sensors  
RU2354991C2 (en)  Method of using three gps frequencies for resolving integral uncertainties of carrier phases  
US6720914B2 (en)  Carrierphasebased relative positioning device  
US8682581B2 (en)  Vehicle navigation using nonGPS LEO signals and onboard sensors  
JP4446569B2 (en)  Carrier phase differential positioning device  
CA2494417C (en)  Estimation and resolution of carrier wave ambiguities in a position navigation system  
US6424914B1 (en)  Fullycoupled vehicle positioning method and system thereof  
US8255160B2 (en)  Integrated mobile terminal navigation  
EP2376936B1 (en)  Methods and systems to increase accuracy in the navigation of single frequency receivers  
US7117417B2 (en)  Method for generating clock corrections for a widearea or global differential GPS system  
US9035826B2 (en)  Satellite differential positioning receiver using multiple baserover antennas  
US7504996B2 (en)  Satellitebased positioning receiver with improved integrity and continuity  
US20030154049A1 (en)  Attitude angle detecting apparatus  
US5831576A (en)  Integrity monitoring of location and velocity coordinates from differential satellite positioning systems signals  
US6781542B2 (en)  Method and system for estimating ionospheric delay using a single frequency or dual frequency GPS signal  
JP3408593B2 (en)  Method and apparatus for predicting the position of a satellite in a navigation system for a satellitebased 
Legal Events
Date  Code  Title  Description 

121  Ep: the epo has been informed by wipo that ep was designated in this application 
Ref document number: 07808560 Country of ref document: EP Kind code of ref document: A1 

NENP  Nonentry into the national phase in: 
Ref country code: DE 

122  Ep: pct app. not ent. europ. phase 
Ref document number: 07808560 Country of ref document: EP Kind code of ref document: A1 