Method and system for rapidly predicting underground soil temperature spatial and temporal distribution

By combining Fourier decomposition and least squares method, meteorological data is used to predict underground soil temperature, which solves the problems of high computational complexity and inaccurate boundary condition assumptions of traditional models, and achieves fast and accurate underground temperature prediction. It is suitable for soil temperature field monitoring and control in agriculture, construction, energy and other fields.

CN120633235APending Publication Date: 2025-09-12SHANDONG HUAYI GREEN ECOLOGICAL DEV CO LTD
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202510939229.X
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2025-07-08
Publication Date
2025-09-12

AI Technical Summary

Technical Problem

Traditional finite-difference numerical and analytical models have high computational complexity and are time-consuming when predicting the spatiotemporal distribution of underground soil temperature, making it difficult to meet the needs of high-frequency monitoring and rapid decision-making. In addition, existing methods lead to distorted predictions when boundary condition assumptions are inaccurate, affecting the accuracy and efficiency of applications such as ground-source heat pump site selection, deep foundation pit insulation, and smart irrigation.

Method used

A method combining Fourier decomposition and least squares method is adopted, and ground air temperature and solar radiation data are used as boundary conditions. Multiple decomposed cosine functions are obtained through Fourier decomposition. Combined with the surface heat flow boundary balance equation and underground temperature analytical heat transfer calculation, rapid prediction of underground temperature spatiotemporal distribution is achieved.

Benefits of technology

Without relying on surface temperature assumptions, accurate underground temperature prediction is achieved using meteorological data, which reduces computational complexity, improves prediction efficiency, and provides the flexibility to adjust computational accuracy and speed. It is suitable for precise calculations under complex boundary conditions.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN120633235A_ABST
    Figure CN120633235A_ABST
Patent Text Reader

Abstract

The invention belongs to the technical field of temperature spatio-temporal distribution, and discloses a method for quickly predicting underground soil temperature spatio-temporal distribution, which comprises the following steps of: acquiring annual hourly air dry-bulb temperature, air relative humidity, total solar radiation intensity irradiated to a horizontal plane, wind speed and earth surface soil humidity of a calculation site; the heat conductivity coefficient, the specific heat capacity, the density and the annual air temperature data decomposition limiting amplitude threshold value of underground soil are calculated; setting storage variables for storing the surface temperature and the temperature at the underground depth z of the soil; calculating annual hourly earth surface convection heat exchange amount, earth surface absorbed solar radiation amount, earth surface and sky long wave radiation amount and earth surface moisture evaporation heat dissipation amount; calculating the annual hourly comprehensive heat exchange coefficient; the annual hourly outdoor air dry-bulb temperature is decomposed into trigonometric function signals according to Fourier. According to the invention, data of outdoor air dry-bulb temperature and horizontal plane solar radiation intensity are obtained from meteorological data, and the underground temperature of any underground position and any time point can be predicted and calculated.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The present invention belongs to, but is not limited to, the technical field of temperature spatiotemporal distribution, and in particular relates to a method and system for rapidly predicting the spatiotemporal distribution of underground soil temperature. Background Art

[0002] Shallow subsurface soil temperature is driven periodically by upper boundary conditions such as solar radiation, air temperature, and rainfall, as well as by the coupling of internal mechanisms such as soil thermal conductivity, phase change, and water migration. This makes it highly spatiotemporally dynamic in agricultural, ecological, and engineering scenarios. While traditional finite-difference numerical models can characterize unsteady heat conduction processes moment by moment, they require fine-grid vertical soil measurements from 0–30 m and require iterative solutions over multiple years. Each increase in resolution or computational time increases the likelihood of convergence and exponentially increases the computational effort, making them difficult to meet the demands of high-frequency monitoring or rapid decision-making.

[0003] To reduce computational complexity, analytical models commonly used in recent years simplify surface temperature to a single-frequency cosine signal, using its mean and annual amplitude as input parameters. This directly provides a closed-form solution for the depth-time distribution. However, actual surface temperature often exhibits multimodal, irregular, or even abrupt segmented characteristics, and in most regions, a representative annual mean and single amplitude are lacking. Replacing the true boundary with a single-frequency cosine assumption inevitably introduces systematic biases, resulting in distorted predictions of surface soil temperature (0–5 m) and inaccurate estimates of thermal inertia in shallow soil layers (>10 m).

[0004] As a result, designers are unable to obtain accurate and efficient underground temperature profiles in industrial applications such as ground-source heat pump site selection, deep foundation pit insulation, smart irrigation, and farmland microclimate control. Numerical models are accurate but time-consuming, making them difficult to integrate into real-time control. Analytical models are fast but lack adaptability, leading to mismatches in equipment capacity selection, energy efficiency assessment, and safety margin design, further increasing system maintenance costs and energy consumption risks. This contradiction is becoming a core technical bottleneck facing industries such as green building, precision agriculture, and geological disaster early warning. Summary of the Invention

[0005] In response to the problems encountered by existing numerical heat transfer models and analytical heat transfer models when calculating the spatiotemporal distribution of underground soil temperature, the present invention aims to propose a new calculation method. This method is based on the existing analytical heat transfer model framework, but uses more easily obtained meteorological data such as ground air temperature and solar radiation as boundary conditions. It comprehensively utilizes Fourier decomposition and least squares method to perform Fourier decomposition on the hourly changing outdoor air temperature to obtain multiple decomposed cosine functions. The surface heat flow boundary balance equation and least squares method are used to obtain the surface cosine function expression. Then, the analytical heat transfer calculation formula for underground temperature is superimposed with the Fourier signal to restore the prediction results of the spatiotemporal distribution of underground temperature under complex meteorological conditions.

[0006] The present invention is achieved by providing a method for rapidly predicting the spatiotemporal distribution of underground soil temperature, characterized in that the method comprises:

[0007] S1: Before the calculation begins, it is necessary to prepare the hourly air dry bulb temperature T of the calculation location throughout the year. air , relative humidity RH, total solar radiation intensity G on the horizontal surface, wind speed u, wetness of surface soil f; and thermal conductivity k of underground soil s Specific heat capacity c p , density ρ, annual air temperature data decomposition limit amplitude threshold A limit ;

[0008] S2: Set initial values ​​for the soil surface temperature and the temperature at depth z below the soil surface, and set storage variables for storing the surface temperature and the temperature at depth z below the soil surface at each moment. The initial value of the soil surface temperature only needs to assume an average value for the entire year.

[0009] S3: Calculate the surface convection heat transfer q hourly throughout the year c , the amount of solar radiation absorbed by the surface q sun , the long-wave radiation q from the earth's surface and the sky sky , ground water evaporation heat dissipation q evap , and the total net heat gain q of the ground at the current moment is calculated as follows:

[0010] q c =h c (T air -T s ) (1)

[0011] h c =5.7+3.8×u (2)

[0012] q sun =αG sun (3)

[0013]

[0014] T sky =T air (∈ sky ) 1 / 4 (5)

[0015] ∈ sky =∈ sky,clear +(1-∈ sky,clear )C (6)

[0016] q evap =0.0168fh c (p ws -p w ) (7)

[0017] q=q c +q sun +q sky +q evap (8);

[0018] Where: h c Represents the convection heat transfer coefficient of the ground, W / m 2 .K;T s represents the surface temperature, ℃; α represents the absorption rate of solar radiation by the ground; G sun Represents the total solar radiation received by the horizontal surface, W / m 2 ; ε represents the long-wave radiation emissivity of the earth's surface; σ represents the Stefan-Boltzmann constant; Tsurface represents the surface temperature, ℃; Tsky represents the sky temperature, ℃; ε sky represents the long-wave emissivity of the sky; ε sky,clear represents the long-wave emissivity of an ideal clear sky; C represents the cloud thickness parameter of the sky; p ws Represents the saturated partial pressure of water vapor in the air near the ground, pa; p w Represents the water vapor partial pressure in the air near the ground, pa;

[0019] S4: Based on the calculation of S3, the net heat gain q of the ground every hour throughout the year, as well as the initial assumed average value of the outdoor air temperature and the surface temperature, calculate the comprehensive heat transfer coefficient h every hour throughout the year q , the calculation formula (9) is as follows:

[0020] h q =q / (T air -T s ) (9)

[0021] S5: The outdoor air dry bulb temperature T air According to Fourier decomposition into several trigonometric function signals, the calculation formula is as follows:

[0022]

[0023] Where: T air,m,i represents the average value of the ith signal after the hourly data of air temperature is decomposed throughout the year, A air,i represents the amplitude of the ith signal after the hourly data of air temperature throughout the year is decomposed, and satisfies A air,i limit ; t0 is the phase difference of signal i; N is the number of decomposed signals; t represents time; t y Represents the total length of the calculation cycle;

[0024] S6: Calculate the hourly surface temperature T throughout the year s According to Fourier decomposition into several trigonometric function signals, the calculation formula is as follows:

[0025]

[0026] Where: T s,m,i represents the average value of the ith signal after the hourly data of the surface temperature throughout the year are decomposed. s,i represents the amplitude of the ith signal after decomposing the hourly data of air temperature throughout the year, and t0 is the phase difference of the signal; but this T s,m,i Value and A s,i The values ​​are initially two unknowns, waiting to be solved in the next step. In the calculation, it is assumed that the phase difference t0 of the signal decomposed by the air temperature is consistent with the phase difference t0 of the signal decomposed by the surface temperature.

[0027] S7: Using the boundary condition that the heat flow from the air into the ground at the surface interface is equal to the heat flow from the ground into the ground, as shown in formula (12), formula (12) is further organized into an unknown number to be solved T air,m,i Value and A air,i The numerical form of equation (13):

[0028]

[0029] Where: w is the frequency of data change throughout the year, w = 2π / ty; λ is the thermal conductivity of the soil, W / mK; a is the thermal diffusivity of the soil, m 2 / s; x is an intermediate calculation variable, x=(w / (2a)) 0.5 ; The subscript i represents the number of signals decomposed from the annual temperature data, ranging from 1 to N;

[0030] S8: Assign the time t in formula (13) a value of 1-8760, thereby expanding one equation (13) into 8760 equations, and using the least squares method to solve the unknown number T​s,m,i Value and A s,i Numeric value;

[0031] S9: The calculation formula for the soil temperature T at the underground position z at time t is finally obtained:

[0032]

[0033] In combination with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present invention are as follows:

[0034] First, the method proposed in the present invention for predicting the temporal and spatial distribution of underground temperature can predict the underground temperature at any location and at any time, provided that the surface temperature is unknown and only requires the outdoor air dry-bulb temperature and horizontal solar radiation intensity data, which are easily obtained from meteorological data.

[0035] The calculation method proposed in the present invention does not need to assume that the outdoor air temperature satisfies the ideal trigonometric function variation law. The method can accurately calculate under the action of boundary conditions with arbitrarily complex changes and has a wide range of applications.

[0036] The calculation method provided by the present invention also has the advantage of high calculation flexibility. During the calculation, the threshold of the Fourier decomposition signal amplitude can be customized. The lower the threshold, the more signals are decomposed, which makes the calculation result more accurate. Conversely, the higher the threshold, the fewer signals are decomposed, which makes the calculation faster. Therefore, the present invention provides a parameter that can easily adjust the calculation accuracy and calculation speed.

[0037] Although the present invention provides a parameter that can easily adjust the calculation accuracy and calculation speed, the calculation method of the present invention can still maintain a high-accuracy calculation result under the premise of fewer Fourier decomposition signals.

[0038] Second, the technical solution of the present invention fills the technical gap in the industry at home and abroad:

[0039] There are currently several methods for predicting the spatiotemporal distribution of underground temperature, but most commonly used are numerical solutions using finite differences or formulas to calculate underground temperature trends when the surface temperature is known to vary according to a trigonometric function. This method requires simplifying the complex surface temperature variations into a trigonometric function, which is obviously inconsistent with the actual situation. Furthermore, while air temperature data is generally readily available, surface temperature data is often unavailable. This makes existing methods difficult to implement or requires numerous assumptions (for example, assuming that the surface temperature follows an ideal cosine function). On the other hand, while the finite difference method can accurately calculate complex outdoor boundary conditions, it is slow and inefficient. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 This is a flow chart of a method for rapidly predicting the temporal and spatial distribution of underground soil temperature provided by an embodiment of the present invention;

[0041] Figure 2 The embodiment of the present invention provides a year-round hourly outdoor temperature and humidity and soil temperature and humidity at different depths underground;

[0042] Figure 3 1 is a diagram showing the verification results of the underground analytical heat transfer model provided by an embodiment of the present invention. (a) The analytical model calculates the soil temperature at 0.3 meters underground; (b) The analytical model calculates the soil temperature at 0.5 meters underground; (c) The analytical model calculates the soil temperature at 1.0 meters underground; and (d) The analytical model calculates the soil temperature at 2.0 meters underground. DETAILED DESCRIPTION

[0043] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0044] This method addresses the challenges of traditional soil temperature prediction models, such as high computational complexity, complex parameter coupling, and poor real-time performance. By separating the surface heat exchange process from the underground conduction process and introducing Fourier signal decomposition, it significantly reduces computational complexity and improves prediction accuracy. First, at the multi-source meteorological data input, an amplitude threshold is applied to the hourly air temperature series throughout the year to filter out high-frequency noise components. Finite sine terms are then used to accurately reconstruct the temperature variation pattern, replacing the time-domain difference calculations used in traditional numerical models and reducing redundant iterations.

[0045] Secondly, the surface energy balance was carefully analyzed, with convective heat transfer, longwave radiation, shortwave absorption, and evaporative heat dissipation independently modeled as physical processes. Specifically, the coupling coefficient fitting method for evaporative heat dissipation was improved to account for the dynamic changes in atmospheric humidity and surface moisture content, thereby maintaining an energy balance error below 5% even under high humidity conditions.

[0046] Thirdly, based on the matching of the heat diffusion equation with the boundary conditions, frequency-domain analysis techniques were introduced to transform the continuity of heat flux at the surface interface into a finite set of algebraic constraints on the amplitude and phase of sinusoidal signals. The least-squares method was used to solve the boundary equations for 8,760 hours, resulting in a set of globally optimal Fourier coefficients that accurately describe the time-varying characteristics of the surface temperature and automatically adapt to different soil thermophysical parameters.

[0047] Subsequently, a hybrid algorithm combining analytical and numerical solutions was employed for the underground conduction process. The relationship between the depth attenuation coefficient and phase delay was theoretically derived, and the reconstructed surface temperature signal was then fed into the analytical heat conduction formula term by term. Given the rapid decay of high-frequency terms, only low-order terms were retained for the conduction calculation, ensuring the continuity of the deep temperature field while significantly reducing calculation time.

[0048] At the industrial application level, this method can be embedded in smart agricultural irrigation systems to achieve real-time early warning of soil freeze-thaw, heat stress, and root zone temperature; it can also be used in ground-source heat pump design to accurately predict underground thermal energy storage and recovery efficiency; in large-scale tunnel, subway construction and pipeline laying, it can provide dynamic assessment of seasonal ground temperature fields during construction, providing data support for material selection and the formulation of anti-freeze measures.

[0049] This invention takes physical processes as its core and frequency domain decomposition as its technical means, breaking the computational bottleneck of traditional time-domain difference models. It not only takes into account both prediction accuracy and real-time performance, but also has good parameter interpretability and engineering scalability. It is suitable for soil temperature field monitoring and control in multiple fields such as agriculture, construction, and energy.

[0050] like Figure 1 As shown, an embodiment of the present invention provides a method for quickly predicting the spatiotemporal distribution of underground soil temperature, the method comprising:

[0051] S1: Before the calculation begins, it is necessary to prepare the hourly air dry bulb temperature T of the calculation location throughout the year. air , relative humidity RH, total solar radiation intensity G on the horizontal surface, wind speed u, wetness of surface soil f; and thermal conductivity k of underground soil s Specific heat capacity c p , density ρ, annual air temperature data decomposition limit amplitude threshold A limit ;

[0052] S2: Set initial values ​​for the soil surface temperature and the temperature at depth z below the soil surface, and set storage variables for storing the surface temperature and the temperature at depth z below the soil surface at each moment. The initial value of the soil surface temperature only needs to assume an average value for the entire year.

[0053] S3: Calculate the surface convection heat transfer q hourly throughout the year c , the amount of solar radiation absorbed by the surface q sun , the long-wave radiation q from the earth's surface and the sky sky , ground water evaporation heat dissipation q evap , and the total net heat gain q of the ground at the current moment is calculated as follows:

[0054] q c =h c (T air -T s ) (1)

[0055] h c =5.7+3.8×u (2)

[0056] q sun =αG sun (3)

[0057]

[0058] T sky =T air (∈ sky ) 1 / 4 (5)

[0059] ∈ sky =∈ sky,clear +(1-∈ sky,clear )C (6)

[0060] q evap =0.0168fh c (p ws -p w ) (7)

[0061] q=q c +q sun +q sky +q evap (8);

[0062] Where: h c Represents the convection heat transfer coefficient of the ground, W / m 2 .K;T s represents the surface temperature, ℃; α represents the absorption rate of solar radiation by the ground; G sun Represents the total solar radiation received by the horizontal surface, W / m 2; ε represents the long-wave radiation emissivity of the earth's surface; σ represents the Stefan-Boltzmann constant; Tsurface represents the surface temperature, ℃; Tsky represents the sky temperature, ℃; ε sky represents the long-wave emissivity of the sky; ε sky,clear represents the long-wave emissivity of an ideal clear sky; C represents the cloud thickness parameter of the sky; p ws Represents the saturated partial pressure of water vapor in the air near the ground, pa; p w Represents the water vapor partial pressure in the air near the ground, pa;

[0063] S4: Based on the calculation of S3, the net heat gain q of the ground every hour throughout the year, as well as the initial assumed average value of the outdoor air temperature and the surface temperature, calculate the comprehensive heat transfer coefficient h every hour throughout the year q , the calculation formula (9) is as follows:

[0064] h q =q / (T air -T s ) (9)

[0065] S5: The outdoor air dry bulb temperature T air According to Fourier decomposition into several trigonometric function signals, the calculation formula is as follows:

[0066]

[0067] Where: T air,m,i represents the average value of the ith signal after the hourly data of air temperature is decomposed throughout the year, A air,i represents the amplitude of the ith signal after the hourly data of air temperature throughout the year is decomposed, and satisfies A air,i limit ; t0 is the phase difference of signal i; N is the number of decomposed signals; t represents time; t y Represents the total length of the calculation cycle;

[0068] S6: Calculate the hourly surface temperature T throughout the year s According to Fourier decomposition into several trigonometric function signals, the calculation formula is as follows:

[0069]

[0070] Where: T s,m,i represents the average value of the ith signal after the hourly data of the surface temperature throughout the year are decomposed. s,i represents the amplitude of the ith signal after decomposing the hourly data of air temperature throughout the year, and t0 is the phase difference of the signal; but this T s,m,i Value and A s,i ​The values ​​are initially two unknowns, waiting to be solved in the next step. In the calculation, it is assumed that the phase difference t0 of the signal decomposed by the air temperature is consistent with the phase difference t0 of the signal decomposed by the surface temperature.

[0071] S7: Using the boundary condition that the heat flow from the air into the ground at the surface interface is equal to the heat flow from the ground into the ground, as shown in formula (12), formula (12) is further organized into an unknown number to be solved T air,m,i Value and A air,i The numerical form of equation (13):

[0072]

[0073] Where: w is the frequency of data change throughout the year, w = 2π / ty; λ is the thermal conductivity of the soil, W / mK; a is the thermal diffusivity of the soil, m 2 / s; x is an intermediate calculation variable, x=(w / (2a)) 0.5 ; The subscript i represents the number of signals decomposed from the annual temperature data, ranging from 1 to N;

[0074] S8: Assign the time t in formula (13) a value of 1-8760, thereby expanding one equation (13) into 8760 equations, and using the least squares method to solve the unknown number T s,m,i Value and A s,i Numeric value;

[0075] S9: The calculation formula for the soil temperature T at the underground position z at time t is finally obtained:

[0076]

[0077] The specific application fields or related products of the present invention are:

[0078] 1. This invention can be used to accurately predict temperature changes throughout the year within the 0-5 meter range underground, thereby better assessing the growing environment of crops or implementing some intervention measures, and providing guidance for the growth of crops or some cash crops;

[0079] 2. The present invention can be used to predict the temperature of shallow underground soil, which helps to assess the potential of local shallow geothermal energy utilization and can be used in the application of geothermal energy buried pipe systems, geothermal energy geothermal water extraction systems, etc.

[0080] The relevant evidence for the technical effects achieved by the embodiments of the present invention is:

[0081] like Figure 2Figure 1 shows a year-round hourly outdoor temperature and humidity graph, along with soil temperature and humidity at different depths, according to an embodiment of the present invention. Hourly outdoor air temperature and humidity were measured in Wuhan throughout 2024. Temperature sensors were also used to measure temperatures at depths of 0.3, 0.5, 1, and 2 meters underground.

[0082] The calculation method proposed in the present invention was used to calculate the temperature values ​​of 0.3 meters, 0.5 meters, 1.0 meters, and 2.0 meters underground in Wuhan in 2024. These calculated values ​​were compared with the measured values ​​one by one. The results are as shown in the embodiment of the present invention. Figure 2 As shown in Table 1, the relevant error analysis shows that the average error of the calculation does not exceed 5%, which proves that the calculation method of the present invention has high calculation accuracy.

[0083] Table 1 Error analysis of analytical model

[0084]

[0085] It should be noted that the embodiments of the present invention can be implemented by hardware, software, or a combination of software and hardware. The hardware portion can be implemented using dedicated logic; the software portion can be stored in a memory and executed by an appropriate instruction execution system, such as a microprocessor or dedicated design hardware. Those skilled in the art will appreciate that the above-mentioned devices and methods can be implemented using computer-executable instructions and / or contained in processor control code, for example, such as a carrier medium such as a disk, CD or DVD-ROM, a programmable memory such as a read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuits such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., can also be implemented by software executed by various types of processors, or can be implemented by a combination of the above-mentioned hardware circuits and software, such as firmware.

[0086] The above description is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with this technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered by the scope of protection of the present invention.

Claims

1. A method for rapidly predicting the temporal and spatial distribution of underground soil temperature, characterized in that: The following steps are involved: S1 collects the hourly external air dry-bulb temperature Tair, air relative humidity RH, total solar radiation intensity G, wind speed u, surface soil wetness f, and underground soil thermal conductivity ks, specific heat capacity cp, and density ρ at the calculation location throughout the year, and sets the air temperature data decomposition amplitude threshold Alimit; S2, set initial values ​​for the surface temperature Ts and the soil temperature T at any depth z underground and establish corresponding hourly storage variables; S3, based on the external meteorological parameters obtained in step S1, calculate the surface convective heat transfer qc, the surface absorbed solar radiation qsun, the surface and sky longwave radiation qsky, and the surface water evaporation heat dissipation qevap, and obtain the surface net heat gain q at that moment; S4, calculate the comprehensive heat transfer coefficient hq using the net heat gain q and the external air temperature Tair; S5, perform Fourier series decomposition on the hourly outside air temperature Tair throughout the year to obtain N trigonometric function components, and the amplitude of each component does not exceed Alimit; S6, performing Fourier series decomposition on the hourly surface temperature Ts throughout the year to obtain N trigonometric function components corresponding to step S5, and setting the phase difference between the two components to be consistent; S7, establish a system of equations based on the boundary condition that the heat flow between the air and the surface is equal to the heat flow conducted from the surface to the subsurface, and take the mean and amplitude of the decomposed components of the surface temperature as unknown quantities; S8, expanding the system of equations at 8,760 moments throughout the year and solving the unknowns in step S7 using the least squares method; S9: Substitute the unknown quantity obtained into the analytical formula for underground heat conduction to obtain the soil temperature T at any depth z and any time t and output it.

2. The method according to claim 1, characterized in that In step S3, the surface convective heat transfer qc is calculated by combining the external air temperature Tair, the surface temperature Ts and the wind speed u with a correction coefficient, and the surface absorbed solar radiation qsun is obtained by combining the total solar radiation intensity G and the surface albedo.

3. The method according to claim 1, characterized in that The Fourier decomposition in step S5 uses the fundamental frequency as the annual frequency change = π divided by the number of seconds in a year, and the number of components N is automatically determined under the condition that all component amplitudes are less than Alimit.

4. The method according to claim 1, wherein In step S7, the underground thermal diffusion coefficient a is calculated based on the underground soil thermal conductivity coefficient ks, specific heat capacity cp, and density ρ.

5. The method according to claim 1, wherein In step S8, the least squares method is solved by using a direct solution algorithm based on a pseudo-inverse matrix to increase the calculation speed.

6. A fast computing system for implementing the method according to any one of claims 1 to 5, characterized in that: The method comprises a processor, a memory and a communication interface, wherein the processor is configured to execute all steps of the method according to any one of claims 1 to 5 and store output results.

7. The fast calculation system according to claim 6, characterized in that: The processor is further configured to automatically generate a soil temperature profile file with a depth of zero to thirty meters and a time resolution of one hour after the user inputs a target depth range.

8. An underground soil temperature prediction device, characterized in that: include: A data acquisition unit is used to obtain external air meteorological parameters and underground soil physical parameters; a data processing unit, communicatively connected to the data acquisition unit, configured to execute the method according to any one of claims 1 to 5 and generate a soil temperature prediction result; The display unit is connected to the data processing unit and is used to display the prediction curve and numerical values ​​in real time.

9. A computer-readable storage medium having instructions stored thereon, which, when executed on a computer, causes the computer to execute all the steps of the method according to any one of claims 1 to 5.

10. A computer program product, characterized in that The invention comprises a plurality of executable instructions, which implement the method according to any one of claims 1 to 5 when the instructions are loaded and executed in a computing device.