Simulation of an NSGA-III Based Fireball Inner-Temperature-Field Reconstructive Method

For the inner-temperature-field reconstruction of a fireball, a detecting method was proposed, using multi-channel visible spectral remote sensing theory. In our proposed method, the reconstructive algorithm based on multi-channel-detection was considered as a multi-objective optimization problem (MOOP), and a fast non-dominated sorting genetic algorithm based on reference-point strategy (NSGA-III) was employed as the solution of this problem. Besides, a so-called ambient pressure operator, based on the unique detecting model, was proposed and employed during the iteration process, for dynamic genetic parameter adjustment. To verify some performance of our proposed method, several numerical reconstructive simulations were carried out, from simple GA to NSGA-III, using several artificial 2-D virtual data, with different kinds of crossover and mutation functions. The simulation results show that, limited to our problems, the NSGA-III can effectively reconstruct different 2-D data, by a fixed crossover rate and a dynamic mutation rate under the proposed ambient pressure operator and an adaptive mutation rate function. The algorithm, limited to the field-reconstruction, also performs well on stability, but still has some deficiencies to be optimized.


I. INTRODUCTION
In order to obtain the inner-temperature distribution of a fireball, several non-contact and contact detection method have been developed. Our work in fund project is using several ultra-high-speed visible spectral cameras to gain multi-channel images of the target projection, and then using the images to reconstruct the inside temperature field, where the 3-D reconstructive algorithm is urgently required after the image acquisition. Our currently works were, besides the development of the optical detection system, proposing a NSGA-III based reconstructive method, simplifying the projection model and using several artificial 2-D virtual data to prove the effectiveness of our algorithm. In this section, the optical detection system, multi-channel-images based The associate editor coordinating the review of this manuscript and approving it for publication was Paul Yoo .

A. THE VISIBLE LIGHT DETECTION SYSTEM FOR FIREBALL
The detection system proposed in our project for fireball inner temperature field, was considered having several high-speed visible spectral cameras. As is shown in figure 1, these cameras shall be implemented around the blast fireball, outside the dangerous distance, to obtain the projections from different directions of the fireball, geometric optical theory followed, and the storage module will record each frame of the projection data, collected in a certain direction. After sampling for a period of time, reprocess the collected images, and a reconstructive field will finally be available. There are two important hypotheses in the detection system, and both play important roles in the entire reconstruction theory.
Firstly, as is shown in figure 2, each camera in our detection system, must have different trigger and storage time, which  could be efficiently reduced by modern and expensive techniques. However, when developing the reconstructive algorithm, assuming that the camera array having a sufficiently high sampling frequency and synchronization hypothesis is quite helpful, by which images from different cameras in similar frames can consider to be captured at a same time.
The second hypothesis is the simplification of the optical transfer function (OTF), which we call it Temperature-Intensity relation, converting the inner temperature to the gathered projection. The actual overall relation, which considered quite complicated, is currently simplified as three sub relations, as is briefly described in figure 3, which are thermal excitation of elements, reduction of the radiation and receiving property of our multi-channel sensors. In the second hypothesis above, shown in figure 3, based on thermal radiation theory, we consider that each position inside a transient fireball must have a determined temperature message, exciting the elements contained in this position and generating radiation with special spectrum [4], [5]. So, we can define formulation (1) to represent this process. I e,(x,y,z) = f excitation (x,y,z) , ε λ . (1) where (x,y,z) represents the temperature of one node inside the fireball area, whose coordinate value are (x, y, z) and emissivity is ε λ . By a certain function f excitation based on black-body radiation law, the emissive thermal radiation I e,(x,y,z) can finally be acquired.
In the next place, considering one of the radiation lines generated, four environmental conditions were considered, that are propagation inside the fireball, propagation through the vapor or smoke, propagation in the atmosphere, and propagation through the optical devices, causing the attenuation of the radiation energy. This attenuation process [6] can be defined as formulation (2).
where H represents the set of different attenuation coefficients. The intensity of radiation after attenuation I a,(x,y,z) can be calculated through the attenuation process f attenutation . Finally, the rest radiation I a will be converted according to the photoelectric conversion properties [7], [8], and f photoreception as formulation (3) can define the relation between the gray level D gray and the rest light intensity I a .
D gray = f photoreception (I a , t). ( where t is the integration time of pixels. From a different point, the gray level of one pixel might correspond, not to the temperature of a single node, but to the temperature of multiple nodes based on geometrical optics as figure 4. This relation plays an significant role in the reconstructive system which means the error generated in one pixel, must mainly come from the value of points related to this certain pixel. As a result, we can achieve expressions (4), describing the gray-level D (m.n) of each pixel in one entire image located at (m, n). So the gray level of the whole image can be described as D i , the set of all D (m.n) in the photographic film of channel i.
where any (x, y, z) ∈ V and V is the corresponding space area of No.(m, n) pixel in one detecting channel, i.e. (m, n) → V .
Summarize from formulation (1) to (4), the temperature of an inner point and the change of pixel gray level can be generally described as a general formulation (5).

B. MOOP BASED RECONSTRUCTIVE METHOD
Based on the images collected in section 1.1, the field reconstructive system can be described as the schematic shown in figure 5. At first, coordinate range of the fireball area should be calculated by surface 3D reconstructive algorithm according to the multi-channel images. The area calculated can further be divided into a fixed grid, with the initial temperature information given to the grid node. Then, by using the relation defined in expressions (5) and (6), multi-channel projections of the artificial temperature field above, can be calculated. Numerical difference between the projections calculated and the projections sampling by the camera array, must exist and could be calculated as formulation (7). Iterative method is chosen to adjust the temperature value of a node inside the inner fireball mesh, aiming at reducing the error calculated by the sum of squared errors (SSE) of all pixels in one image, between the calculated projection and the collected projection. By using this kind of SSE, one SSE can represents the fitness of the current iteration result in a certain detection channel, and the purpose of iteration is to find the best numerical temperature field solution, fitting all channels' SSE well, which is similar to a multi-objective optimization problem (MOOP) [9], [10] and is suitable for many multi-objective evolution algorithms (MOEA) [11], [12]. Finally, the best solution found in the process of iteration, will be the reconstructed result of the fireball inner temperature field.

II. IMPLEMENT OF NSGA-III TO THE RECONSTRUCTIVE ALGORITHM
Since the multi-channel-projection based field reconstructive method has been classified as MOOP, NSGA-III is chosen to solve this problem [13]. There are two significant points in implement of the algorithm. The first is the simplification of the detection model to build a simple mathematical model. The second is the combination of NSGA-III and the simplified model.

A. SIMPLIFIED MATHEMATICAL MODEL
In the primary stage of the reconstructive project, it is difficult to verify the reconstructive method is actually work, without having real detection data. So, a simplified mathematical model was designed to simulate the detection process, with three main concepts.
The first concept is trying 2-D simulations before 3-D simulations. Same to the 3-D situation, 2-D simulation can keep the corresponding relationship between the projection and internal value, only having difference in the dimension of variables in calculation. Guided by this theory, virtual 2-D fireball data and virtual visible spectral camera array was randomly created.
The second concept is the generation of 2-D fireball data sets. The fireball may have several high temperature regions inside the explosion area. So, we randomly choose three centers and three radii, generating several virtual fireball data sets, as is shown in  The third concept is using the principle of pinhole imaging, one of the simple but frequently-used OTF, in calculating projections. In 2-D situation, for one virtual pixel in figure 7, two vectors can be calculated by one virtual focus and the two boundary points of the pixel.
In the 2-D vector space, the two vectors can be the base vectors as they are linearly independent, and all points can be represented by the two vectors as formulation (6), where the virtual focus point can be described as O i = o x o y , one random object point can be described as P = p x p y , and the neighbouring two boundary points of one pixel are P 1,i = x 1,i y 1,i and P 2,i = x 2,i y 2,i .
According to the principle of pinhole imaging, the point is judged to project towards the pixel, when two coefficients k 1,i and k 2,i , derived from formulation (7-9), are both positive.
Based on the concepts above, figure 8 shows a 2-D 10-channel virtual detection system, with 40 pixels, implemented around a random 2-D artificial fireball.

B. APPLICATION OF NSGA-III IN THE RECONSTRUCTIVE ALGORITHM
In the optical model of the detection system, the Temperature-Intensity relation might quite complex, and maybe quite difficult to solve the backward relation. So genetic algorithm (GA) was chosen, only using the forward relation to participate iteration. Different from the simple genetic algorithm, NSGA use non-dominated sorting [14], [15] as the elitist strategy to choose better individuals in iteration, and drive the iteration based on GA theory [16], [18]. From NSGA to NSGA-III, population diversity of individuals has been taken into consideration in NSGA-II, and furthermore, NSGA-III uses a reference point association to effectively overcome this problem. The application of NSGA-III can be described as what is shown in figure 9.
Step1. Load an artificial 2-D fireball and set the virtual camera array by the given parameters, like what is shown in figure 8, which automatically generated in MATLAB by our own designed program. Then use the OTF virtual pinhole-imaging OTF as figure (7) and equations (7)(8)(9), to generate the standard projections set D s .
Step2. Use projections generated to produce a mean back-projection data Ps, and implement a fixed mutation rate p mutation on random location (x, y) of Ps to produce one initial individual P by equations (10).
where i represents the number of one channel, ranging from 1 to n, and mean back-projection Ps can be calculated as below.
Step3. Calculate each channel's projection D prj of each individual P with equations (12).
Compare D prj with standard projections D s in the first step, to get AE (i, j) and SE (i, j) matrix with equations (13), where i represents the number of one channel and j represents the number of one pixel in channel i.
Step4. Also, as is shown in figure 9, generate child-group by crossover operator and mutation operator, two basic GA operators. During the mutation, a linear mutation rate function [19] was used to control the number of mutation nodes, in order to make the algorithm being adaptive [20], [21]. The mutation rate p m could be calculated as equations (14).
where N is the number of objects (i.e. number of detective channels in this paper), p n can be calculated as follows.
where t and t m are the current and maximum generation, respectively; p is a fixed real number, and in this paper p can be given as follows. In order to control the intensity of the mutation, we proposed an ambient pressure operator based on the D-value between D prj (i, j) and D s (i, j) of the individuals in paternal generation. The aim of development of this operator, based on two theory. On the one hand, the final individuals must be similar to the mean of back-projections of each detecting channels, and all individuals in the iteration can gradually approach to the best individual along the curve, constructed by origin point and fitness point of back-projections' mean value, in this high dimension MOOP. On the other hand, though, the mutation position and mutation value should be and can be controlled by the current error between calculated and observed values, though they are generated by the current mutation rate.
The generation of this proposed ambient pressure operator is described as below.
Firstly, the number of nodes under mutation N m is controlled by the basic nodes' number N p and the mutation rate in this generation. The mutation number N m (i, j), in the corresponding region of pixel j of channel i, can be calculated by N p (i, j), the total number of nodes in the corresponding region of the pixel, and p m , the mutation rate in the current generation, as formulation (17).
Next, when considering the construction of the mutation value, the more error, the more mutation value. As equation (18), both the D-value between D prj (i, j) and D s (i, j), and a real number K are used to control the size of the mutation corresponding to pixel j of channel i. Besides, the direction of the mutation is defined as the reverse of the D-value. At last, average distribute the mutation value to N m (i, j) nodes.
The real number K is recommended in this paper as follow. This recommended value can restrict the total mutation value to less than 1, which corresponding to a certain pixel, making individuals can converge steadily along, floating in a small range near the mean back-projection Ps.
Meanwhile, a fixed rate from experience was chosen as follow to control the number of nodes under crossover process.
The results of different P c values were tested as is shown in figure 10. From figure 10(a), the algorithm can converge to the same value with different Pc and other same parameters. From figure 10(b), when set P c = 0.7, the algorithm can converge a little bit faster than other values. So,P c = 0.7 was chosen as a primary recommended crossover rate from experience. It is worth noting that the nodes under both crossover and mutation process are randomly generated in the range of the virtual fireball region, which is quite different from the method given in Jakub's paper [16], where all points in related region are chosen to participate in the iteration, as is shown in figure 11(a), which may lead to both the increase of operation and the rise of error. In the proposed fireball  Step5. After the basic genetic process, by merging both child-group and parent-group into one group, and achieving the SE matrix by step 3, as figure 12, the SSE matrix of any individual, inside the united group, can easily be acquired.
The SSE value of individual m towards channel i, i.e. sse (m, i), can be achieved by formulation (21).
where j is the serial number of pixels corresponding to channel i, and the structure of SSE matrix (i.e. SSE sort ) is as follow.
SSE sort is an M × N matrix, where M is the number of individuals and N represents the number of channels. Any row of (i.e. SSE sort (m, :) as below) represents this individual's fitness to each projection channel. SSE sort (m, :) = sse m1 sse m2 · · · sse mN (23) As a result, non-dominated sorting strategy could be used to elect elite individuals from the united groups, based on the multi-fitness vectors (i.e. matrix SSE sort ).
In the process of electing elite individuals, firstly, the reference-points association is implemented to keep the diversity of the current population [22]- [24]. In order to construct the reference points in the N-dimensional multifitness-vector domain, linear hyperplane method was introduced, whose expression can be described as equation (24).
where N is the dimension of the whole multi-fitness-vector space. The constant a can briefly be defined as 1, because that all the vectors above need to be standardized by multi-dimensional coordinate translation (i.e. equations (24)), and multi-dimensional normalization as equations (25).
where m i is the minimum value of column i in operator matrix T and T is an 1 × N matrix of the each channel's minimum value and UT is an M × N matrix, translated from SSE sort by the row matrix T constructed by m i in all dimension of the fitness. The elements in UT can be recognized as ut ij where i represents the channels and j represents the individuals. By selecting the maximum value of each column in UT to form the normalized parameter matrix, the normalization of matrix UT is realized, and the standardization of SSE sort is finally completed. The process is described as equations (25).
where operator matrix N saves reciprocal values of all channels' maximum SSE value. After dividing the columns of N by the maximum of each channel, the SSE matrix under normalization (i.e. SSE matrix standardized, SSE STD ) can easily achieved. After achieving the SSE STD matrix, in order to keep the diversity of individuals' distribution, try the best to relate at least one individual to each reference point, as is shown in figure 13, by comparing the distance between each reference vector with each fitness vector. The reference vector above is generated by the original point and reference points on the N-dimensional reference hyper-plane which constructed by equation (23).
In the equations above, SSE STD is an M × N matrix, ref i is one of the vector on the reference plane, defined as the equation (28) below.
There are three conditions while counting the number of related individuals as figure 13, and we try to realize that one reference vector only have one related vector. One of the useful solutions is using priority level, where non-dominated individuals is firstly taken into consideration, randomly select one individual from individuals related to a certain reference vector. If the number of individuals remain is not enough, we go back to the non-dominated sets removing in the former step, and randomly pick out several individuals. After that, if the individuals' number is still not enough, we have no choice but picking out individuals in under-dominated sets. Figure 14 shows the Pareto Front of the individuals during the iteration under a 10-channels reconstruction simulation, where the fitness of all individuals were separated with several 2-D or 3-D fitness by projection for visualization. The individuals were forming several niche, which could be observed easily during the iteration, which may be caused by the normalized reference points. Step6. The individuals picked out from the merging paternal group make up the filial elite set. The elite set, when their SSE matrix cannot satisfy our expectation, will be sent back to step 2 as the paternal individuals, to carry out a new round of iteration. The iteration process shall stop automatically only when the error of the individuals can satisfy the requirements or when the number of iterations is out of maximum bound.

III. SIMULATION RESULTS
Based on the simplified mathematical model and the reconstructive strategy above, several numerical experiments have been implemented, using the same projection mathematical model (i.e. pinhole imaging model with the same relative parameters) and same size of virtual fireball data, but under other different parameter conditions. The simulative tests of the proposed reconstructive method contain four different situations.
All of our simulations were carried out on a personal computer with a Intel(R) Core(TM) i7-8700K CPU, Kingston DDR4 3600Hz 32G RAM. The software platform are Windows 10 professional and MATLAB R2019a.
A. COMPARISON BETWEEN GA, NSGA AND NSGA-III The same virtual fireball and the same position parameters of all detecting channels, but different reconstructive theory and different number of the detecting channel were employed, in order to show that NSGA-III works completely better than simple GA and basic NSGA in this special MOOP. The reconstructive results of same artificial data by different method are compared in figure 15 with rising number of channels. From the visual qualitative analysis shown in figure 15, NSGA-III do have the obviously better reconstructive results than basic NSGA and Simple Genetic Algorithm. This may mainly because of the reference points method in NSGA-III, by which the diversity of individuals can stably be maintained and the local optimal solutions during the iteration can effectively be avoided.

B. COMPARISON BETWEEN FIXED RATE METHOD AND PROPOSED METHOD
Iterations with fixed parameters and proposed dynamic control method were tested under the same condition. The mutation rate has effect on both value and number of the individuals as equation (17). When using the same parametric condition and the same target, the mutation rate have a great effect on the convergent speed as is shown in figure 16. The explanation of this phenomenon is inferred as the influence of both the dynamic rate and the proposed operator. The proposed operator has the ability in applying more adjustment to the place with larger error.

C. COMPARISON BETWEEN DIFFERENT NUMBER OF CHANNELS
Based on the fact that NSGA-III has been proved to be the most suitable method in this multi-channel reconstructive problem, the same virtual data and different number of channels were used, to test what may be caused by different number of detective channel. The results of simulations with different number of channels are shown in figure 17, where the same color represents the same value. Also from the visual qualitative analysis, reconstructive results become better with the increasing number of detecting channels, but when over 10 channels, the degree of optimization becomes gradually not obvious.
From all elite individuals reconstructed with different number of channels, the mean of standard error (S.E.) curves shown in figure 18 and standard deviation (S.D.) curves these elite group can easily be saved by the iteration program.  All S.E. and S.D. were generated by a solitary program which compare the current results with the original artificial data as figure 17(a) during the iteration, and have no influence on accuracy in the iteration process.
In figure 19, the S.D. value of elite individuals' S.E. means that, based on the proposed ambient pressure operator,   the iteration process can keep a stable group diversity, and avoid local optimization in some degrees.
In the meanwhile, the mean S.E. curves in figure 15 indicates that the more number of channels, the better  reconstructive results. However, the more number of channels means the higher cost of calculating time. So, considering the artificial data only have three set value(i.e. 84, 169, and 255), the percentage of S.E. to the truth value (i.e. artificial data as figure 17(a)), was introduced to evaluate the reconstructive results with different number of channels as figure 20.
Considering the visible results in figure 14, the S.E. curves in figure 15, and the time cost, 10 or 12 channels are recommended in this proposed reconstructive theory. The first reason is, as figure 14(b) and (c), the reconstructive results figure 17, the reconstructive accuracy rise slowly but pay a high time cost in computation.

D. COMPARISON BETWEEN DIFFERENT TARGETS
The same number of detective channels and different virtual 2-D fireball data were used to prove the fitness of the reconstructive method when dealing with different targets. Over 30 artificial data were used in this simulation, and the results were all basically meet our expectations, five of our reconstructive results is shown as figure 21.
The relative mean value curves of elite individuals' S.E., generated by different artificial fireball data with 10 same detecting channels, are shown in figure 22, from which all targets can be observed converging to a similar S.E. value.

E. ALGORITHM STABILITY
Finally, all parameters were set in the same to conclude the repeatability of the proposed reconstructive method. The result is described in figure 23, in which the same number of channels can finally converge to the same S.E. value. However, the results also indicate that different number of channels may lead to a different converging results, the same as the results in subsection 3.2.

IV. CONCLUSION
The introduction of MOOP can solve the fireball-field reconstructive problem, and the reconstructive results of NSGA-III is much better than those without using NSGA theory. From the guide of the 2-D reconstructive results, 10 channels will primarily be used in the construction of detecting system. In addition, there are four advantages when implement NSGA-III in fireball-field reconstruction.
Firstly, this reconstruction theory we considered is suitable for different size of artificial data and 3-D reconstruction problem. The different size of artificial data may reduce the accuracy of reconstruction, but which can be solved by adding the number of pixels in detecting channel. In the meanwhile, there are two differences between 2-D and 3-D reconstruction, the rising of dimension and the change of OTF. All differences are suitable for the reconstructive theory and can easily be realized when programming. What needs to be stated in advance is the reason why the smallest size of 2-D artificial data (i.e. a 50 × 50 real matrix) was used in simulation experiments above, is reducing the computational period while calculating the high-dimensional MOOP by NSGA-III.
To the second, different from general computed tomography (CT), only forward relative equations are needed in the reconstructive system. Difficulty of operation can be reduced because that backward equations or solutions of the forward relative equations are not needed in the algorithm.
To the third, the limitation of iterated area by an accessible method (e.g. surface reconstructive method in 3-D condition) does have a better reconstructive result than Jakub's in the iteration. This may because the limitation can reduce the risk of correction value being distributed to the nodes outside the correct region while performing GA manipulation.
To the fourth, the use of ambient pressure operator based on the D-value between detecting value and calculated value of each unit in one detecting channel, can limit the mutation rate in a reasonable range as figure 16, keeping both the individuals' similarity to the mean of back-projections, and the diversity of the elite group. This is our first trial in using NSGA-III to solve field reconstruction problem based on many objective optimization which means there are still many works to do in the future. The iteration results and the performance of our algorithm can merely satisfy our primary needs, and whether this method is suitable to other complicated problem should be answered very precisely. Two main weaknesses in our research are as below.
Firstly, from both the reconstruction results and the iteration curves, we believe that the algorithm can be further optimized for better results and smaller error. The niche-preservation approach should be considered and the balance of convergence and diversity should be kept better. We plan to seek the solution to the limitation and weakness in this field-reconstruction problem by the above two significant points.
Secondly, We do agree with the idea that, whether the modified algorithm could achieve better performance than other state-of-the-art algorithms, should be tested on benchmark problems and using some statistical analysis as is mentioned in paper [25], which gives us a great idea to make our work more credible.
Finally, What to do in the future are the regularization of algorithm, the implement of real OTF and the rising of dimension(i.e. from 2-D simulation to 3-D simulation).
XUE BING received the B.S. degree in measurement and control technology and instruments from the North University of China, Shanxi, China, in 2017. He is currently pursuing the M.S. degree in instrument science and technology. His research interests include back-end data processing, multi-objective optimization problem, high-temperature detecting theory, and three-dimensional field reconstructive systems.