Artificial neural networks for retrieving absorption and reduced scattering spectra from frequency-domain diffuse reflectance spectroscopy at short source-detector separation

: Diffuse reflectance spectroscopy (DRS) based on the frequency-domain (FD) technique has been employed to investigate the optical properties of deep tissues such as breast and brain using source to detector separation up to 40 mm. Due to the modeling and system limitations, efficient and precise determination of turbid sample optical properties from the FD diffuse reflectance acquired at a source-detector separation (SDS) of around 1 mm has not been demonstrated. In this study, we revealed that at SDS of 1 mm, acquiring FD diffuse reflectance at multiple frequencies is necessary for alleviating the influence of inevitable measurement uncertainty on the optical property recovery accuracy. Furthermore, we developed artificial neural networks (ANNs) trained by Monte Carlo simulation generated databases that were capable of efficiently determining FD reflectance at multiple frequencies. The ANNs could work in conjunction with a least-square optimization algorithm to rapidly (within 1 second), accurately (within 10%) quantify the sample optical properties from FD reflectance measured at SDS of 1 mm. In addition, we demonstrated that incorporating the steady-state apparatus into the FD DRS system with 1 mm SDS would enable obtaining broadband absorption and reduced scattering spectra of turbid samples in the wavelength range from 650 to 1000 nm.


Introduction
Diffuse Reflectance Spectroscopy (DRS) is a model based technique. It could be applied to study the properties and physiological functions of biological tissues at depths from the submm range to a few cm and has been used in various clinical applications [1][2][3][4]. In a typical DRS measurement scheme, light is injected into the sample at a certain location and the diffusely reflected light is collected by a detector placed at a distance from the light injection location. One of the DRS variants employs an intensity modulated light source to investigate the light amplitude demodulation and phase delay introduced by the turbid sample which can be in turn utilized to quantify the sample's absorption (μ a ) and reduced scattering (μ s ') properties [5,6]. Such techniques are usually categorized as the frequency-domain (FD) DRS method and has been successfully used in conjunction with the standard diffusion equation (SDE) [7] to characterize the optical properties of deep tissues, such as breast and brain, with source-detector separations (SDSs) longer than 10 mm [8][9][10]. Moreover, its miniaturization and system cost reduction have been demonstrated [11,12]; thus, it has a great potential of being widely used in the clinical setting for diagnostic applications and cancer treatment monitoring [12]. In spite of the success of FD DRS in deep tissue applications, to date, only very limited literatures with respect to the application of FD DRS to superficial tissue studies can be found. Kienle et al. theoretically studied the optical property recovery errors caused by the use of FD standard diffusion equation at SDSs from 1.25 to 10.00 mm, and suggested utilizing the Monte Carlo modeling to avoid the errors due to diffusion approximation [13]. Willmann et al. employed a small-volume FD oximetry with SDSs of 3 to 9 mm to determine the oxygen saturation and total hemoglobin content of skin [14]. Papazoglou et al. studied the skin optical property variation at various wound healing stages using FD DRS with SDSs of 4,8,12, and 16 mm [15]. Hornung et al. employed FD DRS at SDS of around 10 mm to study the normal cervical tissue and epithelial lesions [16,17]. Tseng et al. reported that the validity of FD standard diffusion equation would become tenure when the SDS was reduced to less than 5 mm and proposed an alternative measurement scheme for accurate skin optical property recovery at SDS of 3 mm using a two-layer standard diffusion model [18,19]. While it has be known that DRS measurements carried out at short SDSs facilitate superficial tissue investigation, to our knowledge, FD DRS measurements performed at SDSs shorter than 3 mm for superficial tissue investigation have not been shown. This is majorly due to the modeling and system limitations. Up until now, efficient and accurate modeling the FD diffuse reflectance at SDS of around 1 mm is still a challenging task. Besides, at short SDSs, the phase delay introduced by samples would become small and thus the optical property recovery results could be easily affected by the inherent random system phase fluctuation.
In general, FD DRS systems that employ one modulation frequency usually demand measurements carried out at several SDSs to accurately recover sample optical properties [10,13]; this is because that the frequency domain reflectance measured at a certain frequency and at a single SDS could sometimes ambiguously correspond to several (μ a , μ s ') combinations [13]. In contrast, FD DRS systems that employ multiple modulation frequencies could reliably recover the sample optical properties from the measurements performed at a single SDS [20]. Although the multi-frequency FD DRS systems usually require a more complicated configuration or more expensive instruments than the single-frequency counterpart, in terms of optical property recovery accuracy of homogeneous samples, these two apparatuses do not have a significant difference. Nevertheless, for investigating biological tissues of layered structure, one clear advantage for carrying out measurements at a single SDS is that the sampling volume is definite during the measurement. On the other hand, DRS systems that require taking measurements at several SDSs would have different sampling volumes for different SDS pairs and thus would lead to biased optical properties recovery results as well as reduced sensitivity to the tissue inhomogeneity.
In this study, we revealed the construction method of artificial neural networks (ANNs) that can accurately and efficiently compute the amplitude demodulation and phase delay at various modulation frequencies introduced by turbid samples at short SDSs (e.g., 1 mm). The ANNs trained by Monte Carlo simulation generated databases not only had accuracy level comparable to the Monte Carlo method, but also could work with a least-squares minimization routine to determine the optical properties of samples from the FD reflectance measured at short SDSs in less than 1 second. We further investigated the necessity of employing FD reflectances measured at multiple frequencies for accurate sample optical property recovery at 1 mm SDS. Finally, we demonstrated that incorporating the steady-state (SS) apparatus into the FD DRS system with 1 mm SDS would enable obtaining broadband absorption and reduced scattering spectra of turbid samples in the wavelength range from 650 to 1000 nm.

Steady-state and frequency-domain combined DRS system
The SSFD (Steady-State and Frequency-Domain) DRS technique has been employed for deep tissue optical property studies [21,22] and its configuration was described in detail elsewhere [23]. In this study, we constructed a SSFD DRS system consisting of two sub-systems; 1) a bank of frequency modulated laser diodes and avalanche photodiode as shown in Fig. 1(a), and 2) a steady-state light source and spectrometer as illustrated in Fig. 1(b). The source and detector optical fibers were multimode fibers with 400-μm core diameter and 0.22 numerical aperture, and their center-to-center distance was 1 mm. The SS sub-system equipped with a tungsten halogen light source (HL2000, Ocean Optics, FL) and spectrometer (QE65000, Ocean Optics, FL) to collect wavelength dependent reflectance of turbid samples. The FD apparatus employed four laser diodes with wavelengths of 660 (120mW, Thorlabs, NJ), 780 (150mW, Qphotonics, CA), 830 (150mW, Qphotonics, CA), and 980 nm (15 mW, NCKU, Taiwan) which were mounted on laser diode mounts (LDM9T, Thorlabs, NJ) and driven by a laser diode controller (LD3908, Newport, CA). In addition, the laser diode on each mount was sinusoidally modulated by a RF signal generated from a network analyzer (N5230C, Agilent, CA) at frequencies ranging from 50 to 400MHz. The optical power out of the 400 μm diameter optical fiber was 26 mW at 660 nm, 9 mW at 780 nm, 55 mW at 830 nm, and 8 mW at 980 nm. The frequency domain diffuse reflectance was detected by an avalanche photo diode (C5658, Hamamatsu, Japan). The detected signal was transmitted to the network analyzer to determine the phase delay and amplitude demodulation introduced by the samples. For reasons still under investigation, all four of our piezoelectric element based optical switches would interfere the FD light signal in a way that the detected FD reflectance (both amplitude and phase) would resonant over the whole frequency range. Thus, the FD and SS light sources were manually switched to the same source fiber.
In general, the amplitude demodulation decreases parabolically with modulation frequencies and the phase delay increases linearly with modulation frequencies [20]. The measured data were fit to an appropriate photon transport model to obtain absorption and reduced scattering coefficients for each laser wavelength. The broadband sample diffuse reflectance obtained from the SS sub-system was calibrated using a spectrum acquired from the integrating sphere, and the optical properties derived from the FD sub-system were used to scale the calibrated broadband diffuse reflectance. By fitting the reduced scattering coefficient to a scattering power law, the full spectrum reduced scattering coefficients in the region from 600 to 1000nm can be obtained [23]. Once the reduced scattering coefficients at the wavelength region of interest were defined, the absorption coefficient spectrum could subsequently be calculated.

GPU-based Monte Carlo model and artificial neural network
DRS is a model based technique and usually need to incorporate a photon transport model to relate the measured diffuse reflectance to the sample optical properties. In a typical sample optical property recovery procedure, a photon transport model is iteratively used to generate diffuse reflectances for various sample optical property combinations to find a best fit to the measured diffuse reflectance, and thus the optical properties of the sample can be determined. An important prerequisite for ensuring the robustness of this procedure is that the model relating optical properties to diffuse reflectances possesses one-to-one correspondence.
The photon transport model employed in this study was adopted from the GPU-based (Graphic Processing Unit) Monte Carlo model developed by Alerstam et al. [24]. In the code, the diameters of source and detector were both 400 μm. The radius and thickness of the sample were both set to 10 6 cm to simulate the semi-infinite sample geometry. The refractive indices of sample and source/detector were set to 1.4 and 1.45, respectively. Photon packets propagating in the medium followed the Henyey-Greenstein phase function with the anisotropy factor g set to 0.8. In the Monte Carlo model, the sample was assumed to be homogeneous and the sample optical properties were varied from 0.01 to 0.30 mm −1 in a 0.01 mm −1 step for μ a and varied from 0.55 to 2.5 mm −1 in a 0.05 mm −1 step for μ s '. Each simulation ran continuously until a total number of 20,000 photon packets arrived at the detector to ensure the standard deviations of amplitude and phase were within 0.6% and 0.01°, respectively, in the frequency range from 10 to 500 MHz. A typical simulation for SDSs of 1 or 3 mm would take about a few seconds on a computer equipped with an NVIDIA GTX 780-Ti graphic card.
Although the computation speed of a GPU-based Monte Carlo model is usually faster than a conventional Monte Carlo model by about three orders, using such a model for real time (e.g., less than 1 second) sample optical property recovery would require tremendous computing resources. To achieve real time sample optical property recovery at 1 mm SDS with commonly available computers, we constructed artificial neural networks (ANNs), which were trained by databases composed of numerous Monte Carlo simulation results, to efficiently and accurately determine the FD diffuse reflectance. The essential steps for constructing ANN models for a spatially resolved SS DRS system were elaborated in our previous paper [25]; specifically, we demonstrated the use of ANNs to accurately compute the sample optical property spectra from the SS diffuse reflectances measured at SDSs of 1 and 2 mm. In this study, we aimed to obtain two ANNs for the FD sub-system, and each of them could take sample optical properties as input parameters and generate amplitude #259000 demodulation or phase delay at various modulation frequencies (e.g., 100 frequencies) in the range from 10 to 500 MHz. Directly constructing such 2-in-100-out ANNs is usually time consuming or in some cases an impossible task. We discovered that a practical, robust way to produce ANNs of required function was to first preprocess the training databases to reduce the dimension of output parameters. We found that the frequency domain amplitude demodulation and phase delay versus modulation frequencies could be modelled by a second order polynomial y = ax 2 + b and a linear function y = cx, respectively. Here y and x represent amplitude/phase and modulation frequency, respectively. We could fit the databases to the polynomial functions to obtain the three coefficients (a, b, and c), and thus the dimension of the amplitude and phase databases could be reduced to 2-in-2-out and 2-in-1-out, respectively. The total database construction time was around 4-5 hours and a typical successful ANN training would require 1-2 hours. In addition, for the SS sub-system, an ANN that is capable of taking reflectance and reduced scattering coefficient as input parameters and returning the absorption coefficient of the sample is necessary. We generated a 2-in-1-out database using the Monte Carlo method for training such an ANN for the SS subsystem.
The databases were then sent to the MATLAB (Mathworks, MA) training function "trainbr" to construct the ANN models. The network training function "trainbr" utilizes the Bayesian regularization process to minimize a combination of squared errors and weights so as to produce a network that generalizes well. We used one hidden layer for all ANN trainings, and the number of neurons was adjusted between 35 and 45 to optimize the ANN performance. The successfully trained ANNs (one amplitude-ANN and one phase-ANN) for the FD sub-system can work in conjunction with a nonlinear optimization algorithm, such as the "lsqcurvefit" function in MATLAB, to determine μ a and μ s ' of the sample from the FD diffuse reflectance. The nonlinear least-squares solver "lsqcurvefit" was employed in this study to find a best (μ a , μ s ') combination which resulted in a minimum deviation between the ANN-derived FD diffuse reflectance and the measured one. A full wavelength region μ s ' spectrum is obtained by fitting the reduced scattering coefficients derived from the FD subsystem. The μ s ' spectrum and the calibrated steady state reflectance are then sent to the ANN for the SS sub-system to determine the absorption spectrum in the full wavelength region.

Silicone phantoms
By mixing the two-part silicone with different proportions of absorbing (black India ink) and scattering (TiO2 powder) agents, four homogeneous silicone phantoms of different optical properties were created. The benchmark optical properties of these phantoms were determined by inverse adding-doubling method from 650 to 1000 nm [26], and the values at 660 nm and 830 nm are listed in Table 1. Phantom-1 and phantom-2 were used to mimic low absorption samples, such as light skin or deep tissue [23,25,[27][28][29], while phantom-3 and phantom-4 had optical properties similar to those of high absorption samples such as dark skin [25,27,28]. All silicone phantoms had diameter of 7.5 cm and thickness of 4.5 cm.  Phantom-1 and phantom-3 were employed as the calibration phantoms to determine the system response. The process of calibration is elaborated in the following. First, the GPUbased Monte Carlo program was utilized to calculate the theoretical FD reflectances of the calibration phantoms at laser diode wavelengths. Next, the experimental FD reflectances of the calibration phantoms were divided by the theoretical amplitude demodulation and subtracted from the theoretical phase delay to obtain the system response. The experimentally measured FD reflectance of target phantoms was assumed to be influenced by a system response that was the same as the one obtained in the previous step. Thus the experimentally measured reflectance could be calibrated to be independent of the system response. The calibrated FD diffuse reflectance was then sent to the nonlinear optimization routine to recover the sample optical properties.

Performance of the standard diffusion model at SDS of 1mm
To evaluate the performance of the SDE at 1 mm SDS, we calculated the percent deviation of the SDE derived amplitude demodulation from those obtained from the Monte Carlo model at 100 MHz and the results are displayed in Fig. 2(a). The FD SDE used in this study was adopted from the derivation provided by Haskell et al. [7], and the sample refractive index was set to 1.4. It can be seen in the figure that the amplitude error was within 150% for SDS of 1 mm. The highest SDE amplitude errors occur in the region where samples have reduced scattering coefficient lower than 1 mm −1 and absorption coefficient higher than 0.2 mm −1 . For typical biological tissues, our results suggest that the SDE is especially not suitable for predicting the amplitude demodulation introduced by samples at wavelengths higher than 1000 nm where the reduced scattering property is generally lower than 1 mm −1 and the absorption property is higher than 0.1 mm −1 due to strong water absorption above 1000 nm. On the other hand, the difference in degree between the phase predicted by the SDE and the Monte Carlo model at SDS of 1 mm is displayed in Fig. 2(b). It can be seen in Fig. 2(b) that the SDE phase error is in general less than 0.2°. The highest SDE phase errors occur in the region where samples have reduced scattering coefficient lower than 1.3 mm -1 and absorption coefficient lower than 0.1 mm -1 . In summary, the results demonstrated in this subsection show that the SDE is not capable of accurately predicting FD diffuse reflectance at 1 mm SDS. In addition, we analyzed the amplitude and phase errors of SDE at 3 mm SDS and found that for typical biological tissue optical properties, the amplitude and phase errors of SDE were within 20% and 0.05°, respectively. This result verified that the accuracy SDE was improved as the SDS increased. The experimental performance of the SDE will be studied in Section 3.3.

Accuracy of the ANN in predicting frequency-domain reflectance at SDS of 1 mm
The accuracy of the amplitude of FD reflectance predicted by the amplitude-ANN was evaluated by calculating the percent difference between the amplitude values determined by the ANN and the Monte Carlo method. The results for 100 MHz modulation frequency are shown in Fig. 3(a) for SDS of 1 mm. It can be seen in the figure that the ANN was well trained and had amplitude errors less than 2.5%. We also checked the ANN's amplitude prediction errors at 200, 300, and 400 MHz, and their maximum values were all less than 3% (data not shown). Our results imply that using a second order polynomial to approximate the amplitude of FD reflectance at various frequency is a viable method for biological tissue relevant sample optical property combinations. The difference in terms of degree between the 100 MHz FD phase delays predicted by the phase-ANN and the Monte Carlo model at SDS of 1 mm for various sample optical property combinations were calculated and illustrated in Fig. 3(b). We found in Fig. 3(b) that the phase errors were less than 0.01° for the majority of sample optical property combinations, except for some samples with lowest absorption (μ a = 0.02 mm −1 ) at which the phase errors could reach 0.06°. It should be noted that all phase errors plotted in Fig. 3(b) were all less than 1% of their corresponding total phase delay. We expect that the phase errors at the low absorption boundary could be further reduced by training the ANNs with an extended database range covering absorption coefficients lower than 0.01 mm −1 . From the data shown in Figs. 3(a) and 3(b), it can be concluded that the amplitude-ANN and phase-ANN for SDS of 1 mm were successfully trained and had accuracy comparable to the Monte Carlo method. In addition, we also constructed two ANNs for computing FD amplitude and phase at SDS of 3 mm. Their amplitude and phase prediction errors were within 3% and 0.07°, respectively.

Theoretical sample optical property recovery error induced by measurement uncertainty
As mentioned earlier, modeling and system limitations impeded the extensive use of FD DRS at short SDSs. The ANN construction method proposed in this paper enables one to obtain a model capable of instantly generating accurate FD diffuse reflectance at SDS of 1 mm; this provides a solution to the modeling issue. As for the system limitation, it has been reported inevitable measurement uncertainty of the FD apparatus would affect the sample optical property recovery accuracy [13,20,30]. To understand the inherent limitation of FD DRS at SDS of 1 mm, it is important to investigate the influence of the measurement uncertainty on the optical property recovery error.
Measurement uncertainty for the FD method is a result of multiple factors, such as laser diode performance, sample motion artifact, and the apparatus design for resolving phase delay and amplitude demodulation, and thus is system dependent. For example, Pham et al., Boas et al., and Kienle et al. defined the phase and amplitude uncertainties to be (0.5°, 3.0%), (0.08°, 0.3%), and (0.1°, 0.1%), respectively, for their FD DRS systems [13,20,30]. We monitored the long term stability of the four laser diodes in our FD DRS system and found that the amplitude and phase fluctuation were within 1.89% and 0.16° for modulation frequencies from 50 to 400 MHz in a 1-hour time frame. Thus, we defined that any two sets of FD diffuse reflectance that had differences within 1.89% and 0.16° for amplitude and phase, respective, were not distinguishable and therefore not unique. We exhaustively examined the uniqueness of all diffuse reflectance sets in the database, and calculated the largest optical property variation between the indistinguishable diffuse reflectance sets, and the results are depicted in Fig. 4 for SDS of 3 mm. The optical property variations shown in Fig. 4 represent the optical property recovery uncertainty or the theoretically largest optical property recovery error at a certain sample optical property combination introduced by the measurement uncertainties in amplitude and phase. We can observe in Fig. 4(a) that the theoretical recovery errors of μ a can reach nearly 100% in the low μ a , low μ s ' region. On the other hand, the theoretical μ s ' recovery errors are close to 200% for samples have μ s ' lower than 0.7 mm −1 . Nevertheless, in the typical human skin optical property variation range in the wavelength region from 600 to 1000 nm outlined in Fig. 4 [25,27,28], the theoretical optical property recovery errors are reduced to within 40% for μ a and 70% for μ s '. These are the worst error values if the 100 MHz FD reflectance measured with our system were used to recovery sample optical properties at 3 mm SDS. It should be noted that for light complexion subjects, the skin absorption coefficient is typically close to 0.05 mm −1 in the 600-1000 nm range, while for dark complexion subjects, the absorption could reach as high as 0.5 mm −1 in the same wavelength region [28,31]. We performed similar calculations to determine the theoretically largest optical property recovery errors for our system at SDS of 1 mm, and the results are illustrated in Fig. 5. Compared to Fig. 4, the low error region marked in red color is greatly reduced in Fig. 5. This was because the fraction of the phase measurement uncertainty to the total phase delay increased as the SDS decreased from 3 mm to 1 mm. Further, in Fig. 5, the maximum values are located within the outlined area indicating typical human skin optical property variation range. This suggests that FD DRS carried out at single frequency 100 MHz, 1mm SDS is not suitable for skin applications.  Fig. 6 are reduced roughly from 200% to 65%. This result implies that it is advantageous to employ higher modulation frequencies for sample optical property determination. In addition, although in Fig. 6 the maximum values are still located within the typical human skin variation range, they are shifted to the higher absorption and higher reduced scattering region in contrast to Fig. 5.

Benefit of simultaneously using multiple frequencies for sample optical property recovery
It can be seen from the results shown in the previous section that accurately recovering sample optical properties from the FD reflectance measured at one modulation frequency is not feasible at SDSs shorter than 3 mm due to measurement uncertainties. We carried out one run of Monte Carlo simulation to generate the FD reflectance at 1 mm SDS, 100 MHz for the sample optical property combination (μ a = 0.1 mm −1 and μ s ' = 1.2 mm −1 ) where the maximum value in Fig. 5 located. The FD reflectance was then repeatedly added with 1.89% amplitude and 0.16° phase random noises 50 times to produce 50 sets of FD reflectance. We utilized the ANNs discussed in section 3.2.1 and the MATLAB least-square optimization function "lsqcurvefit" to fit the simulated data. The sample optical property recovery percent error for each set of FD reflectance was determined. It was found that the 50 μ a recovery errors had an average value of 1.820% with the standard deviation of 64.352%, and the 50 μ s ' recovery errors were 6.560% in average with the standard deviation of 35.455%. Our simulation results discussed here support the data shown in Fig. 5. In addition, they verified that accurately recovering sample optical properties from the FD reflectance measured at one modulation frequency is not feasible at 1 mm SDS owing to inevitable measurement uncertainties.
We conducted a theoretical investigation to understand whether the measurement uncertainty issue could be alleviated by simultaneous employing FD reflectances measured at multiple frequencies for the sample optical property recovery. Again, we used the Monte Carlo method to obtain one 50 to 400 MHz FD reflectance of a sample with (μ a = 0.1 mm −1 and μ s ' = 1.2 mm −1 ) at 1 mm SDS. The FD reflectance was added with 50 sets of random noise of the same level mentioned earlier at all frequencies to generate 50 sets of test data. In the inverse calculation, the total number of modulation frequencies was varied, but the lowest and the highest modulation frequencies were kept at 50 and 400 MHz, respectively. The average and standard deviation of 50 optical property recovery errors were calculated and the results are listed in Table 2. It can be noted that the standard deviation of optical property recovery errors of the 2-frequency group is much less than that of the 1-frequency group discussed earlier. Further, we can see from the table that the standard deviation of optical property recovery errors gradually reduced as more data measured at different modulation frequencies were employed in the inverse fitting. The merit of employing multiple frequency data for alleviating the measurement uncertainty's influence on the optical properties recovery accuracy at 1 mm SDS can be clearly seen here.  Similar theoretical investigation was carried out for SDS = 3 mm, and it was found that the measurement uncertainty issue only slightly affected the recovery results for SDS = 3mm. When using 2-frequency data at 50 and 400 MHz for data recovery, the maximum optical property recovery errors were 6.0% for μ a and 4.9% μ s '. Our results imply that when carrying out FD DRS at SDSs shorter than 3 mm, using FD reflectances measured at multiple frequencies for sample optical property recovery can alleviate the influence of the measurement uncertainties on the accuracy of results.
Moreover, to investigate the influence of the upper bound of modulation frequency on the recovery accuracy, we kept the total frequency number and the lower bound of modulation frequency at 36 and 50 MHz, respectively, and gradually decreased the upper bound of modulation frequency from 400 to 100 MHz. We found that to keep the standard deviation of recovery errors less than ± 10% for both μ a and μ s ', the upper bound of the modulation frequency should be at least 150 MHz.

FD DRS at 1mm and 3mm SDSs
The experimental performance of the ANNs constructed for 3 mm SDS was evaluated and compared to that of the SDE. Typical 50-150 MHz, 830 nm data fitting outcomes of phantom-2 by using the ANNs and the SDE are illustrated in Fig. 7. System response calibration was attained by using phantom-1 as the calibration standard. The signal level of the data in the frequency range from 150 to 400 MHz reached the noise level of the photo detector and were not included in the analyses. It can be observed from Fig. 7 that the fitting quality using both models were decent, but the amplitude and phase predicted by the two models were different. The differences were 5.6% between the average amplitude values and 0.1 degree between the phase values at 150 MHz shown in Figs. 7(a) and 7(b). This difference is understandable since in the process of calibration, system response was determined from a phantom with known optical properties. The theoretical reflectance of the phantom determined using the ANN or the SDE would be different and thus the system response would be not the same when different theoretical model was employed.
The optical property recovery results for all four laser diodes are listed in Table 3. The maximum optical property recovery errors were 14.74% for μ a and 5.37% for μ s ' for the ANN model, and were 18.61% for μ a and 7.35% for μ s ' for the SDE. Overall, the ANNs and the SDE had comparable performance in recovering deep tissue relevant sample optical properties at SDS of 3 mm.  Further, we fit the 1 mm SDS, 50-400 MHz FD reflectance data of phantom-2 to the corresponding ANNs and SDE. A typical data fitting outcome for the 830 nm laser wavelength are illustrated in Fig. 8. Since the signal attenuation introduced by the sample at SDS of 1 mm was much less than that at 3 mm SDS, the upper bound of the data fitting frequency was extended to 400 MHz. Compared the 150 MHz data shown in Figs. 7 and 8, we noted that the amplitude and phase varied by about ten-fold and 1.5°, respectively, when the SDS was reduced from 3 to 1 mm. The phase uncertainty of 0.16° in our system played a significant role in the 1 mm SDS data recovery and we found that employing multiple frequencies rather than single frequency in the data fitting would produce much better fitting outcomes. The optical property recovery results for all four laser diodes at SDS of 1 mm are listed in Table 4. The maximum recovery error was within 9.97% when using the ANNs while the recovery error could reach 21.94% for the SDE. From Table 4, the advantage of using the ANN method for sample optical property recovery is apparent.  Since the optical properties of superficial tissues generally are more absorbing than deep tissues, it is thus important to evaluate the performance of the 1 mm SDS ANNs for high absorption samples. We employed phantom-4 that had optical properties similar to those of human dark skin for testing the performance of the ANNs and the SDE at SDS of 1 mm. Here, phantom-3 served as the calibration phantom. The optical property recovery results at four laser diode wavelengths are listed in Table 5. It can be observed that the performance difference between the two models are remarkable; the maximum recovery errors were within 5.07% for the ANNs and 65.49% for the SDE. The results presented in this section accentuate the need for an advanced model as well as multi-frequency measurements when using FD DRS at short SDSs for superficial tissue applications. As mentioned in section 2.1, one could use the scattering power law to fit the reduced scattering coefficients derived from the FD measurement conducted at a certain SDS to obtain the sample broadband reduced scattering spectrum. Using the broadband reduced scattering spectrum along with the calibrated steady state reflectance measured at the same SDS, the broadband absorption spectrum of sample can be determined with an appropriate model. In this subsection, we will present the phantom broadband optical property spectra recovery results using our SSFD DRS system. The broadband optical property spectra of phantom-4 recovered using the ANN for the SS DRS apparatus at SDS of 1 mm are depicted in Fig. 9 as solid thick lines. In addition, the dotted lines illustrated in Fig. 9 were obtained from the inverse adding-doubling method [26] and were regarded as the benchmark values. It can be seen that the optical property spectra derived from the ANNs agree well with the benchmark values. The slight deviation (maximum of 5.5% at 1000 nm) of the SSFD DRS derived absorption spectrum from the benchmark values was caused by the small disagreement between the reduced scattering spectrum predicted by the scattering power law and that derived from the spatially resolved DRS.
Besides, we also employed the SDE to calculate the sample optical property spectra at the same SDS, and the results are illustrated in Fig. 9 as solid thin lines. Overall, the spectra derived from the SDE did not follow the same trend with the benchmark spectra and their differences reached 39.7% for μ a and 32.7% for μ s '. Fig. 9. Broadband optical property spectra (a) μ a and (b) μ s ' of phantom-4 determined using the ANN for the SS DRS apparatus (solid thick lines) and the SDE (solid thin lines) at 1 mm SDS. Squares and circles represent the coefficients obtained from the FD DRS method using ANN and SDE, respectively. Benchmark spectra are depicted as dotted lines.

Conclusion
FD DRS is a mature platform and its characteristics and the potential for deep tissue applications have been carefully studied. This platform has been employed by many groups in bedside studies, and its miniaturization has been demonstrated. Due to system and modelling constraints, using FD DRS in superficial tissue studies has rarely been shown. In this study, we found that the key to overcome the influence of measurement uncertainties on the sample optical property recovery accuracy at short SDSs is to perform FD measurements at multiple frequencies. Moreover, we revealed the construction method of ANNs that could relate the sample optical properties and the FD diffuse reflectance at SDS near 1 mm. We constructed ANNs for computing FD reflectance at 1 mm SDS and found them not only very efficient but also had precision level comparable to that of the Monte Carlo method. The ANNs could work in conjunction with a least-square optimization routine to recover the optical property of turbid samples at discrete laser diode wavelengths. In addition, it was shown that broadband absorption and reduced scattering spectra of turbid samples in the wavelength range from 650 to 1000 nm could be derived by incorporating the SS apparatus into the FD DRS system with 1 mm SDS. Currently, the primary use of FD DRS systems is in deep tissue applications, and our findings disclosed in this study provide a practical means to extend their capability in superficial tissue investigations.