# Inverse Transient Radiative Analysis in Two-Dimensional Turbid Media by Particle Swarm Optimizations

Authors: Yatao Ren; Hong Qi; Qin Chen; Liming Ruan
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Three intelligent optimization algorithms, namely, the standard Particle Swarm Optimization (PSO), the Stochastic Particle Swarm Optimization (SPSO), and the hybrid Differential Evolution-Particle Swarm Optimization (DE-PSO), were applied to solve the inverse transient radiation problem in two-dimensional (2D) turbid media irradiated by the short pulse laser. The time-resolved radiative intensity signals simulated by finite volume method (FVM) were served as input for the inverse analysis. The sensitivities of the time-resolved radiation signals to the geometric parameters of the circular inclusions were also investigated. To illustrate the performance of these PSO algorithms, the optical properties, the size, and location of the circular inclusion were retrieved, respectively. The results showed that all these radiative parameters could be estimated accurately, even with noisy data. Compared with the PSO algorithm with inertia weights, the SPSO and DE-PSO algorithm were demonstrated to be more effective and robust, which had the potential to be implemented in 2D transient radiative transfer inverse problems.


---

## Body

## 1. Introduction

The problem of transient radiative transfer in participating media has attracted increasing interest over the last two decades due to the emergence of ultrashort pulse laser and its application to the picosecond level or higher accuracy level of time-resolved techniques [1]. Lasers with pulse durations of pico- to femtoseconds have been widely utilized in the areas of laser light-scattering flame diagnose, the prediction of the temperature distribution in combustion chambers, atmospheric science, remote sensing, nondestructive testing, laser machining, and medical diagnostics, to name a few [2–4]. In recent years, the medical applications of short pulse laser, such as the laser tissue welding [5, 6], the laser tissue ablation [7, 8], the optical tomography in medical imaging [9], and the photodynamic therapy [10], have drawn more and more attentions all over the world. In particular, the short pulse laser of near-infrared spectrum is extensively used in noninvasive techniques to retrieve the characteristics of the semitransparent media such as the biological tissues and turbid media [11–14]. Theoretically speaking, the research emphases of noninvasive reconstruction technology are focused on using the model-based iterative imaging reconstruction techniques, which employ a precise and efficient numerical forward model based on solving complete transient radiative transfer equation (TRTE) and then establish a proper inverse algorithm to retrieve the optical property parameters or internal geometry of the medium by using the measured time-resolved transmittance and reflectance signals. Compared with the traditional measurement methods, its biggest advantages are noninvasive, in situ measurement, containing a lot of time-resolved information about internal composition and characteristic of the participating media.For solving the inverse transient radiation problems in the noninvasive reconstruction, the first step is to solve the direct TRTE model accurately and efficiently. Nowadays, several numerical strategies have been developed to solve the forward model, including the discrete ordinate method (DOM) [15–18], the finite volume method (FVM) [18–21], the integral equation (IE) method [22], the finite-element method FEM [23, 24], and the Monte Carlo method [25]. Each of these methods has its own relative advantages and disadvantages and none of them is superior to others in all aspects. However, many numerical complications may arise from incompatibilities between the radiative transport model and the discrete formulation employed for the fluid dynamics and other heat transfer modes. The FVM for RTE avoids many of these complications and has the advantage of greater compatibility with existing finite volume based heat transfer software. Besides, the FVM also could simplify the problems of complex geometries using the unstructured meshes. Through decades of development, FVM has become a sophisticated technology in mechanics and thermal analysis. There have been several investigations focused on the interaction of ultrafast radiation with participating media using FVM [19, 20]. In these studies, both the homogeneous and inhomogeneous media have been considered. However, most of the studies of inverse radiation problems concentrate on one-dimensional (1D) cases and there are very few researches about multidimensional problems, especially the inverse transient radiative problems.To solve the inverse problems, a wide variety of inverse techniques have been successfully employed in the inverse radiation analyses which can be grouped into two categories, that is, the traditional algorithm based on gradient and the intelligent optimizations [26]. To date, many gradient-based techniques have been employed in the inverse radiation analysis such as the constrained least-squares method [27], the conjugate gradient (CG) method [28], and Levenberg-Marquardt method [29]. However, all these traditional methods depend on the initial value or the derivatives and gradients which are difficult to be solved accurately by numerical simulation in some cases. Furthermore, if the initial value is not chosen properly, the solution may be infeasible in the original domain, or in a worse case, not be convergent. In addition, these conventional methods using derivatives are generally restricted to the conditions of continuity, sensitivity, convexity, monotonicity, and nonlinearity of both objective functions and constraints [30]. Compared with traditional gradient-based methods, the intelligent optimization algorithms have a couple of outstanding characteristics; that is, the ill-posed inverse problem could also be solved, the derivative of the objective function is not necessary, and the a priori information is not needed. A typical characteristic of the intelligent optimization methods is that they can solve the global optimal problem reliably and obtain high quality global solutions with enough computational time, especially for higher problem dimensions [31]. This is primarily due to the fact that the random search technique is used to modify locally examined regions, which gives the algorithm a global optimization capability. The Particle Swarm Optimization (PSO), as an intelligent swarm optimization, first introduced in 1995 by Eberhart and Kennedy [32], has been studied extensively by many researchers in recent years. It is a potential heuristic bionic evolutionary algorithm originally inspired by the locking behavior of birds. Generally speaking, PSO is characterized to be simple in concept, easy to implement, and computationally efficient. Early studies have found the implementation of PSO to be effective and robust in solving problems featuring nonlinearly, nondifferentiability, and high dimensionality. Many modifications have been made to improve the convergence rate of the original PSO algorithm. More recently, our research group has done several studies on the PSO algorithm to solve the one-dimensional inverse radiation problems [33–37]. However, to the best of our knowledge, few researches are presented to solve the multidimensional inverse transient radiation problems by the intelligent optimization methods.The objective of present work is to apply the PSO-based algorithms to solve the inverse transient radiation problem in two-dimensional (2D) participating medium. The optical properties, the size, and location of the circular inclusion were retrieved, respectively. The remainder of the paper is organized as follows. The detailed mathematical formulation and computational steps of the direct model are described in Section2. The theoretical models of the standard PSO, the Stochastic PSO, and the hybrid Differential Evolution and Particle Swarm Optimization (DE-PSO) algorithms are described in Section 3. The solving model for direct problem and the inverse transient radiation analysis by the PSO-based algorithms are examined in Section 4. In Section 5, the main conclusions and perspectives are provided.


## 2. The Direct Model

In the present work, the thermal effect caused by the incident laser was ignored and only the transient radiative transfer was considered in the numerical model. Thus, the transient radiative transfer equation can be expressed as [38](1)1c∂Is,Ω,t∂t+∂Is,Ω,t∂s=-βIs,Ω,t+σs4π∫4πIs,Ω′,tΦΩ′,ΩdΩ′,where s is the length along the direction Ω, radiative intensity I is the function of position s, direction Ω, and time t, β, κa, and σs denotes extinction, absorption, and scattering coefficients, respectively, and β=κa+σs. ΦΩ′,Ω represents the scattering phase function between incoming direction Ω′ and scattering direction Ω. Ω′ is the solid angle in the direction Ω′. In the 2D transient radiative transfer problem caused by collimated irradiation, the radiative intensity within the medium is separated into two parts [35]: (a) Ic, the remnant of the collimated beam after partial extinction, by absorption and scattering, along its path; (b) Id, the diffuse part, which is the result of emission from the boundaries, emission from the medium inside, and the radiation scattered away from the collimated irradiation. Thus, the radiative intensity in the media can be written as(2)I=Ic+Id.For the square pulse laser, the intensity within the participating media can be expressed as [35](3)Ic=qinexp⁡-∫0s0κa+σsds·Hct-s0-Hct-ctp-s0δΩ-Ω0,where H is Heaviside function [35]; δ denotes Dirac function [35]; qin is the heat flux of incident laser at boundary x=0; Ω0 is the incident direction; and s0 represents the geometric distance along the direction Ω0.For the 2D participating media, the discrete equation could be obtained using the FVM model:(4)1c∂Ids,Ω,t∂t+∂Ids,Ω,t∂s=-βIds,Ω,t+Scs,Ω,t+Sds,Ω,t=-βIds,Ω,t+Sts,Ω,t,where(5)Sc=σs4π∫4πIcs,Ω,tΦΩ′,ΩdΩ′,Sd=σs4π∫4πIds,Ω,tΦΩ′,ΩdΩ′.In (4), St=Sc+Sd denotes the total source term. The source terms Sc and Sd result from the collimated radiation Ic and the diffused radiation Id, respectively. For the diffuse reflection boundary, the boundary condition is given as [35](6)Isw,Ω,t=εwIbsw,t+1-εwπ∫n^·Ω′>0Isw,Ω′,tn^·Ω′dΩ′hhhhhhhhhhhhhhhhhhhhhhhhn^·Ω<0,where Isw,Ω,t and Isw,Ω′,t are the outgoing and the incoming intensities at the boundary, respectively; n^ is the outward unit normal vector at the boundary; εw is the emissivity of the boundary. Substituting the length-dimensional time term t*=ct into (4), we can get(7)∂Ids,Ω,t*∂t*+∂Ids,Ω,t*∂s=-βIds,Ω,t*+Sts,Ω,t*.The total source termSt can be depicted as(8)Sts,Ω,t*=σs4π∫4πIcs,Ω,t*ΦΩ0,ΩdΩ′+σs4π∫4πIds,Ω,t*ΦΩ′,ΩdΩ′,where Ω0 denotes the direction of incident radiation. Using the fully implicit scheme, the diffusion intensity Id for the control volume P in the direction sm at time t* can be written as(9)Id,Pmt*=DxmId,xumt*fy+DymId,yumt*fx+SPmDxm+Dym+βΔxΔΩm/Bfxfy,SPm=ΔVΔΩmSt,pmfxfy+CΔVΔΩmBId,Pmt*-Δt*fxfy,St,Pm=Scm+Sdm=σs4πIcΦΩ0,Ω+σs4πId,PmΦΩm,ΩΔΩm,Dxm=∫ΔΩmn·s^mdΩ,where c represents the velocity of light; Δt* is the time step; Dm represents the direction weight; B and C are set as Δt*/(1+βΔt*) and 1/1+βΔt*, respectively. Id,xum and Id,yum denote the m-direction radiative intensities of node P at the upstream boundaries of control volume which are along axis x and axis y, respectively. The time-domain thermal signals, that is, the transmittance signal ρT and reflectance signal ρR, can be expressed as(10)ρTlx,y,t*=1qin2π∫μ>0Idlx,y,μ,t*dμ+qclx,y,Ω0,t*,ρR0,y,t*=-2πqin∫μ<0Id0,y,μ,t*dμ,where μ is the direction cosine and qc denotes the direct transmission flux of the collimated light. The FVM was chosen to solve the equation of TRT. For the sake of simplicity, the details of FVM are available in [20] and are not repeated here.


## 3. The Inverse Model

Commonly, there are three ways to obtain the internal information of media, that is, the continuous wave method, the time-domain method, and the frequency-domain method [39]. All these methods have the similar procedure in the experimental studies and the numerical simulations. Taking the time-domain method as an example, the reconstruction scheme consists of three major parts [40]: (1) a forward model that predicts the detector signals based on the solution of the transient radiative transfer equation; (2) an objective function that provides criterion of the differences between the detected and the predicted data; (3) the reconstruction technique which can minimize the objective function to get new guesses of the estimated parameters such as the optical properties or the geometrical parameters. Based on the new guesses of the optical properties, a new forward calculation is performed to get the corresponding detector predictions. The reconstruction process is completed when the value of the objective function is less than a preset value or the number of iterations exceeds the maximum number of iterations. Foremost, the forward model must be solved precisely enough so that the measurement data obtained by detectors could be simulated correctly. Consequently, the forward model based on the complete transport-theory TRTE should be utilized [40]. The details of the inverse optimization algorithms are shown as follows.


### 3.1. Theoretical Model of Basic PSO Algorithm

The PSO algorithm was first introduced by Eberhart and Kennedy [32] inspired by social behavior simulations of flocking birds. The flocking population is called a swarm, and the individuals are called the particles. The basic idea of PSO can be described simply as follows. Each particle in a swarm represents a potential solution which changes its position and updates its velocity based on its own and its neighbors’ experience. The velocity and position of the ith particle at generation t+1 could be expressed as [32](11)Vit+1=wVit+c1·r1·Pit-Xit+c2·r2·Pgt-Xit,(12)Xit+1=Xit+Vit+1,where c1 and c2 represent two positive constants called acceleration coefficients; r1 and r2 are uniformly distributed random numbers in the interval 0,1; w is the inertia weight factor which is used to control the impact of the previous velocities on the current velocity. In the present work, w is defined as w=1.0-0.6×t/N, where N represents the maximum generation number; Xit=Xi1t,Xi2t⋯XintT denotes the present location of particle i, which represents a potential solution; Vi(t)=[Vi1(t),Vi2(t)⋯Vin(t)]T is the present velocity of particle i, which is based on its own and its neighbors’ flying experience; Pi(t)=[Pi1(t),Pi2(t),…,Pin(t)]T is the “local best” in the tth generation of the swarm. The global best position of the swarm is expressed as Pg.


### 3.2. Theoretical Model of DE-PSO Algorithm

In order to avoid premature convergence of the standard PSO and ensure the overall convergence, a considerable number of modified PSO algorithms were proposed in recent years, like the Stochastic PSO (SPSO), the Multiphase PSO (MPPSO), the Quantum-Behaved PSO (QPSO), the hybrid Ant Colony Optimization and PSO (ACO-PSO), and so forth [33, 36, 41]. All these improved PSO-based algorithms are proved to be adaptive and robust parameter-searching techniques. In the present work, the SPSO and the hybrid Differential Evolution and PSO (DE-PSO) are utilized in the inverse procedure [42]. The DE-PSO is the combination of two stochastic algorithms and is expected to accelerate the convergence process by introducing the DEA which has strong global searching ability to the PSO. The detailed model description of SPSO is available in [41] and will not be repeated here. The theoretical model of DE-PSO is described as follows.The differential evolution algorithm (DEA) is a novel parallel, direct, and stochastic searching method which was proposed to solve the continuous domain problems in 1995 [43]. The basic idea of DEA is to use the differential of two individuals in the current generation, namely, vectors r2 and r3, to combine with a third vector r1 to generate a new particle, which is called the trial vector. Then, the objective function was calculated to decide whether the trial vector is to replace vector k or not. The trial vector T is generated according to [43](13)Tkt=Xr1t+dXr2t-Xr3t,where d is a real and constant factor which controls the amplification of the differential variation Xr2t-Xr3t, and in this paper it was set as 0.5. The integers r1, r2, and r3 are chosen randomly from the interval [1,n] and are different from running index i. Other parameters are defined as same as those of the PSO. Figure 1 illustrates the vector-generation process defined by (13).Figure 1
The schematic diagram of vector-generation process for DEA.In the process of DE-PSO, the differential evolution operator is introduced to the position updating equation, which can be expressed as [42](14)Xit+1=Xit+Vit+1+βεDE-Fmin⁡,where β=dXr2t-Xr3t, which is named as differential evolution operator; Fmin⁡ represents the minimum value of objective function in the present generation; εDE denotes a preset constant value which satisfies εDE≤Fmin⁡.The main function of the third term on the right side of (14) is to avoid local convergence. However, the differential evolution operation may lead to instability of the algorithm, so only when the objective function of the new particle is smaller than that of the last generation, (14) is applied as the position update equation. Otherwise, the position updates according to (12). The implementation of the approach for solving the inverse transient radiation problem by DE-PSO algorithm can be carried out according to the following routine.Step 1.
Initialize the positionXit and the speed Vit of each particle and calculate the objective function. Afterwards, the global best position Pg and every local best position Pit can be determined.Step 2.
Check if the value of the objective function is smaller than the preset parameterε or the maximum number of iterations is reached. If so, go to Step 6; otherwise, go to the next step.Step 3.
Update the velocity and position of each particle according to (11) and (14), respectively. Each velocity and position component is bounded by the maximum and minimum values to which each particle in the space is constrained. The range of the particle position is decided by the physical situation; for example, the single scattering albedo can only be set in the range of 0,1. Then, calculate the objective function of each particle and update the global beat position and local best positions.Step 4.
Calculate the objective function of each newly generated particle. If the value of the objective function is smaller than that of the corresponding particle in the last generation, then the new particle can be maintained. Otherwise, regenerate the particle according to (12).Step 5.
Update the generation fromt to t+1 and go to Step 2.Step 6.
Output the global best position and its corresponding value of objective function.To illustrate the difference between PSO, SPSO, and DE-PSO algorithms more clearly, the flowcharts of the three algorithms are shown in Figure2.Figure 2
The flowcharts of (a) PSO, (b) SPSO, and (c) DE-PSO algorithms.
(a)
(b)
(c)


### 3.3. Inverse Analysis Procedure

The transient radiative inverse problems for estimating the internal radiative properties of the 2D participating media are solved by minimizing the objective function which can be defined as(15)Fa=12∑j=14∫ts1ts2ρj,esta-ρj,meaa2dt,where j represents the different wall of the 2D media, ρj,meaa denotes the measured time-resolved transmittance or reflectance signals, which will be simulated by the forward model using the FVM, and ρj,esta is the estimated time-resolved transmittance and reflectance for an estimated vector a=a0,a1,…,aNT. Moreover, ts1 and ts2 are start and end point of the sampling span, respectively.To demonstrate the effects of measurement errors on the predicted terms, the random errors were considered. In the present paper, the simulated measured signals with random errors were obtained by adding normally distributed errors to exact reflectance and transmittance as follows:(16)ρmea=ρexa+rand·σ,where ρmea and ρexa represent the measured reflectance or transmittance signals with and without noisy data, respectively, and rand denotes a normal distributed random variable with zero mean and unit standard deviation. The standard deviation of the measured values is denoted by σ, which is defined as(17)σ=ρexa·γ%2.576,where γ% is the measured error and 2.576 arises from the fact that 99% of a normally distributed population is contained within ±2.567 standard deviation of the mean value. The absorption and scattering coefficients or the inclusion locations were estimated by minimizing the objective function using SPSO and DE-PSO algorithms. The inverse algorithms were adopted for minimization of the objective function, which is also called the fitness function. The stop criterion of the inverse algorithms is that the number of iterations exceeds the maximum number of iterations or the best fitness function is less than a specified small value.


## 3.1. Theoretical Model of Basic PSO Algorithm

The PSO algorithm was first introduced by Eberhart and Kennedy [32] inspired by social behavior simulations of flocking birds. The flocking population is called a swarm, and the individuals are called the particles. The basic idea of PSO can be described simply as follows. Each particle in a swarm represents a potential solution which changes its position and updates its velocity based on its own and its neighbors’ experience. The velocity and position of the ith particle at generation t+1 could be expressed as [32](11)Vit+1=wVit+c1·r1·Pit-Xit+c2·r2·Pgt-Xit,(12)Xit+1=Xit+Vit+1,where c1 and c2 represent two positive constants called acceleration coefficients; r1 and r2 are uniformly distributed random numbers in the interval 0,1; w is the inertia weight factor which is used to control the impact of the previous velocities on the current velocity. In the present work, w is defined as w=1.0-0.6×t/N, where N represents the maximum generation number; Xit=Xi1t,Xi2t⋯XintT denotes the present location of particle i, which represents a potential solution; Vi(t)=[Vi1(t),Vi2(t)⋯Vin(t)]T is the present velocity of particle i, which is based on its own and its neighbors’ flying experience; Pi(t)=[Pi1(t),Pi2(t),…,Pin(t)]T is the “local best” in the tth generation of the swarm. The global best position of the swarm is expressed as Pg.


## 3.2. Theoretical Model of DE-PSO Algorithm

In order to avoid premature convergence of the standard PSO and ensure the overall convergence, a considerable number of modified PSO algorithms were proposed in recent years, like the Stochastic PSO (SPSO), the Multiphase PSO (MPPSO), the Quantum-Behaved PSO (QPSO), the hybrid Ant Colony Optimization and PSO (ACO-PSO), and so forth [33, 36, 41]. All these improved PSO-based algorithms are proved to be adaptive and robust parameter-searching techniques. In the present work, the SPSO and the hybrid Differential Evolution and PSO (DE-PSO) are utilized in the inverse procedure [42]. The DE-PSO is the combination of two stochastic algorithms and is expected to accelerate the convergence process by introducing the DEA which has strong global searching ability to the PSO. The detailed model description of SPSO is available in [41] and will not be repeated here. The theoretical model of DE-PSO is described as follows.The differential evolution algorithm (DEA) is a novel parallel, direct, and stochastic searching method which was proposed to solve the continuous domain problems in 1995 [43]. The basic idea of DEA is to use the differential of two individuals in the current generation, namely, vectors r2 and r3, to combine with a third vector r1 to generate a new particle, which is called the trial vector. Then, the objective function was calculated to decide whether the trial vector is to replace vector k or not. The trial vector T is generated according to [43](13)Tkt=Xr1t+dXr2t-Xr3t,where d is a real and constant factor which controls the amplification of the differential variation Xr2t-Xr3t, and in this paper it was set as 0.5. The integers r1, r2, and r3 are chosen randomly from the interval [1,n] and are different from running index i. Other parameters are defined as same as those of the PSO. Figure 1 illustrates the vector-generation process defined by (13).Figure 1
The schematic diagram of vector-generation process for DEA.In the process of DE-PSO, the differential evolution operator is introduced to the position updating equation, which can be expressed as [42](14)Xit+1=Xit+Vit+1+βεDE-Fmin⁡,where β=dXr2t-Xr3t, which is named as differential evolution operator; Fmin⁡ represents the minimum value of objective function in the present generation; εDE denotes a preset constant value which satisfies εDE≤Fmin⁡.The main function of the third term on the right side of (14) is to avoid local convergence. However, the differential evolution operation may lead to instability of the algorithm, so only when the objective function of the new particle is smaller than that of the last generation, (14) is applied as the position update equation. Otherwise, the position updates according to (12). The implementation of the approach for solving the inverse transient radiation problem by DE-PSO algorithm can be carried out according to the following routine.Step 1.
Initialize the positionXit and the speed Vit of each particle and calculate the objective function. Afterwards, the global best position Pg and every local best position Pit can be determined.Step 2.
Check if the value of the objective function is smaller than the preset parameterε or the maximum number of iterations is reached. If so, go to Step 6; otherwise, go to the next step.Step 3.
Update the velocity and position of each particle according to (11) and (14), respectively. Each velocity and position component is bounded by the maximum and minimum values to which each particle in the space is constrained. The range of the particle position is decided by the physical situation; for example, the single scattering albedo can only be set in the range of 0,1. Then, calculate the objective function of each particle and update the global beat position and local best positions.Step 4.
Calculate the objective function of each newly generated particle. If the value of the objective function is smaller than that of the corresponding particle in the last generation, then the new particle can be maintained. Otherwise, regenerate the particle according to (12).Step 5.
Update the generation fromt to t+1 and go to Step 2.Step 6.
Output the global best position and its corresponding value of objective function.To illustrate the difference between PSO, SPSO, and DE-PSO algorithms more clearly, the flowcharts of the three algorithms are shown in Figure2.Figure 2
The flowcharts of (a) PSO, (b) SPSO, and (c) DE-PSO algorithms.
(a)
(b)
(c)


## 3.3. Inverse Analysis Procedure

The transient radiative inverse problems for estimating the internal radiative properties of the 2D participating media are solved by minimizing the objective function which can be defined as(15)Fa=12∑j=14∫ts1ts2ρj,esta-ρj,meaa2dt,where j represents the different wall of the 2D media, ρj,meaa denotes the measured time-resolved transmittance or reflectance signals, which will be simulated by the forward model using the FVM, and ρj,esta is the estimated time-resolved transmittance and reflectance for an estimated vector a=a0,a1,…,aNT. Moreover, ts1 and ts2 are start and end point of the sampling span, respectively.To demonstrate the effects of measurement errors on the predicted terms, the random errors were considered. In the present paper, the simulated measured signals with random errors were obtained by adding normally distributed errors to exact reflectance and transmittance as follows:(16)ρmea=ρexa+rand·σ,where ρmea and ρexa represent the measured reflectance or transmittance signals with and without noisy data, respectively, and rand denotes a normal distributed random variable with zero mean and unit standard deviation. The standard deviation of the measured values is denoted by σ, which is defined as(17)σ=ρexa·γ%2.576,where γ% is the measured error and 2.576 arises from the fact that 99% of a normally distributed population is contained within ±2.567 standard deviation of the mean value. The absorption and scattering coefficients or the inclusion locations were estimated by minimizing the objective function using SPSO and DE-PSO algorithms. The inverse algorithms were adopted for minimization of the objective function, which is also called the fitness function. The stop criterion of the inverse algorithms is that the number of iterations exceeds the maximum number of iterations or the best fitness function is less than a specified small value.


## 4. Results and Discussion

### 4.1. Homogeneous Media

Consider a 2D homogeneous semitransparent medium with length 1.0 m × 1.0 m whose left side exposes to a collimated square pulse laser beam in thex-direction with pulse width ctp=0.1m. The absorption coefficient κa and the scattering coefficient σs are set as 0.02 m−1 and 9.98 m−1, respectively. The emissions of the medium and its boundaries are neglected. The media is assumed to be anisotropic scattering with black boundaries and the phase function is defined as(18)ΦΩ′,Ω=1+a0μμ′+ηη′+ξξ′,where a0=1, a0=0, and a0=-1 represent the backward scattering, isotropic scattering, and forward scattering, respectively. μ, η, and ξ and μ′, η′, and ξ′ denote the direction cosines of the incoming direction Ω and scattering direction Ω′, respectively. The control volumes were set as Nx×Ny=500×500 and solid angles were divided into Nθ×Nϕ=20×10 in the FVM. The total observed time span ct was 10 m. As shown in Figure 3, the results of FVM model are compared with the solutions of Discrete Ordinate Method (DOM) [16] and Least Square Finite Element Method (LSFEM) [44]. It can be seen that the results calculated by FVM approximation in the present study are in good agreement with those in [16, 44] which means the FVM model in this paper is proved to be valid.Figure 3
Validation of the FVM solution for transient radiative transfer of a 2D homogeneous rectangle medium. (a) Transmittance and (b) reflectance.
(a)
(b)


### 4.2. Nonhomogeneous Media with Circular Inclusions

The estimation in this paper is expected to be applied in the nondestructive detection of human or small animals. Figure4 shows the model of a rat liver with tumor which is simplified to a rectangle medium with a circular inclusion. To illustrate the influence of the inclusions to the radiative signals, four cases are investigated in this section that is, two circular inclusions with different radiuses, location and optical properties in a rectangular medium. The media size is 1.0m×1.0m which exposed to a pulse laser with pulse width ctp=0.1m. The absorption coefficient and scattering coefficient of the 2D medium are set as κa0=2.0m-1 and σs0=3.0m-1. The geometric locations of the two inclusions are identified by three parameters, namely, z1, z2, and r, among which z1 and z2 represent the coordinate values of the inclusion center and r denotes the radius of the circular inclusions. Table 1 shows the control parameters of the four cases. It is worthy noticing that the four cases are designed on the principle that, compared to Case  1, the other three cases only change one set of parameters, that is, the size, location, or optical properties, for both inclusions. For the case of clarity, the diagrams of the four cases are shown in Figure 5. The control volumes are set as Nx×Ny=100×100 and the time step is taken as Δct=0.01m. The time-resolved reflectance and transmittance signals are shown in Figures 6(a) and 6(b). It is apparent that the size, location, and properties of the inclusions have a significant influence on the radiative signals which can be concluded that these parameters can be retrieved theoretically. All the cases were implemented using FORTRAN code and the developed program was executed on an Intel Xeon E5-2670 PC.Table 1
The control parameters of the four cases.
Parameters
Case  1
Case  2
Case  3
Case  4
Size [m]
r
1
0.1
0.2
0.1
0.1
r
2
0.2
0.1
0.2
0.2
Location [m]
(
z
1,1
,
z
2,1
)
(0.3, 0.45)
(0.3, 0.45)
(0.7, 0.65)
(0.3, 0.45)
(
z
1,2
,
z
2,2
)
(0.75, 0.65)
(0.75, 0.65)
(0.35, 0.25)
(0.75, 0.65)
Optical properties [m−1]
(
κ
a
1
,
σ
s
1
)
(1.0, 7.5)
(1.0, 7.5)
(1.0, 7.5)
(4.0, 5.5)
(
κ
a
2
,
σ
s
2
)
(0.2, 9.8)
(0.2, 9.8)
(0.2, 9.8)
(1.5, 9.0)Figure 4
The simplification of actual models.Figure 5
The physical model contains two inclusions.
(a)
(b)
(c)
(d)Figure 6
The (a) reflectance and (b) transmittance signals in the function ofct.
(a)
(b)


### 4.3. Reconstruction of the Inclusion Parameters

#### 4.3.1. Estimation of the Optical Parameters

To test the reconstruction ability of the three swarm intelligence optimization algorithms, two cases for estimating the optical parameters were considered. For the laser incidence, two different conditions were included (see Figure7): (1) point-irradiated (PI) condition which means one single laser beam irradiates one spot of the media wall; (2) face-irradiated (FI) condition which means the laser beams incident on the whole boundary wall. In this section, the two cases were proceeded under both two laser incident conditions. In Case  1, the inclusion was in the center of the medium which means z1=z2=0.5 m and r=0.25 m. The responses of detectors D1, D2, D3, and D4, which are shown in Figure 8(a), were simulated by the forward model. The properties of the inclusion were estimated and the inverse analysis was taken by using the PSO, SPSO, and DE-PSO algorithms. Meanwhile, the optical properties of the 2D medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true values of the inclusion properties were κa1=1.0m-1 and σs1=9.0m-1, respectively. The searching intervals were set as [0,10 m−1]. The estimated results are listed in Table 2. In Case  2, two inclusions inside the medium were considered as shown in Figure 8(b). The optical properties of both inclusions were estimated by PSO, SPSO, and DE-PSO. In this case, the optical properties of the medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true value of the inclusions which needed to be estimated was κa1,σs1=1.0m-1,9.0m-1 and κa2,σs2=4.5m-1,5.5m-1, respectively. The size and location of the two inclusions were z1,1,z2,1,r1=0.3,0.3,0.2 m and z1,2,z2,2,r2=0.7,0.6,0.15 m. The retrieving results of optical parameters are shown in Tables 3 and 4.Table 2
The retrieving results of the optical properties for Case  1 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Condition
True values or inversion results [m−1]
Relative errors [%]
κ
a
1 = 1.0
σ
s
1 = 9.0
κ
a
1
σ
s
1
PSO
PI
1.0185
9.1372
1.846
1.524
FI
0.9964
8.9671
0.356
0.366
SPSO
PI
0.9880
8.8960
1.200
1.155
FI
0.9979
9.0151
0.206
0.167
DE-PSO
PI
1.0098
8.9138
0.983
0.958
FI
1.0008
8.9956
0.077
0.049Table 3
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO of different generations.
Generations
True value [m−1]
Algorithm
Retrieving results [m−1]
Relative errors [%]
10
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9750, 9.3216)
(2.4962, 3.5738)
(7.4025, 5.9205)
(64.5009, 7.6456)
SPSO
(0.9803, 8.9818)
(1.9718, 0.2021)
(5.4861, 0.3347)
(21.9131, 93.9149)
DE-PSO
(0.9844, 8.8901)
(1.5635, 1.2207)
(3.3674, 3.7704)
(25.1700, 31.4477)
50
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.4496, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9995, 8.9996)
(0.0478, 0.0040)
(3.5513, 4.6093)
(21.0814, 16.1946)
DE-PSO
(1.0021, 9.0036)
(0.2071, 0.0397)
(3.3049, 4.0275)
(26.5584, 26.7718)
100
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.449, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9998, 8.9995)
(0.0234, 0.0054)
(3.5119, 4.8832)
(21.9588, 11.2147)
DE-PSO
(1.0007, 9.0042)
(0.0696, 0.0472)
(3.6868, 4.8677)
(18.0717, 11.4957)Table 4
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Time [s]
True values [m−1]
Retrieving results [m−1]
Relative errors [%]
PSO
25476
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(0.9985, 9.0019)
(0.1486, 0.0207)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(6.9803, 7.5019)
(55.118, 36.398)
SPSO
32775
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0000, 8.9996)
(0.0034, 0.0046)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.4444, 5.4556)
(1.2349, 0.8077)
DE-PSO
49554
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0005, 9.0036)
(0.0544, 0.0397)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.5248, 5.4454)
(0.5508, 0.9925)Figure 7
The incident types of pulse laser.Figure 8
The schematic of pulse laser irradiation and the measurement position of time-domain signals of the two test cases.
(a)
(b)It can be seen from Table2 that κa1 and σs1 can be retrieved accurately under both PI and FI conditions by PSO, SPSO, and DE-PSO when there is only one inclusion. Furthermore, it also demonstrates that the retrieving results under FI condition are more accurate than those under PI condition, the reason of which is that, under FI condition, the reflectance and transmittance signals are stronger and are carrying more information available. Meanwhile, the DE-PSO algorithm shows its priority in reconstructing the inclusion properties. Under both conditions, the relative errors of retrieval results are smaller than those of the SPSO and PSO. Table 3 demonstrates the retrieving results of the optical properties for Case  2 of different generations. As shown in Table 4, the results of SPSO and DE-PSO are more accurate than those of the PSO. Furthermore, under FI condition, for Case  2, the optical properties of both inclusions can be estimated accurately without measurement errors and with the fact that the estimated relative errors of results for inclusion B (see Figure 8) are larger than those of inclusion A. However, the calculation time required of the DE-PSO is more than that of the SPSO and PSO. Thus, it can be concluded that SPSO and DE-PSO can retrieved the optical properties in a 2D medium accurately but with longer calculation time. The value of the objective function of PSO, DE-PSO, and SPSO with two inclusions is shown in Figure 9. It can be seen that the PSO fall into local convergence quickly and stop searching. In the earlier stage of the inverse process the objective function of SPSO declines faster than that of DE-PSO. However, in the later period of the inverse process, the SPSO almost stops searching and the objective function of DE-PSO continues declining which leads to the result that the retrieving results of DE-PSO are more accurate than those of the SPSO.Figure 9
The objective function of PSO, DE-PSO, and SPSO in Case  2.


#### 4.3.2. Estimation of Geometrical Parameters

In the inverse procedure, the long sampling span will reduce the computational efficiency. Hence, it is of great significance to select the proper sampling span. Thus, the sensitivity of radiative signal to the inverse parameters needs to be analyzed. The sensitivity coefficient is one of the most important characteristic parameters in the sensitivity analysis, which is the first derivative of the radiative signals to a certain inverse parameter. The sensitivity coefficient is defined as(19)Smiρj=∂ρj∂mimi=m0=ρjm0+m0Δ-ρjm0-m0Δ2m0Δ,where mi denotes the independent variable which stands for z1, z2, and r; Δ represents a tiny change; ρj is the radiative signals of each boundary.The sensitivities of the inclusion location (z1,z2) and size r for the signals D1–D4 are shown in Figure 10. The parameters are set the same as Case  1 described in Section 4.3.1. It demonstrates that the sensitivity coefficient is relatively larger in the span of ct∈0,3m. So this span can be chosen as the sampling span for inverse analysis. The same results are expected in other cases.Figure 10
The sensitivity of the radiation signals for the geometric parameters of the inclusion.
(a)
(b)
(c)In many conditions, the shapes and optical properties of inclusions in the media are known in the actual application cases of reconstructing the medium’s inside information. For example, a tumor or cyst in the biological tissue is spherical or spherical-like, and the defects in the specimen, such as the air pore, are mainly the shape of a cube, sphere, or ellipsoid. In addition, the optical properties of human and animals tissues, in vitro or in vivo, are studied thoroughly [45]. Consequently, the geometrical feature and optical properties of inclusions could be used as a priori information in these cases, which will reduce the complexity of the reconstruction technique and improve the quality and efficiency of the reconstruction results. Therefore, on the condition that the inclusion shape and the optical properties are known, the inverse problem turns into the retrieval of the size and location of the inclusions.In this section, the inclusion location (z1, z2) and size r of cases (z1 = z2 = 0.25 m, r=0.25 m), (z1 = z2 = 0.5 m, r=0.25 m), and (z1 = z2 = 0.75 m, r=0.25 m) were estimated simultaneously. The size of the 2D medium was set as 1.0m×1.0m. The absorption coefficient κa0 and scattering coefficient σs0 of the medium were 0.2 m−1 and 0.3 m−1, respectively. There was one circular inclusion in the 2D medium with absorption coefficient κa1 of 1.0 m−1 and scattering coefficient σs1 of 9.0 m−1, but the inclusion size and location were unknown, which needed to be retrieved. The FI incidence type and the data of the detectors 1, 2, 3, and 4 were used for inverse analysis by using the PSO, SPSO, and DE-PSO algorithms with measurement errors. The search span of z1, z2, and r was set as 0,1 m. The results are shown in Tables 5 and 6. Figure 11 shows the comparison of the exact and reconstructed reflectance signal when z1=z2=r=0.25 m with 10% measurement error.Table 5
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO without measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2446
0.2457
0.2479
2.163
1.712
0.837
SPSO
0.2499
0.2500
0.2500
0.027
0.007
0.017
DE-PSO
0.2500
0.2500
0.2500
0.015
0.005
0.007
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5127
0.5045
0.2584
2.552
0.893
3.369
SPSO
0.500
0.500
0.2499
0.042
0.002
0.033
DE-PSO
0.4999
0.500
0.2499
0.019
0.006
0.021
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.7371
0.7694
0.2281
1.718
2.591
8.779
SPSO
0.7497
0.7500
0.2500
0.041
0.005
0.105
DE-PSO
0.7498
0.7500
0.2498
0.026
0.005
0.067Table 6
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO with 10% measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2804
0.2669
0.2676
12.18
6.747
7.051
SPSO
0.2509
0.2506
0.2492
0.356
0.243
0.339
DE-PSO
0.2340
0.2412
0.2452
6.395
3.514
1.934
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5993
0.3354
0.1238
19.88
32.91
50.49
SPSO
0.5006
0.5267
0.2453
0.120
5.331
1.894
DE-PSO
0.5150
0.5302
0.2386
2.996
6.045
4.549
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.8371
0.5694
0.2981
11.61
24.08
19.22
SPSO
0.7356
0.7681
0.2265
1.926
2.419
9.418
DE-PSO
0.7237
0.7786
0.2324
3.508
3.821
7.041Figure 11
The comparison of the exact and reconstructed reflectance signals.It can be seen that the size and location can be estimated accurately using PSO, SPSO, and DE-PSO without measurement errors under FI condition. Furthermore, the results of DE-PSO are more accurate than those of the SPSO and PSO. The relative errors are reasonable even with 10% measurement error. However, it is worth noting that the relative errors of DE-PSO with 10% measurement error are bigger than those of the SPSO, which means the SPSO algorithm is more robust than the DE-PSO.


## 4.1. Homogeneous Media

Consider a 2D homogeneous semitransparent medium with length 1.0 m × 1.0 m whose left side exposes to a collimated square pulse laser beam in thex-direction with pulse width ctp=0.1m. The absorption coefficient κa and the scattering coefficient σs are set as 0.02 m−1 and 9.98 m−1, respectively. The emissions of the medium and its boundaries are neglected. The media is assumed to be anisotropic scattering with black boundaries and the phase function is defined as(18)ΦΩ′,Ω=1+a0μμ′+ηη′+ξξ′,where a0=1, a0=0, and a0=-1 represent the backward scattering, isotropic scattering, and forward scattering, respectively. μ, η, and ξ and μ′, η′, and ξ′ denote the direction cosines of the incoming direction Ω and scattering direction Ω′, respectively. The control volumes were set as Nx×Ny=500×500 and solid angles were divided into Nθ×Nϕ=20×10 in the FVM. The total observed time span ct was 10 m. As shown in Figure 3, the results of FVM model are compared with the solutions of Discrete Ordinate Method (DOM) [16] and Least Square Finite Element Method (LSFEM) [44]. It can be seen that the results calculated by FVM approximation in the present study are in good agreement with those in [16, 44] which means the FVM model in this paper is proved to be valid.Figure 3
Validation of the FVM solution for transient radiative transfer of a 2D homogeneous rectangle medium. (a) Transmittance and (b) reflectance.
(a)
(b)


## 4.2. Nonhomogeneous Media with Circular Inclusions

The estimation in this paper is expected to be applied in the nondestructive detection of human or small animals. Figure4 shows the model of a rat liver with tumor which is simplified to a rectangle medium with a circular inclusion. To illustrate the influence of the inclusions to the radiative signals, four cases are investigated in this section that is, two circular inclusions with different radiuses, location and optical properties in a rectangular medium. The media size is 1.0m×1.0m which exposed to a pulse laser with pulse width ctp=0.1m. The absorption coefficient and scattering coefficient of the 2D medium are set as κa0=2.0m-1 and σs0=3.0m-1. The geometric locations of the two inclusions are identified by three parameters, namely, z1, z2, and r, among which z1 and z2 represent the coordinate values of the inclusion center and r denotes the radius of the circular inclusions. Table 1 shows the control parameters of the four cases. It is worthy noticing that the four cases are designed on the principle that, compared to Case  1, the other three cases only change one set of parameters, that is, the size, location, or optical properties, for both inclusions. For the case of clarity, the diagrams of the four cases are shown in Figure 5. The control volumes are set as Nx×Ny=100×100 and the time step is taken as Δct=0.01m. The time-resolved reflectance and transmittance signals are shown in Figures 6(a) and 6(b). It is apparent that the size, location, and properties of the inclusions have a significant influence on the radiative signals which can be concluded that these parameters can be retrieved theoretically. All the cases were implemented using FORTRAN code and the developed program was executed on an Intel Xeon E5-2670 PC.Table 1
The control parameters of the four cases.
Parameters
Case  1
Case  2
Case  3
Case  4
Size [m]
r
1
0.1
0.2
0.1
0.1
r
2
0.2
0.1
0.2
0.2
Location [m]
(
z
1,1
,
z
2,1
)
(0.3, 0.45)
(0.3, 0.45)
(0.7, 0.65)
(0.3, 0.45)
(
z
1,2
,
z
2,2
)
(0.75, 0.65)
(0.75, 0.65)
(0.35, 0.25)
(0.75, 0.65)
Optical properties [m−1]
(
κ
a
1
,
σ
s
1
)
(1.0, 7.5)
(1.0, 7.5)
(1.0, 7.5)
(4.0, 5.5)
(
κ
a
2
,
σ
s
2
)
(0.2, 9.8)
(0.2, 9.8)
(0.2, 9.8)
(1.5, 9.0)Figure 4
The simplification of actual models.Figure 5
The physical model contains two inclusions.
(a)
(b)
(c)
(d)Figure 6
The (a) reflectance and (b) transmittance signals in the function ofct.
(a)
(b)


## 4.3. Reconstruction of the Inclusion Parameters

### 4.3.1. Estimation of the Optical Parameters

To test the reconstruction ability of the three swarm intelligence optimization algorithms, two cases for estimating the optical parameters were considered. For the laser incidence, two different conditions were included (see Figure7): (1) point-irradiated (PI) condition which means one single laser beam irradiates one spot of the media wall; (2) face-irradiated (FI) condition which means the laser beams incident on the whole boundary wall. In this section, the two cases were proceeded under both two laser incident conditions. In Case  1, the inclusion was in the center of the medium which means z1=z2=0.5 m and r=0.25 m. The responses of detectors D1, D2, D3, and D4, which are shown in Figure 8(a), were simulated by the forward model. The properties of the inclusion were estimated and the inverse analysis was taken by using the PSO, SPSO, and DE-PSO algorithms. Meanwhile, the optical properties of the 2D medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true values of the inclusion properties were κa1=1.0m-1 and σs1=9.0m-1, respectively. The searching intervals were set as [0,10 m−1]. The estimated results are listed in Table 2. In Case  2, two inclusions inside the medium were considered as shown in Figure 8(b). The optical properties of both inclusions were estimated by PSO, SPSO, and DE-PSO. In this case, the optical properties of the medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true value of the inclusions which needed to be estimated was κa1,σs1=1.0m-1,9.0m-1 and κa2,σs2=4.5m-1,5.5m-1, respectively. The size and location of the two inclusions were z1,1,z2,1,r1=0.3,0.3,0.2 m and z1,2,z2,2,r2=0.7,0.6,0.15 m. The retrieving results of optical parameters are shown in Tables 3 and 4.Table 2
The retrieving results of the optical properties for Case  1 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Condition
True values or inversion results [m−1]
Relative errors [%]
κ
a
1 = 1.0
σ
s
1 = 9.0
κ
a
1
σ
s
1
PSO
PI
1.0185
9.1372
1.846
1.524
FI
0.9964
8.9671
0.356
0.366
SPSO
PI
0.9880
8.8960
1.200
1.155
FI
0.9979
9.0151
0.206
0.167
DE-PSO
PI
1.0098
8.9138
0.983
0.958
FI
1.0008
8.9956
0.077
0.049Table 3
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO of different generations.
Generations
True value [m−1]
Algorithm
Retrieving results [m−1]
Relative errors [%]
10
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9750, 9.3216)
(2.4962, 3.5738)
(7.4025, 5.9205)
(64.5009, 7.6456)
SPSO
(0.9803, 8.9818)
(1.9718, 0.2021)
(5.4861, 0.3347)
(21.9131, 93.9149)
DE-PSO
(0.9844, 8.8901)
(1.5635, 1.2207)
(3.3674, 3.7704)
(25.1700, 31.4477)
50
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.4496, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9995, 8.9996)
(0.0478, 0.0040)
(3.5513, 4.6093)
(21.0814, 16.1946)
DE-PSO
(1.0021, 9.0036)
(0.2071, 0.0397)
(3.3049, 4.0275)
(26.5584, 26.7718)
100
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.449, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9998, 8.9995)
(0.0234, 0.0054)
(3.5119, 4.8832)
(21.9588, 11.2147)
DE-PSO
(1.0007, 9.0042)
(0.0696, 0.0472)
(3.6868, 4.8677)
(18.0717, 11.4957)Table 4
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Time [s]
True values [m−1]
Retrieving results [m−1]
Relative errors [%]
PSO
25476
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(0.9985, 9.0019)
(0.1486, 0.0207)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(6.9803, 7.5019)
(55.118, 36.398)
SPSO
32775
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0000, 8.9996)
(0.0034, 0.0046)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.4444, 5.4556)
(1.2349, 0.8077)
DE-PSO
49554
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0005, 9.0036)
(0.0544, 0.0397)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.5248, 5.4454)
(0.5508, 0.9925)Figure 7
The incident types of pulse laser.Figure 8
The schematic of pulse laser irradiation and the measurement position of time-domain signals of the two test cases.
(a)
(b)It can be seen from Table2 that κa1 and σs1 can be retrieved accurately under both PI and FI conditions by PSO, SPSO, and DE-PSO when there is only one inclusion. Furthermore, it also demonstrates that the retrieving results under FI condition are more accurate than those under PI condition, the reason of which is that, under FI condition, the reflectance and transmittance signals are stronger and are carrying more information available. Meanwhile, the DE-PSO algorithm shows its priority in reconstructing the inclusion properties. Under both conditions, the relative errors of retrieval results are smaller than those of the SPSO and PSO. Table 3 demonstrates the retrieving results of the optical properties for Case  2 of different generations. As shown in Table 4, the results of SPSO and DE-PSO are more accurate than those of the PSO. Furthermore, under FI condition, for Case  2, the optical properties of both inclusions can be estimated accurately without measurement errors and with the fact that the estimated relative errors of results for inclusion B (see Figure 8) are larger than those of inclusion A. However, the calculation time required of the DE-PSO is more than that of the SPSO and PSO. Thus, it can be concluded that SPSO and DE-PSO can retrieved the optical properties in a 2D medium accurately but with longer calculation time. The value of the objective function of PSO, DE-PSO, and SPSO with two inclusions is shown in Figure 9. It can be seen that the PSO fall into local convergence quickly and stop searching. In the earlier stage of the inverse process the objective function of SPSO declines faster than that of DE-PSO. However, in the later period of the inverse process, the SPSO almost stops searching and the objective function of DE-PSO continues declining which leads to the result that the retrieving results of DE-PSO are more accurate than those of the SPSO.Figure 9
The objective function of PSO, DE-PSO, and SPSO in Case  2.


### 4.3.2. Estimation of Geometrical Parameters

In the inverse procedure, the long sampling span will reduce the computational efficiency. Hence, it is of great significance to select the proper sampling span. Thus, the sensitivity of radiative signal to the inverse parameters needs to be analyzed. The sensitivity coefficient is one of the most important characteristic parameters in the sensitivity analysis, which is the first derivative of the radiative signals to a certain inverse parameter. The sensitivity coefficient is defined as(19)Smiρj=∂ρj∂mimi=m0=ρjm0+m0Δ-ρjm0-m0Δ2m0Δ,where mi denotes the independent variable which stands for z1, z2, and r; Δ represents a tiny change; ρj is the radiative signals of each boundary.The sensitivities of the inclusion location (z1,z2) and size r for the signals D1–D4 are shown in Figure 10. The parameters are set the same as Case  1 described in Section 4.3.1. It demonstrates that the sensitivity coefficient is relatively larger in the span of ct∈0,3m. So this span can be chosen as the sampling span for inverse analysis. The same results are expected in other cases.Figure 10
The sensitivity of the radiation signals for the geometric parameters of the inclusion.
(a)
(b)
(c)In many conditions, the shapes and optical properties of inclusions in the media are known in the actual application cases of reconstructing the medium’s inside information. For example, a tumor or cyst in the biological tissue is spherical or spherical-like, and the defects in the specimen, such as the air pore, are mainly the shape of a cube, sphere, or ellipsoid. In addition, the optical properties of human and animals tissues, in vitro or in vivo, are studied thoroughly [45]. Consequently, the geometrical feature and optical properties of inclusions could be used as a priori information in these cases, which will reduce the complexity of the reconstruction technique and improve the quality and efficiency of the reconstruction results. Therefore, on the condition that the inclusion shape and the optical properties are known, the inverse problem turns into the retrieval of the size and location of the inclusions.In this section, the inclusion location (z1, z2) and size r of cases (z1 = z2 = 0.25 m, r=0.25 m), (z1 = z2 = 0.5 m, r=0.25 m), and (z1 = z2 = 0.75 m, r=0.25 m) were estimated simultaneously. The size of the 2D medium was set as 1.0m×1.0m. The absorption coefficient κa0 and scattering coefficient σs0 of the medium were 0.2 m−1 and 0.3 m−1, respectively. There was one circular inclusion in the 2D medium with absorption coefficient κa1 of 1.0 m−1 and scattering coefficient σs1 of 9.0 m−1, but the inclusion size and location were unknown, which needed to be retrieved. The FI incidence type and the data of the detectors 1, 2, 3, and 4 were used for inverse analysis by using the PSO, SPSO, and DE-PSO algorithms with measurement errors. The search span of z1, z2, and r was set as 0,1 m. The results are shown in Tables 5 and 6. Figure 11 shows the comparison of the exact and reconstructed reflectance signal when z1=z2=r=0.25 m with 10% measurement error.Table 5
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO without measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2446
0.2457
0.2479
2.163
1.712
0.837
SPSO
0.2499
0.2500
0.2500
0.027
0.007
0.017
DE-PSO
0.2500
0.2500
0.2500
0.015
0.005
0.007
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5127
0.5045
0.2584
2.552
0.893
3.369
SPSO
0.500
0.500
0.2499
0.042
0.002
0.033
DE-PSO
0.4999
0.500
0.2499
0.019
0.006
0.021
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.7371
0.7694
0.2281
1.718
2.591
8.779
SPSO
0.7497
0.7500
0.2500
0.041
0.005
0.105
DE-PSO
0.7498
0.7500
0.2498
0.026
0.005
0.067Table 6
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO with 10% measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2804
0.2669
0.2676
12.18
6.747
7.051
SPSO
0.2509
0.2506
0.2492
0.356
0.243
0.339
DE-PSO
0.2340
0.2412
0.2452
6.395
3.514
1.934
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5993
0.3354
0.1238
19.88
32.91
50.49
SPSO
0.5006
0.5267
0.2453
0.120
5.331
1.894
DE-PSO
0.5150
0.5302
0.2386
2.996
6.045
4.549
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.8371
0.5694
0.2981
11.61
24.08
19.22
SPSO
0.7356
0.7681
0.2265
1.926
2.419
9.418
DE-PSO
0.7237
0.7786
0.2324
3.508
3.821
7.041Figure 11
The comparison of the exact and reconstructed reflectance signals.It can be seen that the size and location can be estimated accurately using PSO, SPSO, and DE-PSO without measurement errors under FI condition. Furthermore, the results of DE-PSO are more accurate than those of the SPSO and PSO. The relative errors are reasonable even with 10% measurement error. However, it is worth noting that the relative errors of DE-PSO with 10% measurement error are bigger than those of the SPSO, which means the SPSO algorithm is more robust than the DE-PSO.


## 4.3.1. Estimation of the Optical Parameters

To test the reconstruction ability of the three swarm intelligence optimization algorithms, two cases for estimating the optical parameters were considered. For the laser incidence, two different conditions were included (see Figure7): (1) point-irradiated (PI) condition which means one single laser beam irradiates one spot of the media wall; (2) face-irradiated (FI) condition which means the laser beams incident on the whole boundary wall. In this section, the two cases were proceeded under both two laser incident conditions. In Case  1, the inclusion was in the center of the medium which means z1=z2=0.5 m and r=0.25 m. The responses of detectors D1, D2, D3, and D4, which are shown in Figure 8(a), were simulated by the forward model. The properties of the inclusion were estimated and the inverse analysis was taken by using the PSO, SPSO, and DE-PSO algorithms. Meanwhile, the optical properties of the 2D medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true values of the inclusion properties were κa1=1.0m-1 and σs1=9.0m-1, respectively. The searching intervals were set as [0,10 m−1]. The estimated results are listed in Table 2. In Case  2, two inclusions inside the medium were considered as shown in Figure 8(b). The optical properties of both inclusions were estimated by PSO, SPSO, and DE-PSO. In this case, the optical properties of the medium were taken as κa0=0.20m-1 and σs0=0.30m-1, respectively. The true value of the inclusions which needed to be estimated was κa1,σs1=1.0m-1,9.0m-1 and κa2,σs2=4.5m-1,5.5m-1, respectively. The size and location of the two inclusions were z1,1,z2,1,r1=0.3,0.3,0.2 m and z1,2,z2,2,r2=0.7,0.6,0.15 m. The retrieving results of optical parameters are shown in Tables 3 and 4.Table 2
The retrieving results of the optical properties for Case  1 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Condition
True values or inversion results [m−1]
Relative errors [%]
κ
a
1 = 1.0
σ
s
1 = 9.0
κ
a
1
σ
s
1
PSO
PI
1.0185
9.1372
1.846
1.524
FI
0.9964
8.9671
0.356
0.366
SPSO
PI
0.9880
8.8960
1.200
1.155
FI
0.9979
9.0151
0.206
0.167
DE-PSO
PI
1.0098
8.9138
0.983
0.958
FI
1.0008
8.9956
0.077
0.049Table 3
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO of different generations.
Generations
True value [m−1]
Algorithm
Retrieving results [m−1]
Relative errors [%]
10
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9750, 9.3216)
(2.4962, 3.5738)
(7.4025, 5.9205)
(64.5009, 7.6456)
SPSO
(0.9803, 8.9818)
(1.9718, 0.2021)
(5.4861, 0.3347)
(21.9131, 93.9149)
DE-PSO
(0.9844, 8.8901)
(1.5635, 1.2207)
(3.3674, 3.7704)
(25.1700, 31.4477)
50
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.4496, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9995, 8.9996)
(0.0478, 0.0040)
(3.5513, 4.6093)
(21.0814, 16.1946)
DE-PSO
(1.0021, 9.0036)
(0.2071, 0.0397)
(3.3049, 4.0275)
(26.5584, 26.7718)
100
(
κ
a
1
,
σ
s
1
) = (1.0,9.0) (κa2,σs2) = (4.5,5.5)
PSO
(0.9790, 9.0331)
(2.1025, 0.3679)
(6.449, 7.5222)
(43.3244, 36.7678)
SPSO
(0.9998, 8.9995)
(0.0234, 0.0054)
(3.5119, 4.8832)
(21.9588, 11.2147)
DE-PSO
(1.0007, 9.0042)
(0.0696, 0.0472)
(3.6868, 4.8677)
(18.0717, 11.4957)Table 4
The retrieving results of the optical properties for Case  2 using PSO, SPSO, and DE-PSO without measurement error.
Algorithm
Time [s]
True values [m−1]
Retrieving results [m−1]
Relative errors [%]
PSO
25476
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(0.9985, 9.0019)
(0.1486, 0.0207)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(6.9803, 7.5019)
(55.118, 36.398)
SPSO
32775
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0000, 8.9996)
(0.0034, 0.0046)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.4444, 5.4556)
(1.2349, 0.8077)
DE-PSO
49554
(
κ
a
1
,
σ
s
1
) = (1.0,9.0)
(1.0005, 9.0036)
(0.0544, 0.0397)
(
κ
a
2
,
σ
s
2
) = (4.5,5.5)
(4.5248, 5.4454)
(0.5508, 0.9925)Figure 7
The incident types of pulse laser.Figure 8
The schematic of pulse laser irradiation and the measurement position of time-domain signals of the two test cases.
(a)
(b)It can be seen from Table2 that κa1 and σs1 can be retrieved accurately under both PI and FI conditions by PSO, SPSO, and DE-PSO when there is only one inclusion. Furthermore, it also demonstrates that the retrieving results under FI condition are more accurate than those under PI condition, the reason of which is that, under FI condition, the reflectance and transmittance signals are stronger and are carrying more information available. Meanwhile, the DE-PSO algorithm shows its priority in reconstructing the inclusion properties. Under both conditions, the relative errors of retrieval results are smaller than those of the SPSO and PSO. Table 3 demonstrates the retrieving results of the optical properties for Case  2 of different generations. As shown in Table 4, the results of SPSO and DE-PSO are more accurate than those of the PSO. Furthermore, under FI condition, for Case  2, the optical properties of both inclusions can be estimated accurately without measurement errors and with the fact that the estimated relative errors of results for inclusion B (see Figure 8) are larger than those of inclusion A. However, the calculation time required of the DE-PSO is more than that of the SPSO and PSO. Thus, it can be concluded that SPSO and DE-PSO can retrieved the optical properties in a 2D medium accurately but with longer calculation time. The value of the objective function of PSO, DE-PSO, and SPSO with two inclusions is shown in Figure 9. It can be seen that the PSO fall into local convergence quickly and stop searching. In the earlier stage of the inverse process the objective function of SPSO declines faster than that of DE-PSO. However, in the later period of the inverse process, the SPSO almost stops searching and the objective function of DE-PSO continues declining which leads to the result that the retrieving results of DE-PSO are more accurate than those of the SPSO.Figure 9
The objective function of PSO, DE-PSO, and SPSO in Case  2.


## 4.3.2. Estimation of Geometrical Parameters

In the inverse procedure, the long sampling span will reduce the computational efficiency. Hence, it is of great significance to select the proper sampling span. Thus, the sensitivity of radiative signal to the inverse parameters needs to be analyzed. The sensitivity coefficient is one of the most important characteristic parameters in the sensitivity analysis, which is the first derivative of the radiative signals to a certain inverse parameter. The sensitivity coefficient is defined as(19)Smiρj=∂ρj∂mimi=m0=ρjm0+m0Δ-ρjm0-m0Δ2m0Δ,where mi denotes the independent variable which stands for z1, z2, and r; Δ represents a tiny change; ρj is the radiative signals of each boundary.The sensitivities of the inclusion location (z1,z2) and size r for the signals D1–D4 are shown in Figure 10. The parameters are set the same as Case  1 described in Section 4.3.1. It demonstrates that the sensitivity coefficient is relatively larger in the span of ct∈0,3m. So this span can be chosen as the sampling span for inverse analysis. The same results are expected in other cases.Figure 10
The sensitivity of the radiation signals for the geometric parameters of the inclusion.
(a)
(b)
(c)In many conditions, the shapes and optical properties of inclusions in the media are known in the actual application cases of reconstructing the medium’s inside information. For example, a tumor or cyst in the biological tissue is spherical or spherical-like, and the defects in the specimen, such as the air pore, are mainly the shape of a cube, sphere, or ellipsoid. In addition, the optical properties of human and animals tissues, in vitro or in vivo, are studied thoroughly [45]. Consequently, the geometrical feature and optical properties of inclusions could be used as a priori information in these cases, which will reduce the complexity of the reconstruction technique and improve the quality and efficiency of the reconstruction results. Therefore, on the condition that the inclusion shape and the optical properties are known, the inverse problem turns into the retrieval of the size and location of the inclusions.In this section, the inclusion location (z1, z2) and size r of cases (z1 = z2 = 0.25 m, r=0.25 m), (z1 = z2 = 0.5 m, r=0.25 m), and (z1 = z2 = 0.75 m, r=0.25 m) were estimated simultaneously. The size of the 2D medium was set as 1.0m×1.0m. The absorption coefficient κa0 and scattering coefficient σs0 of the medium were 0.2 m−1 and 0.3 m−1, respectively. There was one circular inclusion in the 2D medium with absorption coefficient κa1 of 1.0 m−1 and scattering coefficient σs1 of 9.0 m−1, but the inclusion size and location were unknown, which needed to be retrieved. The FI incidence type and the data of the detectors 1, 2, 3, and 4 were used for inverse analysis by using the PSO, SPSO, and DE-PSO algorithms with measurement errors. The search span of z1, z2, and r was set as 0,1 m. The results are shown in Tables 5 and 6. Figure 11 shows the comparison of the exact and reconstructed reflectance signal when z1=z2=r=0.25 m with 10% measurement error.Table 5
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO without measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2446
0.2457
0.2479
2.163
1.712
0.837
SPSO
0.2499
0.2500
0.2500
0.027
0.007
0.017
DE-PSO
0.2500
0.2500
0.2500
0.015
0.005
0.007
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5127
0.5045
0.2584
2.552
0.893
3.369
SPSO
0.500
0.500
0.2499
0.042
0.002
0.033
DE-PSO
0.4999
0.500
0.2499
0.019
0.006
0.021
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.7371
0.7694
0.2281
1.718
2.591
8.779
SPSO
0.7497
0.7500
0.2500
0.041
0.005
0.105
DE-PSO
0.7498
0.7500
0.2498
0.026
0.005
0.067Table 6
The inverse results of size and location of circular inclusion using PSO, SPSO, and DE-PSO with 10% measurement error under FI condition.
Algorithm
True values or inverse results [m]
Relative errors [%]
z
1
z
2
r
z
1
z
2
r
True value
0.2500
0.2500
0.2500
—
—
—
PSO
0.2804
0.2669
0.2676
12.18
6.747
7.051
SPSO
0.2509
0.2506
0.2492
0.356
0.243
0.339
DE-PSO
0.2340
0.2412
0.2452
6.395
3.514
1.934
True value
0.5000
0.5000
0.2500
—
—
—
PSO
0.5993
0.3354
0.1238
19.88
32.91
50.49
SPSO
0.5006
0.5267
0.2453
0.120
5.331
1.894
DE-PSO
0.5150
0.5302
0.2386
2.996
6.045
4.549
True value
0.7500
0.7500
0.2500
—
—
—
PSO
0.8371
0.5694
0.2981
11.61
24.08
19.22
SPSO
0.7356
0.7681
0.2265
1.926
2.419
9.418
DE-PSO
0.7237
0.7786
0.2324
3.508
3.821
7.041Figure 11
The comparison of the exact and reconstructed reflectance signals.It can be seen that the size and location can be estimated accurately using PSO, SPSO, and DE-PSO without measurement errors under FI condition. Furthermore, the results of DE-PSO are more accurate than those of the SPSO and PSO. The relative errors are reasonable even with 10% measurement error. However, it is worth noting that the relative errors of DE-PSO with 10% measurement error are bigger than those of the SPSO, which means the SPSO algorithm is more robust than the DE-PSO.


## 5. Conclusions

In the present study, the standard PSO, SPSO, and DE-PSO algorithms were applied to estimate the size, location, and optical parameters of the circular inclusions in a 2D rectangular medium. The conclusion was obtained that the retrieving results of absorption coefficient and scattering coefficient under FI condition are more accurate than those under PI condition, the reason of which is that, under FI condition, the reflectance and transmittance signals are stronger and are carrying more information available. In addition, the DE-PSO algorithm shows its priority in reconstructing the inclusion properties. When only one inclusion is considered, under both laser incident conditions, the results’ relative errors of DE-PSO are smaller than those of the SPSO. However, when estimating the size and location of inclusion media, the SPSO is more robust than DE-PSO. Furthermore, compared with the performance of standard PSO, both the SPSO and DE-PSO are proved to be more accurate and robust, which have the potential to be applied in the field of 2D inverse transient radiative problems.



---

*Source: 680823-2015-03-25.xml*