A PPP positioning method, system, computer device, and medium

Abstract

The present invention belongs to the technical field of PPP positioning, and discloses a PPP positioning method, system, computer device and medium. The method of the present invention takes into account the satellite clock difference constant deviation of the PPP-B2b service. By constructing a PPP mathematical model based on the PPP-B2b service, it fully considers the characteristics of the B2b satellite clock difference correction number constant deviation, estimates the satellite clock difference constant deviation as a parameter and adopts a piecewise constant estimation strategy to separate the constant deviation, so as to weaken its influence on other PPP parameters to be estimated, which can effectively shorten the PPP convergence time and improve the PPP positioning performance. In addition, the method of the present invention is simple and reliable, and convenient for real-time implementation. Compared with the conventional method for processing the B2b signal satellite clock difference correction number constant deviation, the strategy of piecewise constant estimation of the satellite clock difference constant deviation proposed by the present invention does not require external provision of additional deviation products, the algorithm is simple and easy to implement, and it is convenient for real-time application.

Description

Technical Field

[0001] The present invention belongs to the technical field of PPP positioning, and particularly relates to a PPP positioning method, system, computer device and medium. The PPP positioning method takes into account the satellite clock difference constant deviation of the PPP-B2b service. Background Technique

[0002] Precise Point Positioning (PPP) technology can obtain high-precision position, atmospheric delay, clock difference and other information globally, and has been widely used in fields such as precise time service, water vapor monitoring, earthquake monitoring, and precise orbit determination. PPP uses the GNSS observations of a single receiver and precise ephemeris for data processing, and the positioning accuracy can reach the decimeter level or even the centimeter level. The International GNSS Service organization provides accurate satellite orbits and clock products to global users through the Internet to meet the needs of real-time precise positioning. However, if the network communication service is interrupted, users will not be able to receive corrections, thus affecting the PPP performance.

[0003] Since 2020, the BeiDou Navigation Satellite System version 3 (BDS-3) has broadcast precise orbit, clock difference and other correction information to the Asia-Pacific and surrounding regions through the geostationary satellite B2b downlink signal, providing regional PPP services (referred to as PPP-B2b). PPP-B2b only provides regional services. In the satellite clock difference estimation, pseudorange observations are required to provide a relative reference at the beginning of a new satellite arc segment, so there are differences in the clock difference constant deviations included in different arc segments in the satellite clock difference estimation results. This deviation will be absorbed by the pseudorange residuals during PPP solution, thus affecting the PPP convergence speed.

[0004] Currently, some studies estimate and correct the satellite clock difference constant deviation and signal distortion deviation existing in the PPP-B2b clock difference product to the observations together, which greatly improves the PPP convergence speed. However, due to the different satellite clock difference constant deviations in different arc segments, the extracted clock difference constant deviation is not applicable to all arc segments, so it is particularly difficult to apply this method in real-time PPP. In addition, some studies use multiple reference stations to jointly estimate the constant deviation and broadcast the constant deviation through the network for user real-time positioning. Nevertheless, the application of this method is significantly limited in areas with insufficient network coverage in remote regions.

[0005] Therefore, how to effectively reduce the influence of the clock difference constant deviation on real-time PPP positioning is of great value for promoting the popularization of PPP-B2b services and improving precise navigation, positioning and timing services. Summary of the Invention

[0006] The object of the present invention is to propose a PPP positioning method considering the satellite clock difference constant deviation of PPP-B2b service. By constructing a PPP mathematical model based on PPP-B2b service, fully considering the characteristics of the B2b satellite clock difference correction constant deviation, parameter estimation of the clock difference constant deviation is carried out and a piecewise constant estimation strategy is adopted to separate the constant deviation, so as to weaken the influence of the clock difference constant deviation existing in the PPP-B2b clock difference correction on PPP positioning, shorten the PPP convergence time, and improve the PPP positioning performance.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] A PPP positioning method includes the following steps:

[0009] Step 1. Obtain the satellite phase observations, pseudorange observation data, broadcast ephemeris, various corrections broadcast by B2b signals, and antenna phase center products at the station.

[0010] Step 2. Use the broadcast ephemeris to check the data quality of the observations, eliminate gross errors, delete satellite data without satellite ephemeris or with incomplete observations, and delete satellites without B2b signal corrections or with abnormal B2b signal corrections according to the validity period of B2b signal correction information and the user ranging accuracy index; at the same time, judge whether the satellite clock difference corrections of the satellites with normal B2b signal corrections in the current epoch are continuous and mark them.

[0011] Step 3. Error correction;

[0012] Step 4. Consider the satellite clock difference constant deviation in the B2b signal clock difference correction , and establish a PPP mathematical model considering satellite constant deviation, including PPP observation equations and stochastic models;

[0013] Step 5. According to the satellite ephemeris continuity mark in Step 2, adopt a piecewise constant estimation strategy;

[0014] Step 6. Perform Kalman filtering estimation according to the mathematical model established in Step 4 and the parameter estimation strategy determined in Step 5 to obtain the three-dimensional position of the receiver, , tropospheric zenith wet delay, ambiguity, and receiver clock difference parameters.

[0015] In addition, on the basis of the above PPP positioning method considering the satellite clock difference constant deviation of PPP-B2b service, the present invention also proposes a corresponding PPP positioning system, which adopts the following technical solutions:

[0016] A PPP positioning system includes the following modules:

[0017] A data acquisition module, which is used to acquire satellite phase observations, pseudorange observation data, broadcast ephemeris, various corrections broadcast by the B2b signal, and antenna phase center products at the measurement station;

[0018] A preprocessing module, which is used to perform data quality inspection and gross error rejection on the observations using the broadcast ephemeris, delete satellite data without satellite ephemeris or with incomplete observations, and delete satellites without B2b signal corrections or with abnormal corrections according to the validity period of the B2b signal correction information and the user's ranging accuracy index; at the same time, judge whether the satellite clock corrections of the satellites with normal B2b signal corrections in the current epoch are continuous and mark them;

[0019] An error correction module, which is used to perform error correction;

[0020] A PPP mathematical model establishment module, which is used to consider the satellite clock constant deviation in the B2b signal clock correction and establish a PPP mathematical model considering the satellite constant deviation, including PPP observation equations and stochastic models;

[0021] A parameter estimation strategy determination module, which is used to, according to the continuous marking of the satellite ephemeris, adopt a piecewise constant estimation strategy;

[0022] And a parameter estimation module, which is used to perform Kalman filter estimation according to the established mathematical model and the determined parameter estimation strategy to obtain the three-dimensional position of the receiver, , tropospheric zenith wet delay, ambiguity, and receiver clock parameters.

[0023] In addition, based on the above PPP positioning method, the present invention also proposes a computer device, which includes a memory and one or more processors. An executable code is stored in the memory, and when the processor executes the executable code, it is used to implement the steps of the above-mentioned PPP positioning method.

[0024] In addition, based on the above PPP positioning method, the present invention also proposes a computer-readable storage medium, on which a program is stored. When the program is executed by a processor, it is used to implement the steps of the above-mentioned PPP positioning method.

[0025] The present invention has the following advantages:

[0026] As described above, the present invention relates to a PPP positioning method that takes into account the satellite clock difference constant deviation of the PPP-B2b service, effectively weakens the error influence, and improves the positioning performance. Specifically, the present invention constructs a PPP mathematical model based on the PPP-B2b service, fully considers the characteristics of the B2b satellite clock difference correction constant deviation, estimates the clock difference constant deviation parameter and adopts a piecewise constant estimation strategy to separate the constant deviation to weaken its influence on other PPP parameters to be estimated, which can effectively shorten the PPP convergence time and improve the PPP positioning performance. In addition, the method of the present invention is simple and reliable and convenient for real-time implementation. Compared with the conventional B2b signal satellite clock difference correction constant deviation processing method, the strategy of piecewise constant estimation of the clock difference constant deviation proposed by the present invention does not require external provision of additional deviation products, the algorithm is simple and easy to implement, and it is convenient for real-time application. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 It is a flowchart of the PPP positioning method that takes into account the satellite clock difference constant deviation of the PPP-B2b service in the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0028] The present invention will be further described in detail below with reference to the drawings and specific embodiments:

[0029] Embodiment 1

[0030] As Figure 1 shown, in order to improve the PPP-B2b service performance, Embodiment 1 of the present invention describes a PPP positioning method that takes into account the satellite clock difference constant deviation of the PPP-B2b service. The PPP positioning method includes the following steps:

[0031] Step 1. Data acquisition.

[0032] Obtain the satellite phase observations, pseudorange observation data, and products required for data processing (i.e., PPP positioning) at the station, including broadcast ephemeris, various corrections broadcast by the B2b signal, and antenna phase center products, etc.

[0033] Step 2. Data preprocessing.

[0034] Using the broadcast ephemeris, perform data quality inspection and gross error rejection on the observations obtained in Step 1, delete the satellite data without satellite ephemeris or incomplete observations, and delete the satellites without B2b signal corrections or with abnormal B2b signal corrections according to the validity period of the B2b signal correction information and the user ranging accuracy index. The purpose is to provide data reliability guarantee, filter out invalid data in advance, reduce unnecessary computational workload, and prevent divergence of iterative calculations caused by data quality problems.

[0035] Meanwhile, it is determined whether the satellite clock corrections with normal B2b signal corrections in the current epoch are continuous and marked. The determination process of whether the satellite clock corrections with normal B2b signal corrections in the current epoch are continuous is as follows:

[0036] If the satellite clock corrections are not continuously broadcast or the B2b signal clock corrections received by the user exceed the data validity period (e.g., 12 s), it is considered that the satellite clock corrections with normal B2b signal corrections in the current epoch are not continuous.

[0037] Here, the purpose of determining whether the satellite clock corrections are continuous is that within a continuous arc segment, the constant bias included in the satellite clock corrections is consistent, and a constant estimation is adopted during the solution, which determines the solution strategy for the constant bias.

[0038] Step 3. Error correction.

[0039] Based on the clock correction parameters and orbit correction parameters provided by the B2b signal, the broadcast ephemeris clock and orbit are corrected respectively; the pseudorange observations are corrected using the inter-code bias correction parameters provided by the B2b signal.

[0040] Meanwhile, the tropospheric delay, relativistic effect, tide, antenna phase center, and earth rotation error are corrected.

[0041] Step 4. Establishment of PPP mathematical model.

[0042] Taking the dual-frequency ionosphere-free combination PPP positioning as an example, the ionosphere-free combination observations are constructed. Considering that the constant bias of the B2b signal clock product will be absorbed by the pseudorange residuals during the PPP solution, thus affecting the PPP convergence speed, and comprehensively considering the PPP model strength, the PPP observation equation considering the satellite constant bias is established as shown in the following formula (1): The PPP observation equation is as follows:

[0043] (1)

[0044] In the formula, and represent the receiver and satellite numbers respectively; represents the ionosphere-free combination pseudorange observation, represents the ionosphere-free carrier phase observation; represents the geometric distance from the satellite at the signal emission time to the receiver at the reception time.

[0045] represents the receiver clock error after absorbing the receiver pseudorange hardware delay bias; represents the satellite clock error; represents the speed of light; Represents the projection function from the tropospheric zenith direction of the station to the inclined direction between the station and the satellite.

[0046] Represents the wet component delay in the tropospheric zenith direction; Represents the wavelength of the ionosphere-free combination ambiguity.

[0047] Is an ambiguity parameter, representing the ionosphere-free combination ambiguity of the satellite clock error correction constant deviation, satellite phase hardware delay, receiver phase hardware delay, and receiver pseudorange hardware delay for the absorbed B2b signal.

[0048] and Represent the pseudorange and carrier observation noise respectively.

[0049] It can be seen from the expression of formula (1) that the ambiguity parameter Absorbs the constant deviation of the satellite clock error correction for the B2b signal. By introducing the clock error constant deviation parameter in formula (1), it is beneficial to analyze the influence of the constant deviation on the positioning result.

[0050] The present invention constructs a PPP mathematical model based on PPP-B2b service, fully considering the characteristics of the constant deviation of the B2b satellite clock error correction, that is, there is a constant deviation in the satellite clock error correction, but this deviation is the same in continuous arcs and different in different arcs. The satellite clock constant deviation is estimated as a parameter, and a piecewise constant estimation strategy is adopted to separate the constant deviation (introducing a constant deviation parameter for estimation, and piecewise estimation is because the constant deviation is different in different arcs) to weaken its influence on other PPP parameters to be estimated, which can effectively shorten the PPP convergence time and improve the PPP positioning performance.

[0051] Calculate the elevation angle of the satellite according to the satellite position and the approximate position of the station, and determine the corresponding stochastic model according to the satellite elevation angle, as shown in formula (2).

[0052] (2)

[0053] In the formula, Represents the observation variance, Represents the satellite elevation angle. Errors such as the antenna phase offset and change between the satellite and the receiving antenna, the earth's rotation, antenna phase winding, tidal correction, and relativistic effect are corrected using corresponding models.

[0054] Step 5. Determine the parameter estimation strategy.

[0055] Since the clock error correction provided by the B2b signal has a satellite-related clock constant deviation and the deviation is inconsistent in different arcs, according to the satellite ephemeris continuity mark in the data preprocessing in step 2, for Adopt the following piecewise constant estimation strategy.

[0056] When the satellite is available and the satellite clock correction is continuous, the clock constant deviation related to the satellite remains unchanged. It is necessary to first After initializing the satellite constant deviation, Adopt the constant estimation strategy, and regard the satellite constant deviation As a constant for estimation, separate the constant deviation to weaken its influence on the PPP parameters to be estimated.

[0057] When the satellite is unavailable, retain the Value and its variance obtained from the last solution when the satellite was available, and when the satellite signal becomes available again, judge whether the satellite clock correction of the currently available satellites is continuously marked.

[0058] If the satellite clock correction is continuous, use the Value and its variance retained for this satellite to estimate the satellite clock constant deviation; otherwise, if the satellite clock correction is not continuous, Re-initialize, and no longer use the Value and its variance retained by the filter. In addition, since the ambiguity parameter absorbs the satellite clock constant deviation, the ambiguity parameter also needs to be re-estimated.

[0059] It should be noted here that when the satellite observation value is unavailable, the signal-to-noise ratio of the satellite observation value is lower than the preset signal-to-noise ratio threshold, or the elevation angle is lower than the preset satellite cut-off elevation angle, it indicates that the satellite is unavailable.

[0060] According to the characteristics of the satellite clock correction, the present invention adopts a piecewise constant estimation and processing strategy for the satellite clock corrections in different arcs, which is beneficial to separating the constant deviation faster to weaken its influence on other PPP parameters to be estimated, and can effectively shorten the PPP convergence time and improve the PPP positioning performance.

[0061] Step 6. Parameter estimation.

[0062] Perform Kalman filter estimation according to the mathematical model established in step 4 and the parameter estimation strategy determined in step 5 to obtain the three-dimensional position of the receiver, , tropospheric zenith wet delay, ambiguity, and receiver clock parameters.

[0063] Embodiment 2

[0064] This Embodiment 2 describes a PPP positioning system, which is based on the same inventive concept as the PPP positioning method in the above Embodiment 1. Specifically, the PPP positioning system in this embodiment includes the following modules:

[0065] A data acquisition module, configured to acquire satellite phase observations, pseudorange observation data, broadcast ephemeris, various correction numbers broadcast by the B2b signal, and antenna phase center products at a measuring station;

[0066] A preprocessing module, configured to perform data quality inspection and gross error rejection on the observations using the broadcast ephemeris, delete satellite data without satellite ephemeris or with incomplete observations, and delete satellites without B2b signal correction numbers or with abnormal correction numbers according to the validity period of the B2b signal correction information and the user's ranging accuracy index; meanwhile, determine whether the satellite clock correction numbers of satellites with normal B2b signal correction numbers in the current epoch are continuous and mark them;

[0067] An error correction module, configured to perform error correction;

[0068] A PPP mathematical model establishment module, configured to consider the satellite clock constant deviation in the B2b signal clock correction number and establish a PPP mathematical model considering the satellite constant deviation, including a PPP observation equation and a stochastic model;

[0069] A parameter estimation strategy determination module, configured to, according to the continuous marking of the satellite ephemeris, adopt a piecewise constant estimation strategy;

[0070] and a parameter estimation module, configured to perform Kalman filter estimation according to the established mathematical model and the determined parameter estimation strategy to obtain the three-dimensional position of the receiver, , tropospheric zenith wet delay, ambiguity, and receiver clock difference parameters.

[0071] It should be noted that in the PPP positioning system of this embodiment, the implementation processes of the functions and roles of each functional module are specifically detailed in the implementation processes of the corresponding steps in the method in the above-mentioned Embodiment 1, and will not be elaborated here.

[0072] Embodiment 3

[0073] This Embodiment 3 describes a computer device, which includes a memory and one or more processors. Executable code is stored in the memory, and when the processor executes the executable code, it is used to implement the steps of the PPP positioning method in the above-mentioned Embodiment 1.

[0074] In this embodiment, the computer device is any device or apparatus with data processing capabilities, which will not be elaborated here.

[0075] Embodiment 4

[0076] This Embodiment 4 describes a computer-readable storage medium, on which a program is stored, and when the program is executed by a processor, it is used to implement the steps of the PPP positioning method in the above-mentioned Embodiment 1.

[0077] The computer-readable storage medium may be an internal storage unit of any device or apparatus with data processing capabilities, such as a hard disk or memory, or may be an external storage device of any device with data processing capabilities, such as a plug-in hard disk, a Smart Media Card (SMC), an SD card, a Flash Card, etc. equipped on the device.

[0078] Of course, the above description is only a preferred embodiment of the present invention. The present invention is not limited to listing the above embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any person skilled in the art under the teaching of this specification fall within the substantial scope of this specification and should be protected by the present invention.

Claims

1. A PPP positioning method, characterized in that, It includes the following steps: Step 1. Obtain satellite phase observations, pseudorange observation data, broadcast ephemeris, various corrections broadcast by the B2b signal, and antenna phase center products at the station; Step 2. Use the broadcast ephemeris to perform data quality inspection and gross error rejection on the observations, delete satellite data without satellite ephemeris or with incomplete observations, and delete satellites without B2b signal corrections or with abnormal B2b signal corrections according to the validity period of the B2b signal correction information and the user's ranging accuracy index; at the same time, determine whether the satellite clock corrections of satellites with normal B2b signal corrections in the current epoch are continuous and mark them; Step 3. Error correction; Step 4. Consider the satellite clock constant deviation O in the B2b signal clock error correction, s and establish a PPP mathematical model considering the satellite constant deviation, including the PPP observation equation and the stochastic model; Step 5. According to the labeling of whether the satellite ephemeris is continuous in Step 2, perform s a piecewise constant estimation strategy; In step 5, for the satellite clock offset constant deviation O s , the following piecewise constant estimation strategy is adopted: When the satellite is available and the satellite clock error correction is continuous, the clock error constant deviation related to the satellite remains unchanged. First, initialize the satellite constant deviation O s Then, for the satellite constant deviation O s Adopt a constant estimation strategy, treat the satellite constant deviation O s as a constant for estimation, separate the constant deviation, and weaken its influence on the PPP parameters to be estimated; When the satellite is unavailable, retain the O value and its variance obtained from the last solution when the satellite was available, and when the satellite signal becomes available again, determine whether the satellite clock correction of the currently available satellites is continuously marked; s ​ If the satellite clock correction is continuous, use the reserved O value of this satellite and its variance to estimate the satellite clock constant deviation; otherwise, re-initialize O, no longer use the reserved O value and its variance, and re-estimate the ambiguity parameter as well; s value and its variance to estimate the satellite clock constant deviation; otherwise, re-initialize O s and no longer use the reserved O s value and its variance, and re-estimate the ambiguity parameter as well; Step 6. Perform Kalman filtering estimation according to the mathematical model established in Step 4 and the parameter estimation strategy determined in Step 5, and obtain the three-dimensional position of the receiver, O s , the tropospheric zenith wet delay, the ambiguity, and the receiver clock offset parameter.

2. The PPP positioning method according to claim 1, wherein in the said Step 2, if the satellite clock corrections are not continuously broadcast or the B2b signal clock corrections received by the user exceed the data validity period, it is considered that the satellite clock corrections of satellites with normal B2b signal corrections in the current epoch are not continuous.

3. The PPP positioning method according to claim 1, wherein in the said Step 3, the error correction process is as follows: Correct the broadcast ephemeris clock and orbit according to the clock correction parameters and orbit correction parameters provided by the B2b signal; correct the pseudorange observations using the inter-code bias correction parameters provided by the B2b signal; At the same time, correct the tropospheric delay, relativistic effect, tide, antenna phase center, and earth rotation error.

4. The PPP positioning method according to claim 1, wherein In step 4, a PPP observation equation considering the satellite constant bias O s is established as shown in formula (1): where r and s represent the receiver and satellite numbers respectively; represents the ionosphere-free combined pseudorange observation; represents the ionosphere-free carrier phase observation; represents the geometric distance from the satellite at the signal emission time to the receiver at the reception time; represents the receiver clock error after absorbing the pseudorange hardware delay deviation of the receiver; dt s represents the satellite clock error; c represents the speed of light; represents the projection function from the tropospheric zenith direction of the station to the inclined direction between the station and the satellite; D w,r represents the wet component delay in the zenith direction of the troposphere; λ IF represents the wavelength of the ionosphere-free combination ambiguity; is the ambiguity parameter, representing the ionosphere-free combined ambiguity of the satellite clock error correction number constant deviation, satellite phase hardware delay, receiver phase hardware delay, and receiver pseudorange hardware delay for absorbing the B2b signal; and represent the pseudorange and carrier observation noise respectively; calculate the elevation angle of the satellite according to the satellite position and the approximate position of the station, and determine the random model corresponding to the satellite elevation angle according to the satellite elevation angle, as shown in formula (2); where p represents the observation variance and θ represents the satellite elevation angle.

5. The PPP positioning method according to claim 1, wherein in the said Step 5, when the satellite observations are unavailable, the signal-to-noise ratio of the satellite observations is lower than a preset signal-to-noise ratio threshold, or the elevation angle is lower than a preset satellite cutoff elevation angle, it indicates that the satellite is unavailable.

6. A PPP positioning system, characterized in that, It includes the following modules: A data acquisition module for obtaining satellite phase observations, pseudorange observation data, broadcast ephemeris, various corrections broadcast by the B2b signal, and antenna phase center products at the station; A preprocessing module for using the broadcast ephemeris to perform data quality inspection and gross error rejection on the observations, deleting satellite data without satellite ephemeris or with incomplete observations, and deleting satellites without B2b signal corrections or with abnormal corrections according to the validity period of the B2b signal correction information and the user's ranging accuracy index; at the same time, determine whether the satellite clock corrections of satellites with normal B2b signal corrections in the current epoch are continuous and mark them; An error correction module for performing error correction; PPP mathematical model establishment module, which is used to consider the satellite clock constant deviation O in the B2b signal clock error correction number s , establish a PPP mathematical model considering satellite constant deviation, including PPP observation equation and stochastic model; A parameter estimation strategy determination module, which is used to mark according to whether the satellite ephemeris is continuous, and for O s adopt a piecewise constant estimation strategy; For the satellite clock offset constant deviation O s , the following piecewise constant estimation strategy is adopted: When the satellite is available and the satellite clock error correction is continuous, the clock error constant deviation related to the satellite remains unchanged. First, initialize the satellite constant deviation O s Then, for the satellite constant deviation O s Adopt a constant estimation strategy, treat the satellite constant deviation O s as a constant for estimation, separate the constant deviation, and weaken its influence on the PPP parameters to be estimated; When the satellite is unavailable, retain the O value and its variance obtained from the last solution when the satellite was available, and when the satellite signal becomes available again, determine whether the satellite clock correction of the currently available satellites is continuously marked; s ​ If the satellite clock correction is continuous, use the reserved O value of this satellite and its variance to estimate the satellite clock constant deviation; otherwise, re-initialize O, no longer use the reserved O value and its variance, and re-estimate the ambiguity parameter as well; s s s ​​​ and a parameter estimation module, which is used to perform Kalman filter estimation according to the established mathematical model and the determined parameter estimation strategy, and obtain the three-dimensional position of the receiver, O s , tropospheric zenith wet delay, ambiguity, and receiver clock offset parameters.

7. A computer device, comprising a memory and one or more processors, wherein executable code is stored in the memory, characterized in that, When the processor executes the executable code, it implements the steps of the PPP positioning method according to any one of claims 1 to 5.

8. A computer-readable storage medium having a program stored thereon, characterized in that, When this program is executed by the processor, it is used to implement the steps of the PPP positioning method according to any one of claims 1 to 5.

Images