A method for predicting the risk of mold growth in wood storage

Abstract

The present application provides a kind of wood storage mouldy risk prediction method, it is related to wireless communication technical field, it is by receiving the temperature and humidity time series data of each monitoring point of wood stack, grid mapping is carried out according to spatial coordinates, the environmental parameter distribution field covering storage space is constructed;Based on the humidity gradient and temperature gradient of adjacent grid node, the humidity conduction rate field and spore diffusion potential field across grid are calculated, and the mouldy risk propagation tensor is generated by space convolution;The time difference characteristic extraction of sliding window is carried out to mouldy risk propagation tensor along time axis, the spatial cluster area showing cumulative enhancement mode is identified, and is marked as high-risk propagation source point set;The future evolution trend of high-risk propagation source point set and environmental parameter distribution field is coupled and deduced in space-time, and the wood stack mouldy risk grade distribution map of partition layer is generated.The present application improves the spatial resolution of mouldy risk monitoring and the timeliness of early warning, which helps to reduce the incidence and loss degree of wood storage mould.

Description

Technical Field

[0001] This invention relates to the field of risk prediction technology, and more specifically, to a method for predicting the risk of mold growth in timber storage. Background Technology

[0002] With the continuous development of timber processing and logistics warehousing technologies, timber, as an important renewable resource, has seen its storage and preservation technologies gradually become a key focus of industry research. Especially under modern large-scale warehousing management models, environmental monitoring and risk early warning systems for timber stacking spaces are constantly being improved. Traditional timber warehousing management mainly relies on regular manual inspections, adjusting storage temperature and humidity parameters using empirical formulas to reduce biological losses such as mold. However, with the introduction of automation and information technology, real-time temperature and humidity monitoring and local intelligent control methods based on the Internet of Things are gradually becoming widespread, improving the scientific and refined level of timber quality management to a certain extent. Currently, domestic and international research on the risk of mold growth in timber storage mainly focuses on environmental factor threshold determination, early warning models for localized temperature and humidity anomalies, or statistical analysis methods, and some results have already been applied to intelligent warehousing supervision systems.

[0003] However, existing technologies for assessing the risk of mold growth in timber storage suffer from the following shortcomings: First, traditional methods, based on local or single-point temperature and humidity data, struggle to reflect the dynamic evolution of environmental parameters within the entire three-dimensional space of the stack, resulting in limited spatial resolution and difficulty in capturing early microenvironmental changes leading to mold growth. Second, current risk modeling methods often employ rule-based threshold discrimination or empirical risk zone definitions, failing to effectively express the complex nonlinear effects of humidity gradients and temperature changes on spore diffusion and risk zone formation under multi-point coupling conditions. Furthermore, some existing solutions focus on planar trend analysis of historical monitoring data, neglecting the spatiotemporal propagation and risk accumulation processes between high-risk areas within the storage space. They lack in-depth modeling and accurate identification capabilities for risk propagation dynamics, particularly limiting their effectiveness in early warning of mold risk and locating local high-risk sources. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention provides a method for predicting the risk of mold growth in timber storage. This method can, to some extent, solve the problem that the current timber storage environment cannot accurately grasp the dynamic distribution of temperature and humidity in different areas of the stack, leading to untimely identification of the spatial propagation of mold risk and high-risk sources, resulting in insensitive local risk warnings and difficulty in taking timely differentiated prevention and control measures.

[0005] According to one aspect of the present invention, a method for predicting the risk of mold growth in timber storage is provided, comprising:

[0006] Receive temperature and humidity time-series data from each monitoring point of the timber stack, and perform gridded mapping of the temperature and humidity time-series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space.

[0007] Based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, the humidity conduction rate field and spore diffusion potential field across the grid are calculated, and the humidity conduction rate field and the spore diffusion potential field are spatially convolved to generate a mold risk propagation tensor.

[0008] The gradient rate of change of the mold risk propagation tensor is extracted within a sliding window along the time axis. Spatial cluster regions that exhibit a cumulative enhancement pattern in the risk propagation tensor are identified, and the spatial cluster regions are marked as a set of high-risk propagation sources.

[0009] By spatiotemporally coupling the set of high-risk transmission sources with the future time-period evolution trend of the environmental parameter distribution field, a partitioned and layered distribution map of the risk level of mold growth in timber stacks is generated.

[0010] Furthermore, when constructing the environmental parameter distribution field, the timber storage space is discretized into a three-dimensional grid structure based on the three-dimensional geometric boundary data of the timber storage space. The side length of the grid cell in any dimension is taken as the smaller of one-third to one-half of the average spacing between sensors in the air channel between adjacent timber stacks and one-third to one-half of the minimum spacing between sensors inside the stack.

[0011] Based on the shortest spatial distance from each grid node to the geometric boundary of the timber stack, a medium attribute is assigned: those located outside the boundary are assigned air medium attributes, those located inside the boundary and whose shortest distance does not exceed the preset surface layer depth are assigned stack surface medium attributes, and those exceeding the preset surface layer depth are assigned stack interior medium attributes.

[0012] Furthermore, when calculating the humidity conduction rate field, the cross-grid humidity conduction coefficient is determined for each pair of adjacent grid nodes based on their medium property identifier:

[0013] If both nodes are air-based media, the air humidity diffusion coefficient is used.

[0014] If both nodes have surface or internal stacking media properties, then the base value of the wood moisture diffusion coefficient and the density correction factor are used. The product of, where The normalized stacking density of the two nodes is the arithmetic mean. Density attenuation coefficient;

[0015] If the two nodes represent the surface medium properties and the interior medium properties of the stack, respectively, the harmonic mean of the wood moisture diffusion coefficients of the two nodes is multiplied by an interfacial conduction reduction factor ranging from 0.3 to 0.6.

[0016] If the two nodes have air medium properties and stack internal medium properties respectively, the cross-grid humidity conduction coefficient is forced to zero.

[0017] Furthermore, when calculating the spore diffusion potential field, for grid nodes that meet the conditions of spore activation temperature range and relative humidity exceeding the spore activation humidity threshold, the local spore diffusion potential benchmark value is calculated by multiplying the temperature response function and the humidity response function. ;

[0018] The temperature response function adopts a Gaussian form. ,in, This is the temperature value. For the optimal temperature, The standard deviation parameter;

[0019] The humidity response function uses a S-shaped logic function. ,in, This is the relative humidity value. It is the half-saturation constant. The slope parameter;

[0020] For mesh nodes that are assigned surface media properties or internal media properties of the stack, Multiply by the stacking density suppression factor , This is the spore diffusion density sensitivity coefficient. This represents the normalized node density.

[0021] For grid nodes that do not meet the conditions for spore activation, Set it directly to zero.

[0022] Furthermore, when generating the aforementioned mold risk propagation tensor, weight values ​​are used. The normalized 3D Gaussian convolution kernel is used to calculate the mold risk propagation intensity value for each target grid node using the following formula. :

[0023]

[0024] in, This represents the set of all grid nodes within the coverage area of ​​the spatial convolution kernel centered on the target grid node. For set The grid node number in the data. This serves as the medium attribute identifier for the target mesh node. For the first The medium attribute identifier of each grid node; For the first Humidity-spore coupling risk term for each grid node For the first Normalized Gaussian space weights for each grid node; For the first The cross-medium propagation attenuation factor of each grid node to the target grid node;

[0025] The formula for calculating the humidity-spore coupling risk term is as follows:

[0026]

[0027] in, For the first The humidity conduction rate value of each grid node. For the first The comprehensive spore dispersal potential value of each grid node. Reference humidity conduction rate; For reference spore dispersal potential; This is a nonlinear enhancement index for coupling risk, with a value ranging from 1 to 3. When the value is greater than 1, the contribution weight of nodes where nonlinear amplification of high humidity conduction and high spore diffusion coexist is calculated. It is a log-stable term.

[0028] Furthermore, during the spatiotemporal coupling simulation, for each grid node in the three-dimensional spatial grid... In every future time step Calculate the environmental suitability function using the following formula :

[0029]

[0030] in, Let the temperature suitability function be... Let the humidity suitability function be... It is the airflow suppression function;

[0031] The temperature suitability function is:

[0032]

[0033] in, This is a dimensionless temperature sensitivity coefficient. For dimensionless temperature, it is expressed as , These are the lower and upper boundaries of the effective temperature window for mold growth. This is a predicted temperature value;

[0034] The temperature suitability function is:

[0035]

[0036] in, The dimensionless humidity sensitivity coefficient Dimensionless humidity, expressed as , To propagate the activation humidity threshold, This is a humidity forecast value;

[0037] The airflow suppression function is:

[0038]

[0039] in, This is a dimensionless airflow suppression coefficient. The dimensionless airflow intensity is expressed as ,in For reference wind speed, This is the predicted airflow velocity.

[0040] Furthermore, the generation of the timber stack mold risk level distribution map includes:

[0041] After completing the projection of all future time steps, for each grid node Calculate the comprehensive risk level assessment value using the following formula:

[0042]

[0043] in, This serves as the baseline value for the risk intensity of the last time slice of the aforementioned mold risk propagation tensor. To normalize the cumulative transmission risk intensity value, The overall risk level of transmission;

[0044] With all grid nodes The upper quartile and median are used as the thresholds for high-risk and medium-risk classification, respectively. Seed-based extended clustering is used to form risk partitions for adjacent connected grid nodes of the same level. The vertical coordinate range of the three-dimensional spatial grid is divided into three levels: bottom, middle and top. The results of planar risk partitioning and vertical stratification are combined to generate a partitioned and stratified distribution map of the risk level of mold growth in timber stacks.

[0045] According to another aspect of the present invention, a timber storage mold risk prediction system is provided, comprising:

[0046] The field construction module is used to receive temperature and humidity time series data from each monitoring point of the timber stack, and to perform grid-based mapping of the temperature and humidity time series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space.

[0047] The risk propagation module is used to calculate the humidity conduction rate field and spore diffusion potential field across grids based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, and to perform spatial convolution operation on the humidity conduction rate field and the spore diffusion potential field to generate a mold risk propagation tensor.

[0048] The source identification module is used to extract the gradient change rate of the mold risk propagation tensor within a sliding window along the time axis, identify spatial cluster regions in the risk propagation tensor that exhibit a cumulative enhancement pattern, and mark the spatial cluster regions as a set of high-risk propagation sources.

[0049] The zoning simulation module is used to perform spatiotemporal coupling simulation of the future evolution trend of the set of high-risk transmission sources and the distribution field of environmental parameters, and generate a zoning and stratified distribution map of the risk level of mold growth in timber stacks.

[0050] Compared with existing technologies, this invention, through spatial gridding modeling and risk propagation dynamics analysis of temperature and humidity time-series data from various monitoring points in timber stacks, can accurately capture the dynamic distribution of humidity and temperature within the storage space and their impact on mold risk. Simultaneously, by employing risk propagation tensor and sliding window analysis, it effectively identifies spatial clusters of areas where risk accumulation is enhanced, enabling early location of high-risk propagation sources and stratified risk level determination. This not only improves the spatial resolution and timeliness of mold risk monitoring but also helps managers develop more targeted prevention and control measures for high-risk areas, thereby reducing the incidence and extent of mold damage during timber storage. Attached Figure Description

[0051] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings:

[0052] Figure 1 A flowchart of a timber storage mold risk prediction method according to an embodiment of the present invention;

[0053] Figure 2 This is a distribution map of the risk level of mold growth in different zones and layers of timber stacks according to an embodiment of the present invention. Detailed Implementation

[0054] Hereinafter, exemplary embodiments according to the present invention will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of the present invention, and not all embodiments of the present invention. It should be understood that the present invention is not limited to the exemplary embodiments described herein.

[0055] As mentioned in the background section above, existing technologies have two main problems: First, existing monitoring methods for timber storage environments are mostly based on single-point or limited monitoring data, which cannot accurately reflect the dynamic distribution and changes of temperature and humidity in different areas of the stacking space, resulting in a lack of refined spatiotemporal perception of the formation and spread of mold risk; Second, traditional risk warning models mostly use empirical thresholds or static area divisions, which make it difficult to identify spatial clusters and high-risk sources where risks gradually accumulate and increase in a timely manner, and ignore the nonlinear coupling effect of environmental factors such as humidity and temperature, as well as the dynamic spread characteristics of mold risk.

[0056] Figure 1 This is a system block diagram of a timber storage mold risk prediction method according to an embodiment of the present invention. Figure 1 and Figure 2 As shown, the methods for predicting the risk of mold growth in timber storage include:

[0057] S1: Receive temperature and humidity time-series data from each monitoring point in the timber stack, and perform gridded mapping of the temperature and humidity time-series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space.

[0058] The system receives time-series temperature and humidity data from multiple temperature and humidity sensors deployed within a timber storage space. Simultaneously, it acquires three-dimensional geometric boundary data of the timber storage space and density distribution data for each timber stack. The temperature and humidity sensors are deployed as follows: at least one surface sensor is placed on the surface of each timber stack; at least two internal sensors are placed at different depths along the thickness direction within each timber stack; and at least one channel sensor is placed in the air channel between adjacent timber stacks. The sensor placement height in the air channel corresponds to the placement height of the sensors within the stacks. Each temperature and humidity sensor has a unique three-dimensional spatial coordinate identifier, with a predetermined corner point of the timber storage space as the origin, recording the sensor's position within the storage space. The coordinate values ​​in the degree, width, and height directions are used. Time alignment processing is performed on the temperature and humidity time-series data collected by each temperature and humidity sensor. Specifically, a preset fixed time interval is used as a unified time reference. The timestamps of each sensor are mapped to the corresponding sampling times under the unified time reference. If a sensor has missing data at a certain sampling time, a linear time interpolation algorithm or a cubic spline interpolation algorithm is used to interpolate and complete the data based on the temperature and humidity values ​​of the two nearest valid sampling times before and after the missing time. If there are fewer than two valid sampling points before and after the missing time, a nearest neighbor interpolation algorithm is used to complete the data based on the temperature and humidity values ​​of the nearest valid sampling time, thus forming a time series of temperature and humidity values ​​for each monitoring point within a continuous time period.

[0059] Based on the three-dimensional geometric boundary data of the timber storage space and the geometric position data of each timber stack, the timber storage space is discretized into a three-dimensional grid structure composed of multiple regular cuboid grid units. Each grid node has a defined grid center coordinate and grid volume. The size of the grid unit is determined by the minimum spacing of the sensors inside the stack and the average spacing of the sensors in the air channels. The side length of the grid unit in any dimension is the smaller of one-third to one-half of the average spacing of the sensors in the air channels and one-third to one-half of the minimum spacing of the sensors inside the stack. Based on the geometric position and density distribution data of the timber stacks, each grid node is assigned a medium attribute identifier. The medium attribute identifier includes at least air medium, stack surface medium, and stack internal medium. For grid nodes assigned stack surface medium attributes or stack internal medium attributes, the identifier is determined based on the location of the grid node. The stack density distribution data corresponding to the location is used to calculate the stack density value of the grid node through spatial location matching or three-dimensional linear interpolation. Specifically, when a grid node is located outside the geometric boundary of all timber stacks, it is assigned the air medium attribute; when a grid node is located inside the geometric boundary of a timber stack, and the shortest spatial distance from the grid node to the geometric boundary of the timber stack does not exceed a preset surface layer depth, it is assigned the stack surface medium attribute; when a grid node is located inside the geometric boundary of a timber stack, and the shortest spatial distance from the grid node to the geometric boundary of the timber stack exceeds a preset surface layer depth, it is assigned the stack interior medium attribute. The preset surface layer depth is determined based on the deployment depth of the stack surface sensors, and its value ranges from one-third to one-half of the vertical depth distance between the stack surface sensor and the nearest internal sensor.

[0060] For each time slice, the temperature and humidity values ​​of each temperature and humidity sensor at the corresponding time of the time slice are taken as known observation values, and the spatial interpolation algorithm is used to calculate the temperature and humidity values ​​of each grid node at that time.

[0061] The spatial interpolation algorithm specifically comprises: for each grid node to be interpolated, a set of temperature and humidity sensors participating in the interpolation is determined within a preset spatial search radius, centered on the grid node. The spatial search radius is 1.5 to 3 times the average spacing between sensors in the air channel. For grid nodes assigned air medium properties, an inverse distance-weighted interpolation function based on distance weights is used to calculate the temperature and humidity values ​​of the grid node. The distance weights are the negative square power of the spatial Euclidean distance from each participating sensor to the grid node. For grid nodes assigned stack surface medium properties or stack internal medium properties, a modified interpolation function considering both stack density weighting and distance weighting is used to calculate the temperature and humidity values ​​of the grid node. The distance weights are the negative square power of the spatial Euclidean distance from each participating sensor to the grid node, and the stack density weighting factor is expressed in exponential form. ,in This is a dimensionless stacking density sensitivity coefficient, ranging from 0.5 to 2.0. To normalize the density difference, the calculation method is to divide the absolute value of the difference between the stack density value of the grid node and the stack density value of the sensor location involved in the interpolation by the stack density value of the grid node. At the boundary between air-medium grid nodes and surface-medium grid nodes or internal-medium grid nodes within the stack, an interpolation discontinuity constraint is set. This constraint is as follows: for air-medium grid nodes whose shortest spatial distance from the stack geometric boundary does not exceed one grid cell side length, sensors located on the stack surface or inside the stack are excluded from the sensor set involved in the interpolation; for air-medium grid nodes whose shortest spatial distance from the stack geometric boundary does not exceed one grid cell side length, sensors located on the stack surface or inside the stack are excluded from the sensor set involved in the interpolation; for air-medium grid nodes whose shortest spatial distance from the stack geometric boundary does not exceed one grid cell side length, sensors located on the stack surface or inside the stack are excluded from the sensor set involved in the interpolation. For grid nodes of medium on the surface of a stack or medium inside a stack that are no more than one grid cell side length away, sensors located in air channels are excluded from the sensor set participating in the interpolation. If the number of same-medium sensors participating in the interpolation after excluding cross-medium sensors is less than a preset minimum number, the spatial search radius is gradually expanded with the grid node as the center until the number of same-medium sensors participating in the interpolation reaches the preset minimum number, but the expanded search range must not include cross-medium sensors. If the preset minimum number cannot be met after expanding the search radius, the observation value of the same-medium sensor closest to the grid node is used as the interpolation result of the grid node.

[0062] The temperature and humidity values ​​of each grid node calculated for each time slice are stored in chronological order, and a spatial index and a time dimension index for the grid nodes are established. By combining the spatial coordinates, medium attribute identifiers and stacking density values ​​of each grid node, an environmental parameter distribution field covering the timber storage space is constructed.

[0063] The environmental parameter distribution field includes a temperature distribution field composed of the temperature values ​​of each grid node at each time moment and a humidity distribution field composed of the humidity values ​​of each grid node at each time moment, as well as a stacking density distribution field and a medium attribute field stored as spatial static attributes. The temperature distribution field and the humidity distribution field form a four-dimensional data structure containing three spatial dimension indices and one time dimension index. The stacking density distribution field and the medium attribute field are stored in a three-dimensional spatial index manner and are associated with the four-dimensional data structure through the grid node spatial index. The temperature value and humidity value of any grid node under any time index can be accessed through the grid node spatial index and the time dimension index. The stacking density value and the medium attribute identifier can be accessed through the grid node spatial index.

[0064] S2: Based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, calculate the humidity conduction rate field and spore diffusion potential field across the grid, and perform spatial convolution operation on the humidity conduction rate field and the spore diffusion potential field to generate a mold risk propagation tensor.

[0065] Based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, the humidity conduction rate field and spore diffusion potential field across the grid are calculated. The humidity conduction rate field and spore diffusion potential field are then spatially convolved to generate the mold risk propagation tensor.

[0066] The calculation of the humidity conduction rate field begins with the humidity values ​​of each grid node in the environmental parameter distribution field. If the environmental parameter distribution field stores relative humidity as a percentage, the relative humidity is first converted to absolute humidity. The conversion method is as follows: based on the temperature values ​​of the grid nodes... Calculate saturated water vapor density ,in The gas constant for air is 287 joules per kilogram per Kelvin. The temperature term in the exp function is in degrees Celsius; the temperature term in the denominator needs to be converted to Kelvin. The unit is kilograms per cubic meter; the absolute humidity value of this grid node is equal to the saturated water vapor density multiplied by the decimal form of the relative humidity. After completing the humidity value conversion, for any grid node in the three-dimensional grid structure, calculate the humidity gradient between this grid node and its six directly adjacent grid nodes in the three spatial dimensions. The humidity gradient is defined as the difference in absolute humidity values ​​between two adjacent grid nodes divided by the spatial distance between the center coordinates of the two grid nodes, and the dimension of the humidity gradient is kilograms per fourth square meter.

[0067] After determining the humidity gradient, for each pair of adjacent grid nodes, the cross-grid humidity transfer coefficient is determined based on the medium property identifiers of the two grid nodes, including:

[0068] If two adjacent grid nodes both have air as the medium property, the cross-grid humidity conduction coefficient is taken as the air humidity diffusion coefficient, with a value range of [value missing]. to Square meters per second. This value range is derived from the measured data of the diffusion coefficient of water vapor in the air under standard atmospheric pressure and normal temperature conditions. The variation range within the typical temperature range of the storage environment does not exceed ±10%, and can be regarded as a constant.

[0069] If two adjacent grid nodes both have surface medium properties or internal medium properties, the cross-grid humidity conduction coefficient is taken as the product of the wood humidity diffusion coefficient and the stack density correction factor, where the baseline value range for the wood humidity diffusion coefficient is [missing value]. to Square meters per second, this benchmark value is based on experimental measurements of the longitudinal and transverse diffusion coefficients of coniferous and broadleaf woods under moisture content conditions of 10% to 30%, which is about two orders of magnitude lower than the air humidity diffusion coefficient, reflecting the hindering effect of the porous medium of wood on moisture conduction; the stacking density correction factor adopts an exponential decay form. , The value is a dimensionless density attenuation coefficient, ranging from 0.3 to 1.0. This parameter was obtained through fitting experiments on the moisture permeation of multi-layer stacks. A larger value indicates a stronger inhibitory effect of increased density on diffusion. The normalized average density is obtained by dividing the arithmetic mean of the stacking density values ​​of two adjacent grid nodes by the reference density value, which is 500 kg per cubic meter, to ensure that the exponent is a dimensionless value.

[0070] If two adjacent grid nodes are respectively air medium property and stack surface medium property, the cross-grid humidity conduction coefficient is taken as the harmonic mean of the air humidity diffusion coefficient and the wood humidity diffusion coefficient. The range and determination method of the air humidity diffusion coefficient are the same as those of the air humidity diffusion coefficient in the case of the aforementioned air medium grid pair. The reference value range and determination method of the wood humidity diffusion coefficient are the same as those of the wood humidity diffusion coefficient in the case of the aforementioned stack surface medium or stack internal medium grid pair.

[0071] If two adjacent grid nodes represent the surface medium property and the internal medium property of the stack, respectively, the cross-grid humidity conduction coefficient is taken as the harmonic mean of the surface layer wood humidity diffusion coefficient and the internal layer wood humidity diffusion coefficient multiplied by the surface-interface conduction reduction factor. The surface layer and the internal layer use the same wood humidity diffusion coefficient. The interface conduction reduction factor ranges from 0.3 to 0.6. This reduction factor is experimentally determined based on the ratio of the surface-interface humidity penetration rate to the internal penetration rate of homogeneous wood. It reflects the increased conduction resistance caused by the discontinuous fiber arrangement, pore blockage, and density gradient at the interface. This reduction factor allows the surface layer to conduct to the internal layer but does not completely block it, thus preserving the formation path of deep humidity accumulation and internal mold risk.

[0072] If two adjacent grid nodes have air medium properties and stack internal medium properties respectively, the cross-grid humidity conduction coefficient is forcibly set to zero due to the violation of the medium continuity assumption, and is not included in the humidity conduction rate field calculation.

[0073] Based on the aforementioned humidity gradient and cross-grid humidity conduction coefficient, for each pair of adjacent grid nodes, the humidity gradient between the two nodes is multiplied by the corresponding cross-grid humidity conduction coefficient to obtain the humidity conduction flux between the two grid nodes, in kilograms per square meter per second. The sign of the flux indicates the direction of conduction; a positive value indicates that humidity is conducted from the high-humidity node to the low-humidity node. For any target grid node, the humidity conduction flux between the target grid node and its six directly adjacent grid nodes is algebraically summed along the three coordinate axes. Flux entering the target node along the positive direction of the coordinate axis is taken as positive, and flux entering the target node along the negative direction of the coordinate axis is equivalent to taking a negative value along the positive direction of the coordinate axis. This yields the net humidity conduction flux components in the three directions. The Euclidean modulus is then calculated for the three components. The modulus is the humidity conduction rate of the target grid node, in kilograms per square meter per second, characterizing the intensity of humidity accumulation or dissipation at that node. A positive net flux indicates humidity accumulation at the node, and a negative net flux indicates humidity dissipation from the node. The magnitude of the modulus reflects the degree of drastic humidity change. The humidity conduction rate of all grid nodes in the three-dimensional mesh structure at a certain moment is stored in spatial coordinates to form the humidity conduction rate field at that moment. The humidity conduction rate field is a three-dimensional scalar field isomorphic to the three-dimensional mesh structure, and each grid node corresponds to a humidity conduction rate scalar value.

[0074] After constructing the humidity conduction rate field, the spore diffusion potential field is calculated. For each grid node in the three-dimensional mesh structure, it is first determined whether the grid node meets the temperature and humidity threshold conditions for fungal spore activation. The conditions are that the temperature value is within the preset spore activation temperature range and the relative humidity value exceeds the preset spore activation humidity threshold. The spore activation temperature range is 5°C to 35°C, and the spore activation humidity threshold is 70% to 80% relative humidity. For grid nodes that meet the conditions, the local spore diffusion potential baseline value is calculated using the product of the temperature response function and the humidity response function. The temperature response function is a Gaussian function. ,in, This is the temperature value. The optimal temperature is set between 20 and 25 degrees Celsius. The standard deviation parameter is denoted by , ranging from 8 to 12 degrees Celsius; the humidity response function is: Type logic functions ,in, Enter the relative humidity value as a percentage. This is the half-saturation constant, which is the relative humidity value corresponding to a humidity response function of 0.5, ranging from 80% to 85%. The slope parameter ranges from 0.1 to 0.3; the local spore dispersal potential baseline value. The value ranges from 0 to 1. For grid nodes that do not meet the conditions, the local spore diffusion potential baseline value is directly set to zero.

[0075] After obtaining the local spore diffusion potential baseline value, for mesh nodes assigned with either the surface medium property or the internal medium property of the stack, the local spore diffusion potential baseline value is further multiplied by a stacking density suppression factor, which adopts an exponential decay form. , The dimensionless spore diffusion density sensitivity coefficient, ranging from 0.8 to 2.0, was determined based on an experimental correlation between air velocity attenuation within the stack and spore settling rate. A higher value indicates greater resistance to spore movement in a high-density environment. The normalized node density is obtained by dividing the grid node stack density value by a reference density value of 500 kg / m³. For air-medium grid nodes, the stack density suppression factor is constant at 1. The local spore diffusion potential baseline value is multiplied by the stack density suppression factor to obtain the local spore diffusion potential value of the grid node, reflecting the tendency of mold spores at that node to diffuse into adjacent spaces.

[0076] After determining the local spore diffusion potential value of each grid node, the spore diffusion driving force between each grid node and its six directly adjacent grid nodes is calculated for each grid node in the three-dimensional grid structure. The driving force for spore dispersal consists of a concentration-driven component and a thermal convection-driven component. The concentration-driven component is defined as the difference in local spore dispersal potential between two adjacent grid nodes divided by the spatial distance between the center coordinates of the two grid nodes. The thermal convection-driven component is defined as the difference in temperature between two adjacent grid nodes divided by the spatial distance to obtain the temperature gradient. This temperature gradient is then divided by a normalized temperature scale, ranging from 10°C to 20°C, reflecting typical seasonal or diurnal temperature variations in the storage environment. The normalized temperature gradient is then multiplied by a thermal convection weighting coefficient to obtain the thermal convection-driven component. This coefficient is dimensionless, ranging from 0.05 to 0.15, and is determined based on a coupling experiment of microscale spore thermophoresis and macroscale air convection, characterizing the proportion of the influence of the temperature gradient on the spore dispersal direction relative to the spore concentration gradient. The spore dispersal driving force is the algebraic sum of the concentration-driven and thermal convection-driven components. It is strengthened when the temperature gradient is aligned with the spore dispersal potential gradient and weakened when they are in opposite directions.

[0077] Based on the driving force of spore diffusion, for each pair of adjacent grid nodes, the spore diffusion effectiveness coefficient is determined according to the medium property identifiers of the two grid nodes. The determination rules are as follows: when both adjacent grid nodes are air medium properties, the spore diffusion effectiveness coefficient is 1.0; when both adjacent grid nodes are stack surface medium properties, it is 0.4 to 0.6; when both adjacent grid nodes are stack interior medium properties, it is 0.1 to 0.3; when both adjacent grid nodes are air medium properties and stack surface medium properties respectively, it is 0.6 to 0.8; when both adjacent grid nodes are air medium properties and stack interior medium properties respectively, it is 0; when both adjacent grid nodes are stack surface medium properties and stack interior medium properties respectively, it is 0.2 to 0.4. The experimental basis for the above values ​​mainly includes: spore diffusion experiments in porous wood media, air velocity attenuation tests inside stacks, spore sedimentation and attachment efficiency determination, measured data of spore concentration gradient distribution in the inner and outer layers of stacks, and analysis of the influence of porosity and pore size distribution on spore throughput.

[0078] Based on the spore diffusion driving force and the spore diffusion effectiveness coefficient, for each pair of adjacent grid nodes, the spore diffusion driving force between the pair of grid nodes is multiplied by the corresponding spore diffusion effectiveness coefficient to obtain the spore diffusion flux, which is a directed vector component along the direction of the line connecting adjacent nodes. A positive value indicates a flow from a high-potential node to a low-potential node, and a negative value indicates the opposite direction. For any target grid node, the spore diffusion flux between the target grid node and its six directly adjacent grid nodes is extracted. If the spore diffusion flux between adjacent nodes connected along a positive direction of a certain coordinate axis is positive, it indicates that spores flow from the adjacent node to the target node, with a positive contribution; if it is negative, it indicates that spores flow from the target node to the adjacent node, with a negative contribution. The six fluxes are then vector-decomposed and algebraically summed along the three coordinate axes to obtain the net spore diffusion flux components in the three directions. , , The Euclidean modulus of the three components is calculated, and this modulus is multiplied by the local spore dispersal potential of the target grid node to obtain the comprehensive spore dispersal potential, which is a dimensionless value that considers both the spore activation conditions and the propagation ability with adjacent spaces. The comprehensive spore dispersal potential of all grid nodes in the 3D grid structure at a certain moment is stored in spatial coordinates to form the spore dispersal potential field at that moment, which is a 3D scalar field isomorphic to the 3D grid structure.

[0079] After obtaining the humidity conduction rate field and the spore dispersal potential field, a spatial convolution operation is performed on both to generate a mold risk propagation tensor. First, a three-dimensional spatial convolution kernel is defined, with a kernel size of [missing information]. × × One grid cell, The number of elements in the convolution kernel is odd, typically 3 or 5. The central element of the kernel corresponds to the target grid node for which the spread intensity of mold risk is to be calculated, and other elements within the kernel correspond to the neighboring grid nodes surrounding the target grid node. The spatial convolution kernel uses a Gaussian weight distribution, where the weight value of each element within the kernel decreases in a Gaussian manner with respect to the spatial distance from that element to the center of the kernel. , Let be the Euclidean distance from a given element within the convolution kernel to the center of the kernel. The spatial scale parameter, ranging from 1 to 2 grid cell side lengths, is determined based on statistical analysis of the mold spread rate and spatial association distance on the stack surface. It characterizes the effective influence radius of the spatial propagation of mold risk, ensuring that the convolution kernel can capture the spatial coupling effect of local risks without causing false associations of distant unrelated nodes due to excessive kernel size. The weight values ​​of all elements in the convolution kernel are normalized so that the sum of all weight values ​​equals 1, ensuring the numerical stability of the convolution operation.

[0080] Using the aforementioned spatial convolution kernel, for any target mesh node in the 3D mesh structure, a spatial convolution kernel is placed centered on that target mesh node. The kernel's coverage area includes the central target mesh node itself and its surrounding neighboring mesh nodes. The humidity conduction rate and comprehensive spore dispersal potential of all mesh nodes within the kernel's coverage area are extracted. The mold risk propagation intensity value of the target mesh node is then calculated using the following formula. :

[0081]

[0082] in, This represents the set of all grid nodes within the coverage area of ​​the spatial convolution kernel centered on the target grid node. For set The grid node number in the data. This serves as the medium attribute identifier for the target mesh node. For the first The medium attribute identifier of each grid node; For the first Humidity-spore coupling risk term for each grid node For the first Normalized Gaussian space weights for each grid node; For the first Cross-media propagation attenuation factor of a grid node to a target grid node.

[0083] Furthermore, The calculation formula is:

[0084]

[0085] in, For the first The humidity conduction rate value of each grid node. For the first The comprehensive spore dispersal potential value of each grid node. To use the humidity conduction rate as a reference, the maximum value of the humidity conduction rate of all grid nodes at the current moment can be taken, or the reference value obtained by calibration experiment under the preset standard temperature and humidity conditions can be taken to make the humidity conduction rate dimensionless. For reference spore dispersal potential; The nonlinear enhancement index for coupling risk, with a value ranging from 1 to 3, enhances the contribution weight of grid nodes with both high humidity conduction and high spore diffusion to the risk of the target node through exponential transformation. When the value is 1, it degenerates into a linear product. When the value is greater than 1, the contribution of the nonlinear amplification to the highly coupled node; As a logarithmic stable term, it ensures that the value inside the parentheses is always positive during logarithmic operations.

[0086] The calculation formula is:

[0087]

[0088] in, For the first The three-dimensional coordinate vector of each grid node. Let be the three-dimensional coordinate vector of the target mesh node. For spatial scale parameters, To normalize the stable term and prevent the denominator from being zero.

[0089] For the first The cross-medium propagation attenuation factor of the nth grid node to the target grid node, according to the nth The media property combinations of each grid node and the target grid node are discretely selected, and the specific rules are as follows:

[0090] like and The same, meaning that both grid nodes have air medium properties, both have stack surface medium properties, or both have stack interior medium properties. ;

[0091] like and When the properties are air medium properties and stack surface medium properties respectively, or when the properties are stack surface medium properties and air medium properties respectively, The attenuation factor is between 0.5 and 0.7, and is based on the determination of the spore settling and attachment efficiency and reverse release efficiency at the air-wood interface.

[0092] like and When the properties are respectively air medium properties and stack internal medium properties, or stack internal medium properties and air medium properties, the media continuity assumption is violated. =0;

[0093] like and When the media properties are respectively the surface media properties and the internal media properties of the stack, or respectively the internal media properties and the surface media properties of the stack, The attenuation factor is set to 0.3 to 0.5, and is based on the measured data of the spore concentration gradient distribution in the inner and outer layers of the stack.

[0094] For the target mesh node at the center of the convolution kernel, since the medium properties are the same... .

[0095] The mold risk transmission intensity value calculated using the above formula The value is dimensionless and comprehensively reflects the humidity accumulation and spore diffusion conditions of the target grid node itself, as well as the coupling effect of the humidity accumulation and spore diffusion conditions of the surrounding neighborhood space on the target node. The humidity conduction rate value and the comprehensive spore diffusion potential value in the formula are calculated according to the aforementioned definition. The cross-medium propagation attenuation factor is consistent with the aforementioned discrete attenuation rule, avoiding the inconsistency of rules caused by the incomplete definition of medium directionality.

[0096] After calculating the mold risk propagation intensity value of all target grid nodes, the mold risk propagation intensity value of all grid nodes in the three-dimensional grid structure at a certain moment is stored according to spatial coordinates, forming the mold risk propagation tensor at that moment.

[0097] The mold risk propagation tensor is a three-dimensional scalar field isomorphic to the three-dimensional mesh structure, with each mesh node corresponding to a scalar value of mold risk propagation intensity. For multiple consecutive time slices, the mold risk propagation tensor corresponding to each time slice is calculated in chronological order and stored by time dimension index, forming a four-dimensional mold risk propagation tensor data structure containing three spatial dimension indices and one time dimension index. In the four-dimensional mold risk propagation tensor data structure, the mold risk propagation intensity value of any mesh node under any time index can be read through the addressing of the mesh node's spatial index and time dimension index. This mold risk propagation intensity value comprehensively characterizes the degree of risk of the mesh node being in an active mold state at the corresponding time and the ability of the risk to propagate and spread to the surrounding space.

[0098] S3: Extract the temporal difference features of the mold risk propagation tensor within a sliding window along the time axis, identify spatial cluster regions in the risk propagation tensor that exhibit a cumulative enhancement pattern, and mark the spatial cluster regions as a set of high-risk propagation source points.

[0099] Obtain the mold risk propagation tensor generated in step S2. This tensor is a four-dimensional data structure containing three spatial dimension indices and one temporal dimension index. Each element records the mold risk propagation intensity value of the corresponding spatial grid node under the corresponding time slice. Set the sliding window length according to the total number of steps in the temporal dimension of the mold risk propagation tensor. Sliding window length Take one-tenth to one-fifth of the total number of time steps, but not less than 3 time steps; if the total number of time steps is less than 3, set the sliding window length to be equal to the total number of time steps and perform single-window analysis.

[0100] The sliding window is initially positioned at the beginning of the time series, so that the window covers the continuous time series starting from the initial time. For each grid node in the 3D spatial grid, the risk propagation intensity value sequence of that grid node across all time steps covered by the current sliding window is extracted from the mold risk propagation tensor. The sequence is arranged in chronological order. The extracted risk propagation intensity value sequence is then subjected to adjacent time step difference operations to obtain the risk intensity increment sequence. The nth time step in the increment sequence is... The element is the first The time step and the first The difference in risk transmission intensity values ​​at each time step. Three characteristics are calculated based on the incremental risk intensity sequence: the first is the average risk change rate, defined as the arithmetic mean of all elements in the incremental sequence; the second is the proportion of positive increments, defined as the number of elements greater than zero in the incremental sequence divided by the number of non-zero elements in the incremental sequence. If the number of non-zero elements in the incremental sequence is zero, the proportion of positive increments is recorded as zero; the third is the cumulative sum of positive increments, defined as the sum of all elements greater than zero in the incremental sequence.

[0101] After completing the feature calculations for all grid nodes within the current sliding window, each grid node is assessed to determine if it meets the criteria for a cumulative enhancement mode. These criteria include: an average risk change rate greater than a preset positive change threshold, a positive increment percentage greater than 0.6, and a cumulative sum of positive increments greater than 0.05 times the arithmetic mean of the risk propagation intensity values ​​of that grid node within the current window. Grid nodes meeting all of these criteria are marked as candidate nodes for cumulative enhancement within the current time window. The preset positive change threshold is either one times the standard deviation of the risk propagation intensity value sequence for that grid node within the current window, or the larger of 0.02.

[0102] It should be noted that the preset positive change threshold is determined based on the statistical distribution of the time difference of the mold risk propagation intensity value under normal storage environment fluctuations. It is higher than the daily random fluctuation range to filter out risk value fluctuations caused by small disturbances in temperature and humidity, and only retain the continuous upward trend with practical prevention and control significance. The proportion of positive increments is greater than 0.6, which means that the proportion of risk increase time steps in the window exceeds 60% of all non-zero change time steps. It is set based on the statistical observation of the alternation frequency of the rising stage and the fluctuation stage in the actual mold risk evolution process, and can distinguish between occasional short-term increases and continuous high-frequency increases. The cumulative sum of positive increments is greater than 0.05 times the arithmetic mean of the risk propagation intensity value of the grid node in the current window. It is used to exclude micro-fluctuations with too small absolute increments, so that the judgment focuses on nodes with substantial increases in risk intensity.

[0103] Spatial clustering analysis is performed on all cumulative enhancement candidate nodes marked within the current time window. Spatial clustering employs a seed expansion algorithm based on spatial connectivity. The neighborhood is defined as the six directly adjacent grid nodes sharing the same surface as the node, consistent with the spatial neighborhood definition in step S2, ensuring the spatial continuity of the medium property spatial structure and risk propagation path. Starting with any candidate node not assigned to a cluster, it is used as a seed node. The search continues among the six directly adjacent grid nodes of this seed node to determine if any candidate nodes exist. If so, the adjacent candidate nodes are added to the current cluster, and the newly added nodes are used as new seed nodes to continue the search until no unadded candidate nodes remain in the neighborhood of any node in the current cluster. The cluster expansion ends at this point. For each formed cluster, if the number of nodes within the cluster is less than the preset minimum cluster size of 3, the cluster is removed, and the nodes within the cluster are marked as noise points. This process is repeated until all candidate nodes are assigned to a cluster or marked as noise points, completing the spatial cluster region identification within the current time window.

[0104] For each identified spatial cluster region, the average of the average risk change rate of all nodes within the cluster and the sum of the cumulative positive increments are calculated as the cumulative enhancement characteristic index of the cluster. Simultaneously, the spatial cluster region of the current time window is matched with the spatial cluster region of the previous adjacent time window. The matching method is to calculate the node overlap rate between the corresponding spatial clusters of the two windows, defined as the number of overlapping nodes divided by the total number of nodes in the smaller cluster. If the node overlap rate is greater than 0.7, the two spatial clusters are considered to be continuations of the same continuous risk evolution cluster in different time windows. The number of time windows in which this continuous risk evolution cluster persists is recorded. If the number of time windows in which it persists is greater than or equal to 2, and the average of the average risk change rate of the cluster in each time window is positive, and the sum of the cumulative positive increments is greater than the average of the cumulative positive increments of all nodes in the corresponding time window, then the spatial cluster region is confirmed to exhibit a risk accumulation enhancement pattern with temporal stability and persistence.

[0105] The parameters involved in the above-mentioned judgment process can be calibrated in actual applications according to the specific warehouse scale, stacking structure, monitoring time step and risk level management requirements, but the relative relationship and magnitude of each parameter remain unchanged to ensure that the judgment process has physical consistency and engineering portability.

[0106] After completing the cumulative enhancement pattern identification and spatial cluster matching for all time window locations, all confirmed risk cumulative enhancement spatial cluster regions with temporal stability and persistence are marked as a set of high-risk transmission sources.

[0107] The three-dimensional spatial location information of each spatial cluster in the set of high-risk transmission sources, the number of nodes in the cluster, the node number in the cluster, the average risk change rate of the cluster as a whole, the cumulative sum of positive increments of the cluster as a whole, the number of persistent time windows, and the index range of the time windows are recorded and output in a structured manner.

[0108] S4: Spatiotemporally couple the set of high-risk transmission sources with the future time-period evolution trend of the environmental parameter distribution field to generate a partitioned and layered distribution map of the risk level of mold growth in timber stacks.

[0109] Obtain the set of high-risk propagation sources output in step S3. This set includes all spatial clusters of risk accumulation enhancement with temporal stability and persistence. Each spatial cluster records its three-dimensional spatial location information, the number of nodes within the cluster, node numbers within the cluster, the overall average risk change rate of the cluster, the overall cumulative positive increment of the cluster, the number of persistent time windows, and the time window index range. Simultaneously, obtain the future environmental parameter distribution field. The future environmental parameter distribution field is obtained through time series prediction of historical environmental monitoring data. The prediction model adopts an autoregressive moving average model. Specifically:

[0110] For each node in the three-dimensional spatial grid, the measured observation sequences of temperature, humidity, and airflow velocity over historical periods are extracted. Each historical observation sequence of environmental parameter is treated as an independent time series, and the stationarity of the sequence is tested. If the sequence shows obvious non-stationary characteristics, it is differentially processed until the sequence tends to be stationary. Then, the autocorrelation structure and partial autocorrelation structure are analyzed based on the stationary sequence. The autoregression order and moving average order are determined according to the truncation or tailing behavior of autocorrelation and partial autocorrelation, so as to determine a suitable autoregressive moving average model structure for the parameter of the node. After the model structure is determined, the model parameters are estimated and fitted using historical observation data. In the rolling prediction phase, historical observation data is first used as the initial input to predict the environmental parameter value for the first future time step. After obtaining the predicted value for the first future time step of the current node, the predicted value of the current node is spatially weighted and corrected using the predicted values ​​generated by neighboring nodes at the same time step. The weight of the spatial weighting correction is the normalized result of the negative square power of the spatial distance between the current node and the neighboring nodes, so that the prediction result is continuously affected by the changes in the surrounding spatial environment in the gradual deduction. Then, the corrected predicted value is included in the input window, replacing the earliest historical observation value to form an updated input sequence. This is then used to predict the environmental parameter value for the second time step, and so on, rolling forward until the prediction of all future time steps is completed.

[0111] Based on the output of the autoregressive moving average model, a future environmental parameter distribution field with the same three-dimensional spatial grid structure as the mold risk propagation tensor is formed. Each grid node corresponds to a sequence of predicted temperature, humidity, and airflow velocity values ​​covering several future time steps. If there is insufficient historical data to support predictive modeling, a constant assumption method is adopted, assuming that the environmental parameters of each node in the future will maintain the same distribution as at the current time, or that the distribution of environmental parameters will be adjusted according to known environmental control strategies.

[0112] After obtaining the distribution field of future environmental parameters, the risk transmission potential score is calculated for each spatial cluster in the set of high-risk transmission sources. For each spatial cluster, the risk growth rate normalization score is obtained by dividing the cluster's overall average risk change rate by the maximum average risk change rate among all clusters; the total risk normalization score is obtained by dividing the cluster's overall cumulative positive increment by the maximum cumulative positive increment among all clusters; the time stability normalization score is obtained by dividing the number of persistence windows by the maximum number of persistence windows among all clusters; and the spatial scale normalization score is obtained by dividing the number of nodes within the cluster by the maximum number of nodes among all clusters. The risk propagation potential score for the cluster is obtained by weighted summation of these four normalization scores. The weighting coefficients are denoted as follows: , , , ,in The normalized score corresponding to the risk growth rate Corresponding to the normalized score of total risk, Corresponding to the time stability normalized score, Corresponding spatial scale normalized score; weight coefficients satisfy ,and ∈[0.25,0.45], ∈[0.25,0.45], ∈[0.05,0.25], ∈[0.05,0.25], at the same time ≥ , ≥ , ≥ , ≥ This ensures that the contribution of the rate of risk growth and the total amount of risk to the propagation potential is no less than that of temporal stability and spatial scale. For example, it is advisable to... =0.35, =0.35, =0.15, =0.15, but this scheme is not limited to this specific value. In practice, it can be adjusted based on the relevance of historical data or expert experience, provided that the above constraints are met. Risk transmission potential score. The value ranges from 0 to 1, with higher scores indicating a stronger overall ability of the cluster to propagate as a risk source in the future. After scoring all spatial clusters, a reference value for propagation potential is set. Lower limit of propagation potential ,Will Divide by Gaining the potential of regularized propagation , ∈[0.01,1], to prevent zero values ​​after normalization from causing abnormal calculations of the space propagation kernel.

[0113] After obtaining the regularized propagation potential of each source point cluster, the environmental suitability function is calculated for each grid node in the three-dimensional spatial grid at each future time step. It is composed of a temperature suitability component, a humidity suitability component, and an airflow inhibition component. The formula is expressed as:

[0114]

[0115] in, Let the temperature suitability function be... Let the humidity suitability function be... This is the airflow suppression function. The value ranges from 0 to 1. The larger the value, the more favorable the environmental conditions at that time step are for the spread of mold risk.

[0116] Furthermore, the calculation of the temperature suitability component includes:

[0117] For future time steps The grid nodes below Extract the predicted temperature value of this node. Humidity forecast value and airflow velocity prediction The formula for calculating dimensionless temperature is as follows:

[0118]

[0119] in This represents the lower limit of the effective temperature window for mold growth. This represents the upper limit of the effective temperature window for mold growth. This dimensionless temperature is... It is positive when the temperature is between 5 and 35 degrees Celsius, with a maximum value of 1, and negative when it is outside the window.

[0120] Then, the temperature suitability function is calculated, expressed by the formula:

[0121]

[0122] in, This is a dimensionless temperature sensitivity coefficient.

[0123] Furthermore, the dimensionless humidity is calculated using the following formula:

[0124]

[0125] in To determine the activation humidity threshold for propagation, a relative humidity of 70% is used. This is a reference humidity scale; the dimensionless humidity is... When it is 0, it is lower than When it is negative, it is higher than The time is right.

[0126] Then, a humidity suitability function is constructed, expressed by the formula:

[0127]

[0128] in, This is a dimensionless humidity sensitivity coefficient.

[0129] Furthermore, the dimensionless airflow intensity is calculated using the following formula:

[0130]

[0131] in, For reference wind speed, the average airflow speed of the storage space under normal ventilation conditions is taken. If it is not possible to obtain the average airflow speed, then 0.1 meters per second is taken.

[0132] Then, an airflow suppression function is constructed, expressed by the formula:

[0133]

[0134] in, It is a dimensionless airflow suppression coefficient.

[0135] It should be noted that the above , , The values ​​must be determined through calibration or sensitivity analysis based on the specific storage environment. Control the rate at which the temperature suitability function rises from zero to the saturation value of 1. A larger value indicates that the system is more sensitive to the optimal temperature range. Propagation rate data under different temperature gradients were measured experimentally, and nonlinear least squares fitting was used. Obtained by curve; Determines the Sigmoid-type humidity suitability function at the critical humidity The steepness of the nearby transition, A higher threshold value indicates a more significant humidity threshold effect, and its calibration requires data collection spanning [a certain range]. The humidity-propagation activity pairing data before and after are solved using logistic regression or gradient optimization; Quantifying the inhibitory effect of airflow velocity on the propagation process, establishing the wind speed-propagation attenuation relationship through wind tunnel experiments or field observations, and then... and When performing linear regression, the absolute value of the regression slope is... Typical applications , , The values ​​range from 0.5 to 5 to ensure that the fitness function can effectively distinguish differences in environmental conditions without causing numerical oscillations due to oversensitivity.

[0136] After completing the environmental suitability calculations for each node at each future time step, a spatiotemporal coupled simulation is executed to model the time-step risk diffusion behavior of high-risk propagation sources under future environmental conditions. Before the simulation, the baseline propagation characteristic scale needs to be determined. The method of determination is as follows:

[0137] First, time segments with clearly defined risk source locations, continuous risk expansion processes, and complete contemporaneous environmental monitoring data are selected from the historical mold risk propagation tensor. For each segment, boundary nodes are determined centered on the identified high-risk spatial clusters. Grid nodes outside the clusters in different directions are arranged from closest to furthest from the cluster boundary. The risk intensity changes of each node within a continuous time step are extracted. If the risk increase of a node occurs after the risk of the source cluster increases, it is considered a node that may be affected by the propagation of that source cluster. Risk increases caused by localized independent environmental anomalies are excluded by combining changes in airflow direction, temperature, and humidity, thus identifying effective propagation nodes.

[0138] After determining the effective propagation nodes, the risk intensity at the source cluster boundary is used as a reference to calculate the relative risk retention at each distance. The arithmetic mean of the relative risk retention of nodes in the same direction within the same distance range is then taken to reduce the impact of monitoring fluctuations and spatial noise from individual nodes. If there are significant differences in risk attenuation in different directions within the same distance range, the risk intensity variation sequences with distance reduction are organized separately according to the main airflow direction, the reverse airflow direction, and the lateral direction. If the directional differences are not significant, the data from each direction are merged. The organized distance attenuation sequences are smoothed to identify the main segment where the risk intensity continuously decreases from the source cluster outwards. An exponential attenuation model is used to fit this segment, and the spatial distance corresponding to when the relative risk retention decreases to 1 / e of the risk intensity at the source cluster boundary is determined as the candidate attenuation feature length for this segment. 1 / e is approximately 36.8%, representing a typical amplitude attenuation of risk intensity.

[0139] When multiple historical segments have obtained candidate attenuation feature lengths, a data integrity score, a propagation process clarity score, and a fitting stability score are calculated for each segment. All three scores are between 0 and 1, and the larger the value, the more reliable the calibration result of the segment. Data integrity score is determined based on the proportion of missing monitoring records within a segment. A score of 1 is given when the missing proportion is less than 5%, 0.8 is given when the missing proportion is between 5% and 15%, 0.5 is given when the missing proportion is between 15% and 30%, and 0.2 is given when the missing proportion exceeds 30%. Propagation process clarity score is determined based on the number of effective propagation nodes and the continuity of the propagation path. A score of 1 is given when the number of effective propagation nodes exceeds 20 and the propagation path has no obvious interruption, 0.7 is given when the number of effective propagation nodes is between 10 and 20 or the propagation path has a small number of interruptions, and 0.4 is given when the number of effective propagation nodes is less than 10 or the propagation path has an obvious interruption. Fit stability score is determined based on the sum of squared residuals of the exponential decay model. The sum of squared residuals is normalized to between 0 and 1. The smaller the sum of squared residuals, the higher the score. A score of 1 is given when the sum of squared residuals is zero, and a score of 0.2 is given when the sum of squared residuals is greater than 0.5. The remaining cases are determined by linear interpolation. The three scores are multiplied together to obtain the overall weight of the segment. Then, the candidate decay feature lengths of all segments are weighted and averaged according to the overall weight to obtain the final decay feature length.

[0140] If the difference between the maximum and minimum candidate attenuation feature lengths in different directions within the same historical segment does not exceed 20% of the smaller value, the directional differences are considered insignificant, and the arithmetic mean of the candidate values ​​in each direction is taken as the unified candidate value for that segment. If the difference exceeds 20%, the candidate values ​​in each direction are retained. When summarizing multiple segments, for cases where directional candidate values ​​are retained, each direction is first weighted and averaged according to a comprehensive weight to obtain the final value for each direction. If the difference between the maximum and minimum final values ​​in each direction still does not exceed 20% of the smaller value, the arithmetic mean of the final values ​​in each direction is taken as the unified attenuation feature length. If the difference exceeds 20%, the directional attenuation feature length is retained, and the corresponding value is selected in subsequent propagation calculations based on the direction of the target node relative to the source cluster. The obtained unified attenuation feature length or directional attenuation feature length is constrained within the range of 3 to 8 grid cell side lengths; if it exceeds this range, the nearest boundary value is taken. Finally, historical propagation segments that were not involved in the acquisition of the above parameters were replayed for verification. The risk intensity, spatial distance, and environmental conditions of the source cluster were substituted into the propagation model to compare the consistency between the simulated risk impact range and spatial attenuation trend and the actual observation results. When the simulation results are consistent with the actual observations, the obtained attenuation feature length is used as the formal parameter in the spatiotemporal coupling deduction. Otherwise, the selection of effective propagation nodes, the identification of the main attenuation section, and the abnormal data are re-examined, and the model is fitted and verified again after correction.

[0141] Determine the reference propagation characteristic scale Then, a spatiotemporal coupling simulation was performed to simulate the time-step risk diffusion behavior of high-risk propagation sources under future environmental conditions. Specifically:

[0142] For each time step in the future Traverse every spatial cluster in the set of high-risk transmission sources. Risk propagation calculations are performed starting from the spatial location of the cluster. For each target node in the three-dimensional spatial mesh... Calculate the relationship between this node and the current source node cluster. Spatial Euclidean distance The spatial distance is calculated by taking the minimum distance between the target node and all nodes within the cluster. This is based on the regularized propagation potential of the current source node cluster. and the defined reference propagation characteristic scale The effective propagation radius is calculated using the following formula:

[0143]

[0144] in The stretching index, ranging from 0 to 1, controls the shape of the attenuation curve. As the contribution cutoff threshold, take ; It is a constant factor.

[0145] If a directional attenuation feature length exists, the corresponding value is selected according to the direction of the target node relative to the source node cluster. If the contribution of the source cluster to the propagation risk intensity of the target node at this time step is set to zero, it means that the risk propagation has decayed to a negligible level.

[0146] like The basic propagation risk intensity is then calculated using the KWW spatial propagation kernel, expressed by the formula:

[0147]

[0148] The basic propagation risk intensity takes into account the propagation potential of the source cluster. And the stretching index decays with spatial distance.

[0149] After obtaining the basic propagation risk intensity, multiply it by the target node. At the current time step Environmental suitability function This yields the actual transmission risk intensity after environmental modulation.

[0150] If multiple source clusters exist, the actual propagation risk intensity of the target node at the current time step is linearly superimposed by all source clusters to obtain the total propagation risk intensity of the target node at the current time step, expressed by the formula:

[0151]

[0152] After completing the calculation for the current time step for all nodes in the three-dimensional spatial grid, proceed to the next future time step and repeat the above process until the derivation of all future time steps is completed.

[0153] After completing the projection of all future time steps, for each grid node The cumulative propagation risk intensity value of a node over the entire future period is obtained by summing the total propagation risk intensity values ​​received at all future time steps over the time dimension. ,in This represents the total number of future time steps. The risk intensity value of the last time slice in the mold risk propagation tensor for this node is extracted as the risk intensity baseline value. ,Will Divide by all nodes The maximum value yields the normalized cumulative propagation risk intensity value. Then, the comprehensive risk level assessment value is calculated, expressed by the formula:

[0154]

[0155] in, As the weight of the risk intensity benchmark value, Weights are assigned to the cumulative transmission risk intensity value.

[0156] The risk intensity benchmark and cumulative propagation risk intensity values ​​for future periods are determined through historical sample playback: multiple historical time periods are selected, and the starting point of each historical time period is used as the simulated current time. The benchmark risk intensity value for each node at that time and the cumulative propagation risk intensity value for future periods are calculated to form a weighted calibration sample; candidate samples are set. The set is For each candidate The value is used to calculate the comprehensive risk level assessment value for each node across all historical periods. and according to Risk levels were assigned, and the results were compared with subsequent actual mold growth observations. Two evaluation indicators were calculated: the first was the current risk hit rate, i.e., the percentage of historical periods where the measured risk was high and... The first is the proportion of nodes classified as high-risk to the total number of high-risk nodes measured at the start of the historical period; the second is the future risk prediction rate, i.e., the proportion of nodes that are subsequently added or expanded to be high-risk. The proportions have been pre-classified as high-risk or medium-high-risk. For each candidate... The value is calculated by averaging the two indicators over all historical periods. The averages of the two indicators are then combined in equal proportions to form a composite score, and the candidate with the highest composite score is selected. The value is used as the weighting coefficient in practice; if the overall scores are tied for the highest, the candidate value with the higher current risk hit rate is selected; if the current risk hit rates are also the same, the value with the higher current risk hit rate is selected. Candidate values ​​with larger values. For example, when historical sample playback results show... When the value is 0.7, the combined current risk hit rate and future risk prediction rate are optimal. Take 0.7, The value is set to 0.3, but the actual weight should be based on the results of historical sample playback.

[0157] Complete the comprehensive risk level assessment of all nodes in the three-dimensional spatial mesh. After calculation, calculate all nodes. The upper quartile is used as the high-risk threshold, and the median is used as the medium-risk threshold. If a node's... If the risk level is above the upper quartile, it is classified as high-risk; if it is between the median and the upper quartile, it is classified as medium-risk; if it is below the median, it is classified as low-risk.

[0158] After classifying the risk levels of all nodes in the 3D spatial mesh, a seed expansion algorithm based on risk level consistency is used for spatial clustering to identify continuous risk partitions. Starting with any node not yet assigned to a partition, it is used as a seed node. The algorithm checks if any of its six directly adjacent mesh nodes have a neighbor node with the same risk level as the seed node. If a neighbor with the same risk level exists, it is added to the current partition and used as a new seed node to continue expansion. If no neighbor with the same risk level exists, or if all adjacent nodes have been assigned, the expansion of the current partition ends. This process is repeated until all nodes are assigned to a risk partition, completing the planar risk partitioning.

[0159] After completing the planar risk zoning, vertical layering is performed based on the vertical coordinates of the 3D spatial grid. The vertical coordinate range is divided into three equal layers: the bottom layer corresponds to the minimum one-third of the vertical coordinate range, the middle layer corresponds to the middle one-third of the vertical coordinate range, and the top layer corresponds to the maximum one-third of the vertical coordinate range. Risk zoning statistical analysis is performed within each vertical layer, calculating the spatial proportions of high-risk, medium-risk, and low-risk zones within that layer. If the spatial proportion of high-risk zones in the bottom layer is greater than 1.5 times the larger of the high-risk zone proportions in the middle and top layers, the bottom layer is marked as a high-risk layer; if the spatial proportion of low-risk zones in the top layer is greater than 1.5 times the larger of the low-risk zone proportions in the bottom and middle layers, the top layer is marked as a low-risk layer; if neither of these conditions is met, the middle layer is marked as a medium-risk layer; if both the bottom and top layers meet their respective conditions, they are marked as high-risk and low-risk layers respectively, while the middle layer remains marked as a medium-risk layer.

[0160] Based on the combined results of planar risk zoning and vertical stratification analysis, a risk level distribution map of timber stacks with zoning and stratification is generated. This map uses 3D volumetric rendering technology to display the risk level distribution at various spatial locations inside and on the surface of the timber stacks. High-risk areas are marked with red, medium-risk areas with yellow, and low-risk areas with green. The map overlays the location markers of spatial clusters within the high-risk transmission source set, using special symbols or bounding boxes to indicate the 3D spatial extent of each cluster. Each risk zone in the map is labeled with its zone number, vertical stratification name, risk level category, minimum and maximum ranges of the comprehensive risk level assessment value, number of nodes within the zone, and zone volume percentage.

[0161] If it is necessary to demonstrate the dynamic evolution of the risk level distribution, the total propagation risk intensity value distribution data saved in the aforementioned time-step extrapolation process for each future time step is used to generate a separate zonal and stratified risk level distribution map for each future time step. The generation method for each time step distribution map is the same as the aforementioned comprehensive assessment method, but when calculating the comprehensive risk level assessment value, the total propagation risk intensity value of that time step is used. Replace cumulative transmission risk intensity value The baseline value for risk intensity remains the current risk intensity value. The distribution maps for all future time steps are organized chronologically to form a dynamic evolution sequence.

[0162] The final output of the partitioned and layered timber stack mold risk level distribution map and its dynamic evolution sequence provides timber storage managers with a spatial distribution view of risk and a temporal evolution trajectory. This is used to identify areas with severe mold risk threats and take intervention measures such as local ventilation enhancement or humidity control. For areas with rising risk, the monitoring frequency and early warning response are increased, while for areas with controllable risk, the existing management strategy is maintained and resource allocation is optimized.

[0163] like Figure 2 As shown in the figure, the spatial distribution of risk levels of each grid node inside and on the surface of the stack is displayed in a 3D volumetric drawing method under an 8×6×5 spatial grid. Among them, the high-risk areas (comprehensive risk assessment value) Areas located above the upper quartile are marked in red, primarily concentrated on the bottom floors and in poorly ventilated areas; medium-risk areas ( Areas located between the median and the upper quartile are marked in yellow, indicating a spread from high-risk sources to the periphery; low-risk areas ( Areas below the median (marked in green) are mainly distributed in the top layer and well-ventilated outer areas. The map also uses black pulse markers to overlay the locations of three initial mold clusters identified from the high-risk propagation source set, illustrating the multi-source risk diffusion calculation process based on the KWW spatial propagation kernel. This distribution map visually demonstrates that the method of this invention can output spatially precise, zonal, and stratified risk level results, achieving refined risk assessment across three vertical levels (bottom, middle, and top) and multiple planar zones. Compared to traditional single-point time-series statistical methods that only provide the overall risk level of the warehouse or the mold index of isolated monitoring points, this invention ensures that managers can directly identify specific spatial areas in the stack that require priority intervention (as shown in the red high-risk concentration area in the figure), and implement differentiated zonal control strategies accordingly. This avoids resource waste caused by over-control of low-risk areas, significantly improving the accuracy and economy of mold control in large-scale timber stacks.

[0164] Those skilled in the art will understand that the specific operations of each step in the above-described method for predicting the risk of mold growth in timber storage have been referenced above. Figure 1 and Figure 2 The method for predicting the risk of mold growth in timber storage is described in detail here, and therefore, its repeated description will be omitted.

[0165] According to another aspect of the present invention, a timber storage mold risk prediction system is provided, comprising:

[0166] The field construction module is used to receive temperature and humidity time series data from each monitoring point of the timber stack, and to perform grid-based mapping of the temperature and humidity time series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space.

[0167] The risk propagation module is used to calculate the humidity conduction rate field and spore diffusion potential field across grids based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, and to perform spatial convolution operation on the humidity conduction rate field and the spore diffusion potential field to generate a mold risk propagation tensor.

[0168] The source identification module is used to extract the gradient change rate of the mold risk propagation tensor within a sliding window along the time axis, identify spatial cluster regions in the risk propagation tensor that exhibit a cumulative enhancement pattern, and mark the spatial cluster regions as a set of high-risk propagation sources.

[0169] The zoning simulation module is used to perform spatiotemporal coupling simulation of the future evolution trend of the set of high-risk transmission sources and the distribution field of environmental parameters, and generate a zoning and stratified distribution map of the risk level of mold growth in timber stacks.

[0170] In summary, the timber storage mold risk prediction method based on the embodiments of the present invention has been clarified. By performing spatial gridding modeling and risk propagation dynamics analysis on the time-series temperature and humidity data of each monitoring point in the timber stack, it can accurately capture the dynamic distribution of humidity and temperature within the storage space and their impact on mold risk. Simultaneously, by employing risk propagation tensor and sliding window analysis, it effectively identifies spatial clusters of areas where risk accumulation is enhanced, enabling early location of high-risk propagation sources and zoning-based risk level determination. This not only improves the spatial resolution and timeliness of mold risk monitoring but also helps managers develop more targeted prevention and control measures for high-risk areas, thereby reducing the incidence and extent of mold damage during timber storage.

Claims

1. A method for predicting the risk of mold growth in timber storage, characterized in that, include: Receive temperature and humidity time-series data from each monitoring point of the timber stack, and perform gridded mapping of the temperature and humidity time-series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space. Based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, the humidity conduction rate field and spore diffusion potential field across the grid are calculated, and the humidity conduction rate field and the spore diffusion potential field are spatially convolved to generate a mold risk propagation tensor. The gradient rate of change of the mold risk propagation tensor is extracted within a sliding window along the time axis. Spatial cluster regions that exhibit a cumulative enhancement pattern in the risk propagation tensor are identified, and the spatial cluster regions are marked as a set of high-risk propagation sources. By spatiotemporally coupling the set of high-risk transmission sources with the future time-period evolution trend of the environmental parameter distribution field, a partitioned and layered distribution map of the risk level of mold growth in timber stacks is generated.

2. The method for predicting the risk of mold growth in timber storage according to claim 1, characterized in that, When constructing the environmental parameter distribution field, the timber storage space is discretized into a three-dimensional grid structure based on the three-dimensional geometric boundary data of the timber storage space. The side length of the grid cell in any dimension is taken as the smaller of one-third to one-half of the average spacing between sensors in the air channel between adjacent timber stacks and one-third to one-half of the minimum spacing between sensors inside the stack. Based on the shortest spatial distance from each grid node to the geometric boundary of the timber stack, a medium attribute is assigned: those located outside the boundary are assigned air medium attributes, those located inside the boundary and whose shortest distance does not exceed the preset surface layer depth are assigned stack surface medium attributes, and those exceeding the preset surface layer depth are assigned stack interior medium attributes.

3. The method for predicting the risk of mold growth in timber storage according to claim 2, characterized in that, When calculating the humidity conduction velocity field, the cross-grid humidity conduction coefficient is determined for each pair of adjacent grid nodes based on their medium property identifier: If both nodes are air-based media, the air humidity diffusion coefficient is used. If both nodes have surface or internal stacking media properties, then the base value of the wood moisture diffusion coefficient and the density correction factor are used. The product of, where The normalized stacking density of the two nodes is the arithmetic mean. Density attenuation coefficient; If the two nodes represent the surface medium properties and the interior medium properties of the stack, respectively, the harmonic mean of the wood moisture diffusion coefficients of the two nodes is multiplied by an interfacial conduction reduction factor ranging from 0.3 to 0.

6. If the two nodes have air medium properties and stack internal medium properties respectively, the cross-grid humidity conduction coefficient is forced to zero.

4. The method for predicting the risk of mold growth in timber storage according to claim 3, characterized in that, When calculating the spore dispersal potential field, for grid nodes that meet the conditions of spore activation temperature range and relative humidity exceeding the spore activation humidity threshold, the local spore dispersal potential benchmark value is calculated by multiplying the temperature response function and the humidity response function. ; The temperature response function adopts a Gaussian form. ,in, This is the temperature value. For the optimal temperature, The standard deviation parameter; The humidity response function uses a S-shaped logic function. ,in, This is the relative humidity value. It is the half-saturation constant. The slope parameter; For mesh nodes that are assigned surface media properties or internal media properties of the stack, Multiply by the stacking density suppression factor , This is the spore diffusion density sensitivity coefficient. This represents the normalized node density. For grid nodes that do not meet the conditions for spore activation, Set it directly to zero.

5. The method for predicting the risk of mold growth in timber storage according to claim 4, characterized in that, When generating the aforementioned mold risk propagation tensor, weight values ​​are used. The normalized 3D Gaussian convolution kernel is used to calculate the mold risk propagation intensity value for each target grid node using the following formula. : in, This represents the set of all grid nodes within the coverage area of ​​the spatial convolution kernel centered on the target grid node. For set The grid node number in the data. This serves as the medium attribute identifier for the target mesh node. For the first The medium attribute identifier of each grid node; For the first Humidity-spore coupling risk term for each grid node For the first Normalized Gaussian space weights for each grid node; For the first The cross-medium propagation attenuation factor of each grid node to the target grid node; The formula for calculating the humidity-spore coupling risk term is as follows: in, For the first The humidity conduction rate value of each grid node. For the first The comprehensive spore dispersal potential value of each grid node. Reference humidity conduction rate; For reference spore dispersal potential; This is a nonlinear enhancement index for coupling risk, with a value ranging from 1 to 3. When the value is greater than 1, the contribution weight of nodes where nonlinear amplification of high humidity conduction and high spore diffusion coexist is calculated. It is a log-stable term.

6. The method for predicting the risk of mold growth in timber storage according to claim 5, characterized in that, During the spatiotemporal coupling simulation, for each grid node in the three-dimensional spatial grid... In every future time step Calculate the environmental suitability function using the following formula : in, Let the temperature suitability function be... Let the humidity suitability function be... It is the airflow suppression function; The temperature suitability function is: in, This is a dimensionless temperature sensitivity coefficient. For dimensionless temperature, it is expressed as , These are the lower and upper boundaries of the effective temperature window for mold growth. This is a predicted temperature value; The temperature suitability function is: in, The dimensionless humidity sensitivity coefficient Dimensionless humidity, expressed as , To propagate the activation humidity threshold, This is a humidity forecast value; The airflow suppression function is: in, This is a dimensionless airflow suppression coefficient. The dimensionless airflow intensity is expressed as ,in For reference wind speed, This is the predicted airflow velocity.

7. The method for predicting the risk of mold growth in timber storage according to claim 6, characterized in that, The generation of the timber stack mold risk level distribution map includes: After completing the projection of all future time steps, for each grid node Calculate the comprehensive risk level assessment value using the following formula: in, This serves as the baseline value for the risk intensity of the last time slice of the aforementioned mold risk propagation tensor. To normalize the cumulative transmission risk intensity value, The overall risk level of transmission; With all grid nodes The upper quartile and median are used as the thresholds for high-risk and medium-risk classification, respectively. Seed-based extended clustering is used to form risk partitions for adjacent connected grid nodes of the same level. The vertical coordinate range of the three-dimensional spatial grid is divided into three levels: bottom, middle and top. The results of planar risk partitioning and vertical stratification are combined to generate a partitioned and stratified distribution map of the risk level of mold growth in timber stacks.

8. A timber storage mold risk prediction system, characterized in that, include: The field construction module is used to receive temperature and humidity time series data from each monitoring point of the timber stack, and to perform grid-based mapping of the temperature and humidity time series data according to the stack space coordinates to construct an environmental parameter distribution field covering the storage space. The risk propagation module is used to calculate the humidity conduction rate field and spore diffusion potential field across grids based on the humidity gradient and temperature gradient of adjacent grid nodes in the environmental parameter distribution field, and to perform spatial convolution operation on the humidity conduction rate field and the spore diffusion potential field to generate a mold risk propagation tensor. The source identification module is used to extract the gradient change rate of the mold risk propagation tensor within a sliding window along the time axis, identify spatial cluster regions in the risk propagation tensor that exhibit a cumulative enhancement pattern, and mark the spatial cluster regions as a set of high-risk propagation sources. The zoning simulation module is used to perform spatiotemporal coupling simulation of the future evolution trend of the set of high-risk transmission sources and the distribution field of environmental parameters, and generate a zoning and stratified distribution map of the risk level of mold growth in timber stacks.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the timber storage mold risk prediction method according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the steps of the timber storage mold risk prediction method as described in any one of claims 1 to 7.

Images