Inverse calculation method for obtaining multi-optical parameter spectral data of turbid medium in real time

Abstract

The invention discloses an inverse calculation method for obtaining multi-optical parameter spectral data of a turbid medium in real time. Scattered light signal measurement values, namely diffuse reflection rate Rd, diffuse transmission rate Td and forward transmission rate Tf, of the turbid medium sample can be accurately and rapidly obtained through one-time measurement at each measurement wavelength center value lambda. The method is simple and easy to operate, collimation transmittance signals attenuated along with the sample thickness D index do not need to be measured, sample attenuation coefficients do not need to be solved, the thickness D value of the prepared turbid medium sample can be two times or more of the thickness allowed by an existing method, therefore, the difficulty and time for preparing the sample with the small thickness D value can be remarkably reduced, and turbid medium multi-parameter spectral data can be obtained in real time. According to the inverse calculation method, a sample signal measurement system and sample preparation steps can be remarkably simplified, and turbid medium multi-parameter spectral data consisting of an absorption coefficient mu a, a scattering coefficient mu s and an anisotropy coefficient g can be accurately and quickly obtained on the basis of accurately and quickly measuring scattered light signal spectral data.

Description

Technical Field

[0001] The present invention relates to a method for obtaining optical parameter spectrum data of a turbid medium, and in particular to an inverse calculation method for obtaining multi-optical parameter spectrum data of a turbid medium in real time. Background Art

[0002] Turbid medium materials have strong light scattering properties, covering many materials including biological tissues, such as milk, wine, synthetic materials, cell samples, blood, skin, environmental sewage and polluted gases. The measurement and research of characteristic parameter spectra of turbid media has wide applications in the fields of material research, food inspection, environmental monitoring and clinical diagnosis, and is one of the major basic scientific issues in the field of optics. Turbid medium materials have particles or non-uniform refractive index structures (micrometer scale) inside, showing strong light scattering properties, and the differences in the molecular structure of the materials (nanometer scale) cause large differences in absorption characteristics in the spectral region, so their spectra contain rich material characteristic information.

[0003] Accurately measuring multiple parameters (hereinafter referred to as multi-parameters) that can describe the optical properties of turbid medium materials and their spectral data is an implicit light scattering inverse problem, and its solution needs to be based on an accurate description of the light scattering mechanism and energy transfer process in the turbid medium material. At present, it is believed that the model that can accurately describe the interaction between light and turbid medium samples whose scale is much larger than the wavelength is the radiation transfer theory. According to this theory, the propagation of light energy in turbid medium materials can be described by the following steady-state, single-energy, passive radiation transfer equation:

[0004]

[0005] Where L(r,s) is the radiance of light at position r, along the direction of unit vector s, μ a is the material absorption coefficient, μ s is the material scattering coefficient, p(s,s') is the material scattering phase function, and dΩ' is the solid angle integral differential in the direction of the unit vector s'. The multiple optical parameters μ of the above turbid medium material a , μ s and p(s, s') are functions of wavelength λ. Compared with other approximate theoretical models such as Kubelka-Munk approximation and light diffusion model, the radiation transmission theory has a higher accuracy in describing the light transmission process and energy distribution in turbid medium materials. The present invention simulates and calculates the light transmission process and energy distribution in turbid medium materials based on the radiation transmission theory.

[0006] The scattering phase function p(s,s') in the above formula (1) describes the distribution of light energy scattered from the s' direction to the s direction. Its accurate measurement requires measuring the scattered light intensity at multiple scattering angles, which generally requires special instruments and a long time, so it is difficult to achieve within a larger wavelength range. Henyey and Greenstein defined the HG function in 1941, namely the Haney-Greenstein scattering phase function p HG (cosθ) can be used to approximately replace the actual scattering phase function p(s,s') of many different turbid medium materials, that is, p(s,s') is approximately equal to p HG (cosθ), where θ is the polar angle between the scattered light direction s' and the incident light direction s. The form of the HG function is completely determined by the parameter g, and the definition of g is exactly the average value of cosθ calculated based on the light scattering process, which is usually defined as the anisotropy coefficient. Since the HG function is simple in form and easy to analyze, it can approximately replace the actual scattering phase function of turbid medium materials. It has been widely used in many studies and has been proven to approximately represent the angular area average distribution of scattering phase functions of many different types of turbid medium materials. The present invention uses the HG function to replace the scattering phase function p(s,s') of the turbid medium sample material, and adopts a simulation model based on the radiation transmission theory to describe the light transmission process in the turbid medium material and calculate the scattered light signal. Assuming that the refractive index n defined by the Fresnel formula is known, the measured optical parameter of the turbid medium material is μ a , μ s and g. It should be noted that μ a , μ s and g are functions of wavelength λ, so μ a , μ s The spectral data of and g can fully describe the optical properties of turbid medium materials under the framework of radiation transfer theory.

[0007] Obtain three optical parameters of turbid medium samples, namely μ a , μ s The spectral data of g generally need to iteratively solve the inverse light scattering problem based on equation (1) at multiple wavelengths, and finally make the light scattering inverse problem based on equation (1) and the assumed μ at each wavelength. a , μ s The calculated values ​​of scattered light obtained from the values ​​of g are consistent with the corresponding measured values ​​within the experimental error range.

[0008] Except for a few boundary value problems where the sample boundary and the light source-measurement sensor geometry are simple shapes (such as a cuboid or a sphere), boundary value problems based on the radiation transfer equation (1) generally cannot obtain analytical solutions or closed-form solutions, and numerical or statistical solution methods such as the finite element method, adding-doubling method, and discrete coordinate method are required. However, these numerical solution algorithms and programming are relatively complex, and any changes in the sample boundary or the light source-measurement sensor geometry require program modifications, which is a lot of work. In view of these problems, many researchers have turned to developing statistical solution methods with relatively simple algorithms, the most representative of which is the Monte Carlo method for solving boundary value problems in radiation transfer theory.

[0009] The measurement methods of multi-parameters and spectral data of turbid medium materials reported in the literature all require solving the inverse problem of light scattering in multiple steps, including methods of measuring scattered light signals using an integrating sphere (see, for example, C. Chen, J. Q. Lu, H. Ding, K. M. Jacobs, Y. Du, and X. H. Hu, "A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm", Optics Express, vol. 14, pp. 7420-7435 (2006) and other references cited in the article) and methods of measuring scattered light signals without using an integrating sphere (see, for example, X. Liang, M. Li, J. Q. Lu, C. Huang, Y. Feng, Y. Sa, J. Ding, X. H. Hu "Spectrophotometric determination of turbid optical parameters without using an integrating sphere", Applied Optics, vol. 55, pp. 2079-2085 (2016) and other references cited in the article). These methods generally require the preparation of samples with a smaller thickness D; otherwise, the measured sample collimated transmittance T c The value will produce a large error due to the small signal-to-noise ratio. c With the D value, the attenuation coefficient μ of the material can be quickly calculated t =μ a +μ s ; Then based on the measured values ​​of other scattered light signals, such as diffuse reflectance R dand full transmittance T d , other optical parameters such as scattering albedo a=μ are solved by inverse calculation method based on physical relationship s / μ t and g, thus obtaining the three optical parameter values ​​μ a , μ s and g. Since the above method requires the measurement of T c The value varies with the attenuation coefficient μ of the sample t and thickness D decreases exponentially, so the thickness D of the prepared sample must be small enough to obtain a T that meets the signal-to-noise ratio requirement. c When the D value is very small, such as less than 0.5 mm, the difficulty of sample preparation increases, the D value measurement error increases, and the error of the obtained optical parameters increases. Summary of the invention

[0010] The technical problem to be solved by the present invention is to overcome the shortcomings of the prior art and provide an inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time, which can accurately measure the spectral data of turbid medium samples with a thickness that is twice or more times the thickness allowed by the existing method.

[0011] The technical solution adopted by the present invention is: an inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time, comprising the following steps:

[0012] 1) obtaining a turbid medium sample and obtaining physical parameters of a system for measuring multi-parameter spectral measurement data of the turbid medium, wherein the physical parameters include: shape, size and refractive index of the sample and the sample dish, shape, size, distance and orientation of the detector relative to the coordinate origin defined by the reference point and reference direction on the surface of the sample dish and the coordinate axis;

[0013] 2) obtaining the shape, area, power distribution, incident angle, measurement values ​​of the center value λ of the measured wavelength and the bandwidth of the incident light beam, as well as the measurement value of the light intensity of the incident light beam, and causing the incident light beam to be incident on the sample;

[0014] 3) Measure the scattered light signals emitted from the sample and the sample dish at multiple angles relative to the sample dish, and obtain the scattered light signal measurement value at the center value λ of the measurement wavelength, including the diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f ;

[0015] 4) Three optical parameters of turbid medium samples at different wavelengths, namely, absorption coefficient μ a , scattering coefficient μ sThe anisotropy coefficient g is represented by the optical parameter vector point P, which is denoted as P(λ)=(μ a ,μ s ,g), where (μ a ,μ s ,g) is called the sample optical parameter space; the initial preset value of the sample optical parameter vector is set to P0(λ) = (μ a0 ,μ s0 ,g0), where μ a0 is the initial absorption coefficient of the sample, μ s0 is the initial scattering coefficient of sample 2, g0 is the initial anisotropy coefficient of the sample, and calculate the covariance m of P(λ) in the sample optical parameter space P -1 adjacent point, take m P =7, there are six adjacent optical parameter vector points P, which are (μ a0 ±s a ,μ s0 ,g0),(μ a0 ,μ s0 ±s s ,g0),(μ a0 ,μ s0 ,g0±s g ), where s a 、s s and g They are the parameter distances of adjacent points, and their values ​​are s a =0.01μ a0 、s s =0.01μ s0 and g = 0.01g0, or m P Take a value between 5 and 27, sa, ss and sg take values ​​between 1 and 10 times of the corresponding values ​​respectively; m including the initial preset value P(λ) of the optical parameter vector P Different optical parameter vector points are introduced into the Monte Carlo ray tracing method based on radiation transfer theory to obtain m P A group of different scattered light signal calculation values, each group of scattered light signal calculation values ​​includes a diffuse reflectance calculation value R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc ;

[0016] 5) According to the scattered light signal measurement value of the sample: diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f , and the calculated value of the scattered light signal of the sample: the calculated value of diffuse reflectance R dc , Diffuse transmittance calculated value Tdc and the calculated forward transmittance T fc , find the first optimization objective function ρ at the optical parameter vector point P δ (P) and the second optimization objective function ρ 1 / δ (P);

[0017] 6) Determine whether the first optimization objective function ρ δ (P) value is greater than the set objective function threshold ρ th , or the second optimization objective function ρ 1 / δ (P) value is less than the set objective function threshold ρ th , if yes, go to step 7); otherwise go to step 8);

[0018] 7) Search and update the sample optical parameter vector point P, and then return to step 4);

[0019] 8) Save the best position P of all particles at present g , P g That is, the solution of the optical parameter vector of the turbid medium sample at the wavelength center value λ: the absorption coefficient μ a , scattering coefficient μ s and anisotropy coefficient g, complete the solution of the inverse scattering problem with the measurement wavelength center value λ;

[0020] 9) Select a new measurement wavelength center value λ and return to step 2) until the absorption coefficient μ of the turbid medium sample of all wavelength center values ​​within the measurement wavelength range is completed. a , scattering coefficient μ s The inverse calculation of the anisotropy coefficient g is then used to obtain the multi-optical parameter spectral data of the turbid medium.

[0021] The inverse calculation method for obtaining the multi-optical parameter spectral data of turbid medium in real time of the present invention can use a device without integration and summation and without measuring the attenuation coefficient for the required multiple scattered light signal measurement values, so the absorption coefficient μ of the turbid medium sample with a thickness of 2 times or more times the thickness allowed by the existing method can be accurately measured. a , scattering coefficient μ s The inverse calculation method of the present invention can significantly simplify the sample signal measurement system and sample preparation steps, and accurately and quickly measure the spectral data of the scattered light signal, and accurately and in real time obtain the spectral data of the absorption coefficient μ. a , scattering coefficient μ s and anisotropy coefficient g composed of multi-optical parameter spectral data of turbid media. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a schematic diagram of the system structure for measuring multi-parameter spectral measurement data of turbid media in the present invention;

[0023] Figure 2 is a schematic diagram of samples and detectors in an embodiment of the present invention;

[0024] Figure 3 It is a flow chart of the Monte Carlo ray tracing method based on the radiation transfer theory of the present invention;

[0025] Figure 4 is a flow chart of the inverse calculation program of the present invention;

[0026] Figure 5a is the p of the present invention δ (P) and ρ 1 / δ Statistical fluctuation error curve of (P) and δ(P);

[0027] Figure 5b is the Δ(ρ δ ) and Δ(ρ 1 / δ ) and q m The relationship curve diagram of

[0028] Figure 6 The spectral data of the scattered light signal measurement value and the error distribution diagram of the fat emulsion with a concentration of 20% of the present invention;

[0029] Figure 7 It is a spectral data diagram of scattered light signal of the fat emulsion with a concentration of 20% of the present invention.

[0030] in:

[0031] 1: Light source unit 2: Modulator

[0032] 3: Spectrum splitter 4: Incident light intensity detection device

[0033] 5: Light path reflector 6: Sample dish

[0034] 7: Sample 8: Incident beam

[0035] 9: First photodetector 10: Second photodetector

[0036] 11: The third photodetector 12: Signal processing unit

[0037] 13: Control computing unit 14: Signal and control cable

[0038] 15: Diffuse reflected light 16: Diffuse transmitted light

[0039] 17: Forward scattered light 18: Coordinate origin DETAILED DESCRIPTION

[0040] The inverse calculation method for obtaining multi-optical parameter spectral data of a turbid medium in real time according to the present invention is described in detail below in conjunction with the embodiments and the accompanying drawings.

[0041] The inverse calculation method for obtaining the spectral data of multiple optical parameters of a turbid medium in real time of the present invention is a method that does not require an integrating sphere device and does not require the measurement of the attenuation coefficient. It only requires the preparation of a sample with a thickness that is 2 times or more than the thickness allowed by the existing method to accurately measure the spectral data of the scattered light signal of the turbid medium sample and obtain the three optical parameters of the sample, namely the absorption coefficient μ by the inverse calculation method of the present invention. a , scattering coefficient μ s and anisotropy coefficient g and its spectral data.

[0042] For the convenience of expression, the absorption coefficient μ of turbid medium samples at different wavelengths is as follows: a , scattering coefficient μ s The anisotropy coefficient g is represented by the optical parameter vector point P, that is, P(λ) = (μ a ,μ s ,g), where (μ a ,μ s ,g) is defined as the sample optical parameter space.

[0043] The inverse calculation method for obtaining multi-optical parameter spectral data of a turbid medium in real time of the present invention comprises the following steps:

[0044] 1) Obtaining a turbid medium sample and obtaining physical parameters of a system for measuring multi-parameter spectral measurement data of the turbid medium, wherein the physical parameters include: shape, size and refractive index of the sample and the sample dish, shape, size, distance and orientation of the detector relative to the coordinate origin defined by the reference point and reference direction on the surface of the sample dish and the coordinate axis.

[0045] like Figure 1As shown, the system for measuring multi-parameter spectral measurement data of turbid media in the present invention includes a sample dish 6 for placing a sample 7, a light source unit 1 arranged in sequence on the optical path, a light modulator 2 and a spectrometer 3 for modulating and splitting the intensity of an incident light beam, an incident light intensity detection device 4 for obtaining a set light beam area and power distribution, and an optical path reflector 5, wherein the reflected light of the optical path reflector 5 constitutes an incident light beam 8 incident on the sample 7, and a first photodetector 9 for receiving diffuse reflected light 15 from the sample 7 and the sample dish 6, a second photodetector 10 for receiving diffuse transmitted light 16 of the sample 7 passing through the sample dish 6, and a forward scattered light detector 9 for receiving the forward scattered light 16 of the sample 7 passing through the sample dish 6. The third photodetector 11 of the incident light 17, the signal output ends of the incident light intensity detection device 4, the first photodetector 9, the second photodetector 10 and the third photodetector 11 are respectively connected to the input end of the signal processing unit 12 for amplifying, filtering and analog-to-digital conversion of the received signal, the output end of the signal processing unit 12 is connected to the signal input end of the control calculation unit 13 for system control and data processing, the signal output end of the control calculation unit 13 is respectively connected to the light source unit 1 and the control input ends of the optical modulator 2 and the spectrometer 3, the distance and direction of the first photodetector 9, the second photodetector 10 and the third photodetector 11 are all based on the coordinate origin 18 as the reference point.

[0046] 2) obtaining the shape, area, power distribution, incident angle, measurement values ​​of the center value λ of the measured wavelength and the bandwidth of the incident light beam, as well as the measurement value of the light intensity of the incident light beam, and causing the incident light beam to be incident on the sample;

[0047] Among them, the shape and area of ​​the incident light beam are smaller than the area of ​​the incident light cross section of the optical window of the sample dish and the sample, such as the shape is elliptical and the area does not exceed 5 square millimeters; the power distribution of the incident light beam must ensure that the power ratio of the same area between the center of the incident light beam and the edge of the incident light beam is less than 104; the incident angle of the incident light beam is 5 to 60 degrees; by controlling the light modulator 2 and the spectrometer 3, the center value of the wavelength, bandwidth and modulation frequency of the incident light beam are obtained, for example, the measurement wavelength center value λ is adjustable in the range of 200nm to 2500nm or part of its wavelength range, the measurement wavelength bandwidth Δλ is adjustable in the range of 0.5nm to 15nm, and the modulation frequency is adjustable in the range of 1Hz to 100MHz; before the incident light beam enters the sample, the incident light intensity signal value I0 is obtained, and then the light beam is incident on the sample in the sample dish.

[0048] 3) Measure the scattered light signals emitted from the sample and the sample dish at multiple angles relative to the sample dish, and obtain the scattered light signal measurement value at the center value λ of the measurement wavelength, including the diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f.like Figure 2 As shown in the figure, diffuse reflected light is the scattered light emitted from the incident light side of the sample after the incident light is scattered by the sample. The measurement of its light intensity needs to avoid the specular reflected light caused by the sample dish. Specular reflected light is the reflected light generated by the sample dish that can be calculated according to the Fresnel formula. Diffuse transmitted light is the transmitted light emitted from the side of the sample opposite to the incident light, excluding the forward transmitted light. Forward transmitted light is the scattered light whose propagation direction is the same or similar to the incident light direction after the incident light is scattered by the sample. It includes:

[0049] When the wavelength center value is measured as λ, the incident light intensity signal value I0 and the diffuse reflection light intensity I measured by the first photodetector 9 are combined. Rd , the diffuse transmission light intensity I measured by the second photodetector 10 Td and the forward scattered light intensity I measured by the third photodetector 11 Tf After amplification, filtering and demodulation, the following formula is used to calculate the scattered light signal measurement value: Diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f :

[0050]

[0051]

[0052]

[0053] 4) Three optical parameters of turbid medium samples at different wavelengths, namely, absorption coefficient μ a , scattering coefficient μ s The anisotropy coefficient g is represented by the optical parameter vector point P, which is denoted as P(λ)=(μ a ,μ s ,g), where (μ a ,μ s ,g) is called the sample optical parameter space. The initial preset value of the sample optical parameter vector is set to P(λ) = (μ a0 ,μ s0 ,g0), where μ a0 is the initial absorption coefficient of the sample, μ s0 is the initial scattering coefficient of the sample, g0 is the initial anisotropy coefficient of the sample, and the covariance m of P(λ) in the sample optical parameter space is calculated. P -1 adjacent point, take m P =7, there are six adjacent optical parameter vector points P, which are (μ a0 ±s a ,μ s0 ,g0),(μa0 ,μ s0 ±s s ,g0),(μ a0 ,μ s0 ,g0±s g ), where s a 、s s and g are the parameter distances of adjacent points, and their values ​​can be taken as s a =0.01μ a0 、s s =0.01μ s0 and g = 0.01g0, or m P Take a value between 5 and 27, sa, ss and sg take values ​​between 1 and 10 times of the corresponding values ​​respectively; m including the initial preset value P(λ) of the optical parameter vector P Different optical parameter vector points are introduced into the Monte Carlo ray tracing method based on radiation transfer theory to obtain m P A group of different scattered light signal calculation values, each group of scattered light signal calculation values ​​includes a diffuse reflectance calculation value R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc ;

[0054] like Figure 4 As shown, the Monte Carlo ray tracing method based on the radiation transmission theory treats the turbid medium sample as a medium containing randomly distributed light absorption centers and light scattering centers. The concentrations of the light absorption centers and light scattering centers are respectively related to the absorption coefficient and scattering coefficient of the turbid medium, and the random distribution of the light absorption centers and light scattering centers is embodied by the calculation of the total distance and travel distance of the photons using a random distribution function. Before the Monte Carlo method calculation begins, it is necessary to input system parameters including the shape, size and refractive index of the sample and the sample dish, the shape, size, position and orientation of the detector relative to the sample dish, and the incident beam parameters including the incident beam shape, area, power distribution, and incident angle, and the number of photons representing the incident beam N0, as well as the preset values ​​of the optical parameters of the sample. Because the Monte Carlo method is a statistical method, its results will contain statistical errors, so the total number of photons N0 that needs to be tracked and calculated must be large enough to reduce the statistical error in the calculation result to a small enough level. However, if N0 is too large, the calculation time will be too long. Generally, N0 is between 10 to the 6th power and 10 to the 10th power. In the present invention, the selected optimization objective function ρ 1 / δ It can effectively reduce the impact of random fluctuations on the results. N0 can be increased from 10 to the 5th power to 10 to the 7th power, which greatly saves calculation time.

[0055] The Monte Carlo ray tracing method based on radiation transfer theory specifically includes:

[0056] (1) For the i-th incident photon, the absorption coefficient μ of the sample is a and random numbers, calculate the total path length of the incident photon in the sample;

[0057] (2) From the scattering coefficient μ of the sample s and random numbers, calculate the length of the next travel distance of the i-th incident photon in the sample;

[0058] (3) Assume that the i-th incident photon is scattered when it travels along the initial travel direction in the sample to the end position of the next travel distance described in step (2), and calculate the end position of the incident photon based on the travel distance of the incident photon;

[0059] (4) Determine whether the i-th incident photon overflows the sample dish at the end position described in step (3): if yes, proceed to step (10); otherwise, proceed to step (5);

[0060] (5) Calculate the total length of the cumulative travel distance of the i-th incident photon in the sample. If the total length of the cumulative travel distance is greater than the total length of the travel distance, determine whether the i-th incident photon is absorbed by the sample. If so, proceed to step (11); otherwise, proceed to step (6);

[0061] (6) From the scattering coefficient μ of the sample s and random numbers, calculate the length of the next travel distance of the i-th incident photon in the sample;

[0062] (7) Calculate the next travel direction of the i-th incident photon in the sample based on the scattering phase function and anisotropy coefficient g of the sample and the random number;

[0063] (8) Assume that the i-th incident photon is scattered when it travels along the direction of the next section of travel described in step (7) to the end position of the next section of travel distance described in step (6), and calculate the end position of the incident photon based on the direction and distance of travel of the incident photon;

[0064] (9) Return to step (4);

[0065] (10) If the sample dish overflows, determine whether the i-th incident photon is received by the first photodetector used to receive the diffuse reflected light of the sample, the second photodetector used to receive the diffuse transmitted light of the sample, and the third photodetector used to receive the forward scattered light of the sample. If yes, add 1 to the accumulated number of photons of the photodetector that receives the i-th incident photon and proceed to step (11); otherwise, proceed directly to step (11);

[0066] (11) Compare i with the total number of incident photons N0. If i is less than N0, add 1 to i and return to step (1). Otherwise, go to step (12).

[0067] (12) The ratios of the accumulated photon number of the first photodetector, the second photodetector, and the third photodetector to N0 are: diffuse reflectivity measurement value R dc , diffuse transmittance measurement value T dc and the forward transmittance measurement T fc .

[0068] 5) According to the scattered light signal measurement value of the sample: diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f , and the calculated value of the scattered light signal of the sample: the calculated value of diffuse reflectance R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc , find the first optimization objective function ρ at the optical parameter vector point P δ (P) and the second optimization objective function ρ 1 / δ (P); specifically including:

[0069] Calculate the measured value of the scattered light signal of the sample and the m in the sample optical parameter space P The difference function δ(P) between the calculated values ​​of the scattered light signals of the points P with different optical parameter vectors is then calculated according to the m P The difference function δ(P) is used to calculate the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ (P);

[0070] The difference function δ(P) is defined as the sum of the squares of the relative differences between the measured value of the scattered light signal of the sample and the calculated value of the scattered light signal under the condition of the optical parameter vector point P:

[0071]

[0072] For the inverse scattering problem, the difference function δ(P) is usually selected as the objective function to measure the accuracy of the result. If δ(P)≤δ th , δ th represents the threshold value determined by the signal measurement error, then the inverse scattering problem is considered to be solved. In order to quantify the change of the calculated value of the scattered light signal caused by the Monte Carlo ray tracing method in the calculation of the scattered light signal, we define the statistical fluctuation error vector f of the calculated value of the scattered light signal c for:

[0073]

[0074] Where R dc∞ ,T dc∞ ,T fc∞ They are the ideal calculated values ​​of diffuse reflectance, diffuse transmittance and forward transmittance obtained by tracing infinite photons when calculating scattered light signals; therefore, the following formula is available:

[0075]

[0076] Among them, S c (P) represents the accuracy vector of the calculated value of the scattered light signal, o represents the operation of multiplying the vectors on both sides of the symbol one by one, S c∞ (P) is the ratio vector of the ideal calculated value and the measured value of the scattered light signal, expressed as follows:

[0077]

[0078] From formulas (6), (7), and (8), we can obtain:

[0079]

[0080] Among them, δ f (P)=|S c∞ (P)-1| 2 , is the ideal value of δ(P) without statistical fluctuation error, δ f (P) = 0 means that P is an ideal solution to the inverse scattering problem without statistical fluctuation error; thus the first optimization objective function ρ is obtained δ (P) and the second optimization objective function ρ 1 / δ (P):

[0081]

[0082]

[0083] Among them, A P is the common m in the sample optical parameter space centered at point P P A set of optical parameter vector points P';

[0084] In order to quantitatively analyze δ(P), ρ δ (P) and ρ 1 / δ (P), q is the probability of the Poisson distribution of the total number of incident photons N0 and the calculated value of the related scattered light signal in the Monte Carlo ray tracing method, q m is the average value of q, that is, the stable value of the signal calculation value corresponding to q, then f in formula (6) c The weight is (q / q m )-1. pIn the range, since all parameter vectors are close to point P, the δ of all points f (P) have the same value, and all f c Function q m When δ(P) ranges from 0.1% to 1%, the objective function ρ δ (P) and ρ 1 / δ (P) about δ f The distribution of (P) is as follows Figure 5a As shown, the solid lines represent δ(P) and δ(P) -1 , the results clearly show that when the difference function δ(P) is between 0.6% and 1%, the calculated value of the scattered light signal calculated by the Monte Carlo ray tracing method based on the radiation transfer theory has a statistical fluctuation error, which causes the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ (P) value fluctuations, but the second optimization objective function ρ 1 / δ The fluctuation of (P) is significantly smaller than the first optimization objective function ρ δ (P), and the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ The fluctuation of (P) is obviously smaller than the difference function δ(P), that is, the second optimization objective function ρ 1 / δ (P) is more suitable as the objective function of the inverse calculation algorithm described in the present invention; in order to further understand the influence of the statistical fluctuation error of the calculated value of the scattered light signal on the accuracy of the solution of the inverse scattering problem, we calculate ρ δ (P) and δ(P) and ρ 1 / δ (P) and δ(P) -1 The relative difference between them, and the average value of the relative difference is taken as Δ(ρ δ ) and Δ(ρ 1 / δ ):

[0085]

[0086]

[0087] Among them, δ f In formula (12) and formula (13), the sample is uniformly between 0.1% and 1.1%, m δ is δ f (P) total number of samples; δ f The δ corresponding to the sample value of (P) f The average values ​​of the signal items in (P) range from 1.8% to 6.1%. Figure 5b The function Δ(ρ δ ), Δ(ρ 1 / δ ) in m δfWhen set to 1000, q m Therefore, in the above case, the mean q of the probability of the Poisson distribution is m The ratio of the total number of incident photons N0 is equal to R dc The ideal value without statistical fluctuation error, q m The larger the corresponding N0 is, the larger the Figure 5b The results show that under the same N0 value, using ρ 1 / δ (P) instead of ρ δ (P) or δ(P) as the optimization objective function can significantly reduce the statistical fluctuation error of the calculated value of the scattered light signal; under the condition that the error is allowed not to exceed a fixed threshold, the total number of incident photons N0 required to be tracked in the Monte Carlo ray tracing method based on radiation transfer theory can be significantly reduced from the commonly used range of 10 to the power of 6 to 10 to the power of 10 to the power of 10 to the power of 5 to the power of 10 to the power of 67, that is, the calculation time is significantly reduced, thereby realizing real-time acquisition of multi-parameter spectral data of turbid media.

[0088] 6) Determine whether the first optimization objective function ρ δ (P) value is greater than the set objective function threshold ρ th , or the second optimization objective function ρ 1 / δ (P) value is less than the set objective function threshold ρ th , if yes, go to step 7); otherwise go to step 8);

[0089] 7) Search and update the sample optical parameter vector point P, and then return to step 4);

[0090] The search and update of the sample optical parameter vector point P is completed by a particle swarm optimization algorithm. The search is performed by multiple processes simultaneously. The particles with index k∈[1,K] are used to represent different search processes. k is the particle number, K is the total number of particles, and index j∈[1,J], j is the number of iterations of the search update process, J is the total number of iterations, and particle k represents different search trajectories in the sample optical parameter space. Trajectory k moves from position P(k,j) to position P(k,j+1). The update formula of P(k,j+1) is as follows:

[0091] P(k,j+1)=P(k,j)+v(k,j+1) (14)

[0092] v(k,j+1)=χ{w0v(k,j)+w p r(P p (k)-P(k,j))+w g r'(P g -P(k,j))} (15)

[0093] Where v(k,j) is the velocity of the k trajectory or particle at the jth iteration, v(k,j+1) is the velocity of the k particle at the j+1th iteration, r and r' are random numbers uniformly distributed between 0 and 1, and P p (k) is the optimal position of particle k, P g is the optimal position of all particles, w0,w p and w g is the weight of the term contributing to v(k,j+1), and χ is the damping constant that increases the stability of the algorithm.

[0094] 8) Save the best position P of all particles at present g , P g That is, the solution of the optical parameter vector of the turbid medium sample at the wavelength center value λ: the absorption coefficient μ a , scattering coefficient μ s and anisotropy coefficient g, complete the solution of the inverse scattering problem with the measurement wavelength center value λ;

[0095] 9) Select a new measurement wavelength center value λ and return to step 2) until the absorption coefficient μ of the turbid medium sample of all wavelength center values ​​within the measurement wavelength range is completed. a , scattering coefficient μ s The inverse calculation of the anisotropy coefficient g is then used to obtain the multi-optical parameter spectral data of the turbid medium.

[0096] In an embodiment of the present invention:

[0097] The light source unit 1 can be composed of an incoherent continuous light source with a continuous spectrum, such as a xenon lamp, an iodine tungsten lamp, or a halogen lamp, and then an optical path composed of a lens or a reflector generates an incident light beam that meets the measurement requirements; the light source unit 1 can also be composed of a coherent continuous laser combination with adjustable wavelength or a combination of multiple continuous light-emitting diodes with different wavelengths.

[0098] In the embodiment of the present invention, the modulator 2 can realize the intensity modulation of the incident light beam by mechanical, electro-optical or acousto-optic modulation, such as the optical chopper model SR540 produced by Stanford Research Systems, the electro-optical modulator model 84503011002 produced by Qioptiq, and the acousto-optic modulator model 23080-1 produced by Gooch & Housego. The light modulation frequency of the modulator 2 is adjustable in the range of 1 Hz to 100 MHz.

[0099] In an embodiment of the present invention, the spectrometer 3 is used to obtain a monochromatic incident light beam within a set measurement wavelength range, wherein the central wavelength value λ is adjustable within a range of 200nm to 2500nm or within the partial wavelength range, and the bandwidth Δλ is adjustable within a range of 0.5nm to 15nm. The spectrometer 3 can use a prism spectrometer as a main component, such as an equilateral prism of model GCL-03013 produced by Daheng Optoelectronics Co., Ltd.; the spectrometer 3 can also use a diffraction grating device as a main component, such as a diffraction grating of model GR25-1204 produced by Thorlabs Co., Ltd.; the spectrometer 3 can also use an interference filter combination as a main component, such as a visible wavelength range interference filter series of model FKB-VIS-10 and a near-infrared wavelength range interference filter series of model FKB-IR-10 produced by Thorlabs Co., Ltd.

[0100] In an embodiment of the present invention, the incident light intensity detection device 4 passes the incident light beam through an optical path composed of a lens and (or) a reflector and an aperture and other devices to obtain a set light beam area and power distribution, and adjusts the incident light beam interface to a circular or elliptical shape, and its diameter is adjustable in the range of 2 mm to 30 mm; the incident light intensity detection device 4 also includes a spectrometer and a photodetector, which are used to monitor the intensity of the incident light beam before entering the sample device composed of a sample dish 6 and a sample 7; when the incident light beam is incident on the sample device, part of the light beam energy will be reflected by the sample dish 6, and the rest will enter the sample 7 and be absorbed or scattered in the sample 7. The scattered light can be emitted from the sample 7 through the sample dish 6 to form diffusely scattered light.

[0101] In the embodiment of the present invention, the photodetector in the incident light intensity detection device 4 and the first photodetector 9, the second photodetector 10 and the third photodetector 11 for detecting scattered light can be photodiodes or photomultiplier tubes, such as photodiodes of model FDS10X10, FDS100 or FGA21 produced by Thorlabs, or photomultiplier tubes of model R5473-02 or R5473-02 produced by Hamamatsu. The first photodetector 9 collects diffuse reflected light emitted from the sample 7 and transmitted through the sample dish 6, the second photodetector 10 collects diffuse transmitted light emitted from the sample 7 and transmitted through the sample dish 6, and the third photodetector 11 collects forward diffuse transmitted light emitted from the sample 7 and transmitted through the sample dish 6. The scattered light intensity analog electrical signals output by the photodetector in the incident light intensity detection device 4 and the first photodetector 9, the second photodetector 10 and the third photodetector 11 are sent to the signal processing unit 12.

[0102] In the embodiment of the present invention, the signal processing unit 12 is composed of circuits such as signal amplification, filtering and demodulation, and analog-to-digital conversion. In the embodiment of the present invention, the input scattered light intensity analog electrical signal can be amplified and filtered and demodulated by the phase-locked amplifier principle according to the modulation frequency of the incident light beam, such as the phase-locked amplifier model SR830 produced by Stanford Research Systems, and then converted into a scattered light intensity digital electrical signal through an analog-to-digital conversion circuit and output. Or in the embodiment of the present invention, the scattered light intensity analog electrical signal can be input to the amplification, filtering and analog-to-digital conversion circuit, and then the scattered light intensity digital signal is sampled and read multiple times and its average value and mean square error value are calculated, and then output to the control calculation unit 13, and the scattered light signal measurement value is calculated according to equations (2) to (4).

[0103] The control calculation unit 13 of the present invention can be implemented by a computer, including a user interface program, a system control program and an inverse calculation program. The system control program controls the light source unit, the light modulator and the spectrometer and the incident light intensity detection device to obtain an incident light beam that meets the measurement requirements, sample and read the scattered light intensity digital signal, calculate and obtain the measurement value of the scattered light signal, and obtain the system physical parameters and the incident light beam parameters. The flow of the Monte Carlo ray tracing method based on the radiation transfer theory in the inverse calculation program is as follows: Figure 3 As shown, this method is based on the initial preset value μ of the sample optical parameter a0 , μ s0 , g0, and obtain the calculated value R of the scattered light signal dc , T dc With T fc The inverse calculation process is as follows: Figure 4 As shown, the program calculates the first optimization objective function such as ρ according to the scattered light signal calculation value and the scattered light signal measurement value. δ function and the second optimization objective function such as ρ 1 / δ function, and the corresponding threshold ρ th If the optimization objective function does not meet the optimization conditions, the sample optical parameter vector point P is searched and updated, and the above steps are iterated until the objective function meets the optimization conditions; at this time, the optimal value P of the sample optical parameter in the inverse calculation program g That is, the sample optical parameter value μ of the turbid medium sample at the measurement wavelength λ a , μ s and g, stored in the output data file.

[0104] The modulator unit and the beam splitter unit are controlled to change the central wavelength of the incident light beam, and the above measurement and calculation process is repeated cyclically until the multi-parameter spectral data of the turbid medium is obtained.

[0105] In order to verify the inverse calculation method for obtaining real-time multi-optical parameter spectral data of turbid media of the present invention, we used 20% fat emulsion (I141-100ML, Sigma-Aldrich) as the turbid medium sample and obtained the scattered light signal measurement value R of the sample under the condition that the incident light beam wavelength ranged from 520 to 1000nm and the step length was 20nm. d ,T d and T f spectral data. Figure 6 The spectrum and error distribution diagram of the scattered light signal measurement value of the sample with a thickness of D = 102 μm; the system parameters are: the incident angle of the incident light beam is 8.5°, the incident light beam is an ellipse with a major axis of 11.4 mm and a minor axis of 4.8 mm, the first photodetector distance is 34.4 mm, the angle is 38°, the second photodetector distance is 39.3 mm, the angle is 39°, and the third photodetector distance is 63.6 mm, the angle is 8.5°; Figure 6 The figure shows the average value and standard deviation obtained by repeating the measurement three times, and the results show that the scattered light signal measurement data obtained by this method is very stable and accurate.

[0106] Figure 7 ρ 1 / δ (P) is the objective function, and the spectral data of the optical parameter vector point P of the sample with thickness D = 102 μm is obtained based on the particle swarm algorithm. The regional boundary of the sample optical parameter space is set at all wavelengths: μ a Value is 1.00x10 -4 and 2.00mm –1 Between, μ s Values ​​between 10.0 and 200 mm -1 The g value is between 0.10 and 1.00. The parameters in the optical parameter vector point P update formula (14) and formula (15) are: K = 10, J = 100, w0 = 0.900, w p =w g =2.0 and χ=0.730. The initial iterative optical parameter vector point P(k,1) of the kth example is randomly distributed in 9 equally divided intervals of the sample optical parameter space to ensure that the global optimal solution is obtained faster and avoid the interference of the local optimal solution of the non-global optimal solution. We found that under the above conditions, the selection of the initial value of P(k,1) in the 9 equally divided intervals of the sample optical parameter space does not affect the accuracy and convergence of the final solution. After calculating the objective function of P(k,1), iterate in the parameter space according to the definition of formula (14) and formula (15). When ρ 1 / δ (P) is greater than ρ th When the search stops, the P g That is, the optical parameter spectrum data of the turbid medium sample. Threshold ρth Corresponding to the fact that δ(P) in formula (11) is all set to 0.8%, m P 7. The calculation time of the spectral data of the turbid medium sample is about 160 seconds on an ordinary desktop computer equipped with a GeForce RTX 2080 Ti graphics card produced by Nvidia. If a newer GeForce RTX 4090 Ti graphics card produced by Nvidia is used, or multiple graphics cards, such as 5 GeForce RTX 2080 Ti graphics cards, the function of obtaining multi-parameter spectral data of turbid media in real time, that is, no more than 30 seconds per wavelength, can be realized.

Claims

1. An inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time, characterized in that: The steps include: 1) obtaining a turbid medium sample and obtaining physical parameters of a system for measuring multi-parameter spectral measurement data of the turbid medium, wherein the physical parameters include: shape, size and refractive index of the sample and the sample dish, shape, size, distance and orientation of the detector relative to the coordinate origin defined by the reference point and reference direction on the surface of the sample dish and the coordinate axis; 2) obtaining the shape, area, power distribution, incident angle, measurement values ​​of the center value λ of the measured wavelength and the bandwidth of the incident light beam, as well as the measurement value of the light intensity of the incident light beam, and causing the incident light beam to be incident on the sample; 3) Measure the scattered light signals emitted from the sample and the sample dish at multiple angles relative to the sample dish, and obtain the scattered light signal measurement value at the center value λ of the measurement wavelength, including the diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f ; 4) Three optical parameters of turbid medium samples at different wavelengths, namely, absorption coefficient μ a , scattering coefficient μ s The anisotropy coefficient g is represented by the optical parameter vector point P, which is denoted as P(λ)=(μ a ,μ s ,g), where (μ a ,μ s ,g) is called the sample optical parameter space; the initial preset value of the sample optical parameter vector is set to P0(λ) = (μ a0 ,μ s0 ,g0), where μ a0 is the initial absorption coefficient of the sample, μ s0 is the initial scattering coefficient of sample 2, g0 is the initial anisotropy coefficient of the sample, and the covariance m of P(λ) in the sample optical parameter space is calculated. P -1 adjacent point, take m P =7, there are six adjacent optical parameter vector points P, which are (μ a0 ±s a ,μ s0 ,g0),(μ a0 ,μ s0 ±s s ,g0),(μ a0 ,μ s0 ,g0±s g ), where s a 、s s and g They are the parameter distances of adjacent points, and their values ​​are s a =0.01μ a0 、s s =0.01μ s0 and g = 0.01g0, or m P Take a value between 5 and 27, sa, ss and sg take values ​​between 1 and 10 times of the corresponding values ​​respectively; m including the initial preset value P(λ) of the optical parameter vector P Different optical parameter vector points are introduced into the Monte Carlo ray tracing method based on radiation transfer theory to obtain m P A group of different scattered light signal calculation values, each group of scattered light signal calculation values ​​includes a diffuse reflectance calculation value R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc ; 5) According to the scattered light signal measurement value of the sample: diffuse reflectance measurement value R d , diffuse transmittance measurement value T d and the forward transmittance measurement T f , and the calculated value of the scattered light signal of the sample: the calculated value of diffuse reflectance R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc , find the first optimization objective function ρ at the optical parameter vector point P δ (P) and the second optimization objective function ρ 1 / δ (P); 6) Determine whether the first optimization objective function ρ δ (P) value is greater than the set objective function threshold ρ th , or the second optimization objective function ρ 1 / δ (P) value is less than the set objective function threshold ρ th , if yes, go to step 7); otherwise go to step 8); 7) Search and update the sample optical parameter vector point P, and then return to step 4); 8) Save the best position P of all particles at present g , P g That is, the solution of the optical parameter vector of the turbid medium sample at the wavelength center value λ: the absorption coefficient μ a , scattering coefficient μ s and anisotropy coefficient g, complete the solution of the inverse scattering problem with the measurement wavelength center value λ; 9) Select a new measurement wavelength center value λ and return to step 2) until the absorption coefficient μ of the turbid medium sample of all wavelength center values ​​within the measurement wavelength range is completed. a , scattering coefficient μ s The inverse calculation of the anisotropy coefficient g is then used to obtain the multi-optical parameter spectral data of the turbid medium.

2. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 1, characterized in that: The shape and area of ​​the incident light beam described in step 2) are smaller than the area of ​​the incident light cross section of the optical window of the sample dish and the sample; the power distribution of the incident light beam ensures that the power ratio of the same area between the center of the incident light beam and the edge of the incident light beam is less than 104; the incident angle of the incident light beam is 5 to 60 degrees; the center value of the wavelength, bandwidth and modulation frequency of the incident light beam are obtained by controlling the optical modulator and the spectrometer in the multi-parameter spectral measurement data system for measuring turbid media; before the incident light beam enters the sample, the incident light intensity signal value I0 is obtained, and then the light beam is incident on the sample in the sample dish.

3. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 1, characterized in that: In step 3), when the wavelength center value is measured as λ, the incident light intensity signal value I0 and the measured diffuse reflection light intensity I Rd , Diffuse transmission light intensity I Td and the forward scattered light intensity I Tf After amplification, filtering and demodulation, the diffuse reflectance measurement value R is calculated using the following formula: d , diffuse transmittance measurement value T d and the forward transmittance measurement T f :

4. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 1, characterized in that: Step 4) includes the initial preset value P(λ) of the optical parameter vector m P Different optical parameter vector points are introduced into the Monte Carlo ray tracing method based on radiation transfer theory to obtain m P A group of different scattered light signal calculation values, each group of scattered light signal calculation values ​​includes a diffuse reflectance calculation value R dc , Diffuse transmittance calculated value T dc and the calculated forward transmittance T fc , specifically including: (1) For the i-th incident photon, the absorption coefficient μ of the sample is a and random numbers, calculate the total path length of the incident photon in the sample; (2) From the scattering coefficient μ of the sample s and random numbers, calculate the length of the next travel distance of the i-th incident photon in the sample; (3) Assume that the i-th incident photon is scattered when it travels along the initial travel direction in the sample to the end position of the next travel distance described in step (2), and calculate the end position of the incident photon based on the travel distance of the incident photon; (4) Determine whether the i-th incident photon overflows the sample dish at the end position described in step (3): if yes, proceed to step (10); otherwise, proceed to step (5); (5) Calculate the total length of the cumulative travel distance of the i-th incident photon in the sample. If the total length of the cumulative travel distance is greater than the total length of the travel distance, determine whether the i-th incident photon is absorbed by the sample. If so, proceed to step (11); otherwise, proceed to step (6); (6) From the scattering coefficient μ of the sample s and random numbers, calculate the length of the next travel distance of the i-th incident photon in the sample; (7) Calculate the next travel direction of the i-th incident photon in the sample based on the scattering phase function and anisotropy coefficient g of the sample and the random number; (8) Assume that the i-th incident photon is scattered when it travels along the direction of the next section of travel described in step (7) to the end position of the next section of travel distance described in step (6), and calculate the end position of the incident photon based on the direction and distance of travel of the incident photon; (9) Return to step (4); (10) If the sample dish overflows, determine whether the i-th incident photon is received by the first photodetector used to receive the diffuse reflected light of the sample, the second photodetector used to receive the diffuse transmitted light of the sample, and the third photodetector used to receive the forward scattered light of the sample. If yes, add 1 to the accumulated number of photons of the photodetector that receives the i-th incident photon and proceed to step (11); otherwise, proceed directly to step (11); (11) Compare i with the total number of incident photons N0. If i is less than N0, add 1 to i and return to step (1). Otherwise, go to step (12). (12) The ratios of the accumulated photon number of the first photodetector, the second photodetector, and the third photodetector to N0 are: diffuse reflectivity measurement value R dc , diffuse transmittance measurement value T dc and the forward transmittance measurement T fc .

5. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 1, characterized in that: Step 5) includes calculating the measured value of the scattered light signal of the sample and the m in the sample optical parameter space P The difference function δ(P) between the calculated values ​​of the scattered light signals of the points P with different optical parameter vectors is then calculated according to the m P The difference function δ(P) is used to calculate the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ (P).

6. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 5, characterized in that: Specifically include: The difference function δ(P) is defined as the sum of the squares of the relative differences between the measured value of the scattered light signal of the sample and the calculated value of the scattered light signal at the optical parameter vector point P: Define the statistical fluctuation error vector f of the calculated value of the scattered light signal c for: Where R dc∞ ,T dc∞ ,T fc∞ They are the ideal calculated values ​​of diffuse reflectance, diffuse transmittance and forward transmittance obtained by tracing infinite photons when calculating scattered light signals; therefore, the following formula is available: Among them, S c (P) represents the accuracy vector of the calculated value of the scattered light signal, o represents the operation of multiplying the vectors on both sides of the symbol one by one, S c∞ (P) is the ratio vector of the ideal calculated value and the measured value of the scattered light signal, expressed as follows: From formulas (6), (7), and (8), we can obtain: Among them, δ f (P)=|S c∞ (P)-1| 2 , is the ideal value of δ(P) without statistical fluctuation error, δ f (P) = 0 means that P is an ideal solution to the inverse scattering problem without statistical fluctuation error; thus the first optimization objective function ρ is obtained δ (P) and the second optimization objective function ρ 1 / δ (P): Among them, A P is the common m in the sample optical parameter space centered at point P P A set of optical parameter vector points P'; When the difference function δ(P) is between 0.6% and 1%, the calculated value of the scattered light signal calculated by the Monte Carlo ray tracing method based on the radiation transfer theory has a statistical fluctuation error, which causes the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ (P) value fluctuations, but the second optimization objective function ρ 1 / δ The fluctuation of (P) is significantly smaller than the first optimization objective function ρ δ (P), and the first optimization objective function ρ δ (P) and the second optimization objective function ρ 1 / δ The fluctuation of (P) is obviously smaller than the difference function δ(P), that is, the second optimization objective function ρ 1 / δ (P) is more suitable as the objective function of the inverse calculation algorithm; calculate ρ δ (P) and δ(P) and ρ 1 / δ (P) and δ(P) -1 The relative difference between them, and the average value of the relative difference is taken as Δ(ρ δ ) and Δ(ρ 1 / δ ): Among them, δ f In formula (12) and formula (13), the sample is uniformly between 0.1% and 1.1%, m δ is δ f (P) total number of samples; δ f The δ corresponding to the sample value of (P) f The average value of each signal term in (P) ranges from 1.8% to 6.1%. The average value q of the probability of Poisson distribution m The ratio of the total number of incident photons N0 is equal to R dc The ideal value without statistical fluctuation error, q m The larger the value, the larger the corresponding N0; under the same N0 value, use ρ 1 / δ (P) instead of ρ δ (P) or δ(P) as the optimization objective function can significantly reduce the statistical fluctuation error of the calculated value of the scattered light signal; under the condition that the error is allowed to not exceed a fixed threshold, the total number of incident photons N0 required to be tracked in the Monte Carlo ray tracing method based on radiation transfer theory can be significantly reduced from the commonly used 10 to the power of 6 to 10 to the power of 10 to the power of 10 to the power of 5 to the power of 10 to the power of 7, that is, the calculation time is significantly reduced, so that multi-parameter spectral data of turbid media can be obtained in real time.

7. The inverse calculation method for obtaining multi-optical parameter spectral data of turbid media in real time according to claim 1, characterized in that: The search and update of the sample optical parameter vector point P in step 7) is completed by a particle swarm optimization algorithm. The search is performed by multiple processes simultaneously, and particles with index k∈[1,K] are used to represent different search processes. k is the particle number, K is the total number of particles, and index j∈[1,J], j is the number of iterations of the search update process, J is the total number of iterations, and particle k represents different search trajectories in the sample optical parameter space. When trajectory k moves from position P(k,j) to position P(k,j+1), the update formula of P(k,j+1) is as follows: P(k,j+1)=P(k,j)+v(k,j+1) (14) v(k,j+1)=χ{w0v(k,j)+w p r(P p (k)-P(k,j))+w g r'(P g -P(k,j))} (15) Where v(k,j) is the velocity of the k trajectory or particle at the jth iteration, v(k,j+1) is the velocity of the k particle at the j+1th iteration, r and r' are random numbers uniformly distributed between 0 and 1, and P p (k) is the optimal position of particle k, P g is the optimal position of all particles, w0,w p and w g is the weight of the term contributing to v(k,j+1), and χ is the damping constant that increases the stability of the algorithm.

Images