The unknown emissivity of materials is a huge obstacle in multi-wavelength pyrometry (MWP). It leads to a set of ill-posed equations that cannot be directly inverted to obtain the true temperature from a set of multi-wavelength measurements. Constraint optimization algorithms such as the gradient projection (GP) and internal penalty function (IPF) algorithms provide solutions without any emissivity model assumptions but require a narrow fixed emissivity range and an appropriate initial emissivity input value, otherwise, accuracy and computational efficiency are greatly affected. Here, we propose a generalized inverse matrix-exterior penalty function (GIM-EPF) algorithm to realize an efficient and accurate inversion without limiting the emissivity range in advance. First, a set of emissivities is obtained by the generalized inverse matrix method; these emissivities are used as initial values in the exterior penalty function iteration algorithm, from which temperature and spectral emissivity are obtained. Simulation results show that the GIM-EPF algorithm provides results superior to IPF, especially in computational efficiency. The proposed GIM-EPF method is 8 times faster than the IPF method with a 0.56% relative error at the 1800 K true temperature. The GIM-EPF method also allows near real-time diagnosis of rocket exhaust flame temperature. Rocket nozzle temperature experiment results show that the temperatures derived by the GIM-EPF algorithm agree well with the theoretical design temperature and the IPF inversion temperature.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Multi-wavelength pyrometry (MWP) is useful for the measurement of high or ultra-high temperature targets whose temperatures cannot be measured by thermocouples [1–6]. It uses the intensity of the emitted thermal radiation in multiple spectral channels to directly retrieve true temperature and associated material spectral emissivity based on Planck's formula [7–13]. Compared with other non-contact temperature measurement methods, such as Holographic Interferometry , CARS , and TDLAS  method, the MWP method has the advantages of simplicity of instrumentation, minimal influence from external parameters, and direct inversion of Planck’s formula to obtain true surface temperature and surface spectral emissivity. It has been widely used in jet and rocket engine testing, plasma diagnostics, and other applications [17,18]. Data processing in MWP requires the solution of n equations based on Planck's law for n spectral channels; however, these n equations contain n + 1 unknowns, namely the n unknown spectral emissivities and the true temperature . Before 2017, all the data processing methods for MWP assumed beforehand a model for the emissivity. The models including emissivity-wavelength models [20–23] and emissivity-temperature models [24–26]. Obviously, if the model is not consistent with reality, the accuracy of the parameters derived from the measurements will be greatly affected.
In 2017, in order to achieve true temperature inversion without any emissivity assumption, we proposed that the MWP data processing problem be converted into a constraint optimization problem. Gradient Projection (GP) and Internal Penalty Function (IPF) algorithms are standard approaches for solving this kind of problem . The results demonstrate that the IPF algorithm is much better than GP in both accuracy and efficiency. IPF needs as input a reasonable emissivity range. However, the efficiency of the IPF approach is poor, and simulation results for 6 kinds of materials show that the average inversion time of IPF takes about 100 seconds, while the GP method requires more than 250 seconds. This makes it difficult to meet the requirement of real-time temperature measurement. In the present paper we develop and test a new Generalized Inverse Matrix-Exterior Penalty Function (GIM-EPF) data processing algorithm. The GIM-EPF algorithm greatly improves the processing speed for MWP data while ensuring good temperature accuracy. Additionally, the algorithm does not require any advanced knowledge of the emissivity.
2. Principles of algorithm
2.1 Reference temperature model
If MWP has n spectral channels, output signal Vi of the ith channel is defined as follows:
In the case of reference temperature T´ of blackbody, the ith channel output signal Vi´ is given by:
This model only needs to measure output Vi´ of each channel under reference temperature T´, as long as reference temperature T´ is stable, it does not affect calculation results of target true temperature T and spectral emissivity ε(λi, T). The reference temperature model are simple and its condition is easy to meet, and error is less than other models such as brightness temperature model and calibration constant model, so it is widely used in data processing of MWP.
2.2 Combination of data processing and generalized inverse matrix of MWPEq. (4) is known, right is unknown. Order , , , , .
So Eq. (4) is transformed as follows:
abbreviated isEquation (7) is an underdetermined equation, it has infinitely groups of solutions.
Because A is not a square matrix, there is no inverse matrix A−1. In order to solve such problems, n × n matrix inverse is extended to n × m matrix inverse, named generalized inverse matrix. Because solution of Eq. (7) is not unique, we further restrict the generalized inverse matrix, that is, Moore-Penrose inverse matrix is used to guarantee the uniqueness of solution.
The matrix G that satisfies the following properties, which is defined as the Moore-Penrose inverse matrix of matrix A, and recorded A+, it is given as follows:
If matrix G satisfies all four formulas of Eq. (8) at the same time, then A+ is unique. For m × n matrix A, when rank r = n (that is, A is full column rank),
when rank r = m (that is, A is full row rank)
2.4 Generalized inverse matrix algorithm simulation for initial emissivity
According to generalized inverse matrix idea, simulation experiment is performed on six kinds of target materials labeled A-F, which are representative of different change trend of emissivity. The spectral emissivity of A-F targets is shown in Fig. 1 (blue ‘*’ market), true temperature is 1800K, reference temperature of blackbody is 1600K, 8 effective wavelengths are 0.4μm, 0.5μm, 0.6μm, 0.7μm, 0.8μm, 0.9μm, 1.0μm, 1.1μm, respectively.
As can be seen from Fig. 1, trend of six target emissivity inversion results is basically the same as the actual distribution, but inversion values are all larger than true values, and the inversion values of partial emissivity are more than 1, which is not in accordance with the actual situation. It can be seen from Table 1 that the mean absolute error of temperature inversion results of six targets is 61.8K at 1800K, average relative error is 3.44%, and accuracy of temperature inversion results is not high. However, the emissivity inversion results based on generalized inverse matrix can be used as initial point of constrained optimization algorithm, so as to improve further accuracy and efficiency of inversion.
2.3 Combination of data processing and exterior penalty function of MWP
According to Eq. (3), ideally inversion temperature T of every spectral channel should be equal to true temperature if ε(λi, T) is known. Therefore, theoretical deviation of calculated temperature should tend to zero for every channel, so it is represented by the mathematical expression:Eq. (12) can be converted into constrain optimization problem by minimized the deviation and made it infinitely close to zero, that is:
After performing logarithmic operations on both sides of Eq. (3), it is given:
Order , ,
then inversion temperature Ti of every channel and E(Ti) can be expressed as follows:
Since emissivity range is [0, 1], that is 0≤ε(λi, T)≤1, the constraint condition is xi<0.
The combination of Eq. (17) and xi<0, then constraint optimization problem of MWP is defined as follows:
This is a standard constraint optimization problem, Eq. (18) can be resolved by some constraint optimization algorithms such as Internal Penalty Function (IPF) method. However, from the results obtained from generalized inverse in Fig. 1, we can see that partial emissivity value obtained by Eq. (11) has been beyond the range of [0,1], IPF is not applicable. Therefore, exterior penalty function constraint algorithm must be used, emissivity starter value may be out of the [0, 1] range for it.
For inequality constrained optimization problemEq. (18) can be conformed into Eq. (19), it can be calculated using exterior penalty function method to process data of MWP.
Programming in MATLAB to achieve the exterior penalty function: minPF.
Call format: [x, minf] = minPF(f, x0, g, c1, p, var, eps)
Where x is independent variable; minf is the minimum value of objective function; f is objective function; x0 is initial point; g is constraint function; c1 is initial value of penalty factor; p is ratio coefficient of penalty factor; var is argument vector; eps is precision.
2.4 Parameters determination of minPF
Some parameters of minPF need to be determined before utilizing exterior penalty function algorithm.
- (1) Initial point x0
Initial point may be any point for exterior penalty function algorithm, but if the initial point is near the true value, calculation efficiency will be much improved. So calculation result from generalized inverse matrix is as the initial point according to Eq. (11), it would be an excellent initial emissivity, so blind choice of initial emissivity points can be avoided.
- (2) Objective function f
From Eq. (15), we can see that penalty condition of exterior penalty function is gj(x) ≥0, i = 1,2,…k, k ϵ Z, where gi(x) is a constraint function, for Eq. (18), is the constraint function. After many simulation experiment, emissivity is an idea range from 0.4 to 0.9, that is 0.4≤εi≤0.9, then −0.9≤xi≤-0.1. Finally, inequality group, that is is used to match constraint function of Eq. (15). Therefore, two matrices are given as follows:
Then is as constraint condition for data processing problem of MWP.
- (3) Initial value c1 and ratio coefficient p
Penalty factor P (Eq. (20) of exterior penalty function algorithm is controlled by initial value c1 and ratio coefficient p, and it has effect on precision and efficiency of algorithm. To make the unconstrained problem close to the original constraint problem, we should choose the larger possible penalty factor, but in order to reduce the difficulty of solving the unconstrained problem, the smaller penalty factor should be selected, so how to select an appropriate penalty factor is very important. Generally, c1ϵ[0.05, 0.5], p ϵ [2,5],in order to determine the optimal c1 and p, following simulation data are analyzed for 6 targets in Fig. 1.
Firstly, keeping unchanged value for p (when p = 2), taking c1 = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5 separately to calculate temperature for A-F six targets.
After c1 = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5 was selected respectively, but c1 has no effect on inversion results of temperature and emissivity, and average relative error is controlled at about 0.56% for temperature (Table 2), inversion trend of emissivity is also close to the actual value (Table 3), it showed a good inversion result from GIM-EPF algorithm.
Secondly, affection on running time from different c1 is considered by GIM-EPF algorithm when p = 2.
Running time of GIM-EPF algorithm in the 6 targets is statistically analyzed. It is known from Table 4 that when c1 = 0.05 is used, running time is shortest and average time of 6 targets is less than 5s.
Thirdly, affection on minimum of target function minf = Eq. (17) from different c1 is calculated by GIM-EPF algorithm when p = 2.
Fourthly, keeping unchanged value of c1 (when c1 = 0.05), taking p = 2, 3, 4, 5 separately to calculate temperature and emissivity for A-F targets by GIM-EPF.
After p = 2, 3, 4, 5 was selected respectively, but p has no effect on inversion results of temperature and emissivity, and average relative error is controlled at about 0.56% for temperature (Table 6), inversion trend of emissivity is also close to the actual value (Table 7), it showed a good inversion result using GIM-EPF algorithm.
Fifthly, running time from different p is considered by GIM-EPF algorithm when c1 = 0.05.
Running time of GIM-EPF algorithm for 6 targets is statistically analyzed. It is known from Table 8 that running time is shortest and average time of 6 targets is less than 5s when p = 2 is used.
Lastly, affection on minimum of target function minf = Eq. (17) from different p is calculated by GIM-EPF algorithm when c1 = 0.05.
To sum up, when c1 = 0.05, p = 2, both emissivity and temperature inversion results can meet requirements very well with the shortest running time and the highest calculation precision, so initial value of penalty factor c1 and ratio coefficient of penalty factor p of GIM-EFP algorithm parameters are c1 = 0.05, p = 2.
3. Comparison between GIM-EPF algorithm and others
In order to verify accuracy, efficiency and noise resistance of GIM-EPF algorithm, GIM-EPF and other constraint optimization algorithms in , such as GP and IPF, are compared respectively under the condition of no noise and 5% random noise on Vi in Eq. (1).
3.1 Comparison results with no noise
From Fig. 2, it can be seen that accuracy of emissivity inversion of GIM-EPF algorithm is higher than that of GP algorithm, which is very closer to that of IPF algorithm. Even in the case of complex emissivity distribution of F material, it is superior to IPF algorithm, which shows that accuracy of emissivity inversion of GIM-EPF algorithm is high.
It can be seen from Table 10 that temperature inversion accuracy of GIM-EPF algorithm is better than that of GP algorithm, which is equivalent to IPF algorithm, and average relative error is 0.56% for A-F materials at 1800K. Therefore, accurate inversion of temperature can be achieved based on GIM-EPF algorithm.
As can be seen from Fig. 3, inversion efficiency is improved greatly by GIM-EPF algorithm, average running time is 4.0s (Operating environment: Intel(R) Core(TM)i5-4210U CPU@1.70GHz 2.40GHz), and average calculation time of IPF algorithm is 34.6s, that of GP is 561.5s, GIM-IPF efficiency is improved 8.7 times higher than that of IPF, and it is more than 140 times higher than that of GP algorithm. This is just because initial value obtained by using generalized inverse matrix idea is brought into EPF, which greatly reduces searching scope of emissivity, thus greatly improving computing efficiency. So if computer performance is better, on-line and high precision data processing of MWP can be realized by using GIM-EPF algorithm.
3.2 Comparison results with 5% random noise
It can be seen from Fig. 4 that in the case of 5% random noise on Vi in Eq. (1), GIM-EPF algorithm can still obtain a better inversion result of emissivity. The new algorithm has a much higher inversion accuracy than GP and IPF algorithm, especially for the E and F materials with large emissivity-wavelength fluctuation, GIM-EPF algorithm has the best results, it shown that GIM-EPF algorithm is good at anti-noise performance.
It can be seen from Table 11 that temperature inversion accuracy of GIM-EPF algorithm is better than that of GP and IPF algorithm with 5% random noise on Vi in Eq. (1), and average relative error is 0.56% for A-F materials at 1800K, especially for the E and F materials with large emissivity-wavelength fluctuation, GIM-EPF algorithm is superior to IPF, it indicated that new algorithm can be adapt to any material with good performance at anti-noise aspect. Therefore, novel GIM-EPF algorithm can be as the MWP data process method.
The most obvious advantage of GIM-EPF algorithm compared with other algorithms is that the calculation time is greatly reduced. From Fig. 5, it can be seen that in the case of 5% noise, average running time of GIM-EPF algorithm is only 4.34s, which is equivalent to the 4.01s without noise, but that of the IPF algorithm is 34.14s, and the GP algorithm is 241.3s. This proves that initial emissivity value based on generalized inverse matrix greatly improves searching efficiency of exterior penalty function algorithm, and it embodies excellent operational efficiency and good noise resistance ability for proposed GIM-EPF algorithm.
4. Inversion of rocket nozzle temperature by MWP
In order to further verify practicability of proposed GIM-IPF algorithm, Prof. SUN measured surface temperature of a rocket nozzle using 8 wavelengths MWP, and gave spectral channel voltage values at reference temperature 2252K and 12 continuous measuring time points of the rocket nozzle in . Using the proposed GIM-IPF algorithm, inversion temperature of each time point is calculated and compared with IPF algorithm.
The theoretical true temperature of the engine flame near the nozzle exit indicated by the rocket engine designer is 2490.0K. Figure 6 is result of temperature inversion of each measurement time node. It shows that GIM-EPF algorithm has a higher inversion accuracy. Compared with the IPF algorithm, although the new algorithm is slightly lower than IPF at some time points, the maximum relative error is only 0.65%, which reflects the higher temperature inversion accuracy than IPF. As can be seen from Fig. 7, GIM-EPF algorithm has a much shorter operation time than IPF algorithm. The average running time of GIM-EPF algorithm is only 3.4s, while the average running time of IPF algorithm reaches 210s. If hardware performance of a computer is improved, calculation time will be shorter. The new algorithm provides a theoretical basis for accurate on-line temperature measurement.
A generalized inverse matrix and exterior penalty function (GIM-EPF) approach was integrated into the processing of multi–wavelength pyrometry data without requiring any advanced knowledge of the emissivity model and without the need to constrain the range of the emissivity beyond requiring that the ln ε ≤ 1. The proposed GIM-EPF algorithm eliminates the influence of emissivity on the inversion results to determine temperature. The simulation results from six materials show that the average relative error at 1800K temperature is 0.56%, a temperature inversion accuracy almost equal to IPF. A important advantage over the IPF is the inversion speed. The calculation speed of the GIM-EPF algorithm with good anti-noise performance is 8 times faster than that of IPF algorithm and is more than 100 times faster than that of Gradient Projection (GP) algorithm. Data from rocket plume temperature measurement confirms that the GIM-EPF algorithm has good inversion accuracy and efficiency. In summary, the GIM-EPF algorithm greatly improves the processing efficiency of MWP data while ensuring good temperature accuracy. The algorithm does not require any advanced knowledge of the emissivity. The GIM-EPF algorithm thus greatly increases the utility of MWP for measurement of high and ultra-high temperatures.
National Natural Science Foundation of China (NSFC) (61405045, 31470714); Fundamental Research Funds for the Central Universities (2572017DB04); “Double First Class” Talents Project Funding of NEFU 2016; Science and Technology Innovation Talents from Harbin City (RC2014QN009026).
1. I. Rancdarbord, G. Baudin, M. Genetier, D. Ramel, P. Vasseur, J. Legrand, and V. Pina, “Emission of gas and Al2O3 smoke in Gas-Al particle deflagration: experiments and emission modeling for explosive fireballs,” Int. J. Thermophys. 39(3), 152–160 (2018).
2. T. Matsumoto, T. Koizumi, Y. Kawakami, K. Okamoto, and M. Tomita, “Perfect blackbody radiation from a graphene nanostructure with application to high-temperature spectral emissivity measurements,” Opt. Express 21(25), 30964–30974 (2013). [CrossRef] [PubMed]
3. S. V. Onufriev, “Measuring the temperature of substances upon fast Heating with a current pulse,” Bull. Russ. Acad. Sci., Physics 82(4), 372–379 (2018). [CrossRef]
4. A. Bendada and M. Lamontagne, “A new infrared pyrometer for polymer temperature measurement during extrusion moulding,” Infrared Phys. Technol. 46(1–2), 11–15 (2004). [CrossRef]
5. F. Benedic, P. Bruno, and P. Pigeat, “Real-time optical monitoring of thin film growth by in situ pyrometry through multiple layers and effective media approximation modeling,” Appl. Phys. Lett. 90(13), 134104 (2007). [CrossRef]
6. A. Seifter and D. C. Swift, “Pyrometric measurement of the temperature of shocked molybdenum,” Phys. Rev. B 77(13), 761–768 (2007).
7. T. Fu, Z. Wang, and X. Cheng, “Temperature measurements of a diesel fuel combustion with multicolor pyrometry,” J. Heat Transfer 132(5), 051602 (2010). [CrossRef]
8. S. Liu, P. Farahmand, and R. Kovacevic, “Optical monitoring of high power direct diode laser cladding,” Opt. Laser Technol. 64(4), 363–376 (2014). [CrossRef]
10. D. Ng and G. Fralick, “Use of a multiwavelength pyrometer in several elevated temperature aerospace applications,” Rev. Sci. Instrum. 72(2), 1522–1530 (2001). [CrossRef]
11. V. Scharf and A. Katzir, “Four-band fiber-optic radiometry for determining the “true” temperature of gray bodies,” Appl. Phys. Lett. 77(19), 2955–2957 (2000). [CrossRef]
13. D. Svet and S. Sergeev, “Triple-wavelength pyrometer that measures true temperature,” Meas. Tech. 54(11), 1273–1275 (2012). [CrossRef]
14. J. S. Pérez-Huerta, T. Saucedo-Anaya, I. Moreno, D. Ariza-Flores, and B. Saucedo-Orozco, “Digital holographic interferometry applied to the investigation of ignition process,” Opt. Express 25(12), 13190–13198 (2017). [CrossRef] [PubMed]
15. D. Han, A. Satija, J. Kim, Y. Weng, J. Gore, and R. P. Lucht, “Dual-pump vibrational cars measurements of temperature and species concentrations in turbulent premixed flames with CO2 addition,” Combust. Flame 181(3), 239–250 (2000).
16. G. Zhang, J. Liu, Z. Xu, Y. He, and R. Kan, “Characterization of temperature non-uniformity over a premixed CH 4 –air flame based on line-of-sight TDLAS,” Appl. Phys. B 122(1), 3–11 (2016). [CrossRef] [PubMed]
17. M. Liang, B. Sun, X. Sun, and J. Xie, “Development of a new fiber-optic multi-target multispectral pyrometer for achievable true temperature measurement of the solid rocket plume,” Measurement 95, 239–245 (2017). [CrossRef]
18. T. Fu, J. Liu, M. Duan, and S. Li, “Sub-pixel temperature measurements in hypersonic plasma jet environments using high speed multispectral pyrometry,” J. Heat Transfer 140(7), 1601–1609 (2018). [CrossRef]
19. P. Coates, “The least-square approach to multi-wavelength pyrometry,” High Temp. High Press. 20(4), 071601 (1988).
20. G. Gathers, “Analysis of multiwavelength pyrometry using nonlinear chi-square fits and Monte Carlo methods,” Int. J. Thermophys. 13(3), 539–554 (1992). [CrossRef]
21. M. A. Khan, C. Allemand, and T. W. Eagar, “Noncontact temperature measurement. I. Interpolation based techniques,” Rev. Sci. Instrum. 62(2), 392–402 (1991). [CrossRef]
22. H. Madura, M. Kastek, and T. Piatkowski, “Automatic compensation of emissivity in three-wavelength pyrometers,” Infrared Phys. Technol. 51(1), 1–8 (2007). [CrossRef]
23. P. Hagqvist, F. Sikström, A. Christiansson, and B. Lennartson, “Emissivity compensated spectral pyrometry for varying emissivity metallic measurands,” Meas. Sci. Technol. 25(2), 405–412 (2014).
24. X. G. Sun, G. B. Yuan, J. M. Dai, and Z. X. Chu, “Processing method of multi-wavelength pyrometer data for continuous temperature measurements,” Int. J. Thermophys. 26(4), 1255–1261 (2005). [CrossRef]
25. J. Xing, S. Cui, W. Qi, F. Zhang, X. Sun, and W. Sun, “A data processing algorithm for multi-wavelength pyrometry-which does not need to assume the emissivity model in advance,” Measurement 67(5), 92–98 (2015). [CrossRef]
27. J. Xing, B. Peng, Z. Ma, X. Guo, L. Dai, W. Gu, and W. Song, “Directly data processing algorithm for multi-wavelength pyrometer (MWP),” Opt. Express 25(24), 30560–30574 (2017). [CrossRef] [PubMed]