Construction method of high-resolution brain blood oxygen dynamic response 3D map

Abstract

The invention discloses a construction method of a high-resolution brain blood oxygen dynamic response 3D atlas, which comprises the following steps: selecting a rectangular network or a honeycomb network to detect a plurality of brain regions, and arranging a plurality of optical fiber probes; dividing an area below the basic unit into a plurality of local areas, and calculating the path length of the time resolution part; a signal inspection strategy is executed under the wavelength 1 and the wavelength 2, and two-wavelength first-measurement diffused light signals are obtained; obtaining a wavelength 1 re-measured diffused light signal at intervals; executing an element signal conjoint analysis strategy on the basic unit to obtain a wavelength 1 local region absorption coefficient change; a remeasured diffused light signal with the wavelength 2 is obtained; executing an element signal conjoint analysis strategy on the basic unit to obtain a wavelength 2 local region absorption coefficient change; solving an equation set of the same wavelength and the local area to obtain blood oxygen change of the local area; the brain blood oxygen dynamic response 3D atlas is obtained. According to the method, the spatial resolution of the brain blood oxygen dynamics response 3D spectrum is improved, and technical reserve is provided for the 3D brain depth information full spectrum.

Description

Technical Field

[0001] The invention belongs to the technical field of near-infrared spectroscopy, and particularly relates to a method for constructing a high-resolution cerebral hemooxygen dynamics response 3D atlas. Background Art

[0002] With the continuous maturity of near infrared spectroscopy analysis technology, various targeted near infrared instruments have emerged, among which functional near infrared spectroscopy (fNIRS) has attracted much attention as a technology with great potential. fNIRS is a non-invasive technology that can be used to monitor the changes of oxyhemoglobin and deoxyhemoglobin during brain activity. The technology is currently in a stage of continuous improvement and development, and is mainly used in the fields of advanced cognition, brain science, neurology and psychology in natural environments. Brain tissue activity will cause changes in blood oxygen content, which in turn affects the optical properties of blood oxygen, mainly manifested in absorption and scattering effects. The basic principle of fNIRS is to reflect the changes in blood oxygen and blood volume in the brain by monitoring the changes in the concentration of oxyhemoglobin and deoxyhemoglobin in the blood, thereby inferring brain activity. Since red blood cells have good absorption characteristics for light with wavelengths below 650nm, and water absorbs light with wavelengths above 900nm, fNIRS uses near-infrared light with a wavelength range of 650nm to 900nm [Guo Jianghui. Research on high-time resolution fNIRS and its application in lung diseases [D]. University of Electronic Science and Technology of China, 2023. DOI: 10.27005 / d.cnki.gdzku.2023.003016.]. In this band, it has good detection capabilities for oxygenated hemoglobin and deoxygenated hemoglobin. fNIRS systems are mainly divided into three types: the first is continuous wave spectroscopy system (Continuous Wave, CW); the second is frequency domain spectroscopy system (Frequency Domain, FD); the third is time domain spectroscopy system (Time Domain, TD).

[0003] The incident light used in time-domain fNIRS technology is a pulse with a short pulse width, generally in the order of picoseconds (1 picosecond = 10 -12Seconds). Usually, the measurement of light intensity by time-domain fNIRS is to record the time of photons through a high-sensitivity detector or a single-photon counter. Since the distance that light travels in tissue is proportional to the flow time and the diffuse reflection in tissue is random, when light starts from the light source, different light intensities may exist at different times (generally hundreds of picoseconds). When the detector receives a large number of photons, the transmission distance of the photons will produce a certain probability distribution. The more photons there are, the more obvious this distribution is. Correspondingly, the distribution of the light intensity of the outgoing light over time can be obtained. This distribution is called the temporal point spread function (TPSF) in optics. By analyzing the photon flight time encoding information contained in the TPSF, the optical parameters of the tissue can be obtained, thereby obtaining the changes in the HbO and HbR concentrations in the tissue [Torricelli A, Contini D, Pidderi A, et al. Time domain functional NIRS imaging for human brain mapping [J / OL]. NeuroImage, 2014, 85: 28-50.]. The time domain measurement method is applicable to each acquisition channel, that is, each source-detector pair, without the need for multi-channel detection analysis like continuous wave near-infrared spectroscopy. The common measurement methods currently include: fitting method, moment calculation method and gating method (also called time gating or time window method). The principle of the gating method is to use the separation of TPSF signals in different time gates to represent different photon arrival times, as well as the influence of multi-layer tissue characteristics on the flight path of the detection light over time to evaluate the hemodynamic changes in different brain regions [Lange, Frédéric, and Ilias Tachtsidis. "Clinical brain monitoring with time domain NIRS: a review and future perspectives." Applied Sciences 9.8 (2019): 1612.].

[0004] Mean Partial Pathlengths (MPP) is the average flight path length of the diffuse light received in fNIRS measurement in different medium domains, which indicates the sensitivity of the measured optical density to the absorption changes in the layer. In 2000, Steinbrink et al. extended the mean partial path length to the time domain diffuse light measurement and proposed the time-dependent mean partial path length (TMPP) [Steinbrink, J., et al. "Determining changes in NIR absorption using a layered model of the human head." Physics in Medicine & Biology 46.3 (2001): 879.]. TMPP describes the effect of absorption changes in different layers on time-resolved reflectivity and can be used to quantify the effect of multi-layer tissue characteristics on the flight path of the detection light that changes over time.

[0005] In the study of light diffusion in multilayer turbid media, the analytical solution of the radiation transfer equation approximated by the diffusion equation is usually used to analyze the target medium and calculate the theoretical expression of each TMPP. However, as the composition and optical parameters of the turbid medium become more complex, the above solution method is no longer robust. There are a large number of literatures that use Monte Carlo photon transport simulation to calculate TMMP [Vera, Demián A., et al. "Retrieval of chromophore concentration changes in a digital human head model using analytical mean partial pathlengths of photons." Journal of Biomedical Optics 29.2 (2024): 025004-025004.]. Monte Carlo photon transport simulation is a computer simulation technology that has been widely used in bio-optical applications, medical imaging, optoelectronic device design and other fields. The basic idea is to generate a large number of photons randomly, simulate their motion trajectories according to certain physical laws, and calculate the reflection, transmission, absorption and other parameters at different wavelengths on the surface or inside of the object through statistical analysis of these trajectories, so as to infer the information such as the object's tissue structure, composition or morphology. Monte Carlo simulation is used as a stochastic method to solve the radiation transfer equation, which governs the photon migration in general cases, where light can not only be absorbed but also scattered in an anisotropic manner. Therefore, MC simulations are the gold standard for simulating light transmission experiments in multilayer turbid media, because they can accurately control all parameters involved in the photon migration process in turbid media. In the time-domain-based MC simulation, the TMPP of different layers can be obtained through time-resolved measurements and formula derivation.

[0006] A key issue in the application of near-infrared photon migration is how to simultaneously increase the penetration depth of the measurement and reduce the measurement area, which is especially true for detecting the brain. The scalp, skull and meninges weaken the optical changes produced by the cerebral cortex. When detecting some human tissues (such as the breast), local positioning is also required to detect and characterize small deep lesions. Time-resolved diffuse reflectance measurement at zero source-detector distance solves the above problems to a certain extent compared to time-resolved diffuse reflectance measurement at a longer source-detector distance. Its advantages are: (1) increasing the number of photons collected at any arrival time; (2) producing higher contrast in local areas; (3) providing better spatial resolution; (4) easily locating the brain tissue area of ​​interest [Pifferi, Antonio, et al. "Time-resolved functional near-infrared spectroscopy at null source-detector separation." Biomedical Optics. Optica Publishing Group, 2008.].

[0007] The self-detection-mutual-detection combined time-domain fNIRS is based on a bidirectional functional diffuse flight time-resolved fiber brain-computer interface probe and a single-photon avalanche diode with extremely high time resolution. It has both the self-detection function of a single fiber probe and the mutual detection function of dual fiber probes. Through time resolution, it can not only effectively improve the lateral and longitudinal spatial resolution of detection, but also enhance the detection effect of the depth of brain tissue. It greatly improves the resolution of traditional functional near-infrared cerebral hemodynamic response maps, making it possible to obtain high-resolution dynamic cerebral blood flow and cerebral blood oxygen 3D maps.

[0008] Arranging the fiber optic probes of fNIRS in a certain spatial pattern and source-detector distance to detect a wide area of ​​the brain can be applied to single-layer or tomographic detection imaging of brain tissue. Tomographic detection imaging of fNIRS, also known as diffuse optical tomography (DOT), can provide 3D images with higher spatial resolution than conventional 2D detection imaging of fNIRS. DOT measures the diffuse light emitted from the boundary of biological tissue, obtains the temporal and spatial distribution information of the light signal, and establishes a forward model based on the light transmission process. The original data is inverted to reconstruct the optical parameters (such as absorption coefficient and scattering coefficient) at different spatial positions inside the tissue, thereby performing tomographic imaging of the tissue and obtaining information such as hemoglobin content, blood oxygen saturation, and blood volume. Due to the use of diffuse light, the spatial resolution of DOT is relatively low (usually 1 cm or less), and the imaging quality needs to be improved. DOT in continuous wave measurement mode requires fNIRS equipment that supports flexible arrangement of a large number of optical fibers, and its application and analysis are usually difficult [Song Bowen, Zhao Yanyu. Diffuse optical imaging methods and applications (invited) [J]. Laser and Optoelectronics Progress, 2024, 61(08): 11-36.]. Applying self-detection-mutual detection combined time-domain fNIRS to DOT can reduce the number and spacing limitations of optical fiber probes, improve the lateral and depth resolution and detection depth of DOT, and obtain high-resolution brain neural activity and hemodynamic response maps.

[0009] At present, the probe geometric arrangements of fNIRS systems are mostly square and diamond arrangements, which can increase the utilization of transmitting and receiving probes [Ma Pei, Shen Wushuang, Shen Huijuan, et al. Review of functional near-infrared brain imaging systems [J]. Optical Instruments, 2022, 44(05): 1-13.]. However, the probe geometric arrangement dimensions of most commercial multi-channel fNIRS systems are not suitable for exploring or locating functional activations in the human cerebral cortex. These systems provide a spatial measurement grid with a spacing of approximately 30 mm, while the human cerebral cortex is composed of gyri approximately 10-15 mm wide, which can be used to indicate independent functions. The detection signals of each pair of sources and detectors with a spacing of 30 mm have a low discrimination of functional brain areas in the cerebral cortex, and the depth specificity needs to be improved. Therefore, this multi-channel measurement may ignore the activity of some cortical gyri (false negative error) [Yamada, Toru, et al. "Wearable fNIRS device of higher spatial resolution realized by triangular arrangement of dual-purpose optodes." Neural Imaging and Sensing 2023. Vol. 12365. SPIE, 2023.]. The detection depth of time-domain fNIRS is not limited by the spacing of fiber probes. Incorporating the zero-source detector distance and the small-source detector distance into the geometric arrangement of the probe can be used to increase the channel density and spatial resolution and provide accurate local brain area positioning. Therefore, the use of self-detection-mutual detection combined time-domain fNIRS to detect whole-brain functional activities requires the reconstruction of the optical information network structure of fiber probe self-detection and mutual detection. For brain areas with special functions, it is necessary to study the scheme of optimizing the detection grid such as probe position and spacing parameters, that is, brain area-specific probe geometric arrangement.

[0010] In order to overcome the shortcomings of the prior art and construct a more accurate high-resolution cerebral blood oxygen dynamics response 3D atlas, the present invention selects a detection network according to the channel density, detects several brain areas, deploys several fiber optic probes of the fiber optic cerebral blood oxygen detection system in the detection network, selects a suitable detection network according to the calculated channel density of different detection networks, divides the brain area corresponding to the basic unit into local areas, solves the time-resolved partial path length, adopts a signal inspection strategy to obtain a time-resolved reflectivity signal, executes a primitive signal joint analysis strategy for each basic unit, performs curve inversion on the time-resolved reflectivity signals obtained successively, solves the absorption change equation group, obtains the local area cerebral blood oxygen change value, constructs a three-dimensional RGB map of cerebral blood oxygen dynamics response, obtains a high-resolution cerebral blood oxygen dynamics response 3D atlas, obtains a method for multi-channel time-resolved measurement of brain tissue blood oxygen changes from surface to domain, and obtains a high-resolution cerebral blood oxygen dynamics response 3D atlas, which provides a technical reserve for obtaining the goal of a full spectrum of 3D brain depth information. Summary of the invention

[0011] The purpose of the present invention is to provide a method for constructing a high-resolution 3D map of cerebral hemooxygen dynamic response, aiming to utilize the zero-spacing source-detector fiber probe of the fiber optic cerebral blood oxygen detection system and the multi-channel time-resolved measurement from surface to domain to improve the spatial resolution and detection depth of time-domain fNIRS, and obtain a high-resolution 3D map of cerebral hemooxygen dynamic response.

[0012] In order to achieve the above-mentioned object of the invention, the present invention adopts the following technical solutions:

[0013] A method for constructing a high-resolution 3D atlas of cerebral hemodynamic response, characterized in that the method steps include:

[0014] 1) Select n brain A brain region to be detected is used, and a detection network is used to surround the area of ​​the cerebral scalp corresponding to the brain region. The detection network is a rectangular network with a size of MxN with a rectangle as the basic unit or a hexagonal honeycomb network with a perimeter of C with an equilateral triangle as the basic unit. For the rectangular network, the basic unit is set to be a rectangle of M1xN2, where M1 is less than or equal to M, and N2 is less than or equal to N. For the hexagonal honeycomb network, the basic unit is set to be an equilateral triangle with a perimeter of C1, where C1 is less than or equal to C. Each intersection in the detection network is a probe arrangement point for placing the optical fiber probe of the optical fiber brain blood oxygen detection system. The optical path emitted by the optical fiber probe into the brain region, which undergoes absorption and scattering, and is received by another optical fiber probe is defined as a mutual detection channel, and the optical path emitted by the optical fiber probe into the brain region, which undergoes absorption and scattering, and is received by the same optical fiber probe is defined as a self-detection channel. The channel spacing ρ is the spacing between the transmitting optical fiber and the receiving optical fiber;

[0015] 2) The calculation formula for the channel density of the rectangular network and the hexagonal honeycomb network is:

[0016]

[0017] Where: is the number of channels of the ith basic unit, S i (j) is the surface area of ​​the jth brain region covered by the ith basic unit. i When (j) is zero, let S i (j) is 1;

[0018] A number of optical fiber probes are respectively set in the rectangular network and the hexagonal honeycomb network. If the channel density of the rectangular network is greater than or equal to the channel density of the hexagonal honeycomb network, the rectangular network is used as the detection network, otherwise the hexagonal honeycomb network is used as the detection network;

[0019] 3) All fiber optic probe arrangement points are numbered, the time resolution of the single-photon avalanche diode coupled to the fiber optic probe is set to T, and R fiber optic probes are deployed at probe arrangement points with different numbers; if a rectangular network is used, the brain area below each basic unit is taken as a rectangular body area, the length, width and height of the rectangular body area are not greater than three times the side length of the basic unit and not less than twice the side length of the basic unit, and it is divided into n local areas, and the local areas are divided into n-1 measurement areas and other areas, the thickness of each measurement area is not less than 1mm and not more than 1cm, and the remaining area excluding the measurement area is the other area; if a hexagonal honeycomb network is used, the brain area below each basic unit is taken as a regular triangular prism area, the edge length of the regular triangular prism area is not greater than three times the side length of the basic unit and not less than twice the side length of the basic unit, and it is divided into n local areas, and the local areas are divided into n-1 measurement areas and other areas, the thickness of each measurement area is not less than 1mm and not more than 1cm, and the remaining area excluding the measurement area is the other area;

[0020] 4) Seven parameters are set for each local area, including length, width, thickness, refractive index, absorption coefficient, scattering coefficient and anisotropy factor. On this basis, the number of photons is simulated by using the voxel-based time-resolved Monte Carlo simulation method. 8 Using Sobel sampling, we calculate the time-dependent average partial path length L of the jth local area for different channel spacings ρ, wavelengths λ and detection time t for different basic units i. i,j (ρ,λ,t);

[0021] 5) The signal inspection strategy is as follows: select one of the fiber optic probes as the core probe, drive the laser diode coupled to the core probe to emit near-infrared pulse light with a light intensity of S and a main wavelength of λ, and at the same time, let the fiber optic probes in all basic units including the core probe receive the diffuse light signal of the brain area. The time for receiving the signal is set to T r, T r is an integer multiple of T. After the signal is received, all fiber optic probe numbers that receive the diffuse light signal of the brain area are placed in the waiting queue in clockwise order, and the head of the queue is taken out. The fiber optic probe corresponding to the head number is used as the core probe. The above operation is repeated until there is no fiber optic probe number in the queue. If there is a deployed fiber optic probe that does not emit pulse light as a core probe, it is used as a core probe to repeat the above operation until all fiber optic probes have emitted pulse light as core probes. The time T for measuring a diffuse light signal of the brain area is calculated. λ =T r xR, set the main wavelength of the pulse light of the laser diode to λ1, execute the signal inspection strategy once, obtain the first measured brain area diffuse light signal of wavelength λ1, then adjust the main wavelength of the pulse light of the laser diode to λ2, execute the signal inspection strategy once, and obtain the first measured brain area diffuse light signal of wavelength λ2;

[0022] 6) After an interval of ΔT, the signal inspection strategy is repeatedly executed once at wavelength λ1 to obtain the second measured brain diffuse light signal at wavelength λ1, and the brain diffuse light signals of all channels at wavelength λ1 are converted into time-resolved reflectance signals using software;

[0023] 7) Execute the basic unit signal joint analysis strategy for each basic unit to obtain the absorption coefficient change value of the local area of ​​each basic unit at wavelength λ1;

[0024] 8) Repeat the signal inspection strategy once at wavelength λ2 to obtain the second measured brain diffuse light signal at wavelength λ2, use software to convert the brain diffuse light signals of all channels at wavelength λ2 into time-resolved reflectance signals, perform the basic unit signal joint analysis strategy on each basic unit, and obtain the absorption coefficient change value of the local area of ​​each basic unit at wavelength λ2;

[0025] 9) Absorption coefficient change Δu in local area a,j The relationship between (λ) and blood oxygen concentration is expressed as:

[0026] Δu a,j (λ)=ε HbO (λ)Δ[HbO] j +ε HbR (λ)Δ[HbR] j (1)

[0027] Where: ε HbO (λ) is the molar extinction coefficient of oxygenated hemoglobin at wavelength λ, ε HbR (λ) is the molar extinction coefficient of deoxyhemoglobin at wavelength λ, Δ[HbO] j is the change in oxygenated hemoglobin concentration in the jth local area, Δ[HbR]j is the change of oxygenated hemoglobin concentration in the jth local area;

[0028] The change in absorption coefficient obtained in step 7) and step 8) is expressed by equation (1), and the equations for the same wavelength and local area are combined:

[0029] The simultaneous equations are solved to obtain the changes in the concentrations of oxyhemoglobin and deoxyhemoglobin in all local areas of each basic unit.

[0030] 10) constructing a three-dimensional brain image file based on the AAL template or the Brodmann partition template, setting the brain surface, brain region position, and brain region edge, adjusting the three-dimensional image output layout size, background color, surface transparency, brain region color, brain region edge color, and image resolution, and then further dividing the corresponding brain region in the three-dimensional brain image according to the division of the local brain region detected in step 3), so that each local area displays the normalized RGB value of the absorption coefficient change, and finally obtaining a three-dimensional RGB image of the brain hemodynamic response, and dividing the brain hemodynamic response depth map layer by layer according to the three-dimensional RGB image of the brain hemodynamic response to obtain a 3D brain hemodynamic response map.

[0031] Furthermore, ΔT should be greater than or equal to the time for solving the change in the absorption coefficient minus the time for measuring the diffuse light signal in the brain area at wavelength λ1. and wavelength λ2 to measure the diffuse light signal in the brain area The sum.

[0032] Furthermore, the optical fiber cerebral blood oxygen detection system is composed of a driving circuit, a light source selection circuit, and a transmitting part consisting of a pulse laser 1 and a pulse laser 2, an optical fiber coupler, an optical fiber splitter, an MPO optical fiber jumper, an optical fiber probe, an optical fiber combiner, a single photon counter and a driving circuit 2, and a receiving part consisting of a time-correlated single photon counter (TCSPC) and a SPAD, and a PC signal processing end, wherein the optical fiber probe is composed of a housing, a dual-fiber collimator, a reflector, and a focusing lens; the process of the system performing zero-source-detector distance time-resolved local area detection through a self-detection channel is as follows: the driving circuit transmits a control signal to the light source selection circuit and simultaneously transmits a timing signal to the TCSPC, and the light source selection circuit controls the pulse laser 1 or the pulse laser The optical device 2 emits light of a specified wavelength. The optical signal passes through the optical fiber coupler and the optical fiber splitter, and then reaches the transmitting optical fiber connected to the optical fiber probe via the MPO optical fiber jumper. Then, after passing through the collimating lens and the reflecting surface of the reflector, the beam angle is adjusted by the focusing lens and injected into the cerebral cortex. Part of the light signal is absorbed and scattered by the brain tissue, passes through the same focusing lens, reflecting surface and collimating lens, and returns to the MPO optical fiber jumper via the receiving optical fiber connected to the optical fiber probe, and then reaches the single photon counter through the optical fiber combiner. After the driving circuit 2 detects the signal of the single photon counter, it sends a stop signal to the TCSPC, and then the TCSPC sends a single photon signal to the SPAD. After the signal is converted and amplified by the SPAD, it is sent to the PC signal processing end for detection data processing.

[0033] Furthermore, the basic signal joint analysis strategy specifically includes the following steps:

[0034] The following operations are performed on the channel signals in the basic unit: if the basic unit deploys four fiber optic probes, the time-resolved reflectivity signals of 12 mutual detection channels and 4 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 16m≥n; if the basic unit deploys three fiber optic probes, the time-resolved reflectivity signals of 6 mutual detection channels and 3 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 9m≥n; if the basic unit deploys two fiber optic probes, the time-resolved reflectivity signals of 2 mutual detection channels and 2 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 4m≥n; if the basic unit deploys one fiber optic probe, the time-resolved reflectivity signal of 1 self-detection channel of the fiber optic probe is extracted, and the sliding window method is used to divide the time-resolved reflectivity signal into m time windows, with a sliding step size of Δw, and m≥n; the channel spacing ρ, wavelength λ and the time-dependent average partial path length L of the jth local area at the detection time t of the basic unit are recalculated using the formula i,j (ρ,λ,t):

[0035]

[0036] Where: R(K z ,K l ,N l ,γ z ,γ l ,ρ,t) is the time-resolved reflectivity multilayer approximate diffusion equation of the channel, K z ,K l ,γ z ,γ l ,N l The solution of the system of equations obtained by solving the approximate diffusion equation using the eigenfunction method;

[0037] The time-resolved reflectivity signals of the first and second measurements of each channel are subjected to curve inversion and the absorption change equations are solved, and the absorption change values ​​of the obtained linear equations of all channels are solved.

[0038] Further, the curve inversion and absorption change equation solution specifically include the following steps:

[0039] Measure the system response function s(t), take an initial time-resolved reflectivity signal R0(ρ,λ,t), and use convolution to obtain the optimization curve y(t), that is, Calculate the mean squared fitting error:

[0040] ∑(R m (ρ,λ,t)-y(t)) 2 (4)

[0041] Where: R m (ρ,λ,t) is the measured time-resolved reflectivity signal; The sum of squares of the mean square fitting error is minimized by the least squares method, and the number of iterations is set to hundreds of thousands, and the threshold is T. When the sum of squares of the mean square fitting error is less than the threshold T, the fitting time-resolved reflectivity signal R(ρ,λ,t) is obtained for the subsequent theoretical calculation of the absorption coefficient change.

[0042] The width of each time window of the fitted time-resolved reflectivity signal R(ρ,λ,t) is set to Δt, and the first and last time windows are discarded. In order to describe the small independent absorption coefficient changes in the local area, the modified Lambert-Beer law is used to convert the number of detected photons in each time window into:

[0043]

[0044] Where: R g is the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the second measurement, R 0,gis the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the first measurement, Δu a,j (λ) is the change in the absorption coefficient of the jth local area in the second measurement relative to the first measurement, L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local area in the gth time window with the channel spacing ρ of the basic unit i and the wavelength λ, and the calculation formula is:

[0045]

[0046] Where: L i,j (ρ, λ, t) is the channel spacing ρ of the basic unit i, the wavelength λ and the time-dependent average partial path length of the jth local area at the detection time t;

[0047] Transform formula (5) into:

[0048] The linear equation system consisting of the relation (7) corresponding to all time windows is obtained:

[0049]

[0050] Where: L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local region in the gth time window with channel spacing ρ and wavelength λ of basic unit i.

[0051] Further, the absorption change value is solved, specifically comprising the following steps:

[0052] The linear equations of all channels obtained by solving the curve inversion and absorption variation equations are combined into one linear equation:

[0053]

[0054] Where: is the time-dependent average partial path length of the jth local region of the gth time window of the cth channel of the basic unit i, R c,g is the number of photons detected in the g-th time window of the c-th channel for the second measurement, R c,0,g is the number of photons detected in the g-th time window of the c-th channel for the first measurement, ρ c is the spacing of the cth channel, λ is the main wavelength of the near-infrared pulse light;

[0055] The coefficient matrix of formula (9) is obtained:

[0056]

[0057] If the rank of the coefficient matrix is ​​equal to the rank of the augmented matrix, and both are equal to the number of local regions n, then the absorption coefficient change Δu is solved a,j (λ), otherwise, according to the previously obtained absorption coefficient change value, Δu is predicted based on the Swin-Transformer prediction model. a,j (λ) is used for prediction.

[0058] Beneficial effects of the present invention: Compared with the traditional single-channel time-domain fNIRS detection with a large source-detector distance, the present invention adopts the zero-spacing source-detector fiber probe of the fiber optic cerebral blood oxygen detection system to realize the joint analysis of multi-channel light diffuse signals from surface to domain, improves the lateral and longitudinal spatial resolution of the local area of ​​cerebral tissue, and obtains a high-resolution 3D atlas of cerebral blood oxygen dynamic response. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0060] Figure 1 Schematic diagram of the detection network arrangement: (a) Schematic diagram of the hexagonal honeycomb network arrangement with triangles as the basic units; (b) Schematic diagram of the rectangular network arrangement with rectangles as the basic units.

[0061] Figure 2 It is a schematic diagram of the process of obtaining a 3D atlas of cerebral hemooxygen dynamics response according to the present invention.

[0062] Figure 3 Schematic diagram of the positions of the optical fiber probes and the corresponding brain regions placed in the embodiment.

[0063] Figure 4 Schematic diagram of the process of obtaining a 3D atlas of cerebral hemodynamic response through a detection network in an embodiment: (a) a tiled diagram of the positions of the optical fiber probes placed in the embodiment and the corresponding brain regions; (b) a schematic diagram of time-resolved multi-channel local area detection within a basic unit; (c) a 3D atlas of cerebral hemodynamic response obtained by the detection network.

[0064] Figure 5 It is a flowchart for constructing a high-resolution 3D atlas of cerebral hemooxygen dynamic response. DETAILED DESCRIPTION

[0065] In order to deepen the understanding of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. The embodiments are only used to explain the present invention and do not limit the protection scope of the present invention.

[0066] As shown in the figure, the method for constructing a high-resolution cerebral hemodynamic response 3D atlas of this embodiment includes the following steps:

[0067] 1) The dorsolateral prefrontal cortex (Brodmann area 9) and the frontal pole area (Brodmann area 10) were selected as the brain areas to be detected, and the detection network was used to surround the area of ​​the scalp corresponding to the brain area (the widest part of the dorsolateral prefrontal cortex was about 58 mm, and the widest part of the frontal pole area was about 31 mm). The detection network was pre-selected as a rectangular network with a size of 65 mm x 58 mm with a rectangle as the basic unit, and its basic unit was a rectangle of 10 mm x 10 mm, or a hexagonal honeycomb network with a perimeter of 25.6 mm with an equilateral triangle as the basic unit, and its basic unit was an equilateral triangle with a side length of 1.07 mm; each intersection in the detection network was a probe layout point, which was used to place the optical fiber probe of the optical fiber brain blood oxygen detection system;

[0068] 2) The calculation formula for the channel density of the rectangular network and the hexagonal honeycomb network is:

[0069]

[0070] Where: is the number of channels of the ith basic unit, S i (j) is the surface area of ​​the jth brain region covered by the ith basic unit. i When (j) is zero, let S i (j) is 1;

[0071] 29 and 36 optical fiber probes are set in the above rectangular network and hexagonal honeycomb network respectively, and the channel density of the rectangular network and the hexagonal honeycomb network is calculated. The channel density of the rectangular network is greater than that of the hexagonal honeycomb network, so the rectangular network is used as the detection network;

[0072] Define the optical path where the optical fiber probe emits light into the brain area and absorbs and scatters it, and then receives it from another optical fiber probe as the mutual detection channel, and the optical path where the optical fiber probe emits light into the brain area and absorbs and scatters it, and then receives it from the same optical fiber probe as the self-detection channel. The channel spacing ρ is the spacing between the transmitting fiber and the receiving fiber, i.e., 1.414 cm or 1 cm. Number all the fiber probe deployment points, set the time resolution of the single-photon avalanche diode coupled to the fiber probe to 50 ps, ​​and deploy 19 fiber probes at probe deployment points with different numbers.

[0073] 3) The brain area below each basic unit is taken as a rectangular area, the length, width and height of which are all 3 cm, and it is divided into 28 local areas, which are divided into 27 measurement areas and other areas. The thickness of each measurement area is 5 mm, and the remaining area after removing the measurement area is the other area;

[0074] 4) Seven parameters are set for each local area, including length, width, thickness, refractive index, absorption coefficient, scattering coefficient and anisotropy factor. On this basis, the number of photons is simulated by using the voxel-based time-resolved Monte Carlo simulation method. 8 Using Sobel sampling, we calculate the time-dependent average partial path length L of the jth local area for different channel spacings ρ, wavelengths λ and detection time t for different basic units i. i,j (ρ,λ,t);

[0075] 5) The signal inspection strategy is as follows: select one of the fiber optic probes as the core probe, drive the laser diode coupled to the core probe to emit near-infrared pulse light with a light intensity of 40mW and a main wavelength of λ, and at the same time, let the fiber optic probes in all basic units including the core probe receive the diffuse light signal of the brain area. The time for receiving the signal is set to 5ns. After the signal is received, put all the fiber optic probe numbers that receive the diffuse light signal of the brain area into the waiting queue in clockwise order, take out the head of the queue, and let the fiber optic probe with the corresponding head number be the core probe. Repeat the above operation until there is no fiber optic probe number in the queue. If there is a deployed fiber optic probe that does not emit pulse light as a core probe, let it be the core probe and repeat the above operation until all fiber optic probes have emitted pulse light as core probes; the time T for measuring a diffuse light signal of the brain area is calculated. λ The main wavelength of the pulse light of the laser diode is set to 125ns, the main wavelength of the pulse light of the laser diode is set to 735nm, the signal inspection strategy is executed once, and the first measured brain area diffuse light signal with a wavelength of 735nm is obtained, and then the main wavelength of the pulse light of the laser diode is adjusted to 905nm, the signal inspection strategy is executed once, and the first measured brain area diffuse light signal with a wavelength of 905nm is obtained;

[0076] 6) After an interval of ΔT, the signal inspection strategy is repeated once at a wavelength of 735 nm to obtain the second measured brain diffuse light signal at a wavelength of 735 nm, and the brain diffuse light signals of all channels at a wavelength of 735 nm are converted into time-resolved reflectance signals using software;

[0077] 7) Execute the basic unit signal joint analysis strategy for each basic unit to obtain the absorption coefficient change value of the local area of ​​each basic unit at a wavelength of 735nm;

[0078] 8) Repeat the signal inspection strategy once at a wavelength of 905 nm to obtain the second measured brain diffuse light signal at a wavelength of 905 nm, use software to convert the brain diffuse light signals of all channels at a wavelength of 905 nm into time-resolved reflectance signals, perform the basic unit signal joint analysis strategy on each basic unit, and obtain the absorption coefficient change value of the local area of ​​each basic unit at a wavelength of 905 nm;

[0079] 8) Absorption coefficient change Δu in local area a,j The relationship between (λ) and blood oxygen concentration is expressed as:

[0080] Δu a,j (λ)=ε HbO (λ)Δ[HbO] j +ε HbR (λ)Δ[HbR] j (1)

[0081] Where: ε HbO (λ) is the molar extinction coefficient of oxygenated hemoglobin at wavelength λ, ε HbR (λ) is the molar extinction coefficient of deoxyhemoglobin at wavelength λ, Δ[HbO] j is the change in oxygenated hemoglobin concentration in the jth local area, Δ[HbR] j is the change of oxygenated hemoglobin concentration in the jth local area;

[0082] The change in absorption coefficient obtained in step 7) and step 8) is expressed by equation (1), and the equations for the same wavelength and local area are combined:

[0083] The simultaneous equations are solved to obtain the changes in the concentrations of oxyhemoglobin and deoxyhemoglobin in all local areas of each basic unit.

[0084] 9) constructing a three-dimensional brain image file of a Brodmann partition template, setting the brain surface, brain region position, and brain region edge, adjusting the three-dimensional image output layout size, background color, surface transparency, brain region color, brain region edge color, and image resolution, and then further dividing the corresponding brain region in the three-dimensional brain image according to the division of the local brain region detected in step 3), so that each local area displays the normalized RGB value of the absorption coefficient change, and finally obtaining a three-dimensional RGB image of cerebral hemooxygen dynamics response, and dividing the cerebral hemooxygen dynamics response depth map layer by layer according to the three-dimensional RGB image of cerebral hemooxygen dynamics response to obtain a 3D atlas of cerebral hemooxygen dynamics response.

[0085] Specifically, ΔT should be greater than or equal to the time to solve the change in the absorption coefficient minus the time to measure the diffuse light signal of the brain area at wavelength λ1 and wavelength λ2 to measure the diffuse light signal in the brain area According to the software's computing environment, the estimated time for solving the absorption coefficient change is 2-5ms, so ΔT is 6ms.

[0086] Specifically, the above-mentioned optical fiber cerebral blood oxygen detection system is composed of a driving circuit, a light source selection circuit, and a transmitting part consisting of a pulse laser 1 and a pulse laser 2, an optical fiber coupler, an optical fiber splitter, an MPO optical fiber jumper, an optical fiber probe, an optical fiber combiner, a single photon counter and a driving circuit 2, and a receiving part consisting of a time-correlated single photon counter (TCSPC) and a SPAD, and a PC signal processing end, wherein the optical fiber probe is composed of a shell, a dual-fiber collimator, a reflector and a focusing lens; the process of the above-mentioned system performing zero-source-detector distance time-resolved local area detection through a self-detection channel is as follows: the driving circuit transmits a control signal to the light source selection circuit and transmits a timing signal to the TCSPC at the same time, and the light source selection circuit makes the pulse laser 1 or the pulse laser Device 2 emits light of a specified wavelength. The optical signal passes through the optical fiber coupler and the optical fiber splitter, and then reaches the transmitting optical fiber connected to the optical fiber probe via the MPO optical fiber jumper. Then, after passing through the collimating lens and the reflecting surface of the reflector, the beam angle is adjusted by the focusing lens and injected into the cerebral cortex. Part of the light signal is absorbed and scattered by the brain tissue, passes through the same focusing lens, reflecting surface and collimating lens, and returns to the MPO optical fiber jumper via the receiving optical fiber connected to the optical fiber probe, and then reaches the single photon counter through the optical fiber combiner. After detecting the signal of the single photon counter, the driving circuit 2 sends a stop signal to the TCSPC, and then the TCSPC sends a single photon signal to the SPAD. After the signal is converted and amplified by the SPAD, it is sent to the PC signal processing end for detection data processing.

[0087] Specifically, the above-mentioned basic element signal joint analysis strategy specifically includes the following steps:

[0088] The following operations are performed on the channel signals in the basic unit: if the basic unit deploys four fiber optic probes, the time-resolved reflectivity signals of 12 mutual detection channels and 4 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into 3 time windows; if the basic unit deploys three fiber optic probes, the time-resolved reflectivity signals of 6 mutual detection channels and 3 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into 4 time windows; if the basic unit deploys two fiber optic probes, the time-resolved reflectivity signals of 2 mutual detection channels and 2 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into 7 time windows; if the basic unit deploys one fiber optic probe, the time-resolved reflectivity signal of 1 self-detection channel of the fiber optic probe is extracted, and the time-resolved reflectivity signal is divided into 28 time windows using the sliding window method with a sliding step of 100 ps. The channel spacing ρ, wavelength λ and the time-dependent average partial path length L of the jth local area at the detection time t of the basic unit are recalculated using the formula i,j (ρ,λ,t):

[0089]

[0090] Where: R(K z ,K l ,N l ,γ z ,γ l ,ρ,t) is the time-resolved reflectivity multilayer approximate diffusion equation of the channel, K z ,K l ,γ z ,γ l ,N l The solution of the system of equations obtained by solving the approximate diffusion equation using the eigenfunction method;

[0091] The time-resolved reflectivity signals of the first and second measurements of each channel are subjected to curve inversion and the absorption change equations are solved, and the absorption change values ​​of the obtained linear equations of all channels are solved.

[0092] Specifically, the above curve inversion and absorption change equation solution specifically include the following steps:

[0093] Measure the system response function s(t), take an initial time-resolved reflectivity signal R0(ρ,λ,t), and use convolution to obtain the optimization curve y(t), that is, Calculate the mean squared fitting error:

[0094] ∑(R m (ρ,λ,t)-y(t)) 2 (4)

[0095] Where: R m (ρ,λ,t) is the measured time-resolved reflectivity signal; The sum of squares of the mean square fitting error is minimized by the least squares method, and the number of iterations is set to hundreds of thousands, and the threshold is 0.001. When the sum of squares of the mean square fitting error is less than the threshold 0.001, the fitting time-resolved reflectivity signal R(ρ,λ,t) is obtained for the subsequent theoretical calculation of the absorption coefficient change.

[0096] The width of each time window of the fitted time-resolved reflectivity signal R(ρ,λ,t) is set to 500ps, and the first and last time windows are discarded. In order to describe the small independent absorption coefficient changes in the local area, the modified Lambert-Beer law is used to convert the number of detected photons in each time window into:

[0097]

[0098] Where: R g is the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the second measurement, R 0,gis the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the first measurement, Δu a,j (λ) is the change in the absorption coefficient of the jth local area in the second measurement relative to the first measurement, L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local area in the gth time window with the channel spacing ρ of the basic unit i and the wavelength λ, and the calculation formula is:

[0099]

[0100] Where: L i,j (ρ, λ, t) is the channel spacing ρ of the basic unit i, the wavelength λ and the time-dependent average partial path length of the jth local area at the detection time t;

[0101] Transform formula (5) into:

[0102] The linear equation system consisting of the relation (7) corresponding to all time windows is obtained:

[0103]

[0104] Where: L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local region in the gth time window with channel spacing ρ and wavelength λ of basic unit i.

[0105] Specifically, the above absorption change value is solved, specifically comprising the following steps:

[0106] The linear equations of all channels obtained by solving the curve inversion and absorption variation equations are combined into one linear equation:

[0107]

[0108] Where: is the time-dependent average partial path length of the jth local region of the gth time window of the cth channel of the basic unit i, R c,g is the number of photons detected in the g-th time window of the c-th channel for the second measurement, R c,0,g is the number of photons detected in the g-th time window of the c-th channel for the first measurement, ρ c is the spacing of the cth channel, λ is the main wavelength of the near-infrared pulse light;

[0109] The coefficient matrix of formula (9) is obtained:

[0110]

[0111] If the rank of the coefficient matrix is ​​equal to the rank of the augmented matrix, and both are equal to the number of local regions 28, then the absorption coefficient change Δu is solved a,j (λ), otherwise, according to the previously obtained absorption coefficient change value, Δu is predicted based on the Swin-Transformer prediction model. a,j (λ) is used for prediction.

[0112] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and the specification only describe the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A method for constructing a high-resolution 3D atlas of cerebral hemodynamic response, characterized in that: The method steps include: 1) Select n brain A brain region to be detected is used, and a detection network is used to surround the area of ​​the cerebral scalp corresponding to the brain region. The detection network is a rectangular network with a size of MxN with a rectangle as the basic unit or a hexagonal honeycomb network with a perimeter of C with an equilateral triangle as the basic unit. For the rectangular network, the basic unit is set to be a rectangle of M1xN2, where M1 is less than or equal to M, and N2 is less than or equal to N. For the hexagonal honeycomb network, the basic unit is set to be an equilateral triangle with a perimeter of C1, where C1 is less than or equal to C. Each intersection in the detection network is a probe arrangement point for placing the optical fiber probe of the optical fiber brain blood oxygen detection system. The optical path emitted by the optical fiber probe into the brain region, which undergoes absorption and scattering, and is received by another optical fiber probe is defined as a mutual detection channel, and the optical path emitted by the optical fiber probe into the brain region, which undergoes absorption and scattering, and is received by the same optical fiber probe is defined as a self-detection channel. The channel spacing ρ is the spacing between the transmitting optical fiber and the receiving optical fiber; 2) The calculation formula for the channel density of the rectangular network and the hexagonal honeycomb network is: Where: is the number of channels of the ith basic unit, S i (j) is the surface area of ​​the jth brain region covered by the ith basic unit. i When (j) is zero, let S i (j) is 1; A number of optical fiber probes are respectively set in the rectangular network and the hexagonal honeycomb network. If the channel density of the rectangular network is greater than or equal to the channel density of the hexagonal honeycomb network, the rectangular network is used as the detection network, otherwise the hexagonal honeycomb network is used as the detection network; 3) All fiber optic probe arrangement points are numbered, the time resolution of the single-photon avalanche diode coupled to the fiber optic probe is set to T, and R fiber optic probes are deployed at probe arrangement points with different numbers; if a rectangular network is used, the brain area below each basic unit is taken as a rectangular body area, the length, width and height of the rectangular body area are not greater than three times the side length of the basic unit and not less than twice the side length of the basic unit, and it is divided into n local areas, and the local areas are divided into n-1 measurement areas and other areas, the thickness of each measurement area is not less than 1mm and not more than 1cm, and the remaining area excluding the measurement area is the other area; if a hexagonal honeycomb network is used, the brain area below each basic unit is taken as a regular triangular prism area, the edge length of the regular triangular prism area is not greater than three times the side length of the basic unit and not less than twice the side length of the basic unit, and it is divided into n local areas, and the local areas are divided into n-1 measurement areas and other areas, the thickness of each measurement area is not less than 1mm and not more than 1cm, and the remaining area excluding the measurement area is the other area; 4) Seven parameters are set for each local area, including length, width, thickness, refractive index, absorption coefficient, scattering coefficient and anisotropy factor. On this basis, the number of photons is simulated by using the voxel-based time-resolved Monte Carlo simulation method. 8 Using Sobel sampling, we calculate the time-dependent average partial path length L of the jth local area for different channel spacings ρ, wavelengths λ and detection time t for different basic units i. i,j (ρ,λ,t); 5) The signal inspection strategy is as follows: select one of the fiber optic probes as the core probe, drive the laser diode coupled to the core probe to emit near-infrared pulse light with a light intensity of S and a main wavelength of λ, and at the same time, let the fiber optic probes in all basic units including the core probe receive the diffuse light signal of the brain area. The time for receiving the signal is set to T r , T r is an integer multiple of T. After the signal is received, all fiber optic probe numbers that receive diffuse light signals from the brain area are placed in a waiting queue in clockwise order, the head of the queue is taken out, and the fiber optic probe corresponding to the head of the queue is used as the core probe. The above operation is repeated until there is no fiber optic probe number in the queue. If there is a deployed fiber optic probe that does not emit pulse light as a core probe, it is used as a core probe to repeat the above operation until all fiber optic probes have emitted pulse light as core probes. Calculate the time T of measuring a diffuse light signal in the brain area λ =T r xR, set the main wavelength of the pulse light of the laser diode to λ1, execute the signal inspection strategy once, obtain the first measured brain area diffuse light signal of wavelength λ1, then adjust the main wavelength of the pulse light of the laser diode to λ2, execute the signal inspection strategy once, and obtain the first measured brain area diffuse light signal of wavelength λ2; 6) After an interval of ΔT, the signal inspection strategy is repeatedly executed once at wavelength λ1 to obtain the second measured brain diffuse light signal at wavelength λ1, and the brain diffuse light signals of all channels at wavelength λ1 are converted into time-resolved reflectance signals using software; 7) Execute the basic unit signal joint analysis strategy for each basic unit to obtain the absorption coefficient change value of the local area of ​​each basic unit at wavelength λ1; 8) Repeat the signal inspection strategy once at wavelength λ2 to obtain the second measured brain diffuse light signal at wavelength λ2, use software to convert the brain diffuse light signals of all channels at wavelength λ2 into time-resolved reflectance signals, perform the basic unit signal joint analysis strategy on each basic unit, and obtain the absorption coefficient change value of the local area of ​​each basic unit at wavelength λ2; 9) Absorption coefficient change Δu in local area a,j The relationship between (λ) and blood oxygen concentration is expressed as: Thu a,j (λ)=e HbO (λ)Δ[HbO] j +e HbR (λ)Δ[HbR] j (1) Where: ε HbO (λ) is the molar extinction coefficient of oxygenated hemoglobin at wavelength λ, ε HbR (λ) is the molar extinction coefficient of deoxyhemoglobin at wavelength λ, Δ[HbO] j is the change in oxygenated hemoglobin concentration in the jth local area, Δ[HbR] j is the change of oxygenated hemoglobin concentration in the jth local area; The change in absorption coefficient obtained in step 7) and step 8) is expressed by equation (1), and the equations for the same wavelength and local area are combined: The simultaneous equations are solved to obtain the changes in the concentrations of oxyhemoglobin and deoxyhemoglobin in all local areas of each basic unit. 10) constructing a brain three-dimensional image file based on the AAL template or the Brodmann partition template, setting the brain surface, brain region position, and brain region edge, adjusting the three-dimensional image output layout size, background color, surface transparency, brain region color, brain region edge color, and image resolution, and then further dividing the corresponding brain region in the brain three-dimensional image according to the division of the local area of ​​the detected brain region in step 3), so that each local area displays the normalized RGB value of the absorption coefficient change, and finally obtaining a cerebral hemooxygen dynamics response three-dimensional RGB image, according to the cerebral hemooxygen dynamics response three-dimensional RGB image, dividing the cerebral hemooxygen dynamics response depth map layer by layer, and obtaining a cerebral hemooxygen dynamics response 3D atlas.

2. The method for constructing a high-resolution 3D atlas of cerebral hemodynamic response according to claim 1, characterized in that: ΔT should be greater than or equal to the time to solve the change in absorption coefficient minus the time to measure the diffuse light signal of the brain area at wavelength λ1 and wavelength λ2 to measure the diffuse light signal in the brain area The sum.

3. The method for constructing a high-resolution 3D atlas of cerebral hemodynamic response according to claim 1, characterized in that: The optical fiber cerebral blood oxygen detection system is composed of a driving circuit, a light source selection circuit, a transmitting part composed of a pulse laser 1 and a pulse laser 2, an optical fiber coupler, an optical fiber splitter, an MPO optical fiber jumper, an optical fiber probe, an optical fiber combiner, a single photon counter and a driving circuit 2, and a receiving part composed of a time-correlated single photon counter (TCSPC) and a SPAD, and a PC signal processing end, wherein the optical fiber probe is composed of a housing, a dual-fiber collimator, a reflector and a focusing lens; the process of the system performing zero-source-detector distance time-resolved local area detection through a self-detection channel is as follows: the driving circuit transmits a control signal to the light source selection circuit and transmits a timing signal to the TCSPC at the same time, and the light source selection circuit makes the pulse laser 1 or the pulse laser 2 The light of the specified wavelength is emitted, and the optical signal passes through the optical fiber coupler and the optical fiber splitter, and then reaches the transmitting optical fiber connected to the optical fiber probe via the MPO optical fiber jumper, and then passes through the collimating lens and the reflecting surface of the reflector, and the beam angle is adjusted by the focusing lens and injected into the cerebral cortex. Part of the light signal scattered by the brain tissue is absorbed, and after passing through the same focusing lens, reflecting surface and collimating lens, it returns to the MPO optical fiber jumper via the receiving optical fiber connected to the optical fiber probe, and then reaches the single photon counter through the optical fiber combiner. After the driving circuit 2 detects the signal of the single photon counter, it sends a stop signal to the TCSPC, and then the TCSPC sends the single photon signal to the SPAD. The signal is converted and amplified by the SPAD and then sent to the PC signal processing end for detection data processing.

4. The method for constructing a high-resolution 3D atlas of cerebral hemodynamic response according to claim 1, characterized in that: The basic signal joint analysis strategy specifically includes the following steps: The following operations are performed on the channel signals in the basic unit: if the basic unit deploys four fiber optic probes, the time-resolved reflectivity signals of 12 mutual detection channels and 4 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 16m≥n; if the basic unit deploys three fiber optic probes, the time-resolved reflectivity signals of 6 mutual detection channels and 3 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 9m≥n; if the basic unit deploys two fiber optic probes, the time-resolved reflectivity signals of 2 mutual detection channels and 2 self-detection channels between the fiber optic probes are extracted, and each time-resolved reflectivity signal is divided into m time windows, 4m≥n; if the basic unit deploys one fiber optic probe, the time-resolved reflectivity signal of 1 self-detection channel of the fiber optic probe is extracted, and the sliding window method is used to divide the time-resolved reflectivity signal into m time windows, with a sliding step size of Δw, and m≥n; the channel spacing ρ, wavelength λ and the time-dependent average partial path length L of the jth local area at the detection time t of the basic unit are recalculated using the formula i,j (ρ,λ,t): Where: R(K z ,K l ,N l ,γ z ,γ l ,ρ,t) is the time-resolved reflectivity multilayer approximate diffusion equation of the channel, K z ,K l ,γ z ,γ l ,N l The solution of the system of equations obtained by solving the approximate diffusion equation using the eigenfunction method; The time-resolved reflectivity signals of the first and second measurements of each channel are subjected to curve inversion and the absorption change equations are solved, and the absorption change values ​​of the obtained linear equations of all channels are solved.

5. The method for constructing a high-resolution 3D atlas of cerebral hemodynamic response according to claim 2, characterized in that: The curve inversion and absorption change equation solution specifically include the following steps: Measure the system response function s(t), take an initial time-resolved reflectivity signal R0(ρ,λ,t), and use convolution to obtain the optimization curve y(t), that is, Calculate the mean squared fitting error: ∑(R m (ρ,λ,t)-y(t)) 2 (4) Where: R m (ρ,λ,t) is the measured time-resolved reflectivity signal; The sum of squares of the mean square fitting error is minimized by the least squares method, and the number of iterations is set to hundreds of thousands, and the threshold is T. When the sum of squares of the mean square fitting error is less than the threshold T, the fitting time-resolved reflectivity signal R(ρ,λ,t) is obtained for the subsequent theoretical calculation of the absorption coefficient change. The width of each time window of the fitted time-resolved reflectivity signal R(ρ,λ,t) is set to Δt, and the first and last time windows are discarded. In order to describe the small independent absorption coefficient changes in the local area, the modified Lambert-Beer law is used to convert the number of detected photons in each time window into: Where: R g is the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the second measurement, R 0, g is the number of photons detected in the gth time window of the fitted time-resolved reflectivity signal of the first measurement, Δu a,j (λ) is the change in the absorption coefficient of the jth local area in the second measurement relative to the first measurement, L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local area in the gth time window with the channel spacing ρ of the basic unit i and the wavelength λ, and the calculation formula is: Where: L i,j (ρ, λ, t) is the channel spacing ρ of the basic unit i, the wavelength λ and the time-dependent average partial path length of the jth local area at the detection time t; Transform formula (5) into: The linear equation system consisting of the relation (7) corresponding to all time windows is obtained: Where: L g,i,j (ρ,λ) is the time-dependent average partial path length of the jth local region in the gth time window with channel spacing ρ and wavelength λ of basic unit i.

6. The method for constructing a high-resolution 3D atlas of cerebral hemodynamic response according to claim 2, characterized in that: The absorption change value solution specifically comprises the following steps: The linear equations of all channels obtained by solving the curve inversion and absorption variation equations are combined into one linear equation: Where: is the time-dependent average partial path length of the jth local region of the gth time window of the cth channel of the basic unit i, R c,g is the number of photons detected in the g-th time window of the c-th channel for the second measurement, R c,0,g is the number of photons detected in the g-th time window of the c-th channel for the first measurement, ρ c is the spacing of the cth channel, λ is the main wavelength of the near-infrared pulse light; The coefficient matrix of formula (9) is obtained: If the rank of the coefficient matrix is ​​equal to the rank of the augmented matrix, and both are equal to the number of local regions n, then the absorption coefficient change Δu is solved a,j (λ), otherwise, according to the previously obtained absorption coefficient change value, Δu is predicted based on the Swin-Transformer prediction model. a,j (λ) is used for prediction.

Images