1. Introduction

While propagating through a highly scattering medium, photons are scattered multiple times, elastically or inelastically, and have a chance of being absorbed. Absorption and elastic scattering are the physical phenomena modelled in diffuse optics and accounted by the diffusion equation (DE) [1]. Diffuse optics techniques enable depth sensitivity of up to a few centimetres, non-invasively through the diffusive media, such as biological tissues. Raman scattering encodes the information about the molecule’s vibrational state and therefore provides excellent chemical specificity. The combination of diffuse optics and Raman spectroscopy can be applied for retrieving chemically specific information about the molecules at different depths inside the diffusive medium. Some of the fields of potential applications are clinical diagnostics [2–7], pharmaceuticals [8–10] and food analysis [11,12], security control [13–15] and cultural heritage [16–18]. In these areas it is of great importance to extract the chemical specificities, i.e., Raman spectra, of deeper layers when examining layered structures. There are several proposed techniques to separate the optical signals coming from different depths: spatially offset Raman spectroscopy (SORS) [2,19–21], frequency offset Raman spectroscopy (FORS) [22], transmittance Raman spectroscopy (TRS) [23,24] and time-domain or time-resolved diffuse Raman spectroscopy (TD-DIRS) [25–27].

TRS allows whole-volume probing with the light source and detector being on the opposite surfaces of the bounded medium (transmittance configuration) [23]. For the other abovementioned techniques, both the light injection point and Raman signal collection point are placed on the same surface (reflectance configuration) and an independent parameter is varied to provide different depth sensitivities, namely the source-detector distance in SORS, the excitation frequency in FORS, or the photon travelling time in TD-DIRS. More in detail, SORS relies on the increased weight on the signal coming from deeper layers of the medium upon increasing the source-detector distance [19–21,28]. FORS exploits the changes in optical properties with wavelength — e.g. the decreasing absorption and scattering coefficients in the 600-700 nm range in biological tissues — which leads to a significant wavelength-dependent mean propagation depth, being dependent for continuous wave (CW) measurement both on absorption and on scattering [22]. TD-DIRS [25–27] relies on the so-called "time encodes depth" principle, that is, photons that travelled longer in the medium also have a larger mean-propagation depth, at every source-detector distance, even null [29–31]. In this approach, in addition to the depth sensitivity, it is also possible to improve the signal-to-noise ratio by selecting the proper time-gating to reduce the contribution of fluorescence from the background. Finally, the covered pathlength contains quantitative information about the Raman emission.

The attempts to reconstruct the Raman spectra of inhomogeneous layered media known in literature so far used more qualitative rather than model-based quantitative approaches. Time-gating of photons has been applied to discriminate between the layers only on an orientational basis – the signals from the earlier time gates have been ascribed to the surface layer, while the later time gates signals have been accounted for both the superficial and deeper layer spectrum contribution [32]. In SORS different empirical approaches have been proposed for reconstruction or depth identification, based – for instance – on the (scaled) subtraction of signals at short source-detector distances from those at large source-detector distances to cancel out the superficial signals [2,19,21] or on Monte Carlo calibration [28]. After the scaled subtraction to separate the layers in SORS, for data processing, more quantitative analysis relying on the multivariate techniques was presented in literature [9,19,33].

There is more information contained in time-domain DE and it can be used to mathematically rigorously reconstruct Raman spectra of different layers. Here we present a novel model and methods for two-layer Raman spectra reconstruction in diffusive media, using the time-domain approach. The analytical model, based on the DE and further assumptions that will be described in the following section, yet simple, enables mathematically rigorous, elegant and fast reconstruction of Raman spectra of layered structures, where each layer is homogeneous with its own Raman spectrum. In this paper, we will demonstrate how the Raman emission spectra of the two layers can be well separated, testing the performance limitations on simulated data in different scenarios and proving the concept on real measurement data from laboratory phantoms. This is, to the best of our knowledge, the first time-domain reconstruction of Raman spectra of multi-layered structures based on analytical models.

2. Model

2.1 Hypotheses

We consider a two-layer medium, represented by a cylinder (see Fig. 1(a)). The medium is highly scattering, with layers having thicknesses, absorption coefficients, reduced scattering coefficients and refractive indices labelled as in Fig. 1(a). We also assume that both layers satisfy the conditions for application of the Diffusion Approximation valid under the condition that the absorption coefficient is much lower than the scattering coefficient [1,34]. A laser excitation light enters perpendicularly the medium from the top surface of the upper layer, at source position (S), while the detector (D) is placed on the same surface (reflectance configuration), at distance $\rho$ from the source (see Fig. 1(a)). The source excitation is assumed in the theoretical modelling as Dirac delta function both in time (delta pulse) and space (point-like). The time-domain diffusion equation and Green’s function method can be applied to find a solution for fluence at any position and time, and for any kind of source and detector (see Chapter 2 of Ref. [1]). Each layer separately is homogeneous – optical properties within the layers are independent of the position (constant). Our goal is to reconstruct the dependence of the Raman scattering coefficient on wavelength, for both layers, that is, Raman spectra of the two layers. We start from the measurements obtained by the detector – see Fig. 1(b) and (c) for an insight in the temporal curves expected at wavelengths around the characteristic Raman peaks of the two layers. These curves were obtained by simulations, using the theoretical model developed in this work (see Section 2.3) with the addition of random noise (Poisson distribution) and they show: 1) what is the expected temporal profile of the signals coming from upper or lower layer and their dynamic ranges – after certain time has passed, signals have comparable levels; 2) how that temporal profile is affected by noise, especially in the case of Raman signal coming from lower layer, when the thickness of the upper layer is increased – compare $d_1=10\,\mathrm {mm}$ (b) with $d_1=15\,\mathrm {mm}$ (c). In the modelling phase, we assume existence of Raman scatterers superimposed to the background. This superimposition will be called “the whole medium”, having Raman scatterers added to “the background medium”. Moreover, Raman scattering is assumed to be the only source of wavelength shift (inelastic scattering) of the migrating photons. Light sources for Raman photons, therefore, can exist anywhere inside the medium, and their intensities are determined by light fluences at specific positions and times. The key hypothesis to develop a Raman two-layer model in this work is that the optical properties of the background medium are identical at excitation and Raman emission wavelengths, which is known as the heuristic approximation [35,36].

 

Fig. 1. a) Two-layer cylindrical slab geometry with top layer thickness $d_1$ and bottom layer thickness $d_2$; Examples of photon counts detected (in simulations) around peak wavelengths corresponding to top and bottom layer signal for: b) $d_1 = 10\,\mathrm {mm}$; c) $d_1 = 15\,\mathrm {mm}$. The top and bottom layer spectra in simulations are taken as Gaussians, with the same standard deviation, but centred at different wavelengths, which are called peak wavelengths – for the top and bottom layer. For the details about the optical properties, see Section 4.1.

Download Full Size | PDF

2.2 Diffusion equation in a two-layer medium

Assuming an isotropic point source at $\vec {r}_s$ which is a delta impulse in time, the diffusion equation for the light fluence $\Phi$ in the inhomogeneous medium can be written as [1]:

In modelling the real problem, it must be considered that the source is an external beam impacting the medium. For this purpose, the position of the isotropic source $z_{s}$ is usually set to $z_{s}=1/\mu _s'$ to have a source of unitary strength used to model an external pencil beam of unitary strength impacting the medium [1,37]. In this scheme, it is assumed that when dealing with a two-layer medium, the position of the source is always located inside the first layer, thus simulating a pencil beam impacting externally from this layer. The detailed derivation of the solution is given in Chapter 6 of Ref. [1] where the fluence in each layer is described with an eigenfunctions’ expansion. At the interface of the two layers is used the continuity of the photon flux, while the fluence at the interface can have a discontinuity generated by a refractive index mismatch between the layers. Here we will just present the final result and apply it to the reflectance measurement case. The assumption is that the source is at $\vec {r}_s = (0,0,z_s)$ in the first layer. The solutions for the fluences in the two layers are given below [1]. For the first layer, when $0\leq z < d_1$, the solution is

For the second layer, when $d_1 \leq z \leq d_1+d_2$, the solution is

Particularly for cylindrical geometry, the Bessel function of the first kind of order zero, $J_0$, appears. Due to the eigenfunctions’ expansion (going from $l=1$ to $+\infty$ and for a particular $l$, from $n=1$ to $+\infty$), the wavenumbers $K_{\rho l}$, $K_{z1ln}$ $K_{z2ln}$ appear. The coefficient $N_{ln}$ is defined in Appendix I of Ref. [1] and it is a normalizing factor for the eigenfunction expansion that assure conservation of energy of the injected light source. The symbol * refers to complex conjugate of a number. It is worth noting that the complex notation is necessary for the expansion of eigenfunctions used to represent the solution, which, depending on the optical properties of the layers, may require imaginary eigenvalues, as can be verified in Ref. [1]. First, the solutions for $K_{\rho l}$ are found from the boundary condition used at $\rho =L$ (see Eqs. (4) and (5)) that dictates the fluence to vanish at the extreme lateral boundary if this is enough far from the source [1], i.e.

Then, the other two wavenumbers, $K_{z1ln}$ and $K_{z2ln}$, are obtained as roots of the following equations (pure imaginary roots are also accepted)

Finally, to obtain the reflectance, we first express the fluence from 7, substituting $z=0$, and then, assuming the direction of the $z$ axis entering the first layer, and using Fick’s law [1]

2.3 Heuristic Raman two-layer model

A heuristic model describing the Raman signal obtained under the hypothesis that the optical properties of the background medium at excitation ($\lambda$) and emission ($\lambda _e$) wavelengths are the same was developed in the previous literature for a homogeneous medium [35,36,38]. From now on the subscript $e$ will be used to denote quantities at the Raman emission wavelength $\lambda _e$, while quantities without such subscript pertain to excitation wavelength $\lambda$.

We extend here the same approach to a two-layer medium with optical properties of the background medium at excitation ($\lambda$) and emission ($\lambda _e$) wavelengths having the same values. Thus, the optical properties can be indifferently referred to $\lambda$ or $\lambda _e$. We denote by $\mu _{a1}$ and $\mu _{a2}$ the absorption coefficients of first and second layer, $\mu _{s1}'$ and $\mu _{s2}'$ the reduced scattering coefficients of first and second layer, $\mu _{sR1}$ and $\mu _{sR2}$ the Raman scattering coefficients of first and second layer, $n_{1}$ and $n_{2}$ the refractive indices in first and second layer, and finally $d_1$ and $d_2$ the thicknesses of first and second layer. We now apply analogous reasoning to the one done by Martelli et al. in [36], but now for a two-layer geometry. Additionally, let us denote by $\mu _{a1b}$ and $\mu _{a2b}$ the absorption coefficients of background (when there are no Raman scatterers) at excitation wavelength $\lambda$, and by $\mu _{a1be}$ and $\mu _{a2be}$ the same at emission wavelength $\lambda _e$. We analogously define the reduced scattering coefficients for the background at $\lambda$: $\mu _{s1b}'$, $\mu _{s2b}'$, and at $\lambda _e$: $\mu _{s1be}'$, $\mu _{s2be}'$. Before going into details about the modelling, it is important to note that the modelization proposed in this section is possible thanks to the fact that, for the purpose of the energy balance at $\lambda$, the Raman scattering expressed by the Raman scattering coefficient can be modelled as an equivalent effective “absorption” term of the radiative transfer equation [36], i.e., its practical effect at $\lambda$ is to remove light intensity as an absorption term does.

The heuristic approximation we apply here assumes the following: $\mu _{a1b}=\mu _{a1be}$, $\mu _{a2b}=\mu _{a2be}$, $\mu _{s1b}'=\mu _{s1be}'$, $\mu _{s2b}'=\mu _{s2be}'$, $v_1=v_{1e}$, $v_2=v_{2e}$ (same speed of light in both layers at $\lambda$ and $\lambda _e$, or in other words same refractive indices at excitation and Raman emission wavelengths).

Since the Raman scattering coefficient is usually very small ($\mu _{sR1,2}\ll \mu _{s1,2b}$ and usually less than $10^{-6}$ mm$^{-1}$), the probability of multiple Raman scattering is negligible, allowing us to assume zero or one Raman scattering event per photon. A Raman scattering event undergone by a photon can be considered equivalent, in energy balance at $\lambda$, to the absorption of a photon at $\lambda$ and its immediate re-emission at wavelength $\lambda _e$.

Adding Raman scatterers to the background, the absorption coefficient at $\lambda$ is affected: $\mu _{a1}=\mu _{a1b}+\mu _{sR1}$, $\mu _{a2}=\mu _{a2b}+\mu _{sR2}$, but not at $\lambda _e$: $\mu _{a1e}=\mu _{a1be}$, $\mu _{a2e}=\mu _{a2be}$, because of the above argument. Reduced scattering coefficients are the same for the background and the whole medium when Raman scatterers are added: $\mu _{s1}'=\mu _{s1b}'$, $\mu _{s2}'=\mu _{s2b}'$ at $\lambda$, and $\mu _{s1e}'=\mu _{s1be}'$, $\mu _{s2e}'=\mu _{s2be}'$ at $\lambda _e$. Before deriving the equations, it is useful to keep a summary – see Table 1 – of optical parameter values in the heuristic approximation at the two wavelengths considered.

Table 1. Overview of absorption and reduced scattering coefficient values at excitation wavelength $\lambda$ and Raman emission wavelength $\lambda _e$, for two layers and in two cases – “the background medium” and “the whole medium”, i.e., background with addition of Raman scatterers.

View Table

The fluence of Raman photons (at $\lambda _e$) in the whole medium can thus be set equal to the fluence of Raman photons in the background medium

The scaling property of absorption [1] can be used to relate the fluence at excitation wavelength in the whole medium to the fluence in the background medium as

Since the optical parameters at $\lambda$ and $\lambda _e$ are the same as well as the scattering phase function, once a Raman scattering occurs, in theory the wavelength of the photon changes from $\lambda$ to $\lambda _e$, but in practice, there is no effect on its propagation. The photon migration continues in the very same way as if the photon had experienced elastic instead of inelastic scattering. Moreover, in Raman scattering events, absorption (removal) of photons at $\lambda$ is not taken into account by keeping background values of absorption coefficients $\mu _{aib}$ instead of using $\mu _{aib}+\mu _{sRi}$, for $i=1, 2$. Therefore, the fluence $\Phi (\vec {r},\mu _{a1b},\mu _{s1b}',\mu _{a2b},\mu _{s2b}',t,\lambda )$ (first term at the right-hand side of Eq. (15)) includes photons at both wavelengths $\lambda$ and $\lambda _e$. The fluence at the left-hand side of Eq. (15) is given by the contribution of those photons that avoided Raman scattering. The difference between the fluence in the background medium at $\lambda$ (given by photons propagating at $\lambda$ when Raman scattering is absent) and the fluence at $\lambda$ in the whole medium (given by photons still propagating at $\lambda$ since did not have Raman interactions) consequently corresponds to the fluence of Raman emitted photons, i.e.

Combining the previous Eqs. (14), (15) and (16) the fluence at the Raman emission wavelength can be rigorously expressed as

Given the typical low values of $\mu _{sR1}$ and $\mu _{sR2}$ (of the order of 10$^{-6}$ mm$^{-1}$), Eq. (17) can be usually further simplified by expanding the exponential in Taylor polynomial of the first order

Calculating the integrals, since the fluence rate term does not depend on the single variables $t_1$ and $t_2$, we obtain

Using Fick’s law, the same relation of Eq. (20) holds true also for the reflectance from the two-layer medium:

The above formula (20) has the same physical meaning as the previous formula obtained for the homogeneous medium in Ref. [36]. The obtained solution consists of two terms that denote the contribution of the two layers to Raman scattering interactions and subsequent re-emission. Each of these two terms is the product of two factors: the first factor is the fluence at the excitation wavelength, while the second factor is proportional to the time spent by detected photons inside the layers. The presence of fluence at the excitation wavelength expresses the fact that Raman photons generated at the emission wavelength continue their migration in the medium as they were doing at the excitation wavelength, and this is due to the assumption used in this approach. The second factor of each term, being proportional to $\mu _{sRi} v_i \langle t_i \rangle$ ($i=1, 2$), expresses the fact that the probability of generating Raman photons in each layer is proportional to the average pathlength covered by photons in the layer.

An approach similar, but not the same, to that here adopted about the assumption used in the presented model was that of Liebert et al. [39] to calculate, by means of Monte Carlo (MC) simulations, the time-resolved fluorescence in layered turbid media. An explanation of a recent MC code with which this approach can be implemented can be found in [40]. We note that, using the MC approach, the implementation of the inverse problem solution may be more problematic and difficult because of the computation time of the procedure and of the noise affecting the final results.

2.3 Another approach to derive the final formula

Let us consider only the detected reflectances at excitation and Raman emission wavelengths, and optical parameters that affect their values. If we apply the same reasoning from the beginning of Section 2.3 with Raman scatterers added to the background, we can express a small change in the detected reflectance at $\lambda$ as

Since the absorption coefficients in the whole medium consisting of background and added Raman scatterers are $\mu _{ai}=\mu _{aib}+\mu _{sRi}$, for both layers $i=1, 2$, the increments of the background absorptions are $\Delta \mu _{aib}=\mu _{sRi}$. In our heuristic approximation, photon migration is unaffected even if the wavelength is shifted from $\lambda$ to $\lambda _e$. Therefore, a deficit in number of detected photons at $\lambda$ must be equal to the number of Raman emitted photons at $\lambda _e$ reaching the detector, or mathematically

This is exactly the same as 23, when average times $\langle t_1\rangle$ and $\langle t_2\rangle$ are expressed through 21 and 22.

Note that this recipe may be easily applied to further generalize two-layer geometry solution to many-layer geometries. Total derivative of a multivariable function, reflectance as a function of many-layer optical parameters, is calculated by simply applying mathematical rules, no matter how many layers there are.

In accordance with the proposed modelling, the reflectance at emission wavelength $R_e$ can be calculated by using an analytical expression for the reflectance from a layered medium at excitation wavelength. Within the diffusion equation we have several available solutions that can be used for this scope in addition to the time domain solutions for a two-layer medium reported in Refs. [1,41,42]. In particular, we mention the $N$-layered DE solution by Liemert and Kienle reported and used in Refs. [34,43–45]. It is worth to notice that such solution is used in time domain NIRS applications for determining the optical properties of tissue [46]. The solution from [1] was presented in Section 2.2.

2.4 Inverse model

The expression for the Raman reflectance detected at some wavelength $\lambda$ and time $t$, by the detector positioned on the top surface of the two-layer medium can be written as

The weights are expressed as functions of both time and wavelength to cover the more general case, even though Eqs. (27), 21, 22 do not predict the dependence of weights on wavelengths due to the heuristic approximation. Optical parameters found in Eqs. (21) and 22 in fact depend on wavelength, therefore, in principle, weights should be functions of both, wavelength and time. However, in the heuristic approximation, weights are the same on excitation and Raman emission wavelength.

If $t^s$ and $t^e$ are some start and end times of a certain time gate, integration of the Eq. (28) yields

From now on we define the gated measurements: $M_e(\lambda,t^{s}, t^{e})=\int _{t^s}^{t^e} R_e(\lambda,t)\mathrm dt$, as well as the gated sensitivity matrix elements: $w_{i}^{g}(\lambda,t^{s}, t^{e})=\int _{t^s}^{t^e} w_i(\lambda,t)\mathrm dt$, for $i=1,2$. If the number of time gates is $N_g$, then the system of equations of form (28) will consist of $N_g$ equations with 2 unknowns, for each wavelength. If the number of wavelengths is $N_w$, then the inverse problem will consist of finding solutions to $N_w$ systems, each of $N_g$ equations with 2 unknowns. For each wavelength $\lambda _i$, $i\in \left \{1,\,2,\,3,\dots,N_w\right \}$, the system to be solved can be written as:

To retrieve the spectrum $\mathbf {S}_{2\times 1}(\lambda _i)$, one would need at least 2 equations for each $\lambda _i$. In practice, the number of gates $N_g$ is greater than 2. From now on, we will assume $N_g>2$, so we have an overdetermined system. To simplify notation, here we skip explicit writing of dependence on wavelength $\lambda$, but always keep in mind that the system to be solved is defined for a specific $\lambda$. There are many possible solutions for the spectrum $\mathbf {S}$, but none of them will be exact, since the 2 parameters $\mu _{sR1}$ and $\mu _{sR2}$ have only 2 degrees of freedom and in general case it is not possible to simultaneously satisfy $N_g>2$ equations. We will look for the least square solution of this overdetermined system of linear equations, which minimises the norm of vector $\mathbf {M_e} - \mathbf {W_g}\mathbf {S}$ [47]. Formally, a vector $\mathbf {S^*}$, the least square solution of the problem, is found from the following system of normal equations:

Here we highlight the term “normal equations”, since this system is now square, $\mathbf {W_g}^{\mathrm T}\mathbf {W_g}$ is a $2\times 2$ matrix in case of 2-layer media. A retrieval of $\mathbf {S^*}$ is now very simple. However, this has to be done on all wavelengths $\lambda _i$, $i=1,2,\dots,N_w$.

3. Methods

3.1 Reconstruction

The process of reconstructing the unknown spectra from the known measurements is actually the process of solving the inverse problem (see Section 2.4). The measurements represent the time-gated values of the detected signal, i.e., number of photon counts detected in specific time intervals (gates). To solve the inverse problem, we need the theoretically calculated forward model (see Section 2.3). The forward model calculations of the gated sensitivity matrices depend on the gating chosen for our measurement data, as well as on the geometrical and optical parameters given as inputs (thicknesses of layers, absorption and reduced scattering coefficients and refractive indices). All the methods used for forward model and inverse problem computations were implemented in Python.

3.2 Numerical simulations

To test our reconstruction methods, we performed numerical simulations in different scenarios, evaluating the performance by means of different figures of merit (see Section 4.2). In each simulation, starting from theoretical (ground-truth) spectra of top and bottom layer, forward model calculation of the output signal is performed based on the two-layer diffuse Raman model described in Section 2.3. The simulated measurements are then obtained by the convolution of this signal with the instrument response function of the system, gating it (summing the number of photons detected over time, in specific time intervals), and finally adding shot noise (simulated as a random variable with Poisson distribution). Although numerical simulations implemented in this way account for realistic effects, there is still a need to validate the model on real experimental, not simulated data. The experimental setup, materials and methods are explained in Section 3.3 and 3.4.

3.3 Experimental setup

The raw experimental data used here to test the reconstruction approach are the same as used in Ref. [27], where all details on the experimental setup and measurements are given. To perform the temporal Raman signal acquisition, a time resolved single photon counting camera (TCSPC) was used [48]. Technical details about the camera are given in [49]. The system setup is schematically represented in Fig. 2(a). The pulsed laser, with wavelength of $532\,\mathrm {nm}$ and power of $100\,\mathrm {mW}$ was used as an excitation source. It was coupled with a $100\,{\mu }$m-core optical fiber to the probe, shown in Fig. 2(b). The probe was operating without contact with the phantom, having ring illumination (provided by axicon lens) and point collection, see Fig. 2(c). For cleaner excitation, a narrow band-pass filter at $532\,\mathrm {nm}$ was used. On the detection side, lenses L1, L2, L3, L4 were used to stop light from unwanted directions, and a long-pass filter, blocking up to $532\,\mathrm {nm}$ was used to select only the useful signal for the spectral analysis. The probe and the spectrometer were connected with the $1\,\mathrm {mm}$-core fiber, with numerical aperture $\mathrm {NA} = 0.39$. The spectrograph had a $200\,\mathrm {\mu }$m entrance slit, $f/4$ diameter of the entrance pupil and grating of 1200 grooves per millimetre. The spectrally resolved signal was then captured by the TCSPC camera. A histogram of Raman photons arrival times was being acquired for $300\,\mathrm {s}$.

 

Fig. 2. a) Schematic of the setup; b) Non-contact probe and the phantom; c) Ring illumination (green) and the collection point (red). Images a) and b) taken from [27].

Download Full Size | PDF

3.4 Phantom

The phantom designed to imitate a two-layer tissue was a silicone elastomer and marble slab, previously described in Ref. [27]. The top layer, made of silicone was $d_1 = 5\,\mathrm {mm}$ thick, with the absorption coefficient $\mu _{a1} = 0.011\,\mathrm {mm}^{-1}$, the reduced scattering coefficient $\mu _{s1}{\prime } = 1.65\,\mathrm {mm}^{-1}$, and the refractive index $n_1 = 1.41$. The bottom layer, made of marble, $\mathrm {CaCO_3}$ was $d_2 = 17\,\mathrm {mm}$ thick, with the optical parameters $\mu _{a2} = 0.003\,\mathrm {mm}^{-1}$, $\mu _{s2}{\prime } = 1.65\,\mathrm {mm}^{-1}$, and $n_2 = 1.66$. Optical characterization of the phantom was performed using the time domain diffuse optical spectrometer described in Ref. [50]. All the optical parameters are expressed at $532\,\mathrm {nm}$ wavelength, assuming to be constant over the Raman measurement range. The source-detector distance was $\rho =10\,\mathrm {mm}$.

4. Results and discussion

4.1 Examples of reconstruction in simulations

A set of exemplary simulations was performed on a two layer geometry represented in Fig. 1, with $d_1=10\,\mathrm {mm}$ and $d_2=60\,\mathrm {mm}$ (thick enough to appear semi-infinite ), and with the source-detector separation of $\rho =10\,\mathrm {mm}$. The optical parameters were chosen to be: $\mu _{a1}=0.011\,\mathrm {mm}^{-1}$, $\mu _{a2} = 0.003\,\mathrm {mm}^{-1}$, $\mu _{s1}{\prime } = 1.65\,\mathrm {mm}^{-1}$, $\mu _{s2}{\prime } = 1.50\,\mathrm {mm}^{-1}$, and $n_1=n_2 = 1.4$. All simulations in this paper assumed gaussian spectra for both layers, with standard deviation of $20\,\mathrm {cm}^{-1}$, the top layer centred around Raman shift of $1411\,\mathrm {cm}^{-1}$ and the bottom layer centred around $1087\,\mathrm {cm}^{-1}$, thus not affecting each other. Figure 3 shows the theoretical spectra of top and bottom layers (Fig. 3(a)), the reconstruction in an ideal case – no added noise and the instrument response function is Dirac delta (Fig. 3(b), the effect of a broader temporal response of $200\,\mathrm {ps}$ FWHM Gaussian (c), and the impact of shot noise (Fig. 3(d), (e), (f)) with instrument temporal resolution of $100\,\mathrm {ps}$. In the ideal case, the reconstruction is perfect within the numerical precision of the PC. The effect of a wider instrument response function (IRF) in the time domain was shown in Fig. 3(c). It was expected that as the temporal width of the IRF increases, or in other words, it is getting further from ideal Dirac delta function, the coupling between the signals coming from different layers becomes more significant. A worse temporal resolution leads to residual contamination of the upper layer spectrum in the lower layer reconstruction. As the number of photon counts increases, as shown in Fig. 3(d), (e), (f), the signal-to-noise ratio also increases, which results in better quality of reconstruction. A total number of counts in excess of $10^6$ is needed to resolve the two peaks and probably even more is required for complex spectral features.

 

Fig. 3. a) Theoretical spectra of the two layers used in simulations; Examples of reconstructions for different cases: b) ideal case; c) no noise, but IRF $\mathrm {FWHM} = 200\,\mathrm {ps}$; noise included, for various numbers of detected photons: d) $10^4$ counts; e) $10^6$ counts; f) $10^7$ counts.

Download Full Size | PDF

4.2 Effect of different parameters

The field of diffuse Raman spectroscopy is still at an early stage, and lacks a common standard terminology and set of relevant figures-of-merit. In SORS, the enhancement factor [51] is used to express the increased sensitivity to the lower layer for a larger $\rho$. Here, we specifically introduce three figures of merit, defined below, to systematically study the effect of various parameters on the quality of reconstruction, and to quantify the results:

Suppression as a parameter is interesting, especially for the lower layer reconstruction, $S_{LU}$, because it shows how convincing is the differentiation between the lower layer peak (what is expected to be in the reconstructed spectrum) and the upper layer peak (what is considered as undesired contamination). It is good if the area of the peak of the reconstructed lower layer spectrum is at least one order of magnitude greater than the area of the contamination in the reconstruction, originating from the top layer spectrum, i.e. if $S_{LU}>10$. For the upper layer reconstruction, values $S_{UL}$ are usually much higher.

The ensamble of different parameters affecting the reconstruction studied here are:

One could study the effect of other parameters, but we opted to systematically examine and present those which we observed to have a significant impact on the quality of reconstruction. These were mainly the optical properties of the upper layer, which highly influenced the levels of the signals coming from the lower layer, as well as one of the most important parameters, the thickness of the upper layer. Moreover, these parameters describe the superficial layer which physically hinders the detection of Raman photons from depth, and it is indeed of high interest to have a reliable reconstruction of the lower layer spectrum, to have it well decoupled from the upper layer spectrum. Number of detected photon counts and temporal width of the instrument response function are natural parameters to study in a detection system. While the former highly affects the reconstruction quality, the latter can be used to estimate when the coupling between the signals from the two layers becomes such that the two spectra cannot be separated well anymore. Fig. 4 sums up the main results which we will elaborate below.

 

Fig. 4. CNR of the top layer (up) and bottom (down) layer vs: a) top layer absorption coefficient; b) photon counts. Suppression of top layer components in the reconstructed spectra of the bottom layer vs: c) FWHM of the instrument IRF; d) top layer thickness for different photon counts; e) top layer absorption coefficient for different scattering coefficients of the top layer.

Download Full Size | PDF

The absorption of the top layer has negligible effect on the contrast-to-noise ratio of the bottom layer reconstruction. While for increasing absorption of the top layer it is evident that the $CNR$ of the top layer reconstruction is decreasing, it still maintains a respectable level (Fig. 4(a)).

The $CNR$ is significantly affected by the number of detected photons (Fig. 4(b)). The larger the number of photon counts, the higher $CNR$ for both layers.

The effect of the instrument response function can be neglected if its temporal profile is narrower than about $150\,\mathrm {ps}$, which can be seen in Fig. 4(c), looking at the upper/lower and lower/upper suppression.

Lower/upper suppression as function of the thickness of the top layer is plotted in logarithmic scale for different numbers of detected photon counts in Fig. 4(d). The suppression roughly decreases exponentially with the thickness of the top layer, while high enough signal level can extend the limit of thickness of the top layer that permits acceptable reconstruction.

Lower scattering coefficients of the top layer are favourable for bottom layer reconstruction, which can be seen in Fig. 4(e). If the scattering in the top layer is not too high, then the more focused beam of photons reaches the bottom layer, thus facilitating the reconstruction of the bottom layer spectrum with the increasing absorption (up to a certain limit).

A brief conclusion that can be drawn from the data provided in Fig. 4 is that the major factors that limit the reconstruction are thickness of the top layer and number of detected photon counts.

4.3 Raw experimental data

The Raman signals detected from the phantom experiments are resolved in time and wavelength. The temporal histogram had 1300 time windows, each $11.8\,\mathrm {ps}$ wide. The Raman shift was covered from $370\,\mathrm {cm}^{-1}$ to $3253\,\mathrm {cm}^{-1}$.

In Fig. 5(a) we present an illustrative example of the measured Raman reflectance signal from the top layer, bottom layer, as well as the background. The top layer signal was obtained by summing the raw signal contributions at wavelengths around the silicone peak ($2888\,\mathrm {cm}^{-1}$) and subtracting the corresponding background. The bottom layer signal was obtained similarly, by summing the contributions at wavelengths around the marble peak ($1087\,\mathrm {cm}^{-1}$) and subtracting its corresponding background. The background signal here was considered as the background corresponding to the silicone, because the signal from the top layer is dominant. Backgrounds corresponding to silicone and marble were obtained by summing the raw signal contributions at wavelengths close to, but just before the silicone peak, and just after the marble peak, respectively.

 

Fig. 5. Raw data measurements: a) temporal evolution of detected signal intensity from top layer (blue), bottom layer (orange) and background (green); b) gated signal (number of photon counts per gate), spectrally resolved and normalized, represented for 12 different time gates defined in the legend.

Download Full Size | PDF

In Fig. 5(b) we present one possible selection of time gates. The choice was to apply 12 time gates, the first one starting at $212.0\,\mathrm {ps}$ and the last one ending at $1486.4\,\mathrm {ps}$ (as it was done in [27]). The time gating was performed in software, after collecting the raw data. The plot shows the number of counts, normalized from 0 to 1, resolved by wavelengths. It is possible to notice the top layer (silicone) peaks (around $2888\,\mathrm {cm}^{-1}$) at early gates. They gradually disappear after around $800\,\mathrm {ps}$, which was expected considering the plot from Fig. 5(a). The bottom layer (marble) peak becomes noticeable after a few time gates (around $600\,\mathrm {ps}$). Figure 5(b) clearly shows that the later time gates mainly contribute to the reconstruction of the bottom layer spectrum, while the earlier time gates mainly contribute to the reconstruction of the top layer spectrum.

4.4 Reconstruction from phantom measurements

After the application of the suggested gating (see Fig. 5(b)), using the forward model with optical and geometrical parameters defined in Section 3.4, we obtain the gated sensitivity matrix $\mathbf {W_g}$. The inversion method described in 2.4, applied on raw gated measurements (without normalization) and with the knowledge of $\mathbf {W_g}$ produces the reconstructed spectra for the two layers $S_1(\lambda )$ and $S_2(\lambda )$. These two spectra are then normalized, and the comparison with silicone and marble reference spectra is plotted in Fig. 6. The reference spectra of silicone and marble were obtained using CW-SORS and FORS Raman spectrometer [22].We can see that the positions of the peaks of both layers are correctly reconstructed. However, there is a background that affects the intensities in the directly (raw) reconstructed and reference spectra. The effect is more noticeable in the raw reconstructed spectra. We believe it has its roots in the nature of the experiment, where the background fluorescence signal is significantly present, especially when the source emission wavelength is relatively short $532\,\mathrm {nm}$ (at such wavelengths, it is more likely to excite the fluorescence emission of the background). Additionally, the reason for the presence of some artefacts in the reconstructed spectra lays in the nature of the mathematical model, where in order to completely decouple the layers, some physical parameters should be as close to ideal as possible. A precise knowledge of the optical parameters of the medium as well as a very narrow temporal profile of the instrument response function and larger number of photon counts are desired. This is in line with the conclusions drawn from the simulations. The raw reconstructed spectra can be further cleared from some of these effects, such as overwhelming background, by applying methods to remove it and leave only significant peaks in the spectrum.

 

Fig. 6. Direct reconstruction (blue) from raw data, normalized, compared to the reference spectrum (orange) for the: a) top layer; b) bottom layer.

Download Full Size | PDF

The background removal can be performed in many ways, but in this paper we present the results of a simple one, cubic spline interpolation through the points not belonging to any peak. These points are all the points in the spectrum that satisfy the conditions: 1) left derivative (first derivative from the left) is negative, with absolute value less than the threshold; 2) right derivative (first derivative from the right) is positive, with absolute value less than the threshold; 3) the threshold is determined as a fraction of the maximum of the first derivative of the spectrum over the whole spectral range considered, and this fraction was $0.1$ for the first layer and $0.2$ for the second layer spectrum in this particular implementation. In Fig. 7 the background interpolation of top and bottom layer spectra are shown on the left – (a) and (b). The reconstructed spectra, compared to the reference spectra, after the background removal, are shown in Fig. 7 on the right, for top (c), and bottom (d) layer. The results show that the spectra of the two layers are indeed well separated.

 

Fig. 7. Reconstructed spectra: a, b) with background exemplary fit; c, d) after background subtraction.

Download Full Size | PDF

5. Conclusion

For the first time, Raman spectra of two-layer diffusive media have been rigorously reconstructed relying on the Diffusion Equation in the time-domain. The heuristic approach was used to derive a simplified model of diffuse Raman propagation from the solution of the Diffusion Equation for a two-layer geometry. Mathematical equations and methods have been implemented in computer software. The software, developed in Python, is capable of calculating the forward model with optical and geometrical parameters as inputs, then solving the inverse problem, obtaining the desired Raman spectra in around half a minute in total (with CPU Intel i5, 11th generation, $2.4\,\mathrm {GHz}$, $16\,\mathrm {GB}$ RAM). The novel model and methods have been applied in simulations, comparing the obtained reconstruction results with the ground-truth, in different scenarios that tested the effects of various parameters. It has been found that the most critical parameters for the reliable reconstruction of the two-layer spectra are the thickness of the top layer and the number of detected photons. Our model and methods have finally been validated on phantom measurements, reconstructing the two-layer Raman spectra of the phantom slab having silicone and marble as top and bottom layer materials, respectively. The top layer thickness was $5\,\mathrm {mm}$, which is below the safe limit of $10\,\mathrm {mm}$ suggested by simulation studies for total detected signal level of $10^6$ photon counts. As a result, silicone and marble Raman spectra have been well separated, despite the presence of background effects, thus proving the concept and showing the potential for future application of this model and methods.

Suggested directions for further study and improvements would be addition of wavelength dependence of the optical properties in the model, aiming to obtain more accurate gated sensitivity matrices. Another important factor to be optimized is the selection of time gates. The model can be further generalized to many-layer geometries. A key open question is how to infer the optical properties of the materials. In the data presented above, these were measured using the same system used exactly at the excitation wavelength (i.e. removing the filters). In some cases literature data could be available. The effect of a mismatch in the exact knowledge of the optical properties on the reconstructed spectra could be further investigated.

Finally, the goal that should be set is to develop methods for application in vivo, making them reliable and fast enough for safe diagnostics in many fields – from clinics to security.