An inline surface measurement method for membrane mirror fabrication using two-stage trained Zernike polynomials and elitist teaching–learning-based optimization

The accuracy of surface measurement determines the manufacturing quality of membrane mirrors. Thus, an efficient and accurate measuring method is critical in membrane mirror fabrication. This paper formulates this measurement issue as a surface reconstruction problem and employs two-stage trained Zernike polynomials as an inline measuring tool to solve the optical surface measurement problem in the membrane mirror manufacturing process. First, all terms of the Zernike polynomial are generated and projected to a non-circular region as the candidate model pool. The training data are calculated according to the measured values of distance sensors and the geometrical relationship between the ideal surface and the installed sensors. Then the terms are selected by minimizing the cost function each time successively. To avoid the problem of ill-conditioned matrix inversion by the least squares method, the coefficient of each model term is achieved by modified elitist teaching–learning-based optimization. Subsequently, the measurement precision is further improved by a second stage of model refinement. Finally, every point on the membrane surface can be measured according to this model, providing more the subtle feedback information needed for the precise control of membrane mirror fabrication. Experimental results confirm that the proposed method is effective in a membrane mirror manufacturing system driven by negative pressure, and the measurement accuracy can achieve 15 µm.


Introduction
In recent years, the space reflector has been used in various fields such as remote sensing, solar energy concentrators and astronomical applications [1][2][3]. The requirements for space optics are continuously growing to satisfy increasing demand in different applications. For purposes of obtaining information about Earth such as geophysical parameters, meteorological data and reconnaissance, all the currently used remote sensing optical instruments are operated in low Earth orbits (LEOs) to reduce the negative imaging effect caused by cloud. There are microwave sensors with the ability to penetrate through clouds, providing detailed ground information. However, these sensors need to work in high Earth orbits. With the increasing observing distance, the quality of the measurement information would deteriorate if the same optical structure was employed. To acquire the same data quality for the sensors in LEOs, a significantly large antenna with an accurate surface is desirable due to the fact that sensitivity and spatial resolution are determined by the optical surface precision, and the signal-to-noise ratio and signal resolution are positive related to the aperture of the mirror [4]. Similarly, in astronomical fields, a space-borne mirror with a large aperture must get rid of atmospheric turbulence to improve the space telescope's performance. Referring to the size of the James Webb Space Telescope, the future requirement for the diameter of primary mirrors is no less than 10 m [5]. The traditional optical reflector cannot meet this new requirement due to the fact that the solid monolithic lens is too large and heavy for the launch mass and storage size of current launch vehicles. Moreover, it is impossible to fabricate the desired mirror by traditional optical manufacturing methods. To meet this challenge, the membrane optic is rapidly developing. As a new optical element, the membrane mirror has the merits of low aerial density, easy deployability and low cost. These characteristics enable the membrane mirror to break through the constraints of traditional optical manufacturing, providing a suitable and alternative choice for the large-aperture and ultralight mirrors required in space telescopes, and other spacebased optical applications.
At present, electrostatic stretch and pneumatic pressure (inflated or vacuum suction) shaping are the two major approaches for fabricating membrane mirrors, both of which require precise measurements of the membrane surface [6][7][8][9]. In addition, the optical surface measurement accuracy determines the final quality of the membrane mirror, which is the critical performance of the whole optical system. Thus, the surface measurement method plays the most important role in the membrane mirror fabrication process. Different from the traditional mirror, the membrane mirror needs to maintain sufficiently tight surface accuracy during the applied process [10]. Considering that the optical surface varies due to environmental factors, such as temperature, pressure and other disturbances, active adjustment is always required to compensate for surface errors. Due to the vulnerability of membrane materials, the contacting measuring method could not be applied to test the surface of the membrane mirror. In ground-based tests, some metrologies have already been developed. One way is to take a photograph of the concerned mirror surface, and then calculate the related coefficients using special software. NASA and SRS utilize this method to evaluate surface accuracy [11,12]. Moreover, researchers from the University of Arizona and the University of New Mexico adopted interferometer and moiré fringes to measure the mirror surface [13,14]. By setting up an imaging system, the image quality assessment of standard pictures is used to describe the surface error by the Air Force Research Laboratory [15]. Similarly, photogrammetry is adopted to measure the surface of inflatable membrane structures [16]. However, the aforementioned methods all require additional expensive equipment, and are not feasible for surfaces with ultra-large apertures. Comparatively speaking, measuring the mirror surface by setting appropriate sensors in the shaping frame is economical and easily implemented. Similarly, the quantity of sensors in the modeling frame is finite due to budget limitations. Considering that the mirror surface is always changing in the forming process, it is impossible to express the whole surface by only using limited information from a few fixed sensors. In addition, sometimes the expected mirror surface is so complicated that it is hard to represent the surface only using partial information.
To overcome this difficulty, approximation methods are adopted in surface measurement systems. Haber et al used a subspace identification technique to obtain a dynamic description of a thermally actuated deformable mirror [17]. Song et al employed a neural network to solve the aberration correction problem in an optic system [18]. However, the adopted methods are all off-line and not suitable for real-time membrane mirror fabrication. Among various surface fitting algorithms, the polynomial approximation is the most popular one. The desirable properties of Zernike polynomials, such as orthogonality, rotational symmetry, relation to classical Seidel aberrations, and simple representation, have made it the most popular basis function for analyzing optical surfaces [19][20][21]. Ares and Royo studied the fitting performance of Zernike polynomials, and found that low-degree Zernike polynomials are suitable for fitting simple wavefronts [22]. MacMartin et al analyzed the structural interaction of the segmented mirror of a telescope using the Zernike basis, offering guidance for structural optimization of mirrors [23]. It is also effective for optical surface measurement. Liu et al employed Zernike poly nomials to fit the deformed surface of a telescope mirror, as a reference for mirror configuration design [24]. In addition, it is worth noting that since the traditional Zernike polynomials are defined within the unit circle, it would lose the original favorable properties in a non-circular region. Fortunately, He et al proposed an effective method to construct novel Zernike polynomials in non-circular regions, which extends their application [25]. However, none of the aforementioned studies have considered the performance of high-order Zernikes. Alkhaldi et al found that interpolation with higher-degree Zernike polynomials brings better performance [26]. However, Runge's phenomenon will appear in polynomial interpolation while the Zernike order is increasing, resulting in a poor approximating performance. The existing studies show that Zernike polynomials are an effective tool to measure the surface of optical elements, but an appropriate strategy is still called for to select the Zernike terms.
Theoretically, surface measurement or reconstruction by Zernike polynomials is a specific expression of the linearin-the-parameters model, which is widely used in nonlinear system identification. However, the most popular orthogonal least squares (OLS) method would cause significantly large computational complexity in dealing with big data [27]. Further, the model performance would drastically deteriorate if OLS was used to solve an ill-conditioned problem, which is very common in mirror surface reconstruction. As a result, a fast model constructing strategy that is non-sensitive to the ill-conditioned matrix is desired. In this paper, a two-stage subset selection scheme, avoiding matrix inverse operation, is employed to realize the surface reconstruction. To maintain the desired properties of the Zernike terms, Gram-Schmidt orthogonalization is employed to construct novel Zernike terms in non-circular regions, which are used as candidate bases. In order to further improve the model accuracy, coefficients are further optimized by a modified elitist teachinglearning-based optimization (ETLBO) algorithm coming after the model structure is determined. Satisfying surface measurement accuracy and speed would be achieved by this method, providing a feasible means to maximize the performance of Zernike polynomials.

Two-stage model selection scheme
The linear-in-the-parameters model has a proven ability to approximate arbitrary nonlinear functions with arbitrary precision, and the general form of this model is where t N 1, 2, ..., = ; N is the number of training data; y t ( ) denotes the model output and t x( ) denotes the model input vector at time constant t; p i M , 1, 2..., i = represents all the candidate nonlinear bases with certain number parameters v i ; i θ denotes the linear coefficients of different nonlinear bases; and t ( ) ξ is the model residual with zero mean. To facilitate computation, the matrix form of (1) could be written as is an M-dimension column vector; and y y N y 1 , ...
are all N-dimension column vectors. If the structure of the model is fixed, the corresponding linear coefficients could be obtained by minimizing the following cost function where the least squares method is adopted, and the optimal coefficients could be given as P P P y.
The approximating capability of this model is rooted in the characteristics of basis functions. As long as the function used is a complete basis, the expected identification performance could be achieved by a finite linear combination of basis functions. The only problem is how to find a compact model with the desired index in a certain period of time. However, in many applications, the training data set for identification is very large, generating a large candidate basis pool. Therefore, if all candidate bases are used like extreme learning machines, the computational complexity of (4) may become extremely high and it may become impossible to solve due to the ill-conditioned matrix. To deal with this problem, a forward recursive algorithm (FRA) is used to generate a parsimonious model, and the coefficients are obtained by heuristic optimizing algorithms avoiding the operation of matrix inversion.

Forward model construction method
The FRA method is a forward construction for the model which adds candidate basis functions one by one. In each forward step, the chosen basis is the one contributing the most to the cost function among all the candidate bases. Suppose K basis functions were selected in the kth step, the corresponding regression matrix is expressed as Then according to (3) and (4), the cost function in the kth step could be calculated by J P y y y P P P P y.
If the cost function does not reach the expected value by the combination of existing chosen bases, a new basis function should be added in the next step. Suppose the new selected basis is p k 1 + , the new regression matrix becomes p P P, This new selected basis should make the cost function reach the minimum value among all the remaining candidates. In other words, the reduction of the cost function formulated as (7) should be maximized by p k 1 + : The new basis p k 1 + satisfies It is easy to find that if the basis selection principle is realized by (7) and (8), a number of matrix inversions would be involved in the model constructing process. Thus the computation complexity is very high, and it also leads to numerical stability problems. To overcome the above deficiencies, a matrix series is defined as The computation of the net reduction of the cost function is substantially simplified using the following properties: Substituting (9)-(6), the cost function becomes J P y R y.
By applying the properties shown in (10)-(13), the net reduction of cost function in the kth step could be calculated by To be specific, the forward model construction procedure is described as follows. Say at step k, a new base from the candidate pool is checked. The matrix R k 1 + should be calculated by (10), and then the net contribution to the model performance by this new base is given by (15). The base that contributes most to the reduction of the cost function in the candidate pool will be added to the model. This process will not be terminated until the reduction by the best new base is insignificant or the desired performance is reached.

Backward basis reselection
Although the forward model construction method provides an efficient way to generate a compact model, the optimality of the obtained model cannot be guaranteed. For the non-independency of candidate bases, the optimal combination of candidates is hard to acquire by the step-wise method. However, the performance of the model could be improved by a refinement operation. The refinement is to re-evaluate the significance (contribution to the cost function) of all selected bases and the remaining bases in the candidate pool one by one. If an unselected base contributes more to the cost function than a previously selected one in the forward stage, replacement will take place. The terminated condition of the refining procedure is that there is no further reduction of the cost function.
It is supposed that p p , , n To review a previously selected base in the generated regressor, say p q in P n , it is moved to the last position of P n at first, as if it were the last selected base. This process could be implemented by interchanging two adjacent terms p k and p k 1 + until the concerned basis is moved to the nth position. To facilitate further computation, another important property of R k should be noted, which is According to (16), any change of base position does not influence the residual matrix R k . So after a series of interchange operations, the only changed residual matrix is R q , which could be recalculated by Subsequently, the reviewed base is moved to the last position of the selected basis vector, and is denoted p n . And then its contribution to the cost function is compared to the remaining bases in the candidate pool. The contribution of p n can be calculated by (15), and the corresponding change in residual matrix should be noted. The contribution of the unselected base in the forward stage, say i φ , can be recalculated by where If there is a remaining base contribution that is more significant than that of the reviewed base, the reviewed base is replaced by the most significant one left in the candidate pool. As a result, the performance of generated model can be further improved.

Model coefficient estimation method
Using the two-stage construction scheme, the structure of the model can be determined. The estimation of model coefficients is then achieved by a modified ETLBO in this paper. Among various meta-heuristic based methods, teaching-learningbased optimization (TLBO) is one of the more popular new tools, and is proposed by Rao in [28]. It mimics a class of teaching process where the teacher and students share ideas to gain group knowledge. The algorithm turns out to be a powerful tool for solving a number of constrained/unconstrained engineering optimization problems [29][30][31]. The original TLBO has potential to be enhanced by variants modifications for specific applications.
Some studies [32,33] adopted an elitist strategy to increase the convergence speed of the TLBO. However, the number of elitists in these methods is fixed all along the iterative optimization of the problem, due to which an improper selection of the number may lead to premature or slow convergence. To relieve this, a novel modified elitist TLBO is proposed in this paper. The number of elitists inertially changes with the iteration process aiming to intelligently keep the elitists without largely gaining possibilities of being trapped in the premature. Two evolution phases, namely the teaching phase and the learning phase, are employed in the algorithm process. It is assumed that a population of particles is a class of students. The elitist strategy is embedded within the teaching phase.
The students then gain knowledge from the difference DMean i as shown below, where St ij new and St ij old are the jth new and old students of the ith iteration. The knowledge of the old and new students is evaluated and the better ones will be retained in the student population for the next phase.

Learning phase
The learning phase is the second step of TLBO and mimics the class learning of the student by personal interaction. In this section, every student will be given a chance to randomly find a classmate and gain knowledge from this classmate. The detailed implementation of the step is shown below, where St ij and St ik are the jth and kth students selected from the population in the ith iteration. The student St ij updates his  knowledge by learning the deviation between him/herself and another randomly selected kth student St ik . The student with better knowledge performance will be the dominant learning direction and the learning student will update their knowledge accordingly.

Elitist strategy
The original TLBO performs well on the majority of benchmark tests in terms of exploration and exploitation ability [34]. However, it has slowly met some problems due to the average focus on the population and the missing of some important solutions. The elitist strategy aims to accelerate the convergence by maintaining the best performing solutions in each iteration rather than updating all the candidates [35]. These best performing particles act as elitists and are assumed to have a higher possibility of achieving the global optimum. In this paper, N e elitists are reserved for the next generation. The number N e is defined as inertial decreasing from N eMax to zero, shown as follows, The decreasing number of elitists would significantly speed up the convergence speed in the first stage of the iteration process with many elitists, and keep the exploitation ability in the later stage by removing the elitists for more potential trials. The elitists are selected by ascending order of fitness function evaluations and used in calculating the Mean i of TLBO in (21). This ETLBO makes it easy to achieve the desired searching results with low computational complexity. It is suitable for estimating the coefficients of the mirror surface model.
In order to illustrate the performance of the new ETLBO algorithm, it was compared with some popular meta-heuristicbased methods, including weighted PSO [36], PSO-CF [37] and classical DE/rand/1/bin [38] as well as the original TLBO method. The PSO method mimics bird swarms adjusting positions using social and cognitive learning methods, associated with two learning parameters c 1 , c 2 , and weighting factor w. The DE method demonstrates individual interactions by crossover and mutation sections, with crossover rate CR and mutation rate F as the featured parameters. The simple variants used here are the most popular ones in their family. All the tests are implemented in 10 well-known functions as defined in [39] and shown below:

Experiment and results
This paper investigates the surface measurement problem in membrane mirror manufacture. The proposed modeling method is used for mirror surface measurement during the fabrication process. In other words, the mirror surface is reconstructed by identification using information measured by several distance sensors.

Experimental setup
In order to shape the membrane mirror, a negative pressure approach is adopted. The negative pneumatic forming device is a closed loop system, including controller, sensor and actuator, shown in figure 2. A mechanical mold with a precise parabolic surface is utilized as the frame of the desired membrane mirror. The membrane material is clamped on the frame edge and its shape is changed by the negative pressure. To generate vacuumed force, the mechanical frame is designed to have a hollow structure, with the air pipeline assembled on its sub-surface. As the power source, a draught fan generates suction at different amplitudes according to the control signal from the controller. Meanwhile, an air-valve array is also designed to provide an additional control channel for adjusting the power of the negative pressure. An industrial computer with a VME (VersaModule Eurocard) bus is selected as the control unit for its excellent real-time performance and high fault tolerance capability.
Different from traditional surface measuring methods, a sensor array is designed inside the top surface of the mold frame, collecting shape information on the membrane mirror and providing training data to the surface reconstructing algorithm. In this paper, an ultrasonic range finder mic + 25/D/ TC is selected to constitute this sensor array. It could achieve 0.015 mm resolution in the range of 30-350 mm. Further, this sensor is non-sensitive to environmental change since it selfcompensates for temperature and pressure. Due to space and cost limitations, it is not possible to make the sensor array too large. On the other hand, more training data for surface identification could be obtained by increasing the quantity of sensors. To reconcile this contradiction, an electromechanical unit is built to enable the sensor to swing in a small range 2.5 , 2.5 [ ] θ ∈ −°°, with 0.5° angle resolution, as illustrated in figure 3. In the experimental system, the sensor array is composed of 36 range finder units distributed on the top surface of the frame, and the desired surface would not be affected by the swing mechanism of the sensor array.   Both the surface reconstruction and control strategy are implemented in the motion control card, using the digital signal processor (DSP) TMS6414t as the computational core. This DSP could operate at a main frequency of 1 GHz, providing powerful computational ability. To transfer data among the closed loop units, such as the data acquisition card, motion control card and power amplifier, a VME bus with 200 Mbit s −1 transfer speed is used, supplying an extra real-time guarantee. An information flow chart of the manufactured system is shown in figure 4

Simulation discussion and data analysis
A two-stage model selection algorithm is used to measure the surface of the membrane mirror. To evaluate the measuring performance, the negative pressure is kept constant to maintain the expected mirror surface (already formed). The distance information collected by the sensor array is used as training data for the identification algorithm.
As mentioned above, there are 36 sensors in the array, and each sensor can generate 11 training points by swinging. Therefore, there are in total 396 training points for surface reconstruction. Considering the layout of the sensor array, a coordinate transformation should be conducted on the distance data, which is defined as  As the most popular optical analysis tool, Zernike polynomials with degree numbers from 0 to 20 are generated and include 231 terms. The membrane mirror in this paper lies in a rectangular region rather than a unit circular region, which means that the properties of the original Zernike polynomials cannot be maintained. It is necessary to construct a novel set of Zernike terms. According to reference [25], the original Zernike terms could project into a new defined region by Gram-Schmidt orthogonalization. The procedure could be expressed as where Z i denotes the standard Zernike term; P i is the novel constructed Zernike term; a i and b i are coefficients to be determined; and m and n are set to be 40 and 50 respectively, based on the dimensions of the rectangular region. As a result, the obtained novel Zernike terms are still orthogonal in the noncircular region and possess clear physical meanings, which could be used as a candidate basis pool for reconstructing the mirror surface. Next, 196 sampling data are employed for identification, and the remaining 196 points are used to validate the model performance. The x and y coordinates are regarded as model input, while the value in the z direction is taken as the output. Applying a two-stage strategy with the root-mean-square (RMS) as the cost function, 90 terms are selected from 231 candidates, with the stopping criterion set at 0.02 mm. The index number of model terms selected by FRA and the two-stage method are listed in table 2. As seen  in table 2, the model performance is further improved by the refinement operation in the second stage, by replacing some bases selected in the first stage. Because the acquired data contains redundant information, the corresponding information matrix is highly ill-conditioned. As a result, although the OLS method makes use of all the candidate bases, the error of the obtained model is still far beyond the tolerance. These results show that the two-stage strategy could make a compact model with desirable performance. Further, it is non-sensitive to ill-conditioned problems.
After the model structure is determined, the coefficient of each model basis can be optimized using the aforementioned ETLBO. The motivation of introducing an elitist strategy is to avoid premature convergence and speed it up. Although a number descending strategy is adopted in the modified ETLBO in this paper, the initial number of elitists or the maximum elite proportion in the population needs to be discussed. The convergence speed increases with the elite number. However, too large a ratio of elitists in the population will cause premature convergence and a fall in the local minimum trap. Too small an elite number slows convergence, weakens the elite strategy, and means optimal performance cannot be obtained, yet an overwhelming penetration of elitists in the population risks causing premature convergence. In the application of ETLBO in the surface reconstruction problem, different values of elitist number N eMax are tested for the linear-in-the-parameters model of the expected mirror surface obtained by the two-stage method. Moreover, three counterpart algorithms, namely weighted PSO [36], TLBO [28] and DE/rand/1/bin [38], compared in the benchmark tests are also employed for comparison here. Moreover, to analyze the sensibility of multiple param eters in each method, five different parameter settings are tested for each algorithm and 30 different runs are employed. The number of the population is set as 30 and the maximum iteration is set as 400; the simulation results are listed in table 3. According to the testing results, the repeatability of each algorithm is validated and reasonable parameter settings of each method are obtained: the weighted PSO uses c 1 = 1, c 2 = 3, w = 0.9; the parameters F = 0.7, CR = 0.9 are set for DE; and the N eMax of ETLBO is set as 8, which could be used in the following surface reconstruction experiment.
In practice, the coefficients of the surface model are searched using these optimization methods with the aforementioned settings. The number of the population is set as 30 and the maximum iteration is set as 400, while the stopping criterion is set as 0.02 mm. The corresponding results are listed in figure 5 and table 4.
As seen from table 4, the modified ETLBO ranks first in both the RMS and the peak error of the reconstructed surface. The expected surface and reconstruction error are shown in figures 6 and 7 respectively. Using the constructed model, the mirror surface measuring accuracy could reach 15 µm.

Conclusion
In this study, the membrane mirror surface measurement problem is solved from the point of view of system identification. The mirror surface is constructed in the form of a linear-in-the-parameters model, using transformed Zernike polynomials which are suitable for the rectangular region as a basis function. The model terms are selected using a two-stage strategy. Then the corresponding coefficients are optimized by a modified ETLBO algorithm.
From analysis of the simulation results and the experiment conducted on a pneumatic membrane mirror fabrication system, the two-stage strategy is shown to provide an effective way to maximize the approximating performance of Zernike polynomials and generate a parsimonious model. Furthermore, this recursive modeling method behaves well with ill-conditioned problems, while the OLS is out of operation. At the same time, the coefficients estimated by the modified ETLBO also help to improve the model performance in comparison with other heuristic methods. According to simulation analysis, the number of elitists in ETLBO should be selected as 25% of the total population in order to reconcile the contradiction between accuracy and convergence speed. Validated by the testing data, the surface measuring accuracy of the obtained model can reach 15 µm. The proposed surface measuring method is also applicable to other manners of membrane mirror manufacture, as long as enough training data can be acquired for the identification method.