Data processing for simultaneous inversion of emissivity and temperature using improved CABCSMA and target-to-best DE algorithms in multispectral radiation thermometry (MRT)

Kaihua Zhang; Jingzheng Dong; Yanfen Xu; Kun Yu; Kun Yu; Yufang Liu; Yufang Liu

doi:10.1364/OE.500998

1. Introduction

Multispectral radiation thermometry (MRT) plays an important role in the field of radiation measurement by measuring the radiation intensity at different wavelengths emitted by the target object to infer the true temperature of the target, which has the advantages of fast response, no upper measurement limit, and no contact with the measured target [1–5], and has been widely used in aerospace, energy, metallurgy, materials and other fields [6–9]. The core is to construct a set of radiation equations containing N + 1 unknowns based on the data of N spectral channels of a multispectral thermometer [10,11], but the set of radiation equations is underdetermined and difficult to solve because the emissivity is unknown.

Traditional multispectral temperature measurement algorithms use least squares, solved equations, and automatic order finding methods, which use emissivity assumption models assuming known emissivity versus wavelength to obtain true temperature versus emissivity results [12,13]. However, when an unknown material is measured, its emissivity versus wavelength is unknown, so there is some blindness when using the emissivity assumption model. Therefore, when an unknown material is measured, other algorithms or methods are needed to determine the emissivity of the material as a function of wavelength in order to obtain accurate temperature measurements.

In recent years, some new algorithms and methods have been proposed and widely used, such as neural network-based and physical model-based inversion algorithms, which can better cope with the emissivity estimation problem when measuring unknown materials and improve the accuracy and reliability of temperature measurement. Professor Dai proposed the quadratic measurement method [14], which can calculate the temperature and emissivity of two consecutive temperature measurement points at the same time. However, the variation of the spectral emissivity of an object is often inconsistent at different measurement locations, so the quadratic measurement method is not applicable to real-time temperature measurements. Therefore, there is an urgent need to develop a new method for multispectral radiometric temperature measurement without the assumption of an emissivity model.

Since neural networks have the function of autonomous learning and fast search for optimal solutions, which can be used for multispectral radiometric temperature measurement to solve for the true temperature of the target. Yang et al. proposed a combined neural network emissivity model to achieve the identification of continuous spectral emissivity and true temperature of an object based only on the measured luminance temperature data [15]. Xi et al. developed a multispectral temperature measurement method based on radial basis function neural network with high accuracy in steel temperature measurement when the spectral emissivity is unknown [16]. Chen et al. established a multispectral temperature measurement method based on an adaptive emissivity model, and combined neural network and genetic algorithm to obtain high-precision temperature inversion results [17]. However, training a neural network model requires a large number of valid data samples and a large amount of time to process the data samples beforehand.

At present, some algorithms that do not require pre-assumption of emissivity models have been proposed. Wang et al. proposed a method to establish constraint conditions based on the development trend of emissivity wavelength curve [18]. The maximum error obtained from the inversion is less than 1%. Y. Zhang established a multi-objective function based on the mathematical model of the reference brightness temperature [19]. Similarly, the maximum error obtained from the inversion is less than 1%. Xing et al. tried the gradient projection (GP) and the internal penalty function (IPF) constrained optimization algorithm [20], in which the IPF function with better results has a maximum error of 25 K and a computation time of 210 s; Liang et al. used the generalized inverse matrix-external penalty function (GIM-EPF) data processing algorithm to invert the temperature of the rocket nozzle [21] with a maximum relative error of 0.65% and a computation time of 3.4 s; K Yu et al. used the Broyden-Fletcher-Goldfarb-Shanno (BFGS) iterative algorithm to optimize the nonlinear objective function and obtained results with a maximum error of less than 0.5% and a computation time of 0.2 s [22]. Subsequently, Tian et al. proposed two new model optimization algorithms, sequential random coordinate shrinkage (SRCS) and multi-population genetics (MPG) [23], and the results showed that the MPG algorithm has an inversion speed of 0.36 s and a maximum relative error of 0.4%. However, in the above studies, the objective function of multispectral temperature measurement falls into local minima, the error of calculation may still be large, and the speed of calculation is still not satisfactory.

In this paper, a new combined algorithm of artificial bee colony and slime mould algorithm (CABCSMA) and a differential evolution (DE) algorithm using a target to best mutation strategy are proposed to process the data based on the mathematical model of reference temperature, so as to achieve the simultaneous inversion of temperature and emissivity. Simulations of six classical emissivity model materials, tungsten samples and rocket motor nozzle temperatures were performed to verify the accuracy and efficiency of the new algorithms.

2. Principles of MRT

2.1 Reference temperature model

According to Planck's law of radiation, for a multispectral radiation thermometer with n wavelength channels, the output signal V_i of the ith channel is as follows:

(1)$${V_i} = {A_{{\lambda _i}}}\varepsilon ({\lambda _i},T)\frac{1}{{\lambda _i^5({e^{{C_2}/{\lambda _i}T}} - 1)}}(i = 1,2,3,...,n)$$

where A_λi is an identification constant that is only related to wavelength but not temperature, which is affected by the spectral responsivity of the detector, the transmittance of the optical element, the geometry of instrument and the first radiation constant. ε(λ_i, T) is the spectral emissivity of the material at the true temperature T, λ_i is the effective wavelength of the ith spectral channel, and C₂ is the second radiation constant.

When C₂ / (λ_iT) << 1, Eq. (1) can be replaced by the Wien approximation as follows:

(2)$${V_i} = {A_{{\lambda _i}}}\varepsilon ({\lambda _i},T)\lambda _i^{ - 5}{e^{ - \frac{{{C_2}}}{{{\lambda _i}T}}}}(i = 1,2,3,\ldots ,n)$$

At the blackbody reference temperature T´, the output signal V_i´ of the ith spectral channel is:

(3)$$V_i^{\prime} = {A_{{\lambda _i}}}\varepsilon ({\lambda _i},{T^{\prime}})\lambda _i^{ - 5}{e^{ - \frac{{{C_2}}}{{{\lambda _i}T}}}}(i = 1,2,3,\ldots ,n)$$

where ε(λ_i,T´) is the emissivity of the blackbody, which is usually identified as 1.

According to the ratio of Eq. (2) and Eq. (3), the reference temperature model can be established:

(4)$$\frac{{{V_i}}}{{V_i^{\prime}}} = \varepsilon ({\lambda _i},T){e^{\frac{{{C_2}}}{{{\lambda _i}}}(\frac{1}{{{T^{\prime}}}} - \frac{1}{T})}}$$

The ingenuity of the model lies in the use of the blackbody as the reference point, while eliminating the identification constant, greatly simplifying the operation. Secondly, it only needs to measure the voltage output signal of each channel at any reference temperature. The calibration is relatively simple, and no matter what value is selected for the reference temperatures, the calculation results will not be affected. However, the number of unknown parameter is always one more than the number of equations in MRT, so the problem of solving the underdetermined system of equations is transformed into a constrained optimization problem.

2.2 Constrained optimization for MRT

Constrained optimization problems are a class of mathematical optimization problems, which consist of two parts: an objective function and constraints related to variables in the objective function, and the optimization process is to optimize (maximize or minimize) the objective function under the constraints. In this paper, the minimum of the objective function is sought and the constrained optimization problem is described as follows:

(5)$$\left\{ \begin{array}{l} \min f(x)\\ Ax \ge b \end{array} \right.$$

where x is the decision variable, i.e., the variable to be optimized, A is the constraint coefficient matrix, and b is the constraint vector.

From Eq. (4), it can be seen that the calculated temperature of each channel is exactly the same and equal to the real temperature when the emissivity is known. Therefore, the ideal deviation of the calculated temperature of each channel is zero, which leads to the following equation:

(6)$$\frac{1}{n}\sum\limits_{i = 1}^n {\frac{{T_i^2}}{{{E^2}({T_i})}}} = 1$$

where T_i is the calculated temperature of each channel and the expression is:

(7)$${T_i} = \frac{1}{{\frac{1}{{{T^{\prime}}}} + \frac{{{\lambda _i}}}{{{C_2}}}[\ln \varepsilon ({\lambda _i},T) - \ln (\frac{{{V_i}}}{{V_i^{\prime}}})]}}$$

E(T_i) is the average calculated temperature of all spectral channels, and its expression is:

(8)$$E({T_i}) = \frac{1}{n}\sum\limits_{i = 1}^n {{T_i}} $$

Since the emissivity is unknown, Eq. (6) can be transformed into an optimization problem based on the principle of minimization of deviations that the temperature deviation tends infinitely to zero as follows:

(9)$$\min F = \left|{(\frac{1}{n}\sum\limits_{i = 1}^n {\frac{{T_i^2}}{{{E^2}({T_i})}}} ) - 1} \right|\to 0$$

In multispectral radiometry, the values of spectral emissivity belong to [0, 1], i.e., 0 < ε (λ_i, T) < 1, and note x_i= lnε(λ_i,T_i). According to Eq. (5), it can be described as:

(10)$$\left\{ \begin{array}{l} \min F = \left|{(\frac{1}{n}\sum\limits_{i = 1}^n {\frac{{T_i^2}}{{{E^2}({T_i})}}} ) - 1} \right|\\ {x_i} \le 0 \end{array} \right.$$

Equation (10) is a standard constrained optimization problem, where F is the objective function and x_i≤ 0 is the equation constraint. CABCSMA and DE algorithms can be used to solve such problems.

3. Principle of CABCSMA algorithm

3.1 Principle of SMA algorithm

Smile mould algorithm (SMA) is an intelligent optimization algorithm proposed based on simulating the oscillatory predatory behavior of individual slime molds [24]. Slime molds approach food through the smell in the air, and their approximation behavior is expressed by:

(11)$$X({t + 1} )\left\{ {\begin{array}{{c}} {{X_b}(t )+ {v_b} \cdot ({W \cdot {X_A}(t )- {X_B}(t )} ),r < p}\\ {{v_c} \cdot {X_t},r \ge p} \end{array}} \right.$$

where X_b(t) is the current optimal individual position, X_A(t) and X_B(t) are the positions of two randomly selected individuals. v_b and v_c are the control parameters, where v_b${\in} $[-a,a], v_c decreases linearly from 1 to 0. r is a random number between [0,1], p = tanh|S(i)-DF|, where i${\in} $1,2,3,…,n. S(i) is the current individual fitness value, DF is the current1 best fitness value, a = actanh(-(t/t_max) + 1), t is the current number of iterations and t_max is the maximum number of iterations. W is the slime mold quality, which represents the fitness weight, and its equation is:

(12)$$W({SI(i )} )\left\{ {\begin{array}{{c}} {1 + r \cdot \log \left( {\frac{{bF - S(i )}}{{bF - wF}} + 1} \right),i = C}\\ {1 - r \cdot \log \left( {\frac{{bF - S(i )}}{{bF - wF}} + 1} \right),i = O} \end{array}} \right.$$

(13) $$SI(i )= sort(S )$$

where r is a random number of [0,1] indicating the best fitness of the current iteration, wF is the worst fitness value of the current iteration, i = C indicates the individual with the top half of the fitness value, i = O indicates the remaining individuals, and SI(i) is the sorted sequence of fitness values.

The enveloping food phase simulates the contraction pattern of the slime molds venous tissue, and the mathematical formula for updating the position of the slime molds in this phase is shown as:

(14)$$X({t + 1} )= \left\{ {\begin{array}{{c}} {rand \cdot ({UB - LB} )+ LB,r < z}\\ {{X_b}(t )+ {v_b} \cdot ({W \cdot {X_A}(t )- {X_B}(t )} ),r < p}\\ {{v_c} \cdot {X_t},r > p} \end{array}} \right.$$

where UB and LB denote the upper and lower boundaries of the search area, respectively, rand denotes a random number taking values between [0,1], and z is the proportion of randomly distributed smile mold individuals to the total number of smile molds.

3.2 Principle of ABC algorithm

The artificial bee colony algorithm (ABC) [25] is an optimization method proposed based on imitating the behavior of honeybees. The algorithm achieves the processing of complex optimization problems by simulating the honey harvesting behavior of natural bee colonies. For the problem that SMA algorithm is easy to converge early, this paper adds the artificial swarm search strategy, and the basic artificial swarm search strategy mathematical model is described as shown in Eq. (15):

(15)$${Z_{i,j}} = {x_{i,j}} + {\phi _{i,j}}({{x_{i,j}} - {x_{k,j}}} )$$

where Z_i,j is the generated candidate solution, x_i,j is the current individual, x_k,j is a random individual, k and j are random parameters, k${\in} ${0,1,…,M}, j${\in} ${1,2,…,d}, M is a fixed value, d denotes dimensionality, and k is not equal to i, ϕ_i,_j takes the value of random numbers of [-1,1]. From Eq. (15), the candidate solutions are generated by randomly selecting two individuals for the difference operation.

Although the manual swarm search strategy has a great advantage in search capability, its exploitation capability is not strong. Therefore, this paper introduces a new search strategy that adds global optimal position guidance to the strong search capability of the manual swarm so as to improve its exploitation capability, and the mathematical model of the improved strategy is described as follows:

(16)$${Z_{i,j}} = {x_{i,j}} + {\phi _{i,j}}({{x_{i,j}} - {x_{k,j}}} )+ {\mathrm{\Omega }_{i,j}}({{P_{g,j}} - {x_{i,j}}} )$$

where Ω is a random number taking values in [0,1.5] and P_j is the global optimal position.

3.3 Principle of CABCSMA algorithm

In this paper, the individuals generated by the improved artificial swarm strategy are introduced in the SMA iteration process, and the greedy strategy is used to retain the better individuals among them to speed up the algorithm. While the powerful search ability of the artificial swarm search strategy can reduce the influence of local extremum points on SMA, thus improving the ability of the algorithm to jump out of the local optimal solution. The main steps of the CABCSMA algorithm are as follows:

Step 1: Parameters initialization: Set the search upper bound UB, lower bound LB, population size N, and maximum number of iterations.

Step 2: Initialize populations: Generate initial populations randomly in the search space.

Step 3: The fitness of each individual in the population is calculated and ranked, and the best fitness bF and the worst fitness wF are recorded.

Step 4: The candidate solution positions are updated according to Eq. (14).

Step 5: An artificial swarm search strategy is introduced to retain the better individuals according to the greedy strategy.

Step 6: Determine the end condition, if the iteration condition is satisfied, the algorithm terminates, otherwise return to step 3.

Step 7: The optimal value is output and the algorithm ends.

The flow chart of the CABCSMA algorithm is shown in Fig. 1:

Fig. 1. Flowchart of the CABCSMA algorithm.

Target	0.4 µm	0.5 µm	0.6 µm	0.7 µm	0.8 µm	0.9 µm	1.0 µm	1.1 µm
A	0.85	0.75	0.67	0.60	0.55	0.50	0.48	0.45
B	0.45	0.48	0.50	0.55	0.60	0.67	0.75	0.85
C	0.45	0.55	0.65	0.75	0.74	0.65	0.55	0.45
D	0.85	0.75	0.65	0.55	0.54	0.65	0.75	0.80
E	0.85	0.65	0.55	0.65	0.84	0.65	0.55	0.50
F	0.50	0.55	0.65	0.84	0.65	0.55	0.65	0.85

		A	B	C	D	E	F
900 K	Inversion Temperature (K)	899.81	897.10	898.42	899.89	899.90	898.17
	Absolute Error (K)	0.19	2.90	1.58	0.11	0.10	1.83
	Relative Error (%)	0.02	0.32	0.18	0.01	0.01	0.20
1200 K	Inversion Temperature (K)	1199.62	1194.36	1197.68	1199.62	1200.03	1196.82
	Absolute Error (K)	0.38	5.64	2.32	0.38	0.03	3.18
	Relative Error (%)	0.03	0.47	0.19	0.03	0.01	0.26
1500 K	Inversion Temperature (K)	1499.63	1491.77	1495.32	1499.40	1500.16	1495.11
	Absolute Error (K)	0.37	8.23	4.68	0.60	0.16	4.89
	Relative Error (%)	0.02	0.55	0.31	0.04	0.01	0.33
1800K	Inversion Temperature (K)	1799.08	1787.70	1794.18	1799.12	1800.15	1793.01
	Absolute Error (K)	0.92	12.30	5.82	0.88	0.15	6.99
	Relative Error (%)	0.05	0.68	0.32	0.11	0.01	0.39

T = 1800K	A	B	C	D	E	F
Inversion Temperature (K)	1799.04	1787.72	1795.86	1799.55	1801.47	1794.69
Absolute Error (K)	0.96	12.28	4.14	0.45	1.47	5.31
Relative Error (%)	0.05	0.68	0.23	0.03	0.08	0.29

		A	B	C	D	E	F
900 K	Inversion Temperature (K)	901.12	898.15	900.48	902.09	901.79	899.96
	Absolute Error (K)	1.12	1.85	0.48	2.09	1.79	0.04
	Relative Error (%)	0.12	0.01	0.05	0.23	0.20	0.01
1200 K	Inversion Temperature (K)	1201.83	1196.63	1201.07	1203.74	1203.04	1199.80
	Absolute Error (K)	1.83	3.67	1.07	3.74	3.04	0.20
	Relative Error (%)	0.15	0.28	0.09	0.31	0.25	0.02
1500 K	Inversion Temperature (K)	1502.85	1494.99	1501.78	1505.89	1504.63	1499.56
	Absolute Error (K)	2.85	5.01	1.78	5.89	4.63	0.44
	Relative Error (%)	0.19	0.33	0.12	0.39	0.31	0.03
1800K	Inversion Temperature (K)	1803.81	1792.43	1801.75	1807.73	1806.32	1799.26
	Absolute Error (K)	3.81	7.57	1.75	7.73	6.32	0.74
	Relative Error (%)	0.21	0.42	0.10	0.43	0.35	0.04

T = 1800K	A	B	C	D	E	F
Inversion Temperature (K)	1803.76	1793.14	1802.86	1806.59	1806.19	1799.81
Absolute Error (K)	3.76	6.85	2.86	6.59	6.19	0.19
Relative Error (%)	0.21	0.38	0.16	0.37	0.34	0.01

Temperature (K)		Wavelength (µm)
		1.434	1.525	1.635	1.706	1.856	1.934	2.061	2.154
Blackbody	923	864.04	1477.52	2589.98	3583.74	5582.90	7277.94	8630.16	9276.10
Silicon Carbide	973	1237.36	2047.60	3480.65	4694.85	7237.20	9117.03	10504.55	11215.13
	1023	2036.30	3273.43	5391.97	7142.47	10650.73	13210.3	14900.20	15669.23
	1073	3206.05	5015.80	8041.75	10474.5	15179.05	18529.8	20470.80	21258.15

Temperature(K)	Wavelength (µm)
Temperature(K)	1.434	1.525	1.635	1.706	1.856	1.934	2.061	2.154
973	0.8179	0.8152	0.8155	0.8145	0.8375	0.8217	0.8131	0.8193
1023	0.8316	0.8298	0.8295	0.8287	0.8467	0.8335	0.8276	0.8331
1073	0.8349	0.8316	0.8321	0.8305	0.8471	0.8343	0.8288	0.8337

	Sample Temperature (K)	Inversion Temperature (K)	Absolute Error (K)	Relative Error (%)	Calculation Time (s)
DE	973	971.56	1.44	0.15	0.06
	1023	1021.53	1.46	0.14	0.06
	1073	1069.97	3.03	0.28	0.06
CABCSMA	973	971.27	1.73	0.18	0.43
	1023	1018.77	4.23	0.41	0.45
	1073	1070.75	2.25	0.21	0.44

Temperature (K)	Wavelength (µm)
Temperature (K)	0.380	0.467	0.665	0.800	1.000	1.500	1.800	2.000
2200	0.482	0.466	0.431	0.410	0.378	0.290	0.225	0.201
2800	0.473	0.458	0.419	0.399	0.371	0.299	0.254	0.233
3400	0.465	0.450	0.407	0.388	0.363	0.308	0.283	0.265

	Sample Temperature (K)	Inversion Temperature (K)	Absolute Error (K)	Relative Error (%)	Calculation Time (s)
DE	2200	2206.50	6.50	0.30	0.06
	2800	2804.59	4.59	0.16	0.06
	3400	3398.26	1.73	0.05	0.06
CABCSMA	2200	2206.90	6.90	0.31	0.44
	2800	2804.06	4.06	0.14	0.43
	3400	3402.11	2.11	0.06	0.44

Abstract

1. Introduction

2. Principles of MRT

2.1 Reference temperature model

2.2 Constrained optimization for MRT

3. Principle of CABCSMA algorithm

3.1 Principle of SMA algorithm

3.2 Principle of ABC algorithm

3.3 Principle of CABCSMA algorithm

3.4 CABCSMA algorithm simulation

4. Differential evolutionary algorithm based on target-to-best variation strategy

4.1 Principle of DE algorithm

4.2 DE algorithm simulation

5. Comparison of CABCSMA and DE algorithms

6. Experiment

6.1 Simulation of silicon carbide

6.2 Simulation of tungsten

6.3 Simulation of rocket nozzle

7. Effectiveness discussion

8. Conclusions

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (15)

Tables (10)

Equations (16)

Optics Express