Intelligent Adjustment and Monitoring Method and System for Skin Phototherapy Dosage

Abstract

This invention relates to the field of phototherapy technology, and more particularly to a method and system for intelligent adjustment and monitoring of phototherapy dose. An initial dose plan is generated based on the initial physiological state information of the skin area to be treated. The area is divided into multiple spatial sub-regions, and baseline tolerance thresholds and conduction coefficients are recorded. The response intensity of each sub-region is monitored in parallel. The cross-influence between adjacent sub-regions is calculated using a diffusion-response coupling equation. A cumulative matrix of cross-influence is constructed and superimposed on the baseline tolerance threshold to obtain the predicted tolerance threshold and remaining tolerance margin. While maintaining a constant overall cumulative energy density, the differentiated light intensity adjustment amount is calculated based on the remaining tolerance margin and then executed. This achieves intelligent and personalized phototherapy dose adjustment, improving treatment safety and effectiveness.

Description

Technical Field

[0001] This invention relates to the field of skin phototherapy technology, and in particular to a method and system for intelligent adjustment and monitoring of skin phototherapy dosage. Background Technology

[0002] Phototherapy is a common technique for treating skin diseases such as psoriasis and vitiligo using light of specific wavelengths. In current clinical practice, treatment plans typically rely on the physician's experience or are based on macroscopic parameters such as the patient's skin type and minimum erythema dose, setting a single fixed dose and applying uniform light intensity across the entire treatment area. During treatment, medical staff primarily rely on visual observation or patient feedback to judge skin reactions, lacking real-time quantitative monitoring of changes in local physiological states. Dosage adjustments are often based on overall efficacy assessments or reactive responses to adverse reactions.

[0003] Human skin is not a homogeneous tissue; the melanin content, stratum corneum thickness, and blood flow distribution vary significantly across different anatomical locations, leading to substantial differences in local tolerance. Uniform irradiation can easily cause erythema, edema, or even blisters in some areas due to excessive doses, while other areas may not reach the treatment threshold due to insufficient doses, resulting in inconsistent therapeutic effects. This spatial heterogeneity is particularly pronounced in heterogeneous lesions such as psoriatic plaques and vitiligo leukoplakia, where current methods struggle to achieve independent and precise control of each sub-region.

[0004] Complex physiological coupling relationships exist between skin regions. When an area experiences an acute inflammatory response due to excessive light exposure, the released inflammatory mediators and heat stress signals can diffuse to adjacent tissues via intercellular fluid or microcirculation, significantly reducing the radiation tolerance threshold of nearby areas. Current phototherapy systems completely ignore this cross-regional transmission effect, leading to a situation where an overreaction at a single site can indirectly cause a sudden decrease in tolerance in adjacent areas, triggering a chain reaction of damage. This dynamic interaction mechanism is not incorporated into existing dose control models, which is a significant reason for the unintended amplification of damage during phototherapy. Summary of the Invention

[0005] This invention provides a method and system for intelligent adjustment and monitoring of skin phototherapy dosage, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a method for intelligent adjustment and monitoring of skin phototherapy dosage, comprising: An initial dose plan is generated based on the initial physiological state information of the skin area to be treated; the skin area to be treated is divided into multiple spatial sub-regions, and the baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dose plan. Based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using the diffusion-response coupling equation. The cross-influence represents the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. A cross-influence accumulation matrix is ​​constructed to record the cross-influence between each pair of spatial sub-regions. The cross-influence is then superimposed on the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. While maintaining a constant overall cumulative energy density in the treatment area, the differential light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region and then sent to the phototherapy execution device for execution.

[0007] The skin area to be treated is divided into multiple spatial sub-regions. The baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dose regimen, including: An initial dosage plan is generated based on the initial physiological state information of the skin area to be treated. Based on the rate of change of physiological parameters in the initial physiological state information, the skin region to be treated is adaptively divided into grids, and the division density is increased in regions where the rate of change of physiological parameters is higher than a preset change threshold to obtain multiple spatial sub-regions. Standard test excitation is applied to each spatial sub-region and the test response intensity is collected. The baseline tolerance threshold of each spatial sub-region is determined based on the test response intensity. The response time-series correlation data of adjacent spatial sub-regions are collected, and the conduction coefficient is calculated based on the response propagation delay and the response amplitude attenuation rate. Response intensity sensing units are deployed in each spatial sub-region to synchronously collect instantaneous physiological response parameters and perform time integration to obtain the response intensity.

[0008] Based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using the diffusion-response coupling equation, including: The central difference operator is used to calculate the spatial gradient of the response intensity in each spatial sub-region in the spatial dimension, and a composite gradient vector is constructed by combining the time rate of change of the response intensity. A diffusion-reaction coupling equation is constructed, which includes a diffusion term and a reaction term. The diffusion term adopts a second-order partial differential form, and the diffusion coefficient is set as the product of the conduction coefficient and the spatial components of the composite gradient vector. The reaction term adopts a first-order kinetic form, and the reaction rate is proportional to the current response intensity. The length of the future prediction time window is adaptively determined based on the ratio of the current response intensity of each spatial sub-region to the baseline tolerance threshold. The larger the ratio, the shorter the length of the future prediction time window. The diffusion-response coupling equation is numerically solved within the future prediction time window using an implicit time discretization method to obtain the predicted response intensity of each spatial sub-region at the end of the future prediction time window. The difference between the predicted response intensity and the current response intensity is calculated as the response intensity increment. When the response intensity increment exceeds the response threshold, it is determined that an over-response has occurred, and the portion of the response intensity increment that exceeds the response threshold is taken as the cross-influence quantity.

[0009] The implicit time discretization method is used to numerically solve the diffusion-reaction coupling equation within a future prediction time window, including: In the spatial sub-region where the magnitude of the composite gradient vector exceeds a set threshold and its adjacent spatial sub-regions, a local mesh refinement is performed to obtain a refined spatial mesh. An adaptive time step is determined based on the mesh size and diffusion coefficient of the refined spatial mesh to meet the numerical stability condition. For the spatial sub-regions located at the boundary of the treatment area, flux boundary conditions are set based on the spatial gradient of the response intensity and the conduction coefficient at the boundary. For the spatial sub-regions located inside the treatment area, periodic boundary conditions are set. The diffusion-reaction coupling equation is spatially discretized based on the encrypted spatial grid. The diffusion term is discretized in time using an implicit Euler scheme, and the reaction term is discretized in time using an explicit scheme. The solution is iteratively solved within the future prediction time window with the adaptive time step until the change in the predicted response intensity of adjacent time steps is less than the convergence threshold, and the predicted response intensity is obtained.

[0010] The predicted tolerance threshold is obtained by superimposing the cross-influence amount onto the baseline tolerance threshold of each spatial sub-region, and the remaining tolerance margin is calculated, including: A cumulative cross-influence matrix is ​​constructed for each spatial sub-region. The row index and column index of the cumulative cross-influence matrix represent the source spatial sub-region and the affected spatial sub-region, respectively. The matrix elements represent the cross-influence between corresponding pairs of spatial sub-regions. The total cross-influence of each spatial sub-region is obtained by summing all the cross-influence in the column vector of the matrix corresponding to each spatial sub-region when it is considered as an affected spatial sub-region. The predicted tolerance threshold is obtained by subtracting the corresponding total cross-influence from the baseline tolerance threshold of each spatial sub-region. The cumulative response load of each spatial sub-region at the current time is collected. The cumulative response load represents the time integral value of the response intensity of the spatial sub-region from the start of treatment to the current time. The difference between the predicted tolerance threshold and the cumulative response load is calculated as the remaining tolerance margin.

[0011] Under the constraint of maintaining a constant overall cumulative energy density in the treatment area, the differentiated light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region, including: A constrained optimization problem is constructed with the goal of maximizing the sum of the ratios of the adjusted light intensity to the remaining tolerance margin in each spatial sub-region. The remaining tolerance margin in each spatial sub-region is used as the upper bound of the inequality constraint, and the overall cumulative energy density of the treatment area is kept constant as the equality constraint. The Lagrange multiplier method is used to solve the basic adjustment components of each spatial sub-region. Traverse all matrix elements in the cumulative cross-influence matrix. When the value of a matrix element exceeds a preset propagation threshold, extract the index of the source spatial sub-region and the index of the affected spatial sub-region corresponding to the matrix element. Calculate the spatial distance between the source spatial sub-region and the affected spatial sub-region. Calculate the propagation compensation component of the spatial sub-region pair. The propagation compensation component is directly proportional to the conduction coefficient between the source spatial sub-region and the affected spatial sub-region and inversely proportional to the spatial distance. For each sub-region of the influence source space, sum up all propagation compensation components of the influence source space to obtain the total propagation compensation component of the sub-region of the influence source space. The total propagation compensation component of each spatial sub-region is superimposed on the corresponding basic adjustment component in the form of a negative value to obtain the differentiated light intensity adjustment amount of each spatial sub-region, which is then sent to the phototherapy execution device for execution.

[0012] Calculating the propagation compensation components of this spatial sub-region pair includes: The basic compensation coefficient is obtained by multiplying the reciprocal of the spatial distance by the transmission coefficient. The excess amount of the matrix element value exceeding the preset propagation threshold is calculated. The excess amount is mapped to a propagation risk factor through the Sigmoid function, such that the larger the excess amount, the closer the propagation risk factor is to the saturation upper limit. The propagation compensation component of the spatial sub-region pair is obtained by multiplying the basic compensation coefficient by the propagation risk factor. A propagation path graph is constructed based on the cumulative cross-influence matrix. The nodes of the propagation path graph represent spatial sub-regions, and the directed edges represent propagation relationships where the cross-influence exceeds the preset propagation threshold. A depth-first search algorithm is used to traverse all propagation paths originating from each influence source spatial sub-region and record the propagation path length. When the propagation path length exceeds the preset hop count threshold, a cascade compensation coefficient proportional to the propagation path length is calculated. The cascade compensation coefficient is then superimposed on the propagation compensation component of the first affected spatial sub-region on the propagation path corresponding to the influence source spatial sub-region.

[0013] A second aspect of the present invention provides a skin phototherapy dose intelligent adjustment and monitoring system, comprising: An initial dose unit is used to generate an initial dose plan based on the initial physiological state information of the skin area to be treated. The partition monitoring unit is used to divide the skin area to be treated into multiple spatial sub-regions, record the baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions, and monitor the response intensity of each spatial sub-region in parallel during the execution of the initial treatment dose plan. The cross-influence unit is used to calculate the cross-influence between adjacent spatial sub-regions based on the response intensity gradient of each spatial sub-region using the diffusion-reaction coupling equation. The cross-influence value characterizes the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. The margin calculation unit is used to construct a cross-influence accumulation matrix to record the cross-influence between each pair of spatial sub-regions, and to superimpose the cross-influence to the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. The adjustment execution unit is used to calculate the differentiated light intensity adjustment amount of each spatial sub-region based on the remaining tolerance margin of each spatial sub-region while maintaining the overall cumulative energy density of the treatment area, and then send it to the phototherapy execution device for execution.

[0014] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0015] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0016] This invention divides the skin area to be treated into multiple spatial sub-regions and records the baseline tolerance threshold and conduction coefficient, achieving refined spatial monitoring. The response intensity of each sub-region is monitored in parallel, capturing the immediate response of the local skin to light exposure. Based on the diffusion-response coupling equation, the cross-influence between adjacent sub-regions is calculated, quantifying the degree to which an excessive response in a sub-region inhibits the future tolerance threshold of the surrounding area, effectively avoiding secondary damage caused by energy diffusion. A cumulative cross-influence matrix is ​​constructed to record the interaction between pairs of sub-regions, which is then superimposed on the baseline tolerance threshold to obtain the predicted tolerance threshold and remaining tolerance margin, making dose adjustment proactive and significantly reducing the risk of burns or overtreatment.

[0017] Under the constraint of maintaining a constant overall cumulative energy density in the treatment area, the method calculates the differential light intensity adjustment based on the remaining tolerance margin of each sub-region, thereby achieving optimized energy redistribution. This method ensures that the total irradiation dose is consistent with the initial plan while flexibly adjusting the local intensity—reducing the intensity in highly sensitive areas and appropriately increasing it in low-sensitivity areas, thus improving light uniformity while avoiding adverse reactions such as erythema and blisters.

[0018] Real-time monitoring and intelligent adjustment form a closed-loop feedback loop, enabling the phototherapy process to adapt to dynamic changes in skin condition. Initial physiological information serves as a baseline, and subsequent actual responses continuously adjust dosage parameters to avoid delayed damage caused by decreased skin tolerance. At the same time, precise energy distribution shortens the total number of treatment sessions. Attached Figure Description

[0019] Figure 1 A flowchart illustrating a method for intelligent adjustment and monitoring of skin phototherapy dosage; Figure 2 This is a flowchart for calculating the cross-influence of spatial sub-regions based on the diffusion-reaction coupling equation. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0022] Figure 1 This is a flowchart illustrating the intelligent adjustment and monitoring method for skin phototherapy dosage according to an embodiment of the present invention.

[0023] Intelligent adjustment and monitoring methods for skin phototherapy dosage include: An initial dosage plan is generated based on the initial physiological state information of the skin area to be treated. The skin area to be treated is divided into multiple spatial sub-regions. The baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dose scheme. Based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using the diffusion-response coupling equation. The cross-influence represents the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. A cross-influence accumulation matrix is ​​constructed to record the cross-influence between each pair of spatial sub-regions. The cross-influence is then superimposed onto the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. While maintaining a constant overall cumulative energy density in the treatment area, the differential light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region and then sent to the phototherapy execution device for execution.

[0024] In one optional implementation, the skin region to be treated is divided into multiple spatial sub-regions, and the baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dosage protocol, including: An initial dosage plan is generated based on the initial physiological state information of the skin area to be treated. Based on the rate of change of physiological parameters in the initial physiological state information, the skin region to be treated is adaptively divided into grids, and the division density is increased in regions where the rate of change of physiological parameters is higher than a preset change threshold to obtain multiple spatial sub-regions. Standard test excitation is applied to each spatial sub-region and the test response intensity is collected. The baseline tolerance threshold of each spatial sub-region is determined based on the test response intensity. Collect response time-series correlation data of adjacent spatial sub-regions, and calculate the conduction coefficient based on the response propagation delay and response amplitude attenuation rate; Response intensity sensing units are deployed in each spatial sub-region to synchronously collect instantaneous physiological response parameters and perform time integration to obtain the response intensity.

[0025] For example, for the skin area to be treated, physiological parameters such as surface temperature distribution, blood perfusion index, melanin concentration, and stratum corneum thickness are collected using a multispectral imaging device. These physiological parameters constitute the basic dataset of initial physiological state information, used for the generation of subsequent dosage protocols. The initial dosage protocol is generated based on the patient's skin type classification, lesion type, and treatment goals, combined with historical treatment records in the clinical database, to determine the initial light intensity, wavelength combination, and expected irradiation duration. This initial dosage protocol serves as the baseline parameters for treatment, providing a reference baseline for subsequent dynamic adjustments.

[0026] The spatial rate of change of each physiological parameter in the initial physiological state information is calculated. Specifically, for temperature parameters, the temperature gradient between adjacent measurement points is calculated; for blood perfusion index, its spatial distribution non-uniformity is calculated; and for melanin concentration, its spatial distribution standard deviation is extracted. These rates of change are compared with preset thresholds. When the rate of change of a physiological parameter in a certain region exceeds the threshold, it indicates that there is a significant difference in the physiological state of that region, requiring more refined monitoring and control. In these high-rate-change regions, the mesh density is increased, reducing the size of individual spatial sub-regions to 50% to 70% of the original planned size, thereby achieving higher spatial resolution. In regions with low rates of change of physiological parameters, larger sub-region sizes are maintained to reduce computational burden. Through this adaptive partitioning strategy, the skin area to be treated is ultimately divided into dozens to hundreds of spatial sub-regions, each typically ranging in area from 0.5 square centimeters to 2 square centimeters.

[0027] The baseline tolerance threshold for each spatial sub-region was determined using a standard test stimulus. The standard test stimulus employed probing light at a dose below the therapeutic level, with an intensity set to 20% to 30% of the expected therapeutic intensity, for a duration of 5 to 10 seconds. During the application of the standard test stimulus, the test response intensity was acquired in real time by multi-parameter sensors deployed in each spatial sub-region. The test response intensity included changes in local temperature rise, blood flow velocity, and tissue reflectance spectrum. These test response intensities were compared to a preset safe response range. The stimulus intensity at which the test response intensity reached 80% of the upper limit of the safe response range was defined as the baseline tolerance threshold for that spatial sub-region. This baseline tolerance threshold reflects the maximum light intensity that the sub-region can safely tolerate under its current physiological state, providing a safe boundary for dose adjustments during subsequent treatment.

[0028] The conduction coefficient between adjacent spatial sub-regions reflects the spatial propagation characteristics of physiological responses. To obtain the conduction coefficient, a pulsed test stimulus is applied to a specific spatial sub-region, while the response time-series data of that sub-region and its adjacent sub-regions are monitored. The response time-series data includes temperature values, blood perfusion index, and tissue metabolic marker concentrations at each time point. The response propagation delay is determined by analyzing the time delay of the response curves of adjacent sub-regions. The response propagation delay is defined as the time difference between the peak response intensity of adjacent sub-regions, typically between 0.5 and 3 seconds. Simultaneously, the response amplitude decay rate is calculated, which is the ratio of the peak response of the adjacent sub-region to the peak response of the stimulated sub-region. The conduction coefficient is calculated by multiplying the reciprocal of the response propagation delay by the response amplitude decay rate; this coefficient quantifies the efficiency of the diffusion of physiological responses from one sub-region to adjacent sub-regions. Regions with higher conduction coefficients indicate a stronger physiological coupling between the two sub-regions, and an overreaction in one sub-region is more likely to affect the tolerance of adjacent sub-regions.

[0029] During the execution of the initial treatment dosage protocol, response intensity sensing units are deployed in each spatial sub-region for parallel monitoring. These units include miniature thermocouples, optical blood flow sensors, and spectral reflectance probes, which synchronously acquire transient physiological response parameters at a fixed sampling frequency. The sampling frequency is set to 10 to 50 times per second to ensure the capture of dynamic changes in the physiological response. Transient physiological response parameters include transient temperature, transient blood flow velocity, and transient tissue optical properties. To convert these transient parameters into cumulative response intensity, each transient physiological response parameter is integrated over time. The calculation window for the time integration extends from the start of treatment to the current time, and the integration result reflects the cumulative physiological load on that spatial sub-region throughout the treatment process. For example, the temperature response intensity is obtained by integrating the portion of the temperature exceeding the baseline temperature over time; this integral value characterizes the cumulative effect of heat in the tissue. The blood flow response intensity is represented by the integral of the change in blood flow velocity over time, reflecting the cumulative response of the vascular system to light stimulation. This time integration method transforms instantaneous physiological fluctuations into quantifiable response intensity indicators, providing a data foundation for subsequent calculations of cross-influence.

[0030] The response intensity data for each spatial sub-region is stored in time-series format, forming a spatiotemporal distribution matrix of response intensity. The row indices of this matrix correspond to different spatial sub-regions, and the column indices correspond to different time sampling points. By calculating the spatial gradient of this matrix, the non-uniform distribution pattern of response intensity in space can be identified, providing input data for calculating cross-influence. Simultaneously, a conduction coefficient matrix records the conduction coefficient values ​​between all pairs of adjacent spatial sub-regions. This matrix is ​​symmetric, with diagonal elements being zero and off-diagonal elements representing the conduction coefficients between corresponding sub-region pairs. Baseline tolerance thresholds are stored in vector form, with each element of the vector corresponding to the baseline tolerance threshold of a spatial sub-region. These data structures collectively constitute the input parameter set for the subsequent dynamic dose adjustment algorithm, ensuring that the adjustment strategy fully considers the individual differences of each spatial sub-region and the mutual influence between adjacent sub-regions.

[0031] In practical applications, the preset change thresholds for adaptive mesh generation are adjusted according to the specific treatment type. For the treatment of pigmentary diseases, the threshold for the rate of change in melanin concentration is set at 15% to 25% per millimeter; for the treatment of vascular diseases, the threshold for the rate of change in blood perfusion index is set at 10% to 20% per millimeter. These thresholds are set based on statistical analysis of a large amount of clinical data to ensure that the divided spatial sub-regions accurately reflect the physiological heterogeneity of the skin region. The selection of standard test stimulus parameters needs to strike a balance between safety and sensitivity. Too low a stimulus intensity may not elicit a measurable response, while too high a stimulus intensity may place an unnecessary burden on the tissue. The optimal test stimulus parameters are determined through multiple preliminary experiments, ensuring that the testing process accurately assesses the baseline tolerance threshold without interfering with subsequent treatments.

[0032] In one optional implementation, based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using a diffusion-response coupling equation, including: The central difference operator is used to calculate the spatial gradient of the response intensity in each spatial sub-region in the spatial dimension, and a composite gradient vector is constructed by combining the time rate of change of the response intensity. A diffusion-reaction coupling equation is constructed, which includes a diffusion term and a reaction term. The diffusion term adopts a second-order partial differential form, and the diffusion coefficient is set as the product of the conduction coefficient and the spatial components of the composite gradient vector. The reaction term adopts a first-order kinetic form, and the reaction rate is proportional to the current response intensity. The length of the future prediction time window is adaptively determined based on the ratio of the current response intensity of each spatial sub-region to the baseline tolerance threshold. The larger the ratio, the shorter the length of the future prediction time window. The diffusion-response coupling equation is numerically solved within the future prediction time window using an implicit time discretization method to obtain the predicted response intensity of each spatial sub-region at the end of the future prediction time window. The difference between the predicted response intensity and the current response intensity is calculated as the response intensity increment. When the response intensity increment exceeds the response threshold, it is determined that an over-response has occurred, and the portion of the response intensity increment that exceeds the response threshold is taken as the cross-influence quantity.

[0033] Combination Figure 2The flowchart illustrating the calculation of cross-influence between spatial sub-regions based on the diffusion-response coupling equation is presented below. After obtaining the real-time response intensity data of each spatial sub-region, it is necessary to quantify the mutual influence relationship between adjacent sub-regions. This influence relationship is not a simple linear superposition, but a complex dynamic process formed by the combined effects of multiple physiological mechanisms such as heat conduction within skin tissue, diffusion of inflammatory factors, and cell signal transduction. To accurately characterize this process, the diffusion-response coupling equation is used to quantitatively calculate the cross-influence between adjacent spatial sub-regions.

[0034] First, the spatial gradient of the response intensity of each spatial sub-region is calculated. Assume the skin region to be treated is divided into... The grid structure, where each grid represents a spatial sub-region, is located at the . Line 1 The spatial subregion of a column is denoted as , at time The response intensity is denoted as The spatial gradients of the response intensity in the horizontal and vertical directions of this spatial sub-region are calculated using the central difference operator. The horizontal spatial gradient is obtained by dividing the difference in response intensity between the left and right adjacent sub-regions by the spatial step size, as shown in the following formula: ,in This represents the physical distance between the center points of adjacent sub-regions in the horizontal direction. The spatial gradient in the vertical direction is calculated in a similar way, with the formula as follows: ,in This represents the physical distance between the center points of adjacent sub-regions in the vertical direction. For spatial sub-regions located at the boundary, since there is no adjacent sub-region on one side, a one-sided difference operator is used for calculation. For example, the horizontal gradient of the sub-region located on the left boundary is calculated using... calculate.

[0035] Based on the obtained spatial gradient, the time rate of change of the response intensity is further calculated. This is done using the current time step. Compared to the previous moment The difference in response intensity is obtained by dividing the time interval; the calculation formula is as follows: ,in This represents the time interval for monitoring and sampling. A composite gradient vector is constructed by combining the two components of the spatial gradient with the rate of change over time. The first two components of this composite gradient vector reflect the spatial non-uniformity of the response intensity distribution, while the third component reflects the rate of evolution of the response intensity over time. Together, they describe the dynamic characteristics of the response state in this spatial sub-region.

[0036] The diffusion-reaction coupling equation comprises two core components: a diffusion term and a reaction term. The diffusion term, expressed in second-order partial differential form, describes the spatial propagation of the response intensity, and its mathematical expression is: ,in This is the diffusion coefficient. This diffusion coefficient is not a fixed constant, but rather varies depending on the spatial sub-region. The diffusion coefficient and the spatial components of the composite gradient vector between the diffusion coefficient and its neighboring sub-regions are dynamically determined. Specifically, the diffusion coefficient is set as the diffusion coefficient. The product of the magnitudes of the spatial components of the composite gradient vector, i.e. This configuration allows the diffusion coefficient to reflect the steepness of the local response intensity gradient; a larger gradient indicates a more significant difference in response between adjacent regions and a stronger diffusion effect.

[0037] The reaction term uses first-order kinetics to describe the self-reinforcing or decaying process of the response intensity, and its mathematical expression is: ,in This is the reaction rate coefficient. This coefficient is proportional to the current response intensity, and is specifically calculated as follows: ,in This is a proportionality constant, usually obtained experimentally based on the physiological characteristics of skin tissue. When the response intensity is high, the response rate coefficient is large, indicating that the response state in this area has a strong self-amplification tendency, which is consistent with the positive feedback mechanism of the inflammatory response of skin tissue after phototherapy stimulation.

[0038] Combining the diffusion and reaction terms yields the complete diffusion-reaction coupling equation. ,in Indicates from the current moment The initial relative time variable. This equation describes the evolution of the response intensity over a future time period and needs to be numerically solved within a specific time window to predict the future response state of each spatial sub-region.

[0039] The length of the future prediction time window is not a fixed value, but is adaptively determined based on the current response state of each spatial sub-region. The current response intensity is calculated. Compared with the baseline tolerance threshold of this sub-region ratio This ratio reflects the saturation level of the current response intensity relative to its tolerance capacity. A larger ratio indicates that the sub-region is closer to its tolerance limit, and the higher the risk of over-response in the future. To promptly capture this high-risk state, an inverse proportional relationship is used to determine the length of the future prediction time window. ,in The baseline forecast time window length is typically set between 30 and 120 seconds. To adjust the index, a value typically ranges from 0.5 to 1.5. When When it is large, Smaller values ​​allow forecasts to focus on the near future, improving sensitivity to warnings of impending overreactions; when When smaller, It is relatively large, allowing for the evaluation of response evolution over a longer time range.

[0040] After determining the length of the future prediction time window, an implicit time discretization method is used to numerically solve the diffusion-reaction coupling equation. The implicit method exhibits better numerical stability compared to the explicit method, and is particularly suitable for cases with large diffusion coefficients or fast reaction rates. The future prediction time window is defined as follows. Divided into There are 1 time steps, each time step having a length of 1. In the first At each time step, the time derivative is discretized using a backward Euler scheme, and the spatial second derivative is discretized using a central difference scheme. For the spatial sub-region... The discretized equation is in the form of ,in Indicates the first At the end of each time step, the spatial sub-region The response intensity is determined. This discrete equation is established for all spatial sub-regions, forming a large linear system of equations. Iterative solvers, such as the conjugate gradient method or multigrid method, are used to solve this system, ultimately yielding the end time of the future prediction window for each spatial sub-region. Predicted response intensity .

[0041] After obtaining the predicted response intensity, the difference between it and the current response intensity is calculated as the response intensity increment. This increment reflects the expected change in the response intensity of this spatial sub-region within the future prediction time window, due to the combined effects of diffusion and reaction mechanisms. A response threshold is pre-set to determine whether overresponse has occurred. This threshold is typically determined based on clinical safety standards and represents the acceptable upper limit of the increase in response intensity. At that time, determine the spatial sub-region There is a risk of overresponse within the future forecast time window. This includes the portion of the response intensity increment that exceeds the response threshold. As the cross-influence of this spatial sub-region on its adjacent sub-regions This cross-influence quantity characterizes the spatial sub-region. The degree to which the overreaction state inhibits the tolerance threshold of adjacent sub-regions within a future time window through a diffusion mechanism. The larger the value, the more significant the negative impact of the sub-region on the surrounding region, and this needs to be given special consideration in subsequent dose adjustments.

[0042] Through the above calculation process, the cross-influence amount is calculated for each spatial sub-region, ultimately forming a cross-influence distribution map describing the internal interaction relationships of the entire skin area to be treated. This distribution map provides a quantitative basis for subsequent construction of the cross-influence accumulation matrix and for differentiated dose adjustment, ensuring that the phototherapy process achieves therapeutic effects while minimizing local over-response and its spread to surrounding areas.

[0043] In one optional implementation, an implicit time discretization method is used to numerically solve the diffusion-reaction coupling equation within a future prediction time window, including: In the spatial sub-region where the magnitude of the composite gradient vector exceeds a set threshold and its adjacent spatial sub-regions, a local mesh refinement is performed to obtain a refined spatial mesh. An adaptive time step is determined based on the mesh size and diffusion coefficient of the refined spatial mesh to meet the numerical stability condition. For the spatial sub-regions located at the boundary of the treatment area, flux boundary conditions are set based on the spatial gradient of the response intensity and the conduction coefficient at the boundary. For the spatial sub-regions located inside the treatment area, periodic boundary conditions are set. The diffusion-reaction coupling equation is spatially discretized based on the encrypted spatial grid. The diffusion term is discretized in time using an implicit Euler scheme, and the reaction term is discretized in time using an explicit scheme. The solution is iteratively solved within the future prediction time window with the adaptive time step until the change in the predicted response intensity of adjacent time steps is less than the convergence threshold, and the predicted response intensity is obtained.

[0044] For example, to improve the solution accuracy and computational efficiency of the diffusion-reaction coupling equation within the future prediction time window, an implicit time discretization method combined with spatial adaptive mesh refinement technology is used for numerical solution. The core of this method lies in refining the local mesh for regions where the response intensity gradient changes drastically, while dynamically adjusting the time step according to the mesh size to ensure numerical stability, and balancing computational accuracy and efficiency through an implicit-explicit hybrid discretization scheme.

[0045] Before performing numerical solutions, the spatial sub-regions that require mesh refinement are first identified. For each spatial sub-region... Calculate its composite gradient vector model The modulus value is compared with a preset gradient threshold. Comparison. When When this occurs, the spatial sub-region and its adjacent spatial sub-regions are determined to be high-gradient regions, requiring local mesh refinement. Gradient threshold. The setting is based on the overall response intensity distribution characteristics of the treatment area, typically taking 1.5 to 2 times the median gradient magnitude of all spatial sub-regions. For identified high-gradient regions, the original spatial sub-regions are further subdivided into smaller, more refined grid units. Specifically, if the horizontal dimension of the original spatial sub-region is... Vertical dimension is During the encryption process, the sub-region is divided into two equal parts along both the horizontal and vertical directions. One, of which The encryption factor is typically 2 or 3. The horizontal dimension of the encrypted grid cell is... Vertical dimension is For spatial sub-regions not identified as high-gradient regions, the original mesh size is kept unchanged, thus forming a non-uniform, refined spatial mesh.

[0046] After obtaining the refined spatial mesh, an adaptive time step needs to be determined based on the mesh size and diffusion coefficient to satisfy numerical stability conditions. For the numerical solution of the diffusion-reaction coupling equation, the time step selection must satisfy the Courant-Friedrichs-Lewy stability condition. For the refined mesh element, the time step... Must meet ,in This represents the maximum diffusion coefficient across all refined mesh cells. For the unrefined original mesh cells, the time step is... Must meet To ensure consistency of time steps across the entire solution domain, the minimum time step value of all grid cells is taken as the global adaptive time step. In practical calculations, to increase the numerical stability margin, the calculated time step is usually multiplied by a safety factor. The safety factor ranges from 0.5 to 0.8, and the final adaptive time step used is... .

[0047] Setting boundary conditions is crucial for solving the diffusion-response coupling equation. For a spatial sub-region located at the boundary of the treatment region, flux boundary conditions are used to describe the conduction behavior of the response intensity at the boundary. Specifically, for a boundary grid cell, its normal flux is determined by the spatial gradient of the response intensity at the boundary and the conduction coefficient. Assuming a boundary grid cell is located on the right boundary of the treatment region, and its normal is horizontal to the right, then the flux at this boundary... It can be represented as ,in The conduction coefficient at the boundary, The spatial gradient of the response intensity at the boundary along the horizontal direction is represented. This spatial gradient is numerically approximated by dividing the difference in response intensity between the boundary grid cell and its adjacent internal grid cells by the grid spacing. For other boundaries of the treatment region (left, top, and bottom), flux boundary conditions are set in a similar manner, simply by aligning the normal direction with the corresponding spatial gradient direction. For spatial sub-regions located within the treatment region, since they are surrounded by other spatial sub-regions and have no direct interaction with the external environment, periodic boundary conditions are set. Periodic boundary conditions mean that during the numerical solution process, the response intensity transmission between internal grid cells and their adjacent grid cells follows the internal evolution law of the diffusion-reaction coupling equation and is not directly affected by the external boundary.

[0048] After completing mesh refinement and boundary condition settings, the diffusion-reaction coupling equations were spatially and temporally discretized. Spatial discretization employed the finite difference method, discretizing the second-order spatial derivative in the diffusion term based on the refined spatial mesh. For each mesh element... The second spatial derivative of its response intensity It can be approximated as ,in The horizontal dimension of this grid cell. Represents grid cells The response intensity. The second spatial derivative in the vertical direction. Discretization is performed in a similar manner. Time discretization employs a hybrid implicit-explicit scheme, using an implicit Euler scheme for the diffusion term to improve numerical stability and an explicit scheme for the reaction term to reduce computational complexity. Let the current time step be... The next step is... Then the grid cell At time step response intensity Satisfying discrete equations , where superscript Indicates the current time step. This indicates the next time step. The discrete equations are established for all grid cells, forming a system of linear equations, which are solved using an iterative solver (such as the conjugate gradient method or the multigrid method).

[0049] The iterative solution process starts from the current time step and uses an adaptive time step. Future prediction time window The process proceeds incrementally. At each time step, the linear equations are solved to obtain the response intensity of all grid cells at that time step. To determine whether the iteration has converged, the change in predicted response intensity between adjacent time steps is calculated. For each grid cell... Its time step With time step The change in response intensity is Calculate the maximum value of the change in response intensity for all mesh elements. The maximum value is compared with the preset convergence threshold. Comparison. When When convergence is reached, the time progression process is terminated. Convergence threshold. The threshold value is set based on the required treatment accuracy, typically ranging from 0.1% to 0.5% of the average baseline tolerance threshold. If the iteration fails to converge, the solution continues for the next time step until the end of the future prediction time window is reached or the convergence condition is met. After convergence, the response intensity of each grid cell at the final time step is used as the predicted response intensity. For spatial sub-regions that have undergone grid densification, the predicted response intensities of the densified grid cells need to be spatially averaged or interpolated and mapped back to the response intensities of the original spatial sub-regions to obtain the predicted response intensities of each spatial sub-region. It is used for subsequent calculation of cross-effects and dosage adjustment decisions.

[0050] In one optional implementation, the cross-influence amount is superimposed onto the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and the remaining tolerance margin is calculated, including: A cumulative cross-influence matrix is ​​constructed for each spatial sub-region. The row index and column index of the cumulative cross-influence matrix represent the source spatial sub-region and the affected spatial sub-region, respectively. The matrix elements represent the cross-influence between corresponding pairs of spatial sub-regions. The total cross-influence of each spatial sub-region is obtained by summing all the cross-influence in the column vector of the matrix corresponding to each spatial sub-region when it is considered as an affected spatial sub-region. The predicted tolerance threshold is obtained by subtracting the corresponding total cross-influence from the baseline tolerance threshold of each spatial sub-region. The cumulative response load of each spatial sub-region at the current time is collected. The cumulative response load represents the time integral value of the response intensity of the spatial sub-region from the start of treatment to the current time. The difference between the predicted tolerance threshold and the cumulative response load is calculated as the remaining tolerance margin.

[0051] For example, after obtaining the cross-influence quantities between each spatial sub-region, these quantities need to be systematically integrated into the tolerance threshold prediction system to provide an accurate safety boundary for subsequent dose adjustments. To achieve this goal, a cross-influence quantity accumulation matrix is ​​constructed as the core data structure, which can completely record the interaction relationships between all pairs of spatial sub-regions within the treatment area.

[0052] Assuming the skin area to be treated is divided into For each spatial sub-region, the cumulative matrix of cross-influence is... for A 3D square matrix. Row indices in the matrix. Indicates the number of the affected source space sub-region, column index. Indicates the number of the affected spatial sub-region, matrix element Represents a spatial subregion For spatial sub-regions The resulting cross-influence. In the actual construction process, all spatial sub-regions are first globally numbered. The numbering rule can adopt row priority or column priority, for example, for those ordered by... The area divided by the grid can be located at the first Line 1 The spatial sub-regions of the column are numbered as follows After numbering, iterate through all spatial sub-region pairs and fill the corresponding positions in the matrix with the calculated cross-influence values. Note that the diagonal elements of the matrix... This represents the influence of a spatial sub-region on itself. In this method, this value is usually set to zero because the baseline tolerance threshold already includes the tolerance characteristics of the region itself.

[0053] matrix After construction is complete, for each spatial sub-region Extract the column vector corresponding to this sub-region when it is considered as the affected object. This column vector contains all other spatial subregion pairs. The cross-influence quantity. Summing all elements in the column vector yields the spatial sub-region. Total cross-influence The calculation formula is: The total cross-influence reflects the combined inhibitory effect of the response states of all other locations within the treatment area on the future tolerance of that spatial sub-region. In practical calculations, a sparse matrix storage method can be used to improve computational efficiency, because the cross-influence between spatial sub-regions that are far apart is usually close to zero and does not require storage or computation.

[0054] After obtaining the total cross-influence of each spatial sub-region, it is superimposed on the baseline tolerance threshold to obtain the predicted tolerance threshold. Specifically, for each spatial sub-region... Its predicted tolerance threshold By using baseline tolerance threshold Subtract the total cross-effect from the middle To obtain, that is The physical meaning of this subtraction operation is that excessive response in adjacent regions consumes part of the tolerance reserve of that region, leading to a reduction in the upper limit of its actual tolerable response intensity. The unit of the predicted tolerance threshold is consistent with the response intensity, usually using a normalized dimension or a quantitative value of a physiological response index. In some cases, if the total cross-influence is too large, causing the predicted tolerance threshold to be negative, it needs to be truncated to zero or a set minimum safety threshold to avoid anomalies in subsequent calculations.

[0055] To assess the remaining safety margin of each spatial sub-region at the current moment, the concept of cumulative response load needs to be introduced. Defined as a sub-region of this space from the start of treatment. Up to the current moment The integral of the response intensity with respect to time is expressed mathematically as follows: ,in For integration variables, Represents a spatial subregion At any moment The response intensity. In the discretized implementation, this integral is approximated by numerical integration methods, such as using the trapezoidal rule to divide the time interval into... If there are 10 sampling points, the cumulative response load can be expressed as: ,in For the first Each sampling time point. The physical meaning of cumulative response load is the total amount of tolerance capacity that has been consumed in this spatial sub-region; the larger the value, the heavier the treatment load that the region has already borne.

[0056] After obtaining the predicted tolerance threshold and cumulative response load, the remaining tolerance margin is calculated. As the difference between the two, that is The remaining tolerance margin characterizes how much additional response load this spatial sub-region can withstand under its current state without exceeding the safety threshold. A positive value indicates that the region still has a safety margin available and can continue to receive phototherapy; a value close to zero indicates that the region is nearing its tolerance limit and phototherapy needs to be reduced or stopped; a negative value indicates that the region has exceeded the safety range and phototherapy must be stopped immediately, and protective measures may be necessary.

[0057] In practical applications, to improve the robustness of predictions, a time-sliding window update can be applied to the cumulative cross-influence matrix. Specifically, a time window length is set. Only the cross-influence calculated within the specified time window is retained; historical influence values ​​exceeding the window are gradually reduced or discarded. This approach is based on physiological observations: skin tissue has a limited memory effect on short-term stimuli, and the impact of premature responses on current tolerance diminishes over time. The selection of the time window needs to be determined based on the specific type of phototherapy and skin characteristics; for ultraviolet phototherapy, a typical time window length is several minutes to tens of minutes.

[0058] Furthermore, when constructing the cumulative cross-influence matrix, a spatial distance attenuation factor can be introduced to weighted attenuate the cross-influence between spatially distant sub-regions. Let the spatial sub-regions be... and The Euclidean distance between them is The corrected cross-influence is then... ,in This is the distance attenuation coefficient, the value of which depends on the thermal conductivity and light scattering properties of the skin tissue. This correction more accurately reflects physiological phenomena because the interaction between distant regions is usually weaker.

[0059] When calculating the cumulative response load, a time-weighted mechanism can be introduced, making the contribution of recent response strengths to the cumulative load greater, while the contribution of earlier response strengths gradually decreases. Specifically, this is achieved by introducing an exponentially decaying weight function into the integral calculation. ,in The time decay factor is the corrected cumulative response load. This weighting method can more accurately reflect the recovery characteristics of skin tissue because the tissue gradually recovers some of its tolerance after a brief response.

[0060] The calculated remaining tolerance margin will serve as the core basis for subsequent dose adjustments. For spatial sub-regions with a large tolerance margin, the light intensity can be appropriately increased to improve treatment efficiency; for spatial sub-regions with a small tolerance margin, the light intensity needs to be reduced to ensure safety. Through this dynamic adjustment mechanism based on real-time monitoring and prediction, the treatment effect can be maximized while ensuring treatment safety, achieving personalized and precise phototherapy.

[0061] In one optional implementation, while maintaining a constant overall cumulative energy density in the treatment area, the differentiated light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region, including: A constrained optimization problem is constructed with the goal of maximizing the sum of the ratios of the adjusted light intensity to the remaining tolerance margin in each spatial sub-region. The remaining tolerance margin in each spatial sub-region is used as the upper bound of the inequality constraint, and the overall cumulative energy density of the treatment area is kept constant as the equality constraint. The Lagrange multiplier method is used to solve the basic adjustment components of each spatial sub-region. Traverse all matrix elements in the cumulative cross-influence matrix. When the value of a matrix element exceeds a preset propagation threshold, extract the index of the source spatial sub-region and the index of the affected spatial sub-region corresponding to the matrix element. Calculate the spatial distance between the source spatial sub-region and the affected spatial sub-region. Calculate the propagation compensation component of the spatial sub-region pair. The propagation compensation component is directly proportional to the conduction coefficient between the source spatial sub-region and the affected spatial sub-region and inversely proportional to the spatial distance. For each sub-region of the influence source space, sum up all propagation compensation components of the influence source space to obtain the total propagation compensation component of the sub-region of the influence source space. The total propagation compensation component of each spatial sub-region is superimposed on the corresponding basic adjustment component in the form of a negative value to obtain the differentiated light intensity adjustment amount of each spatial sub-region, which is then sent to the phototherapy execution device for execution.

[0062] For example, after obtaining the remaining tolerance margin of each spatial sub-region, it is necessary to differentiate the light intensity of each spatial sub-region while maintaining a constant overall cumulative energy density of the treatment area. This process is achieved by constructing a constrained optimization problem, the core objective of which is to maximize treatment efficiency while satisfying safety constraints.

[0063] Let the area of ​​skin to be treated be divided into The spatial sub-region, the first The current illumination intensity of each spatial sub-region is The remaining tolerance margin is Definition of the first The light intensity adjustment for each spatial sub-region is: The adjusted light intensity is The objective function is designed to maximize the sum of the ratios of adjusted illumination intensity to remaining tolerance margin in each spatial sub-region. This objective function is expressed as: The physical meaning of this objective function is that for spatial sub-regions with a large remaining tolerance margin, the light intensity can be appropriately increased to enhance the therapeutic effect; while for spatial sub-regions with a small remaining tolerance margin, the light intensity needs to be reduced to avoid overstimulation. The optimization problem needs to satisfy two types of constraints. The first type is the inequality constraint, which requires that the adjusted light intensity must not exceed the safe upper limit allowed by the remaining tolerance margin of each spatial sub-region, expressed as: ;in For the first The safety utilization factor for each spatial sub-region is typically set between 0.6 and 0.8, used to maintain a safety margin between the remaining tolerance margin and the actual adjustment amount. The second type is the equality constraint, which requires the overall cumulative energy density of the treatment area to remain constant, i.e., the sum of the products of the light intensity and area of ​​each spatial sub-region remains unchanged before and after adjustment. Let the first... The area of ​​each spatial subregion is Then the equality constraint is expressed as: Adjusted light intensity Substituting the equality constraints, we get: This constraint ensures that adjustments to the light intensity are redistributed only among the spatial sub-regions, without altering the total energy received by the treatment area.

[0064] The Lagrange multiplier method is used to solve the above constrained optimization problem. The Lagrange function is constructed as follows: ; in These are the Lagrange multipliers corresponding to the equality constraints. For the first The Lagrange multipliers corresponding to the inequality constraints of each spatial subregion. For the Lagrange function with respect to... Taking the partial derivative and setting it to zero, we get: In practical solutions, inequality constraints are initially ignored, and only equality constraints are considered to solve the basic adjustment components. From the above formula, we can obtain: Substitute this expression into the equality constraints. The Lagrange multipliers can be obtained. The value of . Furthermore, to maximize the optimization objective, adjustments need to be made based on the differences in residual tolerance margins across different spatial sub-regions. Define the normalized residual tolerance margin: The basic adjustment component is calculated as follows: ;in The total area of ​​the treatment region. This formula indicates that spatial sub-regions with a residual tolerance margin higher than the average level will receive a positive adjustment component, i.e., increased light intensity; while spatial sub-regions with a residual tolerance margin lower than the average level will receive a negative adjustment component, i.e., decreased light intensity.

[0065] After obtaining the basic adjustment components, it is necessary to further consider the cross-influence between spatial sub-regions. This involves traversing the cumulative cross-influence matrix. All matrix elements ,in To affect the source space sub-region index, Index the affected spatial sub-regions. When Exceeding the preset propagation threshold When, it indicates a spatial sub-region For spatial sub-regions There are significant cross-influences, requiring analysis of the source spatial sub-regions. The light intensity is further suppressed to prevent its excessive response from further affecting adjacent areas.

[0066] For satisfying Spatial subregion pairs Calculate the spatial sub-region of the influence source. With the affected spatial sub-region Spatial distance between This distance can be calculated using the Euclidean distance between the centers of the spatial sub-regions. The propagation compensation component of this spatial sub-region pair is defined as: ;in The propagation compensation intensity coefficient typically ranges from 0.1 to 0.3. For spatial sub-regions and The conduction coefficient between them; This is the distance regularization constant, used to avoid when... The value becomes unstable when approaching zero, and is typically taken as 0.1 times the characteristic size of the spatial sub-region. The negative sign of the propagation compensation component indicates suppression of the illumination intensity of the influencing source spatial sub-region. This formula reflects the physical law that the propagation compensation component is directly proportional to the cross-influence and conduction coefficient, and inversely proportional to the spatial distance.

[0067] For each sub-region of the influence source space Summing all propagation compensation components that are considered as the source of influence yields the total propagation compensation component for this spatial sub-region: ; The summation process iterates through all target space sub-regions affected by the source space sub-region, accumulating the total propagation compensation. Since each individual propagation compensation component is negative, the total propagation compensation component is also negative, representing the overall suppression of the illumination intensity of the source space sub-region.

[0068] The total propagation compensation component of each spatial sub-region is superimposed on the corresponding base adjustment component in the form of a negative value to obtain the differential illumination intensity adjustment amount of each spatial sub-region: ;because Since the values ​​are already negative, simply adding them together achieves the suppression and correction of the basic adjustment component. In practical applications, to ensure that the adjusted light intensity does not exceed the safe upper limit, the final adjustment amount needs to be truncated. ; This truncation operation ensures the adjusted light intensity It is always a non-negative value and does not exceed the safety limit allowed by the remaining tolerance margin.

[0069] After completing the calculation of the differentiated illumination intensity adjustment for all spatial sub-regions, the adjustment amount will be... The data is sent to the phototherapy execution device. Based on the received adjustment, the phototherapy execution device adjusts the light source output power or irradiation time for each spatial sub-region in real time to achieve differentiated light intensity distribution. During the adjustment process, the response intensity of each spatial sub-region is continuously monitored, and dynamic fine-tuning is performed based on real-time feedback to ensure the safety and effectiveness of the treatment process.

[0070] In one optional implementation, calculating the propagation compensation component of the spatial sub-region pair includes: The basic compensation coefficient is obtained by multiplying the reciprocal of the spatial distance by the transmission coefficient. The excess amount of the matrix element value exceeding the preset propagation threshold is calculated. The excess amount is mapped to a propagation risk factor through the Sigmoid function, such that the larger the excess amount, the closer the propagation risk factor is to the saturation upper limit. The propagation compensation component of the spatial sub-region pair is obtained by multiplying the basic compensation coefficient by the propagation risk factor. A propagation path graph is constructed based on the cumulative cross-influence matrix. The nodes of the propagation path graph represent spatial sub-regions, and the directed edges represent propagation relationships where the cross-influence exceeds the preset propagation threshold. A depth-first search algorithm is used to traverse all propagation paths originating from each influence source spatial sub-region and record the propagation path length. When the propagation path length exceeds the preset hop count threshold, a cascade compensation coefficient proportional to the propagation path length is calculated. The cascade compensation coefficient is then superimposed on the propagation compensation component of the first affected spatial sub-region on the propagation path corresponding to the influence source spatial sub-region.

[0071] For example, after determining that there are significant cross-influences between pairs of spatial sub-regions, it is necessary to further quantify the specific impact of this propagation effect on the illumination intensity adjustment strategy. The calculation of the propagation compensation component aims to transform the interactions between spatial sub-regions into actionable dose adjustment parameters, thereby achieving more precise inter-regional coordination in subsequent optimization processes.

[0072] For the cumulative matrix of cross-influence quantities that satisfies Conditional spatial subregion pairs First, the basic compensation coefficient is calculated. This coefficient comprehensively considers two dimensions: spatial distance and conduction characteristics. The reciprocal of the coefficient reflects the degree of proximity; the closer the distance, the stronger the interaction. The conduction coefficient... This characterizes the energy transfer efficiency of physiological tissue in that direction. Multiplying the two yields the basic compensation coefficient. The calculation formula is: The distance regularization constant is introduced into the denominator. This is to avoid numerical singularities when two spatial sub-regions are closely adjacent. The basic compensation coefficient establishes a direct mapping relationship between the spatial topology and the compensation strength.

[0073] The next step is to quantify the contribution of the degree of excess of cross-influence to the propagation risk. Define the excess amount. The difference between the matrix element values ​​and the preset propagation threshold, i.e. This excess reflects the degree to which the actual cross-effect deviates from the safety threshold. To transform this linear excess into a risk assessment indicator with saturation characteristics, a sigmoid function is used for nonlinear mapping. (Propagation Risk Factor) The calculation expression is: ;in The steepness parameter of the sigmoid function controls the sensitivity of the risk factor to the increase of the excess. When the excess is small, the risk factor exhibits an approximately linear growth; as the excess continues to increase, the risk factor gradually approaches the saturation upper limit of 1. This saturation characteristic ensures that even under extreme cross-influence conditions, the compensation component will not grow indefinitely, leading to system instability.

[0074] Multiply the basic compensation coefficient by the transmission risk factor and introduce the transmission compensation intensity coefficient. Perform global scale adjustment to obtain spatial sub-region pairs Propagation compensation component : This propagation compensation component characterizes the effect of spatial sub-regions. For spatial sub-regions Excessive influence requires spatial sub-regions The additional compensation correction applied during the adjustment of illumination intensity. This correction will be used as a constraint or penalty term in subsequent optimization calculations.

[0075] However, the above calculations only consider the first-order propagation effect between directly adjacent spatial sub-regions. In actual skin phototherapy, the energy response may form a cascading propagation chain along multiple spatial sub-regions. That is, a high-response region not only affects its directly adjacent regions but may also indirectly affect more distant regions through intermediate regions. To capture this multi-hop propagation phenomenon, it is necessary to construct a propagation path map and perform topological analysis.

[0076] The propagation path graph is a directed graph structure whose node set corresponds to all spatial sub-regions and directed edges. The conditions for existence are This means the cross-influence exceeds a preset propagation threshold. The direction of the edges indicates the direction of influence transmission, from the source of influence to the affected area. This graph structure transforms the cross-influence relationships, originally stored in matrix form, into a topological representation that facilitates pathfinding.

[0077] A depth-first search algorithm is employed, starting from each possible spatial sub-region of influence sources, and systematically traversing all reachable propagation paths. Specifically, for each spatial sub-region... Starting from this point, the algorithm recursively visits all satisfying... neighboring nodes and continue from node Proceed to explore the next level of propagation. Record the current path length during the traversal. The length is defined as the number of edges traversed on the path, i.e., the propagation hop count.

[0078] When the length of a certain propagation path exceeds a preset hop count threshold At this point, the path is considered to constitute a significant cascading propagation effect. Therefore, additional cascading compensation needs to be applied to the first affected spatial sub-region along this path. The cascading compensation coefficient is then defined. It is proportional to the propagation path length, and the calculation formula is: ;in This is the cascade compensation ratio coefficient. Indicates the sub-region of the influence source space The first affected spatial sub-region on a certain propagation path from which the influence originates. The design logic of this coefficient is that the longer the propagation path, the wider the range of influence of the source, and the stronger its potential cumulative effect. Therefore, stronger inhibitory compensation needs to be applied at the source.

[0079] The cascaded compensation coefficients are superimposed on the propagation compensation components of the corresponding spatial sub-region pairs to obtain the corrected propagation compensation components. It should be noted that for the same sub-region of the source space... There may be multiple different propagation paths, each corresponding to a different first affected spatial sub-region. In actual calculations, it is necessary to calculate the cascade compensation coefficient for all paths that meet the conditions, and then accumulate it into the corresponding propagation compensation component.

[0080] After completing the propagation compensation component calculation for all spatial sub-region pairs, for each influence source spatial sub-region It is necessary to summarize the total propagation compensation effect it produces on all affected spatial sub-regions. Total propagation compensation component By satisfying all spatial subregion The summation of the corrected propagation compensation components yields: This total propagation compensation component will be used as a spatial sub-region. Additional constraints in differentiated light intensity adjustment ensure that when increasing the light intensity of a region, its potential negative impact on the surrounding area is fully considered, thereby avoiding a chain reaction caused by local overcompensation.

[0081] In practical applications, constructing and traversing the propagation path graph requires efficient data structures. An adjacency list can be used to store the directed graph, where each node maintains a list of outgoing edges, recording all target nodes that meet the threshold condition and their corresponding cross-influence values. During the depth-first search process, an array of visit markers needs to be maintained to prevent infinite loops in the presence of cycles. When a node is detected as having been visited in the current path, the exploration of that branch is terminated, and the process backtracks to the previous level.

[0082] For treatment areas with complex topologies, propagation paths may exhibit a tree-like or network-like distribution. A tree-like structure indicates that the influence spreads outward from a single source, while a network structure indicates that multiple sources of influence interact with each other. In a network structure, a spatial sub-region may simultaneously serve as an intermediate node for multiple propagation paths. In this case, it is necessary to accumulate all the cascaded compensation coefficients received by that node to fully reflect its pivotal role in the propagation network.

[0083] The final calculation result of the propagation compensation component will be compared with the base adjustment component calculated based on the residual tolerance margin. This process integrates various factors to create a differentiated light intensity adjustment that comprehensively considers both local tolerance and spatial propagation effects. The integration method can employ weighted summation or constrained optimization, depending on the system's trade-off strategy between safety and treatment efficiency. Through this multi-layered compensation mechanism, more refined and safer dose allocation can be achieved while maintaining a constant overall energy density in the treatment area.

[0084] A second aspect of the present invention provides a skin phototherapy dose intelligent adjustment and monitoring system, comprising: An initial dose unit is used to generate an initial dose plan based on the initial physiological state information of the skin area to be treated. The partition monitoring unit is used to divide the skin area to be treated into multiple spatial sub-regions, record the baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions, and monitor the response intensity of each spatial sub-region in parallel during the execution of the initial treatment dose plan. The cross-influence unit is used to calculate the cross-influence between adjacent spatial sub-regions based on the response intensity gradient of each spatial sub-region using the diffusion-reaction coupling equation. The cross-influence value characterizes the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. The margin calculation unit is used to construct a cross-influence accumulation matrix to record the cross-influence between each pair of spatial sub-regions, and to superimpose the cross-influence to the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. The adjustment execution unit is used to calculate the differentiated light intensity adjustment amount of each spatial sub-region based on the remaining tolerance margin of each spatial sub-region while maintaining the overall cumulative energy density of the treatment area, and then send it to the phototherapy execution device for execution.

[0085] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0086] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0087] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0088] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for intelligent adjustment and monitoring of skin phototherapy dosage, characterized in that, include: An initial dosage plan is generated based on the initial physiological state information of the skin area to be treated. The skin area to be treated is divided into multiple spatial sub-regions. The baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dose scheme. Based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using the diffusion-response coupling equation. The cross-influence represents the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. A cross-influence accumulation matrix is ​​constructed to record the cross-influence between each pair of spatial sub-regions. The cross-influence is then superimposed on the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. While maintaining a constant overall cumulative energy density in the treatment area, the differential light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region and then sent to the phototherapy execution device for execution.

2. The method according to claim 1, characterized in that, The skin area to be treated is divided into multiple spatial sub-regions. The baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions are recorded. The response intensity of each spatial sub-region is monitored in parallel during the execution of the initial treatment dose regimen, including: Based on the rate of change of physiological parameters in the initial physiological state information, the skin region to be treated is adaptively divided into grids, and the division density is increased in regions where the rate of change of physiological parameters is higher than a preset change threshold to obtain multiple spatial sub-regions. Test excitation is applied to each spatial sub-region and test response intensity is collected. The baseline tolerance threshold of each spatial sub-region is determined based on the test response intensity. Response time-series correlation data of adjacent spatial sub-regions are collected, and the conduction coefficient is calculated based on the response propagation delay and response amplitude attenuation rate. Instantaneous physiological response parameters are collected and integrated over time to obtain the response intensity.

3. The method according to claim 1, characterized in that, Based on the response intensity gradient of each spatial sub-region, the cross-influence between adjacent spatial sub-regions is calculated using the diffusion-response coupling equation, including: The central difference operator is used to calculate the spatial gradient of the response intensity in each spatial sub-region in the spatial dimension, and a composite gradient vector is constructed by combining the time rate of change of the response intensity. A diffusion-reaction coupling equation is constructed, which includes a diffusion term and a reaction term. The diffusion term adopts a second-order partial differential form, and the diffusion coefficient is set as the product of the conduction coefficient and the spatial components of the composite gradient vector. The reaction term adopts a first-order kinetic form, and the reaction rate is proportional to the current response intensity. The length of the future prediction time window is adaptively determined based on the ratio of the current response intensity of each spatial sub-region to the baseline tolerance threshold. The larger the ratio, the shorter the length of the future prediction time window. The diffusion-response coupling equation is numerically solved within the future prediction time window using an implicit time discretization method to obtain the predicted response intensity of each spatial sub-region at the end of the future prediction time window. The difference between the predicted response intensity and the current response intensity is calculated as the response intensity increment. When the response intensity increment exceeds the response threshold, it is determined that an over-response has occurred, and the portion of the response intensity increment that exceeds the response threshold is taken as the cross-influence quantity.

4. The method according to claim 3, characterized in that, The implicit time discretization method is used to numerically solve the diffusion-reaction coupling equation within a future prediction time window, including: In the spatial sub-region where the magnitude of the composite gradient vector exceeds a set threshold and its adjacent spatial sub-regions, a local mesh refinement is performed to obtain a refined spatial mesh. An adaptive time step is determined based on the mesh size and diffusion coefficient of the refined spatial mesh. For the spatial sub-regions located at the boundary of the treatment area, flux boundary conditions are set based on the spatial gradient of the response intensity and the conduction coefficient at the boundary. For the spatial sub-regions located inside the treatment area, periodic boundary conditions are set. The diffusion-reaction coupling equation is spatially discretized based on the encrypted spatial grid. The diffusion term is discretized in time using an implicit Euler scheme, and the reaction term is discretized in time using an explicit scheme. The solution is iteratively solved within the future prediction time window with the adaptive time step until the change in the predicted response intensity of adjacent time steps is less than the convergence threshold, and the predicted response intensity is obtained.

5. The method according to claim 1, characterized in that, The predicted tolerance threshold is obtained by superimposing the cross-influence amount onto the baseline tolerance threshold of each spatial sub-region, and the remaining tolerance margin is calculated, including: The row and column indices of the cumulative cross-influence matrix represent the source and affected spatial sub-regions, respectively. The matrix elements represent the cross-influence between corresponding pairs of spatial sub-regions. The total cross-influence of a spatial sub-region is obtained by summing all the cross-influences in the column vector of the matrix corresponding to each spatial sub-region when it is considered an affected spatial sub-region. The predicted tolerance threshold is obtained by subtracting the corresponding total cross-influence from the baseline tolerance threshold of each spatial sub-region. The cumulative response load of each spatial sub-region at the current time is collected. The cumulative response load represents the time integral value of the response intensity of the spatial sub-region from the start of treatment to the current time. The difference between the predicted tolerance threshold and the cumulative response load is calculated as the remaining tolerance margin.

6. The method according to claim 1, characterized in that, Under the constraint of maintaining a constant overall cumulative energy density in the treatment area, the differentiated light intensity adjustment amount for each spatial sub-region is obtained by solving based on the remaining tolerance margin of each spatial sub-region, including: A constrained optimization problem is constructed with the goal of maximizing the sum of the ratios of the adjusted light intensity to the remaining tolerance margin in each spatial sub-region. The remaining tolerance margin in each spatial sub-region is used as the upper bound of the inequality constraint, and the overall cumulative energy density of the treatment area is kept constant as the equality constraint. The Lagrange multiplier method is used to solve the basic adjustment components of each spatial sub-region. Traverse all matrix elements in the cumulative cross-influence matrix. When the value of a matrix element exceeds a preset propagation threshold, extract the index of the source spatial sub-region and the index of the affected spatial sub-region corresponding to the matrix element. Calculate the spatial distance between the source spatial sub-region and the affected spatial sub-region. Calculate the propagation compensation component of the spatial sub-region pair. The propagation compensation component is directly proportional to the conduction coefficient between the source spatial sub-region and the affected spatial sub-region and inversely proportional to the spatial distance. For each sub-region of the influence source space, sum up all propagation compensation components of the influence source space to obtain the total propagation compensation component of the sub-region of the influence source space. The total propagation compensation component of each spatial sub-region is superimposed on the corresponding basic adjustment component in the form of a negative value to obtain the differentiated illumination intensity adjustment amount of each spatial sub-region.

7. The method according to claim 6, characterized in that, Calculating the propagation compensation components of this spatial sub-region pair includes: The basic compensation coefficient is obtained by multiplying the reciprocal of the spatial distance by the transmission coefficient. The excess amount of the matrix element value exceeding the preset propagation threshold is calculated. The excess amount is mapped to the propagation risk factor. The propagation compensation component of the spatial sub-region pair is obtained by multiplying the basic compensation coefficient by the propagation risk factor. A propagation path graph is constructed based on the cumulative cross-influence matrix. The nodes of the propagation path graph represent spatial sub-regions, and the directed edges represent propagation relationships where the cross-influence exceeds the preset propagation threshold. A depth-first search algorithm is used to traverse all propagation paths originating from each influence source spatial sub-region and record the propagation path length. When the propagation path length exceeds the preset hop count threshold, a cascade compensation coefficient proportional to the propagation path length is calculated. The cascade compensation coefficient is then superimposed on the propagation compensation component of the first affected spatial sub-region on the propagation path corresponding to the influence source spatial sub-region.

8. A skin phototherapy dosage intelligent adjustment and monitoring system, used to implement the method as described in any one of claims 1-7, characterized in that, include: An initial dose unit is used to generate an initial dose plan based on the initial physiological state information of the skin area to be treated. The partition monitoring unit is used to divide the skin area to be treated into multiple spatial sub-regions, record the baseline tolerance threshold of each spatial sub-region and the conduction coefficient between adjacent sub-regions, and monitor the response intensity of each spatial sub-region in parallel during the execution of the initial treatment dose plan. The cross-influence unit is used to calculate the cross-influence between adjacent spatial sub-regions based on the response intensity gradient of each spatial sub-region using the diffusion-reaction coupling equation. The cross-influence value characterizes the degree to which the excessive response of a certain spatial sub-region inhibits the tolerance threshold of adjacent spatial sub-regions within a future time window. The margin calculation unit is used to construct a cross-influence accumulation matrix to record the cross-influence between each pair of spatial sub-regions, and to superimpose the cross-influence to the baseline tolerance threshold of each spatial sub-region to obtain the predicted tolerance threshold and calculate the remaining tolerance margin. The adjustment execution unit is used to calculate the differentiated light intensity adjustment amount of each spatial sub-region based on the remaining tolerance margin of each spatial sub-region while maintaining the overall cumulative energy density of the treatment area, and then send it to the phototherapy execution device for execution.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images