## Abstract

In the acquisition stage of many space applications, such as the Taiji program, the spot center of weak laser light must be accurately determined. Under weak light conditions, the precision of most traditional positioning methods is greatly affected. In this paper, we present a high-precision laser spot center positioning method based on the theoretical analysis of influence factors of precision. It is shown through experimental study that the method’s precision can fulfill the requirement of the Taiji program.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

The discovery of gravitational waves by the Laser Interferometer Gravitational-Wave Observatory (LIGO) [1] promoted the Chinese to propose a space-based gravitational wave detection program called the Taiji program [2,3] to receive gravitational wave signals from more sources. For the Taiji program, three satellites will be used to form an equilateral triangle laser link constellation, whose arm length is $ 3 \times {10^6}\;{\rm km} $. With such a long distance, the transmitting laser beam can only be detected by the receiving satellite with the help of the navigation data. Therefore, a constellation acquisition process is necessary before the science data detection. Like NASA’s Laser Interferometer Space Antenna (LISA) [4–7], the Taiji program proposes to use star trackers and CCD detectors as well as four-quadrant photodetectors (QPD) to realize laser signal acquisition. Even after star trackers are used, the position uncertainty cone of the remote satellite is still larger than the size of the laser beam. Therefore, a scanning maneuver [8] must be performed to cover the whole uncertainty cone. While the sending satellite is actively scanning, the receiving satellite remains in its reference attitude and stares into the corresponding reference direction based on the navigation data. As the field of view (FOV) of the CCD is far larger than the QPD, capturing the signal on the CCD first rather than directly on the QPD will greatly reduce the operating time. At a certain time, the receiving satellite detects the laser spot on its CCD with an offset from the expected reference position due to the ground-provided reference attitude errors. The offset must be computed by accurately determining the center of the receiving laser spot on the CCD surface and corrected by the attitude actuators. Because the residence time of one scanning spot is short, the laser spot centering process should be a real-time process. To make the beam enter the FOV of the QPD, the angular uncertainty induced by the positioning error should not exceed 0.06 μrad [6]. For a CCD with $ 640 \times 512\;{\rm pixel} $, the requirement of the laser spot center positioning precision should be no less than 0.1 pixel. On the other hand, the receiving laser intensity is nearly 100 pW at the CCD surface due to the long propagation distance. Therefore, a high-precision, real-time method applicable to the weak light center positioning method is essential for the Taiji program.

Many algorithms of laser spot center location have been proposed for various areas. The traditional centroid method is used in many areas, such as Hartmann–Shack wavefront sensor [9–11] and laser range sensor [12]. The traditional centroid method is simple and practical for real-time applications, but it usually has low accuracy. For the large 3D surface profile measuring system of laser scanning, Yaoquan Yang proposed the Hough transfer spot centering method [13]. It is effective for applications with laser spot intensity unevenly distributed. However, it can hardly be used in real-time applications because of the large amount of calculation. Rather than traditional differential power sensing method, an alternative positioning formula with the laser intensity of each phase of quadrant detectors is proposed by Song Cui to improve the measurement accuracy [14,15]. However, the photocurrent intensity is only 10 pA magnitude in the Taiji/LISA program. Equivalent current noise of a low noise detector is usually 1 nA magnitude [16]. Therefore, methods with quadrant detectors can hardly be used for weak light conditions. For the Large Sky Area Multi-Object Fiber Spectroscope Telescope (LAMSOT), Lili Wang [17] used the Gaussian fitting centering method, which can reach the precision of 0.01 pixel for a Gaussian spot with little distortion. It is a good choice for applications with a perfect Gaussian beam [18,19]. However, in the Taiji program, because of the long transmitting distance, the receiving beam is a flat-top beam clipped by a telescope rather than a Gaussian beam. With a lens system in front of the CCD, the intensity distribution of the laser spot on the detector subjects to Fraunhofer diffraction distribution. As a result, the Gaussian fitting method can not work well. By contrast, the traditional centroid method is more applicable to the Taiji program if the precision can be improved. Therefore, in this paper, we introduce what we believe is an improved centroid method, which can not only adapt to the real-time requirement but also achieve high accuracy under weak light conditions.

The paper has four sections. In Section 2, we focus on the reasons that limit the precision of the traditional centroid method under weak light conditions. Based on the theoretical analysis, our improved centroid method is presented. In Section 3, we carry out an experiment with a 100 pW receiving laser. We first verify the validity of the theoretical results. Then the performance of the improved centroid method is presented. Finally, we summarize the conclusion in Section 4.

## 2. THEORETICAL ANALYSIS

The traditional centroid method is one of the most efficient methods in the calculation of a laser spot center position. Because the algorithm is computationally simple, it can cope with the real-time requirement of the Taiji program. However, it is sensitive to noises and usually has low precision especially under weak light conditions. In this section, we first reveal how the noises influence the precision.

The centroid position in the horizontal direction (the same form in the vertical direction) on a CCD detector is defined as

where $ {x_i} $ denotes the position of the $ i $th pixel, $ {N_i} $ stands for the detected photo events, and $ n $ is the total pixel number involved in calculation. For the traditional centroid method, $ n $ is the total pixel number of the CCD detector. As $ {N_i} $ contains not only photoelectrons inspired by the incident light ($ {N_{s\!i}} $) but also noise electrons ($ {N_{oi}} $), the centroid position can be rewritten asOne of the greatest contributors to $ {N_{oi}} $ is the background noise of the CCD, including the dark current noise, the charge transfer noise, and the readout noise. To guarantee a stable temperature environment, no cooling mechanism will be used for the CCD in the Taiji program. As a result, the dark current noise will dominate the noise budget. The number of dark current noise electrons approximately obeys Gaussian distribution. We suppose that $ {N_{oi}} \sim N(\mu ,{\sigma ^2}) $, where $\mu$ and $ {\sigma ^2} $ stand for the expectation and the variance of the noise electrons, respectively.

Because $ {N_{oi}} $ is independent of each other for $ i = 1,2, \ldots ,n $, $ E(x) $ can be written as

To simplify the last term in Eq. (4), we define the image area involved in calculation as the computational domain (for the traditional centroid method, it is the whole image). We also suppose that the column number of the computational domain is $ {n_1} $ ($ {n_1} $ is an odd number for simplification), the row number is $ {n_2} $, and the center of the area is $ x_0^\prime $ ($ x_0^\prime $ is an integer). Then,

On the other hand, the noises will also affect the positioning precision by causing a fluctuation of center position around $ E(x) $. After dividing the small computational domain, the number of photoelectrons is far larger than the number of noise electrons. The formula from Eq. (2) can be rewritten as

We can estimate the fluctuation intensity by the root mean square value of $ \Delta x $. Because the number of the noise electrons of each pixel is independent of each other and obeys the Gaussian distribution, there are many similarities between $ {N_{oi}} $ and the Brown motion process [20]. With $ {N_s} \gg \mu$, we rewrite the expression of $ \Delta x $ as

whereWith $ \tau = 1/n $, we know that $ N_{oi}^\prime \sim N(0,\tau ) $, which has the same form with the increment of the standard Brown motion process $ \Delta B(i\tau ) $. Then, we define $ X(k\tau ) $ as

We also introduce $ \Delta X(i\tau ) $ as the increment term, so

Even in the small domain, the pixel number involved in the calculation is still large enough to approximately regard the differential value of $ B(i\tau ) $ and $ X(i\tau ) $ as its difference value. Therefore, we can rewrite Eq. (8) as $ (u = i\tau ) $, so

It can be seen that $ \Delta x $ is a Brown motion process with an interference term that is an Ito process [21,22]. With the statistical properties of the Ito process, we can easily obtain the expectation and the variance of $ \Delta x $, so

Based on the relationship in the vertical direction $ ({x_1} = {x_{\sqrt n + 1}} = {x_{2\sqrt n + 1}} = \cdots ) $, the expression of $ {x_i} $ can be given by

where $ \lfloor\frac{{i - 0.1}}{{\sqrt n }}\rfloor $ is the rounded down number of $ \frac{{i - 0.1}}{{\sqrt n }} $. Here we suppose that $ \sqrt n $ is an integer. Because of the rounding function, we still can hardly get an analytical result. Therefore, we introduce $ {X_i} $ as the sum of the $ i $th line of $ \Delta x $. Then, $ D({X_1}) $ can be given byIf we consider only the highest order term, the approximate expression of $ D({X_1}) $ is written as

Similarly, when the other lines are considered separately, we can also obtain

Because $ {X_1},{X_2}, \cdots {X_{\sqrt n }} $ are independent of each other, we can get

As a result, the root mean square value of $ \Delta x $ can be denoted as

We can conclude that the fluctuation intensity is approximately proportional to the pixel number $ n $ involved in calculation. Therefore, dividing a small computational domain around the real centroid position will also reduce the position fluctuation. Based on the following analysis and the actual working condition of the Taiji program, here is an improved centroid method with three steps:

- 1. Divide a rough computational domain.
As concluded above, a smaller computational domain will greatly reduce the pixel number involved and increase the positioning precision. First, we can obtain the rough domain center by searching for the pixel $ ({x_1},{y_1}) $ with the maximum gray value $ m $. Then, we search for the pixel $ ({x_2},{y_2}) $ which has a number of electrons nearest to $ m/2 $. The distance between $ ({x_1},{y_1}) $ and $ ({x_2},{y_2}) $ is $ r $. To cover nearly all the photoelectrons, the computational domain is defined as a square area whose center is $ ({x_1},{y_1}) $ with side length of $ d = 4r $.

- 2. Find the rough spot center location.
In the computational domain defined above, we use the centroid method and get the rough spot center $ ({x_{{\rm measure}}},{y_{{\rm measure}}}) $.

- 3. Obtain the precision spot center.
Because the domain center should be an integer, we define $ ({x_3},{y_3}) $, which is the nearest integer point of $ ({x_{{\rm measure}}},{y_{{\rm measure}}}) $, as the new domain center. Then, we use the centroid method again and achieve a more accurate center position.

Because we can hardly know the real spot center position, the fixed measurement offset can only be estimated by the maximum possible value of $ \Delta {x_0} $. The reason why we use the centroid method twice is that $ \Delta {x_0} $ has the potential to become smaller at the second time. Then we theoretically verify the feasibility of the method with the real parameters in the Taiji program. The diameter of the laser spot is set to about 7 pixels with the help of the lens system. As a result, the pixel number of the computational domain is approximately $ n = 15 \times 15 \;{\rm pixels} $. The property of noise electrons of the SH640 camera is $\mu= 1.5 \times {10^5}{e^ - },\sigma = 4.5 \times {10^3}{e^ - }{@25^ \circ }{\rm C} $. In part 1, because of the quantization error, the fluctuation of noise electrons, and the influence of diffraction, $ ({x_1},{y_1}) $ may be near but deviate from the real spot center. As the intensity distribution subjects to Fraunhofer diffraction distribution, the intensity changes a lot around $ ({x_1},{y_1}) $. On the other hand, within the pixels around the center, the number of photoelectrons is far larger than the number of noise electrons. Therefore, the real center position can only exist between the two pixels adjacent to the pixel with the maximum gray value. The maximum possible value of $ \Delta {x_0} $ is two pixels. For a 100 pW receiving laser, we can get $ \Delta \le 0.164\;{\rm pixel} $, $ \sqrt {\overline {\Delta {x^2}} } \approx 0.0016\;{\rm pixel} $. The precision of center positioning can be estimated by the sum of $ \Delta $ and $ \sqrt {\overline {\Delta {x^2}} } $. Then, the offset between $ {x_{\rm measure}} $ and the real position is smaller than 0.17 pixel. Because the redivided area center $ {x_3} $ is an integer, as shown in Fig. 1, the maximum value of $ \Delta {x_0} $ reduces to 0.67 pixel in part 3. As a result, the maximum value of $ \Delta $ reduces to 0.055 pixel. Therefore, the theoretical precision of the improved centroid method can be estimated as $ \Delta + \sqrt {\overline {\Delta {x^\prime}^2} } \approx 0.057\;{\rm pixel} $. As the Taiji program calls for 0.1 pixel positioning precision, we believe our method has the potential to fulfill the requirement. Although we can repeat step 3 and finally make $ \Delta {x_0} $ approximately equal to 0.5 pixel, it will not have a great influence on the precision. In the next section, we will present an experiment to prove the rationality and feasibility of the method.

## 3. EXPERIMENTAL RESULTS

The experimental system is shown in Fig. 2. It includes a 1064 nm Nd:YAG laser with an output of about 8 mW, an attenuation system, lenses, and a TEKWIN SH640 camera. The lens group can greatly reduce the spot size and the influence of the laser jitter. The attenuation system consists of three neutral density attenuators whose maximum attenuation rate is 1000 times. Therefore, we can obtain a laser spot of 100 pW on the CCD surface. Figures 3(a) and 3(b) are the collected images without and with the lens group, respectively, for a 100 pW input laser. The diameter of the laser spot in the Fig. 2(b) is about seven pixels.

To verify the feasibility of our improved centroid method, we have to estimate the value of $ \Delta $ and $ \sqrt {\overline {\Delta {x^\prime}^2} } $, respectively.

#### A. Verify the Validity of the Expression of the Fixed Measurement Offset

Because the real spot center can hardly be known, we cannot obtain the accurate value of $ \Delta $ directly. However, we can measure the variation of $ \Delta $ from experiments. Then we can verify the validity of the formula in Eq. (6) by the relationship between $ \Delta $ and $ \Delta {x_0},n $ with the experimental results.

- 1. The relationship between $ \Delta $ and $ \Delta {x_0} $
From Eq. (6) we know that the value of $ \Delta $ is directly proportional to $ \Delta {x_0} $. The theoretical relationship between them is shown in Fig. 4(a). As a result, the variation of $ \Delta $ is also proportional to the variation of $ \Delta {x_0} $. In Fig. 2(b), we suppose that the pixel with the maximum gray value $ ({x_{{\rm max}}},{y_{{\rm max}}}) $ is the reference position and define it as the center of the computational domain. To include nearly all the pixels with photoelectrons, the side length of the domain should be large enough. Here, it is set to 200 pixels. We use the centroid method and get the reference spot center $ ({x_{{\rm ref}}},{y_{{\rm ref}}}) $. To simplify the analysis, only the center in the x-direction is changed because changes in the y direction will achieve the same conclusion. Then we just change the center of the domain to $ ({x_{i0}},{y_{{\rm ref}}}) $ and recalculate the spot center as $ ({X_i},{Y_i}) $. The change of $ \Delta $ can be denoted as $ d\Delta (i) = {X_i} - {x_{\rm ref}} $, and the change of $ \Delta {x_0} $ can be denoted as $ d\Delta {x_0}(i) = {x_{i0}} - {x_{{\rm max}}} $. The relationship between $ d\Delta {x_0} $ and $ d\Delta $ is shown in Fig. 4(b). The results verify that the fixed measurement error is directly proportional to the offset of the domain center.

- 2. The relationship between $ \Delta $ and $ n $

Then, we try to find the relationship between $ \Delta $ and $ n $. Theoretically, the relationship curve is shown in Fig. 5(a) with $\mu = 1.5 \times {10^{(5)}} $ and $ \Delta {x_0} = 2\;{\rm pixel} $. When $ n $ is small, $ \Delta $ is approximately proportional to $ n $. And when $ n $ is big enough, $ \Delta $ approaches a constant. Similarly, we can verify this conclusion with the variation of $ \Delta $. To include as many photoelectrons as possible, we also suppose that the pixel with the maximum gray value $ ({x_{{\rm max}}},{y_{{\rm max}}}) $ is the center of the computational domain in the Fig. 2(b). First, we use the centroid method with $ n = 15 \times 15\;{\rm pixel} $ and get the reference spot center $ (x_{{\rm ref}}^\prime,y_{{\rm ref}}^\prime) $. We gradually increase the side length of the domain $ \sqrt {n(i)} $ and calculate the center position $ (X_i^\prime,Y_i^\prime) $. $ {\Delta _n}(i) $ is defined as the difference between $ X_i^\prime $ and $ x_{{\rm ref}}^\prime $. The difference between $ {\Delta _n}(i) $ and the real fixed measurement value $ \Delta (i) $ is a constant value. Then the relationship between $ {\Delta _n} $ and $ n $ is shown in Fig. 5(b). The difference between Figs. 5(a) and 5(b) is mainly due to different value of $\mu$ and $ \Delta {x_0} $ in the simulation and the experiment as well as the above constant value. As the accurate value of $\mu$ and $ \Delta {x_0} $ is unknown in experiment, only the trend of the relationship curve is important. The two curves in Fig. 5 have the same trend. Therefore, the relationship curve is consistent with the theoretical conclusion.

Therefore, the validity of the formula in Eq. (6) is verified by the experiments. The results in the last section can be directly used. For a 100 pW laser, the maximum value of the fixed measurement offset can be estimated as $ {\Delta _{100\,{\rm pW},{\rm max}}} = 0.055\;{\rm pixel} $ with the help of our improved centroid method.

#### B. Verify the Precision of the Improved Centroid Method

As concluded above, the positioning precision is influenced not only by the fixed measurement offset but also the location fluctuation. Therefore, we also need to estimate the fluctuation intensity. We capture 500 frames of sequential images with the laser intensity of 500 pW, 200 pW, and 100 pW, respectively. Because they will come to the same conclusion, only the experiment with the 100 pW receiving light is presented. Figure 2(b) is one of the images. The laser spot center location $ ({x_{{\rm image},i}},{y_{{\rm image},i}}) $ can be obtained for the image $ i(i = 1,2, \cdots ,500) $ with the help of the traditional method and our improved centroid method, respectively. Similarly, only the center in the x direction is considered. The results are shown in Figs. 6 and 7. The average value of the center position is expressed as

Therefore, the average values of the center position calculated by the two positioning methods are $ {\overline x _{{\rm traditional}}} = 326.939\;{\rm pixel} $ and $ {\overline x _{{\rm improved}}} = 121.825\;{\rm pixel} $, respectively. The two values are quite different from each other due to the influence of the fixed measurement offset.

Then, the fluctuation intensity can be estimated by the root mean square value of the 500 experiment results and expressed as

In addition to the experiment results above, the location fluctuation of the two methods is $ \Delta {x_{rms,\!{\rm traditional}}} = 0.0437\;{\rm pixel} $ and $ \Delta {x_{rms,\!{\rm improved}}} = 0.0032\;{\rm pixel} $, respectively. It can be seen that the improved centroid method also can reduce the location fluctuation, which coincides with the theoretical analysis. Finally, the positioning precision is $ {\Delta _{100\;{\rm pW},{\rm max}}} + \Delta {x_{rms,\!{\rm improved}}} = 0.0582\;{\rm pixel} $. The precision is less than 0.1 pixel. The small offset between the experimental result and the theoretical result is mainly generated from stray light and an approximation error. Therefore, the feasibility of the improved centroid method under weak light conditions is verified by the experiment.

## 4. CONCLUSION

Because the acquisition phase of the Taiji program requires high positioning precision under weak light conditions in a short time, traditional positioning methods cannot be directly used. In this paper, we present an improved centroid method and verify its feasibility under the practical work conditions of the Taiji. We first analyze the factors that influence the precision of the traditional centroid method. Based on the theoretical analysis, we conclude that the positioning precision is affected not only by the fixed measurement offset but also the position fluctuation. Both increase with the pixel number involved in the calculation. Therefore, the core concept of our improved centroid method is dividing a small computational domain around the real center position. Then, we do an experiment to verify the feasibility of the method. For a 100 pW receiving laser, the positioning precision reaches 0.0582 pixel with the help of the improved centroid method. The results fulfill the requirement of the Taiji and verify the feasibility of the method. This method also has the potential to be used in other programs, which call for high positioning precision under weak light conditions.

## Funding

Chinese Academy of Sciences (XDA1502070304, XDA1502070902, XDB23030000).

## Disclosures

The authors declare no conflicts of interest.

## REFERENCES

**1. **B. P. Abbott, S. Jawahar, N. A. Lockerbie, and K. V. Tokmakov, “Gw150914: first results from the search for binary black hole coalescence with advanced LIGO,” Phys. Rev. D **93**, 122003 (2016). [CrossRef]

**2. **W. Hu and Y. Wu, “The Taiji program in space for gravitational wave physics and the nature of gravity,” Nat. Sci. Rev. **4**, 685–686 (2017). [CrossRef]

**3. **G. Jin, “Program in space detection of gravitational wave in Chinese academy of sciences,” in *11th International LISA Symposium* (2017), Vol. 840.

**4. **P. G. Maghami, T. T. Hyde, and J. Kim, “An acquisition control for the laser interferometer space antenna,” Class. Quantum Grav. **22**, S421–S428 (2017). [CrossRef]

**5. **F. Cirillo and P. F. Gath, “Control system design for the constellation acquisition phase of the LISA mission,” J. Phys. Conf. Ser. **154**, 012014 (2009).

**6. **T. di Laurea, “Controller design for the acquisition phase of the LISA mission using a Kalman filter,” Ph.D. thesis (University of PISA, 2007).

**7. **T. T. Hyde, P. G. Maghami, and S. M. Merkowitz, “Pointing acquisition and performance for the laser interferometry space antenna mission,” Class. Quantum Grav. **21**, S635–S640 (2004). [CrossRef]

**8. **F. Ales, P. Gath, U. Johannn, and C. Braxmaier, “Line of sight alignment algorithms for future gravity missions,” in *AIAA Guidance, Navigation, and Control Conference (GNC)* (2015), p. 0094.

**9. **P. M. Prieto and F. Vargas-Martin, “Analysis of the performance of the Hartmann–Shack sensor in the human eye,” J. Opt. Soc. Am. A **17**, 1388–1398 (2000). [CrossRef]

**10. **G. Cao and X. Yu, “Accuracy analysis of a Hartmann–Shack wavefront sensor operated with a faint object,” Opt. Eng. **33**, 2331–2335 (1994). [CrossRef]

**11. **K. L. Baker and M. M. Moallem, “Iteratively weighted centroiding for Shack–Hartmann wave-front sensors,” Opt. Express **15**, 5147–5159 (2007). [CrossRef]

**12. **L. Liu, “Research on spot subdivided location method in laser triangulation measurement,” Comput. Meas. Control **16**, 1396–1398 (2008).

**13. **Y. Yang, “An algorithm to raise the locating precision of laser spot center based on Hough transform,” Acta Opt. Sin. **19**, 040431 (1999).

**14. **S. Cui and Y. C. Soh, “Improved measurement accuracy of the quadrant detector through improvement of linearity index,” Appl. Phys. Lett. **96**, 081102 (2010). [CrossRef]

**15. **S. Cui and Y. C. Soh, “Analysis and improvement of Laguerre–Gaussian beam position estimation using quadrant detectors,” Opt. Lett. **36**, 1692–1694 (2011). [CrossRef]

**16. **F. Cervantes, J. Livas, R. Silverberg, E. Buchanan, and R. Stebbins, “Characterization of photoreceivers for LISA,” Class. Quantum Grav. **28**, 094010 (2011). [CrossRef]

**17. **L. Wang, “Laser spot center location algorithm based on Gaussian fitting,” Appl. Opt. **33**, 985–990 (2012).

**18. **W. Li, “Iteration algorithm of surface fitting in the detection of light-spot position,” Opt. Tech. **30**, 33–35 (2004).

**19. **H. Chen, Z.-H. Yang, P. Guo, Y. C. Zhang, and S.-Y. Chen, “Research of the high precision laser spot center location algorithm,” Trans. Beijing Inst. Technol. **36**, 181–185 (2016).

**20. **B. Hao, “A hundred years of the theory of Brownian motion,” Trans. Beijing Inst. Technol. **36**, 181–185 (2016).

**21. **T. Komatsu, “On the pathwise uniqueness of solutions of one-dimensional stochastic differential equations of jump type,” Proc. Japan Acad. Ser. A Math. Sci. **58**, 353–356 (1982). [CrossRef]

**22. **Q. Xia, “The stochastic orderings of Itô process and its applications,” Ph.D. thesis (Huazhong University of Science and Technology, 2011).