Recent developments on an interferometric multilateration measurement system for large volume coordinate metrology

A mobile multilateration measurement system developed at the Physikalisch-Technische Bundesanstalt (PTB) around 2010 has been thoroughly investigated and refined to gain better performance with smaller uncertainties even when applied to the calibration of large complex workpieces. The mathematical background of multilateration and the propagation of uncertainties for the algorithms involved is explained in detail. Using the example of simple 1D and 2D measuring tasks, the influence of certain parameters characterizing the setup of the measurement system on the overall uncertainty is quantified. A strategy is developed to incorporate multi-stylus measurements which are often inevitable when workpieces feature complex shapes. The findings are verified on a large involute gear which is 2 m in diameter. All measurements are performed on PTB’s large coordinate measuring machine with a working range of 5 m×4 m×2 m .


Introduction
In many industrial sectors, high quality requirements are placed on the production of mechanical components which often feature complex geometries. For example, large gears need to be manufactured very precisely to allow the highly efficient operation of turbines and to prevent damage leading to expensive outages in the long term. Quality assurance makes * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. great demands on manufacturing metrology. Large coordinate measuring machines (CMMs) are suitable for measuring large components in principle. However, they have reached their limits concerning the attainable accuracy, since conventional calibrations using measurement standards (e.g. according to ISO 10360-2 [1]) do not cover the whole measurement volume. This means the measurement volume of a large CMM can only be calibrated partially using the available material standards. Moreover, most of the commercially used CMMs do not conform to Abbe's principle [2], which makes them prone to errors. This especially applies to large measuring objects. Conventional calibration methods, such as length measurement deviations according to ISO 10360-2, are not task-specific. Therefore, it is difficult to draw conclusions about the measurement uncertainty to be expected for a real component, or even almost impossible for complex geometries. The procedure described in this article can solve this problem by establishing real task-specific calibrations.
Early works on interferometry and trilateration as, e.g. [3] tried to determine the volumetric accuracy of machine tools by establishing an error map using a laser interferometer and a ball bar.
Around 2010, the Physikalisch-Technische Bundesanstalt (PTB) developed a mobile multilateration (MLAT) system for 3D measurement applications (M3D3) [4,5]. It is based on LaserTracers (LTs) [6] as an interferometric reference system. Applying a special measurement procedure, it allows the 3D accuracy of the end effector of a positioning machine to be improved. In principle, the solution is suitable for any type of kinematic chain (serial, parallel) [7]. It mainly aims to improve the accuracy of conventional CMMs but can also be used to turn machine tools into qualified 3D measuring devices by determining task-specific error maps.
An LT is a self-tracking laser interferometer following a retroreflector moving in space and continuously executing relative length measurements. These length measurements satisfy Abbe's principle and are therefore well traceable to the metre, the SI unit of length, by calibrating the frequency of the laser light.
This paper summarizes the application of M3D3 as it calibrates large measurement standards in detail. The mathematical background and the current research to optimize the M3D3 method are included. In section 2, the principle of the method and its theoretical background are explained, followed by procedures to validate the basic function and performance in section 3. Strategies for the optimal application of the method in terms of the arrangement of the LTs and the distribution of the measuring points are presented in section 4. Section 5 introduces results of the application of the method to the measurement of a large gear standard.

Principle of the M3D3 method
The M3D3 method is a three stage procedure [4]. The first step is the measurement itself. It is used for generating the sample of 3D points to be corrected. In principle, any point sensor, tactile or optical, can be used. For this article, the method was applied to tactile measurements on a large CMM ( figure 1(a)).
In the second step, which is the task-specific error mapping with the MLAT reference system, the workpiece is removed from the measurement volume. The stylus is then replaced by a retroreflector mounted approximately at the same position as the stylus tip centre during the tactile measurement (i.e. the same stylus tip offset is effective). At least four LTs are placed into the measurement volume. The CMM replays the recorded measuring points of the first step, while the LTs follow the reflector ( figure 1(b)), measuring and recording lengths for each measuring point. Applying the MLAT evaluation [2,4,8], the 3D coordinates of the retroreflector positions can be calculated based on the measured lengths from the LTs. These points are used to correct the initial probing points.
In the final step, the measured quantities are re-evaluated with the corrected points using the original CMM evaluation software, i.e. the part program. Figure 2 illustrates the principle of the local error correction using the example of a circle measurement. The initial CMM measurement of the workpiece results in the measuring points (red filled circles, step 1) as indicated by the CMM. The positions of the CMM's probe may be slightly different when replaying the measuring points for the LT measurements. This is due to the reproducibility of the CMM positioning (in the range of micrometres) and, even more, due to a different offset of the stylus and the retroreflector (in the range of millimetres; red filled square in the inset (2)). The outcome of the MLAT evaluation is a further set of points (represented by a blue square in the inset as an example (2')). These points are used to define the local error vectors (black arrows) as the difference between the CMM's points (red filled square) and the MLAT points (blue square). To obtain the corrected measuring points (blue circles (1')), the local error vectors (black arrows) have to be added to the original measuring points (red filled circles).
It should be noted that this approach (determining the error vector for point 1 at a slightly different position 2) works due to the continuity of the error vector field.
Besides correcting measuring points recorded by a CMM, the M3D3 method has the big advantage of providing traceable measuring values based on the calibration of the frequency of the laser.

Mathematical description of the MLAT evaluation
To locate a point in three-dimensional space, at least three distances to known points (i.e. LTs) must be given. However, since the positions of the LTs are in general unknown as well, a fourth LT is needed in order to calculate both the measuring points and the LT positions [4,7].
Let P i = (x i , y i , z i ) for i = 1, . . . , n denote the measuring points, and let T j = (x 0j , y 0j , z 0j ) for j = 1, . . . , m be the positions of LT j . Since the LTs cannot measure absolute lengths, but only lengths relative to an arbitrary starting point, the total length between T j and the measuring point P i consists of the measured relative length l ij and the unknown dead path l 0j . More precisely, for the m · n measured lengths l ij between the m LTs and the n measuring points P i , the equations: are obtained. Here, s j is a scale factor of LT j , which reflects the uncertainty of the laser frequency, and can vary around 1 within the limits of the calibration uncertainty of the frequency. Moreover, w ij is the residual between the measured distance l ij and the value predicted by the model applied to the estimated parameters. Equation (1) defines a set of m · n equations, and there are m more equations for the scale factors s j . The number of unknowns is given by 5 m + 3n (n point coordinates x i , y i , z i  Red ■ (2): probe position while replaying the measuring points for the LT measurement. They are slightly different due to the limited positioning accuracy of the CMM and different offsets between the stylus and the retroreflector. Blue □ (2'): points from the LT measurement (evaluated by MLAT), resulting in local error vectors (difference between points 2 and 2'). Blue (1'): corrected measuring points, adding the local error vectors to points 1. and five unknowns x 0j , y 0j , z 0j , l 0j , s j of each LT). However, in order to determine the coordinates of the measuring points and the LT positions, a coordinate system has to be fixed, which reduces the number of unknowns to 5 m + 3 n − 6. One necessary condition rendering the given equation system solvable is therefore given by mn + m ⩾ 5m + 3n − 6 or (m − 3)n ⩾ 4m − 6. From this relation, it can easily be seen that at least m = 4 LTs are needed to apply the MLAT algorithm. For exactly four LTs, 10 or more measuring points are needed to generate a solvable equation system. Typical point numbers of 50 to several hundreds lead to an over-determined system of equations with a large degree of freedom ν = (m − 3)n − 4m + 6. To minimize the discrepancies between the equations, a weighted least-squares fit is applied.
Let β be a vector containing the coordinates x i , y i , z i of the measuring points, as well as the unknowns x 0j , y 0j , z 0j , l 0j , s j of each LT, and let the function ϕ be defined by: The measured lengths from the LTs are combined in a vector L i+n( j−1) = l ij . With these definitions, equation (1) can be written as: where w are the residuals of the length measurements. The uncertainties for the length measurements with the LTs are given by: where a and b are constants, describing length-independent and length-dependent uncertainty contributions. Since no correlations between the individual length measurements are taken into account, the weight matrix Q used for the leastsquares fit is a diagonal matrix with the entries q i+n( j−1) = 1/u 2 ij on the diagonal. To ensure the traceability of the measurement, the scale factors s j must be fitted against 1. To take this into account, the function ϕ is extended with ϕ nm+j (β) = s j − 1 for j = 1, . . . , m. Accordingly, the vector L is extended by L nm+j = 0, and the weight matrix is extended by the diagonal elements q nm+j = 1/u(s j ) 2 , where u(s j ) is the standard uncertainty of the calibration of the lasers' wavelengths. The minimization problem then reads: Since this is a non-linear minimization problem, it is approximated by a sequence of linear problems (Gauss-Newton algorithm [9]). Therefore, for any given point β, the function ϕ is linearized by ϕ(β + dβ) = ϕ(β) + Jdβ + . . ., where J = Dϕ(β) is the Jacobian matrix of ϕ at the point β. Withŵ = L − ϕ(β), the linearized problem reads: which implies the equation J t QJdβ = J t Qŵ for dβ.
In order to find a unique solution to this minimization problem, a coordinate system needs to be fixed. One possibility is to use the initial values P 0 for the measuring points reported by the CMM as a reference, i.e. the coordinate system is chosen such that the sum of the squared distances ∥P i − P 0 i ∥ 2 takes a minimum. Let where × denotes the cross product, then the solution must additionally fulfil the constraints ψ(β) = 0. This leads to: with G = Dψ(β), and hence to the constraint ψ(β) = −Gdβ for the calculation of dβ. The constrained system can then be solved by using Lagrange multipliers λ, i.e. the equation: has to be solved. If the iteration is started with a vector β 0 containing the initial measuring points x 0 i , y 0 j , z 0 j , then ψ(β 0 ) = 0 and also ψ(β) = 0 for allβ calculated by iterations. It follows that the solution for dβ is given by: This calculation has to be iterated withβ replaced byβ + dβ until the weighted sum of the squared discrepancies reaches its minimum, where in each step, the matrices J, G and Q must be updated according to the actual value ofβ.

MLAT uncertainty and residuals
To calculate the uncertainties of the single points from the multilateration procedure (in the following referred to as 'MLAT uncertainty'), the covariance matrix C β of the solution β can be used. This covariance matrix is given by the first block of the matrix in equation (10), i.e.
The diagonal entries of C β are the variances of the elements of β, i.e. of the coordinates of the measuring points and LT positions, the dead path lengths and the scale factors. If a measurement value is a function f of β, also the off-diagonal elements of C β must be taken into account for the calculation of the uncertainty of f(β). For a linear function f(β) = Aβ with a certain matrix A, the covariance matrix of f(β) is given by AC β A t . For non-linear functions, in some cases a linearization can be applied as known from the propagation law of uncertainties. For example, let h ij (β) = ∥P i − P j ∥ be the distance between two points, and let B = Dh ij (β) be the Jacobian matrix of h ij at the point β. The propagation of uncertainty then leads to u 2 hij(β) = BC β B t for the uncertainty of this distance. For more complex functions f, the covariance matrix C β can also be used to generate random values in order to perform a Monte Carlo simulation for the calculation of u f(β) .
With the residuals from the fit w = L − ϕ(β), the chi-square χ 2 = w t Qw and the reduced chi-square [10]. If the model fits to the data and the estimates for the uncertainties are reasonable, the Birge ratio should be close to one. (However, a value of the Birge ratio close to one does not imply that the model is correct.) A value significantly larger than one indicates either a poor model, or that the uncertainties for the length measurements were assumed to be too small. A value much smaller than one could mean that the estimates for the uncertainties are too large.
It can be shown that the covariance matrix of the residuals is given by V w = Q −1 − JC β J t . The standardized residualsw ij related to the length measurements are hence given by: with u ij being the uncertainties of the length measurements as defined in equation (4), and H i+n( j−1) , the i + n(j − 1)-th diagonal element of the matrix H = JC β J t . If the model fits to the data, the standardized residuals should follow a normal distribution. A distribution of the residuals that deviates strongly from a normal distribution is also a hint that the assumptions made were not correct.

Stabilization points
In many cases, the measuring points are not distributed in a way that would be desirable for a good MLAT configuration (i.e. rather inappropriate: a more or less coplanar distribution of points). It might thus be necessary to measure some additional points to make the mathematical problem solvable and numerically more stable. These stabilization points should be distributed in space around the workpiece, especially covering the directions the workpiece does not extend to. They serve mainly to determine the positions of the LTs more precisely, which in turn, leads to smaller uncertainties of the measuring points. In sections 4.2 and 4.3, strategies to optimize the distribution and number of stabilization points will be presented.

Further uncertainty consideration
So far, only the MLAT uncertainty has been considered (see section 2.2). For the application of the method described here, however, the so-called task-specific measurement uncertainty [11] is ultimately of interest. This is the uncertainty of the characteristic which is to be measured on the workpiece (e.g. a diameter, form or location parameter). The local error correction (see figure 2) corrects for systematic errors related to the CMM, i.e. residual rigid-body geometric errors, static, nonrigid-body errors and the hysteresis of guideways. The MLAT uncertainty must be combined with the other relevant uncertainty contributions from the CMM, the tactile probing process, the environmental conditions and the workpiece to be calibrated itself. It must consequently be propagated into the task-specific measurement uncertainties. This can be done classically through sensitivity analysis. Due to the complexity of the measurements on CMMs, this will only be possible for a limited number of simple measuring tasks. Therefore, for more complex tasks, a Monte-Carlobased method according to the Guide to the Expression of Uncertainty in Measurement, Supplement 1 [12], and ISO 15530-4 [13], e.g. the Virtual Coordinate Measuring Machine (VCMM) by PTB [14], is suitable.
In this article, unless otherwise stated, task-specific uncertainty means only that part of the uncertainty that depends on the MLAT uncertainty.

Metrological challenges
There are several influences affecting the quality of the M3D3 result and thus the MLAT uncertainty. The most important parameters that have to be chosen carefully when setting up the measurement system are described in the following.
The LTs have a fixed and a length-dependent contribution to the measurement uncertainty of the length measurements of 0.2 µm + 0.3 · 10 −6 · L i (L i measured distance). This uncertainty contains, apart from the uncertainty of the frequency calibration, contributions from the roundness of the sphere, the controlling of the orientation as well as from the environmental data, i.e. temperature, humidity and pressure, that is used to correct the laser's wavelength. Hence, the distances to the measuring points should be as short as possible to ensure small measurement uncertainties.
Also, the angles between the laser beams hitting the retroreflector at a measuring point have to be taken into account to allow good intersection conditions and the minimization of the measurement uncertainty. Ideally, the LTs should face the measuring point from all spatial directions, forming an equilateral tetrahedron [15]. Since in reality this is not feasible (e.g. due to the limited acceptance angle of the retroreflector), the LT positions have to be chosen carefully. Simulations for various LT positions showed that good intersection conditions can be achieved if the angles between the outer LTs (that is, between the left and the right, or between the upper and the lower LT, respectively) are larger than 90 • . Angles larger than 120 • provide no further improvement with regard to the MLAT uncertainty. The waiting time and averaging of the length measurements of the LTs further influence the measurement uncertainty. At each measuring point, the retroreflector halts for a waiting time of several seconds to ensure that no vibrations disturb the measurement. The recorded length is an average over several length measuring values to account for short-time turbulence effects.
The necessary waiting time can make the M3D3 measurement very time-consuming. Since no discontinuity in the local error vectors is to be expected, for closely spaced measuring points, it is possible to thin out the measuring points. Subsequently, calculated local error vectors can be interpolated for the measuring points in between.
To minimise the influence of the environment, the measurements were performed in a measuring room with a very well controlled climate concerning the stability of temperature and humidity: the temperature is temporally stable within ±0.2 K during the entire measuring time, spatially over the whole measuring volume within 0.1 K, the relative humidity is controlled within ±3%.

Validation by reference measurements
One established method for testing the performance of measuring systems is to measure calibrated standards. In [4] and [5], several experiments suitable for the MLAT systems are described. To cross-check after updating the LT hardware and the extensive revision of the measuring software, some of these experiments were repeated on PTB's large CMM. This also served as as starting point for the optimization process.

Line measurements
Four straight lines were generated with ten measuring points each. The lines extended over almost the entire measuring volume of PTB's large CMM. In one measurement process, the CMM moved the reflector to these points to record the lengths between the points. The lengths returned by the CMM  were compared to reference lengths which were recorded simultaneously by a single LT with the tracking motors switched off once the laser beam had been aligned properly. To turn off the tracking has two aims: First, it serves as a cross-check if the reflector is really moved on a straight line as intended. Furthermore, it prevents small readjusting movements of the LT that lead to random noise on the measured length values. The reference LT thus served as a common laser interferometer. Afterwards, with the same points, an M3D3 measurement with four additional LTs was carried out to determine the corrected lengths for comparison with the reference lengths.
The reference LT was placed near a corner in the measurement volume, and the four lines were defined such that they intersect in the centre of the reference LT, as can be seen in For the MLAT measurement setup, 32 stabilization points were defined and arranged equidistantly on the edges of a cuboid frame that spans nearly the entire volume of the measuring points. They were measured in front of the actual measuring points of the four lines. The four LTs for the MLAT measurements were arranged according to the spatial boundaries of the CMM. The outer LTs were placed as far apart as possible in the y direction about 300 mm above the ground. Two LTs were positioned between the outer LTs: one at a height of about 300 mm and one on a stable tripod at a height of approximately 1600 mm.
The length deviations of the CMM for the test lengths were calculated as differences of the distances returned by the CMM and the reference lengths between two corresponding points. As can be seen in figure 4(a), for all lines these deviations are within the specification of the CMM, expressed as the maximum permissible error (MPE = 3.3 µm + 2.5 · 10 −6 · L) with an absolute size of up to 6 µm.
The error bars of the CMM length deviations indicate the expanded measurement uncertainties (k = 2) of the reference LT for the test lengths. The expanded uncertainty U Tj for a test length between point P 0 and point P j of a line, as can be seen in figure 4, is then where U 0 denotes the uncertainty for the laser length between the reference LT and the nearest measuring point of the considered line. The M3D3 length deviations are calculated as differences of the lengths between two measuring points of the MLAT fit and the corresponding reference lengths. These deviations are smaller than 0.8 µm at the largest test length, as can be seen for all lines in figure 4(b). This demonstrates the length correction capabilities of the M3D3 method over a large measuring range. The task-specific uncertainties of the M3D3 fit as length uncertainties between point P 0 and point P j (indicated as error bars in figure 4(b)) were obtained from the covariance matrix for the corresponding points, as shown in section 2.2.

Calibrated hole plate
To examine the influence of the positions of four LTs on the uncertainties of the M3D3 correction, a calibrated Zerodur ® hole plate standard was used. This standard is a twodimensional object with a side length of about 600 mm, featuring 44 holes on its edges, see figure 5. The calibration uncertainty of the distance L between the centre points of any two holes is specified by: The tactile and MLAT measurements were performed on PTB's large CMM. For this test, the numerical error compensation of the CMM was disabled in order to achieve larger correction effects in the further investigations. The hole plate standard was mounted horizontally while being tensionless. For the results of the tactile measurement in comparison to the M3D3 corrected results, see figure 5. Disabling the numerical error correction should help to validate the fundamental principle of the method on a CMM that is quite too good to need any correction with the numerical error correction enabled, see figure 6 for comparison.

Optimization of the method
The main objective of the optimization is to reduce the MLAT uncertainties. These uncertainties depend on the number and the relative positions of the LTs to the measuring points as well as on the number and arrangement of the stabilization points. If the tactile measuring points of the object to be corrected are available, the uncertainties can be determined a priori by simulation. For this purpose, the LT positions and stabilization points are defined, and the lengths between the LTs and the measuring points as well as the added stabilization points are calculated. The multilateration can be conducted with these artificial lengths. The MLAT uncertainties are comparable in both cases (measured or calculated lengths), as they depend almost exclusively on the laser lengths and the intersection conditions. The LT positions and the arrangement of the stabilization points can thus be optimized for small uncertainties prior to the actual MLAT measurement.

Position of the LTs
The influence of the LT positions was investigated on the basis of a configuration for the measurement of the hole plate. To determine optimal positions, simulations were performed based on the CMM coordinates of the tactilely measured points of the hole plate and a set of 32 stabilization points arranged on a cuboid frame around the plate.  The MLAT contribution to the uncertainties of the M3D3 corrected distances between any two holes was obtained by Monte Carlo simulation with 1000 iterations. In each iteration, the MLAT calculation was carried out on simulated laser distances perturbed by random values, drawn according to the LT's uncertainty specification. The centres of the holes were determined by least-squares circles, and the mean values of forward and backward measurements were calculated. Then the 946 distances between any pair of the 44 averaged hole centres were calculated. The length uncertainties were determined as the standard deviation of the simulated distances. To illustrate this, these 946 length uncertainties are plotted over the distances between corresponding holes (see figure 7(b)).
As a starting point, a symmetrical arrangement for the four LTs in the xy orientation was chosen. To achieve small task-specific uncertainties, an optimization of the LT positions should aim at minimizing the x and y components of the measuring point uncertainties. Hence, concerning the intersection conditions (see also section 2.5), a z position of LT 1 and LT 2 at the level of the plate plane is ideal. LT 3 and LT 4 were set at the heights z = 300 mm and z = 1200 mm, respectively.
Varying the LT positions in the simulation setup, the following observations could be made: an upward or downward shift of LT 1 and LT 2 leads to larger uncertainties. In general, larger distances of the LTs to the measuring points lead to larger uncertainties, since the length-dependent part of the measurement uncertainty of the LTs has more impact. The best results were obtained for small average distances of LT 1 and LT 2 to the measuring points. Considering the restriction that the angle α must be smaller than 180 • to guarantee the sight of the laser beams on the reflector, a trade-off had to be found. The less relevant z components of the uncertainties could be reduced by positioning LT 3 and LT 4 closer to the measuring points or by setting the upper LT 4 higher for a larger intersection angle in z. However, this increased the x and y components of the resulting uncertainties. For the optimized 'virtual' setup in figure 7(a), uncertainties even smaller than the uncertainties of calibration U(L) of the plate could be achieved, as can be seen in figure 7(b) (red point cloud compared to blue line).
A setup for a real measurement has some restrictions concerning the spatial conditions in the measuring volume of the CMM and the equipment available for positioning the LTs and adjusting the retroreflector. Taking this into account, a setup for the four LTs was found for the verification measurements close to the optimal positions. LT 1 and LT 2 had to be mounted below the plate at a height of z = 427 mm. LT 3 had a height of z = 288 mm and LT 4 of z = 1181 mm (see black circles in figure 7(a)).
The same stabilization points as in the simulation setup were measured before and after the measuring points. By applying the M3D3 method, a good correction could be achieved, see figure 5(b). The standard deviation σ from the calibration values was only 0.3 µm and the maximum deviation, 0.5 µm. The standard deviation σ was calculated using: where ∆x i and ∆y i are the x and y components of the difference between the calibrated and the M3D3 corrected centre of hole i after the mathematical alignment, and n is the number of holes. The black point cloud (+) in figure 7(b) shows the taskspecific length uncertainties of the M3D3 correction from the Monte Carlo simulation. The main reason for the overall larger uncertainties compared to those of the optimal setup is the fact that LT 1 and LT 2 are not positioned at the height of the plate's plane.
Additionally, a second measurement with worse LT positions (i.e. LT 1 and LT 2 positioned about 800 mm farther away from the plate) was carried out. As expected, this resulted in slightly larger deviations of the M3D3 corrections from the calibration values (standard deviation: 0.4 µm and maximum deviation: 0.8 µm) and, above all, in significantly larger taskspecific length uncertainties (up to a factor of 3 for the larger distances between two holes).

Distribution of the stabilization points
The influence of the distribution of stabilization points was examined based on the measurement of the hole plate standard with the LT positions (black circles) in figure 7(a), the corresponding correction result ( figure 5(b)) and the task-specific length uncertainties (black +) in figure 7(b), respectively. In this basic setup, stabilization points were arranged at 32 locations on a cuboid frame, each measured twice (in reversed order), resulting in 64 distance measurements. To investigate the effect of distribution, subsets of locations were selected and measured several times each in the simulation to maintain an overall number of distance measurements of 64 for comparability.
The arrangement of stabilization points on a frame around the measuring points turned out to be a good choice for common measuring tasks and can be recommended as a default. For the plate measurement considered here, the uncertainties could be further reduced by a more specific arrangement. A subset of the 16 stabilization points near the four LTs (see the setup (red diamonds) and the associated length uncertainties (red circles) in figure 8) has been found from a series of arrangements. The choice of the subset allowed the uncertainties to be slightly reduced as can be seen by comparing the basic uncertainty (black + in figure 7(b)) to the red circles in figure 8(b).
One set, mirroring this arrangement, turned out to be the worst selection from the considered series of arrangements in terms of uncertainties and the quality of the correction (see black squares in figure 8(a)). Here, the stabilization points have a large distance to the LTs. The length uncertainties with this subset are massively increased compared to the basic arrangement on the frame, as can be seen (black crosses) in figure 8(b). In this case, the standard deviation from the calibration values as a measure of the quality of the correction is 0.5 µm, which is also significantly higher than that of the basic arrangement.

Number of the stabilization points
Investigations into the number of stabilization points were carried out by simulations based on the setup regarding line measurements (see section 3.1). In this case, the multilateration problem is solvable even without any stabilization points, since the measuring points on the four lines are sufficiently distributed within the test volume. Nevertheless, the usage of stabilization points is strongly recommended, as even a small set of stabilization points leads to a significant reduction of the uncertainties (e.g. compare figure 4(b) with 32 and figure 9 without any stabilization points).
The basis for further investigations is the cuboid volume that spans from the minimal to the maximal values of the coordinates of all measuring points. The uncertainties were calculated by simulation with a certain number of stabilization points on a frame around the volume and with the same number of points randomly distributed within the volume, respectively. The stabilization points were equidistantly arranged on the edges of the frame. A minimum number of eight stabilization points in the corners of the frame were used. Successively, up to nine further points were added on each of the 12 edges between the corners, which results in an overall number of 116 stabilization points. For each number of points, 100 random sets were generated for evaluation.  In summary, it can be said that an increase in the number of stabilization points steadily leads to smaller uncertainties. This applies both to the frame and, on average, to the randomly distributed stabilization points, at which the frame configuration for all numbers of points yields smaller uncertainties than the randomly arranged stabilization points.
However, while eight stabilization points already reduce the uncertainty significantly from (in case of the frame) about 4 µm for zero stabilization points to below 1.3 µm (correspondent to a factor of about 3), the effect of many more stabilization points is rather small. The maximum number of 116 stabilization points only leads to a further reduction in uncertainty of about five percent.
Random configurations with fewer points have a much greater spread of achieved uncertainties than those with more points. This is obvious, as the chances of drawing an optimal configuration increase with the number of points.
The idea might come up to try to optimize the total measuring time for an acceptable measurement uncertainty by using a certain number of stabilization points. This cannot be done in general, but for a certain measuring task. But since in most measuring tasks there are much more measuring points than stabilization points, there will not be much scope for optimization of the measuring time.
More important might be the question of whether just eight stabilization points can be considered satisfactory, since the benefit of more points is only a few percent. However, experience has proved that it is advantageous to use at least a number of stabilization points so that the MLAT problem is solvable even without any real measuring points. Good indications of the stability of the setup can be obtained by measuring before and after the actual measuring points. A total of 32 points has turned out to be a good compromise in terms of measurement time and achievable uncertainty.
In conclusion, the frame configuration of the stabilization points proves to be superior compared to randomly distributed stabilization points.

Number of LTs
The usage of additional LTs in terms of cost and effort seems to be disproportionate to its benefit. If an additional LT is available, it may be useful to include it in the system in order to have the possibility for additional checks due to the redundancy. Nevertheless, to complete the investigation, the influence of additional LTs on the uncertainties was studied by performing simulations. As a basis, the simulation setup for the hole plate standard in figure 7(a) was used: two LTs are placed at height of the hole plate, one above and one below. Again, the uncertainties were determined using Monte Carlo simulations with 1000 iterations.
To fulfil the aim of minimizing the x and y components of the measuring point uncertainties, the optimum information and the best intersection conditions are provided by LT positions at the height of the plate. Nevertheless, it should be noted that for solvability of the M3D3 algorithm, LTs at a different height, i.e. above and below the plate, are also needed. See also section 4.1. Therefore, the best simulation results were achieved when additional LTs were added at this z position. As can be expected, simulations showed that a successive increase in the number of LTs in this optimal z position continuously reduced the uncertainties. As an example, a setup with seven LTs (five of them at the height of the plate, one below and one above) reduces the uncertainties by about 0.1 µm (corresponding to a reduction of about 25%) for the maximum lengths, compared to the basic setup with four LTs.

Application and verification of the method
The feasibility and performance of the method was verified through measurements of a large involute gear ring standard of about 2 m in diameter with three sets of three external and internal teeth with different helix angles (figure 11) [16]. For the verification measurements, only one set of internal and external teeth was selected. The reference values of this standard are well known from former calibrations and interlaboratory comparisons.
Due to the size and shape of the gear ring, this measurement requires the use of a multi-stylus configuration, i.e. a stylus in the −z direction for the determination of the workpiece coordinate system and a stylus in the x direction to measure the gears (see figure 12). This leads to certain challenges in the M3D3 evaluation that will be addressed in the following section.

Tasks with multiple stylus offsets
The multilateration approach for 3D coordinate acquisition basically works with a single retroreflector and, in the best case, without beam interruption in one measuring cycle. For tactile measurement, this would mean the use of only one stylus per measuring task. However, for measuring tasks on real workpieces, e.g. the large gear ring standard used, this is rarely feasible in practice. When measuring complex parts, it is often necessary to determine the workpiece coordinate system and features to be measured with different styli, having different offsets. The respective retroreflectors have to be mounted afterwards with the same offsets, which can only be realized with limited accuracy. This results in an uncertainty contribution, as the local error correction refers to a slightly different position as is desired for correcting the tactile measured point. Additionally, with different offsets involved, it leads to slightly shifted M3D3 point clouds when evaluating the offset-related measurements separately. Furthermore, the relative orientation of the point clouds cannot be determined correctly. However, this is essential as it contains orientation information caused by rotation errors, the magnitude of which depends on the stylus offset.
In principle, the use of multiple related styli within a test plan should be avoided wherever possible. This holds independently of the application of the M3D3 method, especially if the highest accuracy is required. If this is not possible for practical reasons, the stylus offsets should be as short as possible, and for the multilateration measurement, the offset of the reflector should be as close as possible to the tactile probing tip. The probing system qualification should be undertaken as closely as possible to the place of use. Any difference will contribute to the measurement uncertainty, especially due to the rotational errors of the last kinematic axis. To illustrate this effect and to estimate the magnitude of the uncertainty contribution, a multi-stylus measurement on the setup as shown in figure 13 has been simulated. Differences in the range of approximately 10 mm for the stylus offset (probing tip vs. reflector) and linear rotational errors in the range of 20 µrad m −1 have been assumed. The resulting simulated error vectors are shown in figure 14.
The application of the M3D3 correction procedure requires the following steps (assuming two stylus offsets): (1) Perform tactile measurement as usual.
(2) Define one common point P align to be measured by multilateration using each reflector offset l = 1, 2 (resulting in measured points P align l ). This is the point where the interferometric measurements will be aligned. It should be as near as possible to the point where the alignment of the tactile styli has been undertaken, so that the same rotation errors are effective.   (4) Solve the multilateration problem for stylus one using the measured points for alignment. (5) Solve the multilateration problem for stylus two using the LT positions from (4) for alignment. (6) Shift the point cloud from (5) by the offset between (4) and (5) in order to make P align 1 = P align 2 . (7) In order to obtain the correction vectors in one common coordinate system, solve the multilateration problem for the complete data set (l = 1, 2). A new dead path must be set at the transition between the two reflectors, as the information for the incremental length measurement is lost at this point. The resulting point cloud is consistent in itself. It may be slightly rotated relative to the machine coordinate system, but this has no influence on the final measurement result. (8) Apply the local correction vectors to the originally measured points. (9) As a last step, re-import the points that have been corrected as detailed above into the part program and re-evaluate with all steps for the workpiece.
The reference of the point clouds with different reflector offsets is established via the fixed LT positions. It is obvious that care must be taken to ensure that they do not change during the measuring process, e.g. due to thermal drift or collisions.
Following the steps as proposed above to evaluate the multilateration problem with two stylus offsets reduces the residual errors considerably, as can be seen in figure 15 for the example of the large gear ring.

Practical application
The measurement was performed on PTB's large CMM with the numerical error compensation of the guideways enabled. For the tactile measurement of the gear ring standard, two probes were used as stated above.
The M3D3 measurement was carried out using five LTs. They were positioned according to the presented research results to find the best LT positions for a measuring task (see section 4.1) within the restrictions given by the available equipment (e.g. stands for the LTs). Four LTs were set up at the same height (z ≈ 300 mm) on a slightly curved line. The outer LTs were positioned as far apart as possible to enclose a large angle to the measuring points. The fifth LT was set up in the middle at a height of z ≈ 1600 mm (see also figure 13).
The stabilization points were chosen to form a cuboid volume of 32 points enveloping the whole gear ring standard (see sections 4.2 and 4.3). The stabilization point configuration is also shown in figure 13 (red diamonds).
Since two different probes were used for the determination of the workpiece position and the measurement of the gears, all measuring points (including the stabilization points) were measured with two reflector offsets similar to the offsets of the styli used for the tactile measurement. To obtain the M3D3 corrections, the data (two point clouds) was evaluated following the procedure as suggested in section 5.1.
The results (profile and helix slope deviations of all gears) of the different evaluations were compared to the calibration values of the large gear ring standard (from 2015, [16]), see figure 16 as well as table 1. The contribution of the M3D3 procedure to the task-specific uncertainties was calculated by means of Monte Carlo simulations with 200 iterations. The measured lengths of the LTs were slightly varied by normally distributed random numbers in the size of the uncertainties of the LTs before the multilateration was performed. The evaluation of the profile and helix deviations for every iteration was carried out using the original part program.
The M3D3 corrections for the measurement of the profile and helix slope deviations are quite small but, nevertheless, with very promising small task-specific uncertainties. It should be emphasized again that another advantage of the M3D3 method is the traceability of the measuring results that is ensured by the calibration of the frequency of the laser.

Summary and outlook
An in-depth study of PTB's M3D3 system gives insight into the main uncertainty contributors and helps to find the    ) and helix ( f Hβ ) slope deviation of the M3D3 measurement (with five LTs) of the large gear ring standard (external and internal teeth). U cal is the uncertainty of the calibration, U M represents the MLAT contribution to the task-specific measurement uncertainty.
Profile slope deviation f Hα : Helix slope deviation appropriate setup for any specific application used to improve the performance of CMMs or machine tools as well as ensuring traceability of the measurement. The presented studies quantifying the influence of LT positions as well as the number and arrangement of stabilization points allow an a priori uncertainty estimation. The results can directly be implemented in industrial precision engineering processes, which will benefit from enhanced metrology solutions. For multi-stylus applications, it is recommended that the presented procedure is followed step by step as described in order to minimize systematic errors. The work is embedded in a comprehensive research project on the establishment of a Competence Center for Wind Energy at PTB [17] and complements previous publications in the field. The coordinate metrology part of this project aims to improve the performance of large CMMs by means of various technologies. Ongoing studies focus on the evaluation of a complete task-specific measurement uncertainty that, besides the M3D3 contribution as described here, includes other relevant uncertainty contributions from the CMM, the tactile probing process, the environmental conditions and the workpiece itself by using the VCMM.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.