A Review on Zernike Coefficient-Solving Algorithms (CSAs) Used for Integrated Optomechanical Analysis (IOA)

Abstract

An integrated optomechanical analysis (IOA) can predict the response of an optomechanical system to temperature, gravity, vibrations, and other local loadings; thus, the normal operation of instruments under special conditions is guaranteed. Zernike polynomials are the most popular for fitting the IOA-derived mechanical deformation data. By solving the Zernike coefficients of all deformed optical surfaces, the relationship between aberrations and deformations can be further revealed. The process of IOA is summarized in this article. The principles of four primary Zernike coefficient-solving algorithms (CSAs) were expounded, and the corresponding applications are reviewed in detail, including the least squares method, the Gram–Schmidt orthogonalized method, the Householder transformation, and singular value decomposition (SVD). Artificial neural networks (ANNs) trained for solving a similar overdetermined set of equations are also discussed; an innovative Zernike CSA based on a one-dimensional convolutional neural network (1D-CNN) was proposed, emphasizing its potential for Zernike CSA. The feasibility of the neural network method was verified by conducting experiments on the primary mirror of the front reflection system of a space camera. This review can provide references for the precise optimization of IOA.

1. Introduction

Temperature [1], gravity [2], vibration [3], and other factors [4] seriously degrade the optical performance of an instrument and affect normal operations. An integrated optomechanical analysis (IOA) can comprehensively evaluate the impact of external environmental factors at the design stage, and it guides the design optimization of optomechanical structures using the analysis’s results, thereby significantly improving the adaptability of special optical instruments to the environment, such as aerial cameras, star sensors, space telescopes, etc. [5,6,7]. Currently, there are two main common simulation analysis methods: (1) in the traditional simulation analysis, the optomechanical system is simulated and optimized in each design process; (2) in the IOA, the overall optimization and design of the optomechanical system are carried out.

The overall design process of optical instruments primarily includes the following: optical design, optomechanical design, mechanical analysis, thermal analysis, modal analysis, and other links [8]. In traditional simulation analyses, each link is primarily used to satisfy its own indicators and relies on commercial optical, mechanical, and thermal analysis software independently. This analysis mode cannot achieve the balance optimization of optical, mechanical, thermal, and other factors, and the input conditions are limited to the internals of each link, which is challenging for forming mutual restrictions between links; moreover, enhancing the overall performance of the optomechanical system is challenging. Accordingly, Jacob Miller proposed the concept of integrated optomechanical analyses in 1981. During the design process, the response of an optical system and optomechanical structure to external environment changes is considered comprehensively, and the analysis data can also be seamlessly transferred between different designs and analysis software to achieve the overall optimization of the optomechanical system [9]. All current space optical instruments, such as the European Solar Telescope, Hale Optical Probe, Hubble Space Telescope, James Webb Space Telescope, and Keck Optical Probe, require an IOA prior to launch to verify the adaptability of optomechanical systems relative to vacuum and ultra-low temperature conditions [10,11,12,13].

The most important step in IOA is to use polynomials to fit the deformation data of optical surfaces. Zernike polynomials have orthogonality over a unit circle, and their base terms comprise radial and azimuthal variables suitable for the descriptions of circular apertures. Therefore, the Zernike polynomials are the most commonly used form in IOA. Moreover, the Zernike coefficients of deformed optical surfaces can be utilized by ray-tracing software to further analyze the relationship between deformations and aberrations, providing an efficient method of data transfer from mechanics to optics [14]. Accordingly, the accuracy of IOA is heavily dependent on the Zernike coefficient-solving algorithms (CSAs).

In the following sections, the process of IOA and its typical applications are first presented in Section 2. Zernike polynomials and the corresponding CSAs are reviewed in detail in Section 3, including the least squares method, the Gram–Schmidt orthogonalized method, the Householder transformation, and singular value decomposition (SVD). The neural networks trained for solving a similar overdetermined set of equations are also discussed in Section 4, and the potential of using neural networks to solve Zernike CSAs is analyzed. Section 5 contains the summaries on the characteristics of Zernike CSAs and suggestions for future developments.

2. Typical Applications of IOA

Thermal [15] and modal [16] analyses are two typical applications of IOA. Thermal analyses apply thermal loads to a finite element structural model of an optomechanical system based on specific temperature conditions and then calculate the temperature field under the thermal loads. When temperature fields are applied to the structural model, the thermal strain builds up and causes the displacement of finite element nodes. In order to evaluate the optical performance, the nodal displacement data of the optical surface must be converted into sag data. The Zernike coefficients of all optical surfaces are obtained by fitting the sag data with Zernike polynomials. By importing the coefficients into ray-tracing software, it is possible to quantitatively evaluate the effect of temperature changes on optical performance and, thus, provide a basis for the athermalization design of optomechanical systems.

Figure 1 is an IOA flow diagram of a satellite-borne system. According to the optical design and the satellite’s structure, the structural model of the optomechanical system is designed. Temperature analyses are carried out by using orbit parameters and extremal orbit temperatures. The thermoplastic analysis of the structural model and the temperature analysis of the optical surface together obtain the optical surface’s deformation, and the rigid body displacement of the surface is also calculated by thermoplastic analyses. The surface’s deformation is fitted by polynomials, and the fitting results and the rigid body displacement are optically analyzed together. According to the results of the optical analysis, the structure and temperature control design are optimized. A closed-loop process from the design stage to the analysis stage with respect to design optimizations is formed via IOA.

For ground-based telescopes, ambient temperature changes and axial temperature gradient changes caused by direct sunlight are two key thermal factors. The Nanjing Institute of Astronomical Optics and Technology (NIAOT) evaluated the effect of ambient temperatures on the image quality of an astronomical telescope by IOA [17]. The linear interpolation technique derived from shape functions in the finite element theory was applied as an interface for the thermal model and structural model. Zernike polynomials and power series expansion were compared and used for fitting a paraboloid primary mirror with a diameter of 900 mm. Concurrently, Cho et al. [18] performed thermal analyses for the Thirty Meter Telescope’s (TMT) structure. The telescope structural parts were modeled for various thermal conditions, including air convections, conduction, heat flux loadings, and radiation. The thermal responses of the TMT’s structure were successfully predicted for daytime and nighttime periods. The CAD model of the telescope’s structure and finite element model at a Zenith angle of 30 degrees are shown in Figure 2. Deep space solar observatories need to observe the local sun, and its main mirror should not only adapt to the temperature difference between day and night periods but also avoid the influence of heat sources on the axis with respect to optical performances. In response to this complex thermal environment, Li et al. [19] used parabolic polynomials and Zernike polynomials to fit the thermal deformation of the main mirror under the supported conditions. According to the thermal analysis, the supporting system of the main reflector was improved by using a six-point support, which increased the thermal control accuracies.

Compared with ground-based optical systems, space optical systems have a broader range of ambient temperature variations, which poses a new challenge with respect to the temperature adaptability of structures. Due to Mercury’s radiant temperature range of 100 K–700 K, the thermo-elastic analysis of high-spatial-resolution imaging cameras (HRIC) is essential for evaluating the effect of thermo-elastic deformations on optical performances under extreme thermal environmental conditions [20]. Dipasquale et al. used Zernike polynomials to decompose surface deformations induced by the piston/tilt of the two mirrors. A deformed camera with the contour plot of the displacements is shown in Figure 3. The ray-tracing results showed that the movement of the mirrors introduced a slight optical aberration; however, the pointing error can be neglected.

In addition to the above large temperature difference conditions, the stability and uniformity of the ambient temperature were strictly demanded by the optical instruments used for ultra-high-precision surface measurements. Zhang et al. [21] carried out IOA on the interferometer lens, and the results revealed that the temperature difference should be less than ±0.015 °C along the mirror axis direction. In addition, Tan et al. proposed an annular-type thermal control system for eliminating the excessive heat of a large-aperture primary mirror. Moreover, thermal analyses were conducted by targeting natural and forced convection. The comparative results revealed that the thermal stability of the primary mirror with the control system significantly improved [22].

Modal analyses can predict the response of an optomechanical structure to external dynamic loads, and this is primarily divided into static analysis and vibration analysis; moreover, the analysis process is similar to thermal analysis. Static analyses are used to calculate the response of an optomechanical structure under gravity loads, forced displacements, and inertial loads [23]. For example, space optical instruments are subjected to large loads during the take-off and landing phases of a spacecraft; thus, the displacement and stress responses need to be calculated under weightlessness and overweight conditions, respectively. In contrast, the vibration environment primarily causes relative position changes among the optical components in the optical instrument, including the translation, rotation, and surface shape change in the optical components along the direction of the three axes, and ultimately affects optical performances [24]. In order to find out the reasons why the wavefront error of the space telescope did not achieve the goal, focusing on gravity, flexure bonding, and thermal expansion, Bo-Kai Huang et al. [25] analyzed the structure model of the primary mirror assembly. A new assumption of the bonding effect with zero gravity was proposed and verified by reverse engineering and the synthetic wavefront error method. Lin et al. [26] established a mathematical model based on the theory of the transfer function and line spread function to describe the influence of a carrier aircraft’s vibration on the modulation transfer function (MTF) of an aspherical aerial camera. The effectiveness of the mathematical model was validated by the results of an experimental analysis.

When large ground-based optical systems are designed and validated, it is necessary to consider the influence of modal factors, among which wind-induced vibrations are typical factors. Wind vibration analyses were performed with the baseline configuration of the Giant Magellan Telescope (GMT) structure made by Simpson Gumpertz and Heger Inc. [27]. With the wind pressure and velocity data recorded at the 8 m Gemini South Telescope as the dynamic load on GMT structure, the random response and optical errors were determined. Moreover, the areas that significantly contributed to optical errors were identified by using spectral response curves and mode shapes. The finite element model of the baseline design of the GMT telescope’s structure is shown in Figure 4.

External loads, such as ambient temperature change and vibration, will eventually cause a sag change on the optical surface, leading to the degradation of optical performance. The data interface is an important part of IOA as a bridge between optical surface deformation and optical performance. The polynomials are most commonly used for fitting the sag change data on the optical surfaces in IOA, such as Zernike polynomials, annular Zernike, X-Y, Fourier–Legendre, and aspheric polynomials. By solving the fitting coefficients, wavefront variations can be quantitatively analyzed, and the relationship between deformation data and optical aberration can be established. Compared to other forms, Zernike polynomials became the most commonly used data interface due to the advantages of orthogonality over the continuous unit circle, rotation invariance, and completeness property [28]. In addition, due to the correspondence between Zernike coefficients and aberrations, commercial ray-tracking software can perform optical analysis based on Zernike coefficients without applying conversions. Accordingly, the accuracy of Zernike CSAs directly determines the accuracy of IOA results.

3. Zernike Coefficient-Solving Algorithms (CSAs)

The main Zernike CSAs include the least squares method [29], the Gram–Schmidt orthogonalized method [30], the Householder transformation [31], and singular value decomposition (SVD) [32].

3.1. Least Squares Method

Zernike polynomial fitting is used to describe the deformation of the optical surfaces of arbitrary shapes by using linear combinations of infinite bases. The discretized optical surface deformation is expressed as follows:

$$\Delta s= {\displaystyle \sum _{i=1}^N{a}_i}{z}_i(r,\theta ),$$

(1)

where ∆s is the sag change of each node, N represents the number of Zernike polynomial terms, z(r,θ) denotes the Zernike polynomials in polar form, and a is the corresponding coefficients of the Zernike polynomials. The overall optical surface deformation can be expressed as follows:

$$\left\{\begin{array}{c}\Delta s_1=a_1{z}_1(r_1,{\theta }_1)+a_2{z}_2(r_1,{\theta }_1)+\cdots +a_n{z}_n(r_1,{\theta }_1)\\ \Delta s_2=a_1{z}_1(r_2,{\theta }_2)+a_2{z}_2(r_2,{\theta }_2)+\cdots +a_n{z}_n(r_2,{\theta }_2)\\ \cdots \\ \Delta s_m=a_1{z}_1(r_m,{\theta }_m)+a_2{z}_2(r_m,{\theta }_m)+\cdots +a_n{z}_n(r_m,{\theta }_m)\end{array}\right.$$

(2)

Equation (2) can be transformed into

$$Z\times A=\Delta S,$$

(3)

where Z is an m × n matrix comprising Zernike polynomials, m represents the number of sampling nodes, n represents the number of terms, A is the column vector composed of Zernike coefficients, and ∆S is the column vector comprising surface sag changes. When m > n, Equation (3) is an overdetermined equation [33]. Z is the coefficient matrix of the overdetermined equation.

The Zernike polynomial coefficients can be calculated by solving the overdetermined linear equation. The overdetermined equation has no exact solution. The least squares solution of the overdetermined equation is regarded as an approximate solution, and solution A satisfies ${‖ZA-\Delta S‖}_2^2=\mathit{min}$. ${‖ZA-\Delta S‖}_2^2$ can be decomposed by:

$${‖ZA-\Delta S‖}_2^2={(ZA-\Delta S)}^T(ZA-\Delta S)=A^T{Z}^{T}ZA-2A^T{Z}^T\Delta S+\Delta S^T\Delta S,$$

(4)

$$\frac{\partial {‖ZA-\Delta S‖}_2^2}{\partial A}=2Z^{T}ZA-2Z^T\Delta S=0.$$

(5)

In this case, the minimum value is obtained for Equation (4). Equation (5) can be transformed into:

$$Z^{T}ZA=Z^T\Delta S.$$

(6)

Equation (6) is known as the normal equation [34]. According to this condition, Zernike coefficients are solved by:

$$A={(Z^{T}Z)}^{-1}{Z}^T\Delta S.$$

(7)

Wang et al. [35] analyzed the influence of temperature on the MTF of an infrared membrane diffraction optical system, and the least squares method was adopted to solve the Zernike coefficients, thus establishing the thickness distribution model of the membrane diffractive optical components. The established model was helpful for estimating the diffraction efficiency and the MTF values for various fields of view. SAGUARO is an open source software of IOA developed by the Large Optics Fabrication and Testing group at the University of Arizona, in which the least squares method is used to fit Zernike polynomials [36]. Banyal et al. [37,38] carried out IOA on solar telescope mirrors and used SAGUARO to fit the optical surface; the fitting result is shown in Figure 5. Song et al. [39] examined the effects of an irregular refractive index on optical performances and adopted the least squares method to determine the best gradient equation. The experimental results showed that the influence of the irregular refractive index effect on optical performances is one-third to one-seventh that of the thermoelastic deformation effect.

3.2. Gram–Schmidt Orthogonalized Method

Since IOA needs to discretize the optical surface, Zernike polynomials lose their orthogonality on the continuous unit circle. In this case, the normal equation of the overdetermined equation is ill-conditioned when the scale of the coefficient matrix is large, and the small change in the coefficient matrix leads to a drastic change in the solution vector, and the solution obtained by the least squares method becomes unstable. The orthogonality of Zernike polynomials can be improved by using the Gram–Schmidt orthogonalized method. Using the Gram–Schmidt orthogonalized method, the coefficient matrix is orthogonalized to realize the transformation of the overdetermined equation [40], and the least squares method is continually used to solve the overdetermined equation. Ill-conditioned problems due to large coefficient matrices are avoided.

Let ${\overrightarrow{\alpha }}_m=(z_1(r_m,{\theta }_m),z_2(r_m,{\theta }_m),\cdots ,z_n(r_m,{\theta }_m))$. The linearly independent vector $({\overrightarrow{\alpha }}_1,{\overrightarrow{\alpha }}_2,\cdots ,{\overrightarrow{\alpha }}_m)$ is transformed into orthogonal vector $({\overrightarrow{\beta }}_1,{\overrightarrow{\beta }}_2,\cdots ,{\overrightarrow{\beta }}_m)$ by using the Gram–Schmidt orthogonalized method. The orthogonalization process can be expressed as follows:

$$\left\{\begin{array}{l}{\overrightarrow{\beta }}_1={\overrightarrow{\alpha }}_1\\ {\overrightarrow{\beta }}_2={\overrightarrow{\alpha }}_2-\frac{〈{\overrightarrow{\alpha }}_2\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}{\overrightarrow{\beta }}_1\\ \cdots \\ {\overrightarrow{\beta }}_m={\overrightarrow{\alpha }}_m-\frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}{\overrightarrow{\beta }}_1-\cdots -\frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_{m-1}〉}{〈{\overrightarrow{\beta }}_{m-1}\cdot {\overrightarrow{\beta }}_{m-1}〉}{\overrightarrow{\beta }}_{m-1}\end{array}\right.$$

(8)

Equation (8) can then be expressed as follows:

$$\left\{\begin{array}{l}{\overrightarrow{\alpha }}_1={\overrightarrow{\beta }}_1\\ {\overrightarrow{\alpha }}_2={\overrightarrow{\beta }}_2+\frac{〈{\overrightarrow{\alpha }}_2\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}{\overrightarrow{\beta }}_1\\ \cdots \\ {\overrightarrow{\alpha }}_m={\overrightarrow{\beta }}_m+\frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}{\overrightarrow{\beta }}_1+\cdots +\frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_{m-1}〉}{〈{\overrightarrow{\beta }}_{m-1}\cdot {\overrightarrow{\beta }}_{m-1}〉}{\overrightarrow{\beta }}_{m-1}\end{array}\right.$$

(9)

Let $Z={\left[{\overrightarrow{\alpha }}_1,{\overrightarrow{\alpha }}_2\cdots {\overrightarrow{\alpha }}_m\right]}^T$ and $B={\left[{\overrightarrow{\beta }}_1,{\overrightarrow{\beta }}_2\cdots {\overrightarrow{\beta }}_m\right]}^T$. The orthogonalization from Z to B is Z = V × B. According to Equation (9), V is defined as follows:

$$\mathrm{V}=\left[\begin{array}{cccc}1& & & \\ \frac{〈{\overrightarrow{\alpha }}_2\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}& 1& & \\ ⋮& & \ddots & \\ \frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_1〉}{〈{\overrightarrow{\beta }}_1\cdot {\overrightarrow{\beta }}_1〉}& \cdots & \frac{〈{\overrightarrow{\alpha }}_m\cdot {\overrightarrow{\beta }}_{m-1}〉}{〈{\overrightarrow{\beta }}_{m-1}\cdot {\overrightarrow{\beta }}_{m-1}〉}& 1\end{array}\right]$$

(10)

Equation (3) can be transformed into:

$$VB\times A=\Delta S.$$

(11)

Equation (11) is the overdetermined equation after orthogonalization, and a stable solution can be obtained by solving Equation (11) using the least squares method. In addition, the Gram–Schmidt orthogonalized method can be simplified, and the covariance matrix of the coefficient matrix can be used for linear transformations to prevent the orthogonalization process and reduce the cost of computing resources [41].

Swantner et al. [42] generalized the application of Zernike polynomials with respect to annular sector apertures and verified the inference procedure by using the Gram–Schmidt orthogonalized method. Upton et al. [43] optimized the Gram–Schmidt orthogonalized method and extended the method to apertures of an arbitrary shape. Taking the hexagonal aperture as an example, the feasibility of the improved Gram–Schmidt orthogonalized method was verified. A comparison of the fitting results of Zernike polynomials in terms 7–9 is shown in Figure 6. U denotes the fitting result of the circular Zernike, and V denotes the fitting result of the orthonormal hexagonal Zernike.

The ability of Zernike polynomials to characterize information on spherical cap is limited by the aperture angle. Zheng et al. [44] derived Zernike-like functions for all types of spherical caps and used the Gram–Schmidt orthogonalized method for calculations. The spherical cap’s homeomorphism maps to the unit circle, as shown in Figure 7.

The original intention of using the Gram–Schmidt orthogonalized method is to avoid the ill-conditioned problem of the least squares method. Although the Gram–Schmidt orthogonalized method solves ill-conditioned problems effectively, it introduces new unstable factors. When the least squares method is ill-conditioned, the orthogonal basis function system constructed by using the Gram–Schmidt orthogonalized method is correlated, which also affects the stability of the solution [45]. In addition, when the orthogonal basis functions appear to be correlated, the solution obtained before is still stable without interferences.

3.3. Householder Transformation

The Householder transformation method can solve the overdetermined equation directly without constructing the normal equation in order to avoid ill-conditioned problems in the process of solving equations. The coefficient matrix of the overdetermined equation is orthogonalized and triangulated by Householder transformations [46]. Based on the rules of Householder transformations, Z can be decomposed by the following:

$$Z=Q\times \left[\begin{array}{l}R\\ O\end{array}\right],$$

(12)

where Q is an m × m orthogonal matrix, R is an n × n upper triangular matrix, and O is a (m − n) × n null matrix. The least squares solution of Equation (3) is expressed as follows.

$${‖ZA-\Delta S‖}_2^2={‖Q\left[\begin{array}{l}R\\ O\end{array}\right]A-\Delta S‖}_2^2={‖Q^{T}Q\left[\begin{array}{l}R\\ O\end{array}\right]A-Q^T\Delta S‖}_2^2={‖\left[\begin{array}{l}R\\ O\end{array}\right]A-Q^T\Delta S‖}_2^2$$

(13)

Let $Q^T\Delta S=\left[\begin{array}{l}b\\ c\end{array}\right]$. Equation (13) can then be expressed as follows:

$${‖ZA-\Delta S‖}_2^2={‖\left[\begin{array}{l}R\\ O\end{array}\right]A-\left[\begin{array}{l}b\\ c\end{array}\right]‖}_2^2={‖RA-b‖}_2^2+{‖c‖}_2^2$$

(14)

When ${‖RA-b‖}_2^2=0$, the minimum value is obtained for Equation (14). According to this condition, Zernike coefficients are solved by the following:

$$A=R^{-1}b$$

(15)

Wang et al. [47] proposed a new method for solving active optic correction forces and simulated a 400 mm test mirror in which the Householder transformation was used to solve Zernike coefficients. The comparative results are shown in Figure 8, and the results revealed that the corrective method is effective. Cao et al. [48] described the human eye’s wavefront aberration by using Zernike polynomials and used the least squares method and Householder transformation to calculate the eye aberrations of different pupils. The calculation results showed that the accuracy of the proposed algorithm is comparable to that of the least squares method. Xiao et al. [49] established an integrated simulation method and chose the Householder transformation as a Zernike CSA. The curve of optical distortion with temperature was obtained by analyzing the objective lens.

3.4. Singular Value Decomposition (SVD)

Although the Householder transform can prevent the ill-conditioned problems of overdetermined equations, the computational complexity increases exponentially with the increase in coefficient matrix size. The SVD can enhance the efficiency and accuracy of the solution process while avoiding ill-conditioned problems [50]. Based on the rules of SVD, Z can be decomposed by:

$$Z=U{\left[\begin{array}{c}\sum \\ 0\end{array}\right]}_{m\times n}{V}^T,$$

(16)

where U is an m × m orthogonal matrix, the column vectors of which are the feature vectors of AA^T; V^T is an n × n orthogonal matrix, the column vectors of which are the feature vectors of A^TA; Σ is a diagonal matrix composed of singular values of A.

The least squares solution of Equation (3) is expressed as follows:

$${‖ZA-\Delta S‖}_2^2={‖U\left[\begin{array}{c}\Sigma \\ 0\end{array}\right]V^{T}A-\Delta S‖}_2^2={‖\left[\begin{array}{c}\Sigma \\ 0\end{array}\right]V^{T}A-U^T\Delta S‖}_2^2$$

(17)

The first n column vectors of U are used to form a new matrix, U_{m × n}, and U is defined as follows:

$$U=\left[U_{m\times n},U_{m\times (m-n)}\right]$$

(18)

Equation (17) can then be expressed as follows:

$${‖ZA-\Delta S‖}_2^2={‖\sum V^{T}A-U_{m\times n}^T\Delta S-U_{m\times (m-n)}^T\Delta S‖}_2^2$$

(19)

According to the matrix’s decomposition, the right side of Equation (19) can be transformed into:

$${‖\sum V^{T}A-U_{m\times n}^T\Delta S‖}_2^2+{‖U_{m\times (m-n)}^T\Delta S‖}_2^2\ge {‖U_{m\times (m-n)}^T\Delta S‖}_2^2.$$

(20)

When ${‖\sum V^{T}A-U_{m\times n}^T\Delta S‖}_2^2$, the minimum value is obtained for Equation (19). According to this condition, Zernike coefficients are solved by:

$$A={(\sum V^T)}^{-1}{U}_{m\times n}^T\Delta S=V{\sum }^{-1}{U}_{m\times n}^T\Delta S.$$

(21)

When the overdetermined equation is ill conditioned, the SVD obtains smaller singular values and dramatically amplifies the variances in the solution, which makes the calculation result unreliable [51]. The truncated singular value decomposition (TSVD) can cut off the smaller singular values, which can easily cause instability, in order to obtain a reliable solution [52]. The selection of truncation parameters is the main factor affecting the optimal solution. There are two methods for selecting the truncation parameters: generalized cross validation [53] and L curve [54]. The reconstruction of three-dimensional flame temperature fields in a furnace is an overdetermined problem. Huang et al. [55] used TSVD and the L curve to reconstruct temperature fields, and the reconstructed temperature distribution is shown in Figure 9. Xie et al. [56] compared Tikhonov regularizations and TSVD in the same problem. TSVD exhibited no higher accuracies due to the truncation parameter selection problem.

Currently, Zernike CSAs are broadly used. However, traditional algorithms have their own defects [57]. For reflective optical systems, they generally exhibit large aperture characteristics. In order to ensure the density of finite element nodes, the total amount of node deformation data significantly increased, and the computational complexity of traditional algorithms in the fitting process significantly increased. Although the optical surface size of the transmission optical system is generally relatively small, the number of optical surfaces is significantly more than that of the reflection optical system, which is also limited by the solution’s efficiency. At the same time, the fitting accuracy is related to the number of terms used by Zernike polynomials. It is generally believed that, if more terms are used, the fitting accuracy will be higher. However, when the number of terms exceeds 36, the traditional algorithm will have serious abnormal problems, and the selection of truncated parameters will be more challenging, reducing the accuracy of solving Zernike coefficients. Therefore, enhancing the fitting accuracy and efficiency of traditional Zernike CSAs is critically important.

4. The Potential for ANN Use in IOA

According to Equation (2), when the polar coordinates of the finite element nodes are inserted, fitting sag data essentially solve the overdetermined equation sets defined by Zernike polynomials [58]. With the development of artificial intelligence algorithms, many artificial neural networks (ANNs) have been proposed and used to solve overdetermined equation sets, such as the backpropagation neural network (BPNN) [59], convolutional neural network(CNN) [60], and recursive neural network (RNN) [61]. The training objective of artificial neural networks involves finding a mapping relationship between the input vector and output vector [62]. Therefore, an artificial neural network can be regarded as a function approximation tool. This is similar to the solution of an overdetermined equation. According to the structural similarity between ANNs and overdetermined equations [63], many fields began applying ANNs to solve their own overdetermination problem. Zhang et al. [64] presented a neural network model based on ordinary differential equations (ODE) to compute generalized and restricted singular value decompositions; the neural network’s structure for tracking the CS component is shown in Figure 10. The stability of the neural network was proven by numerical experiments, and it was applied to solve the overdetermined equation; the calculation error of the random matrix ranges from 10⁻¹⁵ to 10⁻¹⁷. The results show that the neural network can effectively solve the ill-conditioned problem of the overdetermined equation.

A Zhang neural network (ZNN) is a special RNN used to solve the time-varying overdetermined system of linear equations. Zhang et al. [65] designed two new error functions for ZNN and verified the effectiveness using computer simulations. Zhang et al. [66] proposed a varying-parameter convergent-differential neural network (VP-CDNN); the neural topology of the VP-CDNN model is shown in Figure 11 [67]. The algorithm can obtain a least squares solution with a super-exponential rate. VP-CDNN has more speed, accuracy, and robustness than RNN and zeroing neural networks. In addition, the VP-CDNN is also used to solve time-varying convex quadratic-programming problems [68].

Based on these studies, using neural networks to solve the overdetermined equation is feasible. Distinct from the above cases, due to the fact that Zernike polynomials comprise one-dimensional data, it is extremely appropriate to use a one-dimensional convolutional neural network (1D-CNN) to solve Zernike coefficients. Figure 12 is a model of a one-dimensional convolutional neural network [69]. Both one-dimensional convolutional neural networks and two-dimensional convolutional neural networks have strong local sensing abilities and parameters. The only difference is that the convolutional kernel of 1D-CNN is one-dimensional, which can extract features from one-dimensional sequences; thus, it can only perform vector operations without complex matrix operations, which significantly reduces memory consumption [70]. According to this characteristic, using 1D-CNN to solve Zernike coefficients has a significantly higher computational speed than the traditional algorithm, particularly in the case of a large number of finite element nodes to be calculated.

Training a 1D-CNN model equivalent to the overdetermined equation is the first step of solving Zernike coefficients. A 1D-CNN model usually includes an input layer, a 1D convolutional layer, pooling layer, a fully connected layer, and an output layer [71]. The Zernike polynomials of each node are taken as the input of the neural network, and the sag change is taken as the output. The number of neurons in the input layer is set according to the number of Zernike items. A convolutional layer is used to extract the features of Zernike polynomial data. The pooling layer is used to reduce the data dimension of a fully connected layer input, and it saves more detail when reducing the amount of computation. A full connection is used to learn the relationship between high-level features obtained from the pooling layer and the output. In the 1D-CNN model, the activation functions used to connect each layer primarily include sigmoid and ReLU. The function’s image is shown in Figure 13 [72]. The sigmoid function can introduce nonlinearity to 1D-CNN and enhance the learning ability of the model. The ReLU function can alleviate the problem of vanishing gradients in the neural network. The dropout technique was introduced to prevent overfitting [73,74]. By relying on the powerful prediction performance of one-dimensional convolution [75] and using the data set obtained by finite element analysis to train 1D-CNN, a model equivalent to the overdetermined equation can be obtained.

The second step was to adjust the structure of 1D-CNN according to the training results. This is the most critical step, as it affects the accuracy and training efficiency of the model and further affects the final accuracy of the solution. In addition, another method for optimizing the neural network model is to change the gradient descent algorithm. Currently, common gradient descent algorithms include adaptive moment estimation (Adam), the random gradient descent algorithm (SGD), and the adaptive learning rate optimization algorithm (Adagrad). Finally, after the optimized 1D-CNN model is obtained, the identity matrix is imported into the model for the prediction of results. As a result of training, the 1D-CNN model is equivalent to the overdetermined equation. The predicted results are shown in Equation (22). Since the training error of 1D-CNN is close to 0, the predicted result can be regarded as the Zernike coefficient.

$$I\times A=A$$

(22)

In order to verify the reliability of the above method, we conducted a test experiment. The data set for training 1D-CNN was obtained by applying finite element analysis on the primary mirror of the front reflection system of the space camera. It is a common condition that gravity direction is parallel to optical axis when the space camera is working. Therefore, the gravity response of the primary mirror when gravity direction was parallel to the optical axis was analyzed. The radius of the primary mirror was 645mm. The rigid body displacement of the primary mirror was separately input into the optical analysis in IOA; thus, the optical surface deformation that removes the rigid body displacement was fitted. Zernike polynomials select the 34-term polynomials that remove the tip/tilt/piston terms [76]. Based on the 34-term Zernike polynomial, the input of 1D-CNN is 34 terms, and the output is a node sag change. The finite element node data were divided into a training set, test set, and verification set according to the proportions of 70%:15%:15%, respectively. The results of the training are shown in Figure 14. Figure 14a shows the error histogram, and Figure 14b shows the regression distribution.

According to the error and regression distribution of the model, the model can be regarded as equivalent to the overdetermined equation constructed by Zernike polynomials. Taking the identity matrix as the input, each Zernike coefficient was obtained. The optical surface deformation was obtained by substituting Zernike coefficients into the overdetermined equation. The fitting result of the optical surface deformation is shown in Figure 15.

To verify the stability of the neural network method, a small random disturbance was applied to the data set, and the model was trained again. The verification result is shown in Figure 16. The results show that the error of the training results with a small random disturbance exhibited the same order of magnitude, and the error distribution was consistent. The accuracy of the neural network model was not affected by random errors, and the neural network method exhibited stability. The neural network method avoided the ill-conditioned problem.

The fitting residual was obtained by subtracting the fitted sag change from the actual sag change. The fitting residual of the neural network method is shown in Figure 17. The root mean square of the residual was 17.02 nm. The average value of the ratio between the absolute value of fitting residual and the actual sag change is used as the criterion for fitting accuracies. The smaller the average ratio, the higher the fitting accuracy. The error of the neural network method was 0.864%. As a comparison, the least squares method was used to solve and fit the sag change. The error of the least squares method was 1.083%. After adding small random disturbance, the error of the neural network method was 0.872%, and the error of the least square method increased to 2.933%. Due to the robustness of neural networks, neural network methods have a strong resistance to small perturbations [77,78,79], which shows that the neural network method can effectively solve the ill-conditioned problem of the overdetermined equation. The accuracy and reliability of the proposed neural network method are verified by the above testing experiments.

5. Discussion

Various past studies have shown that IOA plays an obvious role in the design of optical instruments [80,81,82], and Zernike CSAs have an important impact on the accuracy of IOA. Current Zernike CSAs have different application scenarios according to their characteristics. The least squares method is the simplest Zernike CSA. The method exhibits fast calculation speeds and can maintain calculation accuracies when the number of Zernike polynomial terms is small. However, when the number of terms is large, the solution is unstable due to ill-conditioned problems. Therefore, the least squares method is mainly applicable to the IOA of small optical instruments and optical devices that do not need to calculate more Zernike terms. The Gram–Schmidt orthogonalized method improves the Zernike polynomial and solves the ill-conditioned problem of the least squares method, but it introduces new unstable factors. In non-circular aperture optical instruments, Zernike polynomials completely lose orthogonality. For the non-circular optical instrument, using the Gram–Schmidt orthogonalized method is an effective method for orthogonalizing and then solving equations. The Householder transformation method avoids the construction of normal equations and directly solves them. It is not constrained by the number of Zernike polynomial terms and exhibits good stability and accuracy. Its computational complexity increases exponentially with the increase in data size, and it is suitable for IOAs with lesser finite element mesh. The SVD method has high computational accuracy in theory, but its application is limited by the selection of truncation parameters; thus, it is more suitable for projects with high precision requirements.

For some large-scale optical instruments, increased fine finite element meshing is usually required [76]. Banyal et al. proved that a higher order Zernike polynomial also has a certain influence on the fitting of optical surfaces [37]. Compared with the traditional algorithm, the neural network algorithm exhibits good accuracy and stability; in this case, the number of Zernike polynomial terms and grid division is no longer restricted. In addition, the working environment of some optical instruments is complex and requires multiple IOAs according to different working conditions [83,84,85]. A neural network algorithm has higher computing speeds and is more suitable for the IOAs of such optical instruments. Since the Zernike coefficient is obtained by the neural network algorithm via the direct fitting of the node data of the finite element network, the orthogonality of the Zernike term is not changed. In contrast, for the IOAs of some small optical instruments, the traditional method considers both solving efficiency and solving accuracy; thus, the traditional algorithm is the better choice.

In the test experiment, the neural network algorithm was used to achieve the surface fitting operations of the primary mirror in the front reflection system of the space camera. The RMS of the fitting residual reached 17.02 nm, and the regression evaluation index R of the neural network reached 0.9998. Compared to the least squares method, the neural network method improves the solution’s accuracy and stability. However, the neural network algorithm also has some shortcomings. Overfitting is one of the main factors affecting the fitting accuracy [86,87]. In addition, the local optimization caused by the gradient disappearance of the neural network also affects the final result [88]. Methods for avoiding these problems and for further improving fitting accuracies are the focus of future research studies.

6. Conclusions

Zernike CSAs, as an important link in the IOA that influences analysis results, are indispensable in the IOA of various optical systems in space and on Earth. The calculation accuracy of Zernike CSAs is an important factor that affects the results of IOA, and the optimization of Zernike CSAs can obtain more accurate analysis results. Currently, Zernike CSAs all have certain limitations that impede both efficiency and precision for projects requiring increasingly accurate simulations. The errors of Zernike CSAs lead to failures in obtaining accurate MTF and aberration analysis results. In this review, we summarized and analyzed traditional Zernike CSAs and combined current neural network applications when solving overdetermined equations in other fields. An innovative Zernike coefficient-solving algorithm based on a one-dimensional convolutional neural network was proposed. The neural network algorithm considered both accuracy and efficiency and increased the number of Zernike polynomial terms that can be calculated. More Zernike items can be used to analyze the higher-order aberration of optical devices. In the future, large-scale optical devices need to analyze the impact of high-order aberrations, and space projects also face more complex external environments. As a method to improve the accuracy of IOA, the neural network method is an effective solution. As hardware technology evolves and computing power advances, we expect that neural network algorithms will play an important role in IOAs.