A road long-period icing prediction method based on meteorological parameter inversion

Abstract

The application discloses a kind of road long-period icing prediction methods based on meteorological parameter inversion, comprising: collecting target road section meteorological data, to temperature radiation temperature difference correction term, air cooling correction term and evaporation cooling correction term obtain road surface temperature inversion value;By precipitation data initialization virtual water film thickness and update hour by hour according to evaporation rate, zero when cutting evaporation cooling correction term;By past N meteorological data calculation history road surface temperature inversion value, according to the similarity weighted generation seasonal reference curve and map to future prediction period with each historical day weather;With current road surface temperature inversion value as starting point, fusion forecast road surface temperature inversion value and reference curve obtain final predicted road surface temperature;When it does not exceed 0 ℃ and there is moisture source on road surface, determine icing risk.The application does not need to be laid out road surface contact type sensor, also does not need long-term measured data offline calibration, can realize long-period road surface temperature prediction in observation blind area section, expand icing early warning coverage range.

Description

Technical Field

[0001] This invention relates to the field of traffic meteorological monitoring and intelligent traffic safety technology, specifically to a method for predicting long-term road icing based on meteorological parameter inversion. Background Technology

[0002] Road icing early warning is a core component of ensuring winter traffic safety. Road surface temperature is a key physical quantity determining whether road icing will occur, and accurately predicting the spatiotemporal distribution of road surface temperature is crucial for traffic management departments to formulate anti-icing and snow removal strategies and issue traffic warnings. However, existing road icing prediction technologies face multiple challenges in large-scale engineering applications, urgently requiring breakthroughs in traditional technological approaches.

[0003] Currently, road surface temperature prediction mainly relies on two existing technical solutions. The first is the hardware monitoring solution, which directly acquires road surface temperature data by burying contact road surface sensors; the second is the sensorless prediction solution, which uses standard atmospheric monitoring data in road sections lacking sensors and indirectly estimates the road surface temperature using simple statistical models or empirical formulas, in order to expand the coverage of icing warnings.

[0004] However, the aforementioned existing technologies all have significant limitations in engineering applications. Regarding hardware dependence and construction costs, the long distances along highways mean that the cost of large-scale deployment and maintenance of hardware sensors is extremely high, leaving most ordinary road sections in meteorological observation blind spots. This hardware dependence fundamentally restricts the coverage of icing warnings. In terms of physical logic completeness, existing sensorless prediction schemes often use simple statistical models or empirical formulas, ignoring radiation anomalies under atypical meteorological conditions and the significant impact of latent heat of evaporation after precipitation on road surface energy balance. Specifically, radiation behavior under high humidity sunny days and dry cloudy days deviates significantly from conventional assumptions. Existing models typically treat the influence of cloud cover as a proxy variable for humidity, leading to systematic biases under these atypical meteorological scenarios. Furthermore, due to the lack of real-time evolution assessment of road surface wetness, the models are prone to generating false evaporative cooling prediction errors under dry road surface conditions, resulting in false alarms or missed warnings for icing. Regarding parameter calibration and adaptive capabilities, the radiation, convection, and evaporation coefficients in the model are significantly affected by road material and micro-topography. To achieve high-precision prediction, existing technologies must rely on long-term historical measured data for offline calibration. For newly built roads or road sections lacking sensors, the model lacks online self-evolution capabilities, causing the prediction accuracy to rapidly decline with road surface aging and seasonal changes. Furthermore, the parameter calibration process itself relies on road surface temperature sensor data, creating a logical contradiction between the need for measured data for calibration and sensorless application scenarios. Summary of the Invention

[0005] This invention proposes a long-term road icing prediction method based on meteorological parameter inversion, which solves the problem that existing technologies cannot achieve long-term prediction of road surface temperature without relying on the target road section's own sensors.

[0006] To address the aforementioned technical problems, this invention provides a method for predicting long-term road icing based on meteorological parameter inversion, comprising the following steps: Step S1: Collect meteorological observation data of the target road section at the current moment, apply radiation temperature difference correction term, wind cooling correction term and evaporative cooling correction term to the air temperature based on the meteorological observation data to obtain the current road surface temperature inversion value; initialize the virtual water film thickness based on the precipitation data in the meteorological observation data, and update the virtual water film thickness hourly based on the evaporation rate; cut off the evaporative cooling correction term when the virtual water film thickness drops to zero. Step S2: Collect hourly meteorological observation data for the past N days, calculate the historical road surface temperature inversion value for each historical day, obtain meteorological forecast data, and generate a seasonal baseline curve by weighting the historical road surface temperature inversion value based on the similarity between the forecast weather characteristics and the weather characteristics of each historical day. Map the seasonal baseline curve to the future prediction period according to the daily cycle to obtain the prediction period baseline curve. Step S3: Calculate the predicted road surface temperature inversion value for each future time based on the meteorological forecast data. Using the current road surface temperature inversion value as the prediction starting point, fuse the predicted road surface temperature inversion value with the baseline curve of the prediction period to obtain the final predicted road surface temperature. Step S4: When the final predicted road surface temperature does not exceed 0℃ and there is a source of moisture on the road surface, it is determined that the target road section has a risk of icing.

[0007] Preferably, the expressions for the current road surface temperature inversion value, the air-cooling correction term, and the evaporative cooling correction term in step S1 are as follows: ; ; ; In the formula, This is the current road surface temperature inversion value; Temperature from meteorological observation data; This is a correction term for radiation temperature difference; For air-cooled correction; This is a correction term for evaporative cooling. This is the correction factor for air cooling; The coefficient of performance is the evaporative cooling rate. Wind speed in meteorological observation data; Relative humidity in meteorological observation data; This is the road surface wettability coefficient.

[0008] Preferably, the radiation temperature difference correction term is calculated separately for the daytime and nighttime periods, including the following steps: The sunrise and sunset times are calculated based on the latitude and longitude information of the target road segment. If the current time is between the sunrise and sunset times, it is determined to be the daytime stage; otherwise, it is determined to be the nighttime stage. Calculate the solar altitude angle based on the latitude and longitude information of the target road segment and the current time. And determine the solar elevation angle modulation factor. ,when hour ,when hour ; During the daytime phase, the expression for the radiation temperature difference correction term is: ; During the nighttime phase, the expression for the radiation temperature difference correction term is: ; In the formula, The solar altitude angle; Solar elevation angle modulation factor; This is the radiation correction factor; The diurnal radiation coefficient; This is a dynamic compensation item for atypical daytime weather conditions; Operator for atmospheric window transparency correction.

[0009] Preferably, the daytime atypical meteorological dynamic compensation term and the atmospheric window transparency correction operator are valued in the following manner: During the daytime phase, the weather phenomenon is coded as sunny and When it is greater than 70%, set Positive compensation parameters for daytime conditions under sunny and high humidity conditions ,and ; During the daytime phase, the weather phenomenon is coded as cloudy and When less than 30%, set Negative compensation parameters for daytime conditions under cloudy and low humidity conditions ,and ; During the nighttime phase, the weather phenomenon is coded as clear and When less than 30%, set The negative compensation parameter for nighttime under clear and dry nighttime conditions. ,and ; During the nighttime phase and When the weather phenomenon is greater than 85% or coded as light fog, set Nighttime positive compensation parameters under conditions of high humidity or light fog. ,and ; If none of the above four conditions are met, and All are set to zero.

[0010] Preferably, the virtual water film thickness and the road surface wetting coefficient are determined in the following manner: When precipitation data in the meteorological observation data is identified as precipitation, the virtual water film thickness... Take initial value The road surface wettability coefficient Take 1; After the precipitation ends, the Magnus approximation formula is used to calculate the first... The saturated vapor pressure corresponding to the temperature at that moment: ; Calculate the first according to Dalton's law of evaporation. Evaporation rate of road surface water film at any given time: ; Update the virtual water film thickness according to the integral decay equation: ; The mapping relationship between the road surface wettability coefficient and the virtual water film thickness is as follows: ; In the formula, This is the elapsed time since the end of the precipitation. The time corresponding to the current moment value; The virtual water film thickness; This is the constant for the maximum water film thickness on the road surface; For the first The saturated vapor pressure corresponding to the temperature at a given moment; For the first Temperature at any given moment; For the first The rate of evaporation of the water film on the road surface at any given time; The coefficient of evaporation due to natural convection; The mass transfer coefficient is the wind speed. For the first Wind speed at any given moment; For the first Relative humidity at any given time.

[0011] Preferably, the radiation correction coefficient, the diurnal radiation coefficient, the wind-cooling correction coefficient, the evaporative cooling coefficient, the daytime atypical meteorological dynamic compensation term, and the atmospheric window transparency correction operator are subjected to unsupervised adaptive calibration based on physical boundary constraints of natural meteorological events, including the following steps: Identify natural calibration events from the meteorological observation data, where the road surface temperature is constrained by physical laws to a known value; Based on the aforementioned natural calibration events, pseudo-true values ​​of road surface temperature are extracted; The residual is calculated by comparing the current road surface temperature inversion value with the road surface temperature pseudo value. The parameters are updated with the goal of minimizing the residual using the projection gradient descent method based on a sliding time window. Each calibration only updates the parameters corresponding to the natural calibration event. The sliding time window consists of the most recent preset number of the same type of natural calibration events, and the projected gradient descent method projects the candidate parameter values ​​to the interval defined by the lower and upper bounds of the updated parameters after each iteration.

[0012] Preferably, the natural calibration events include the following five types: Evaporative cooling parameter calibration event: When the virtual water film thickness is greater than zero, the road surface wetting coefficient is not less than the preset wetting threshold, and the wind speed is greater than 2m / s after the precipitation ends, the pseudo-true value of the road surface temperature is taken as the wet-bulb temperature calculated according to the Stull empirical formula, and the evaporative cooling coefficient is updated accordingly. Nighttime radiation parameter calibration event: When the nighttime phase is in progress, the weather phenomenon is coded as sunny, the relative humidity is less than the preset humidity threshold, there is no precipitation, the virtual water film thickness is zero and the wind speed is less than 2m / s, and the meteorological observation record indicates the occurrence of dew or frost, the pseudo-true value of the road surface temperature is taken as the dew point temperature when the air temperature is higher than 0℃ and as the frost point temperature when the air temperature is not higher than 0℃, and the diurnal radiation coefficient is updated accordingly. Wind-cooled parameter calibration events: no precipitation, the virtual water film thickness is zero, the cloud cover in the meteorological observation data remains at the same cloud cover level within a preset time period, the cloud cover level is divided according to a preset cloud cover threshold and the wind speed appears in a stage not less than a preset range within a preset interval, the difference between the current road surface temperature inversion value at the moment of minimum wind speed in the stage and the reference road surface temperature at the same moment is used as the bias correction amount, the sum of the reference road surface temperature at each moment in the stage and the bias correction amount is used as the pseudo value of the road surface temperature at the corresponding moment, the difference between the pseudo value of the road surface temperature and the air temperature is used to perform a linear regression on the wind speed, and the negative value of the regression slope is used to update the wind-cooled correction coefficient accordingly, the reference road surface temperature is the median of the measured road surface temperature of the reference road section under the same solar altitude angle and similar meteorological conditions, and the climate zone type, road surface material and micro-topographic features of the reference road section are the same as those of the target road section; Daytime radiation parameter calibration event: When the daytime phase is in progress, the weather phenomenon is coded as sunny, there is no precipitation, the virtual water film thickness is zero, the wind speed is less than 2m / s and the temperature is higher than 0℃, and the road surface snow or ice begins to melt, the pseudo-true value of the road surface temperature is set to 0℃, and the radiation correction coefficient is updated accordingly. Extreme freezing constraint event: When snow is detected on the road surface and the air temperature rises from no higher than 0°C to higher than 0°C, the pseudo-true value of the road surface temperature is set to 0°C and continues until the snow melt ends. The corresponding daytime atypical meteorological dynamic compensation item or atmospheric window transparency correction operator triggered in the current period is updated, and the extreme freezing constraint event is executed before the daytime radiation parameter calibration event.

[0013] Preferably, the similarity in step S2 is cosine similarity, and the historical road surface temperature inversion value is weighted using meteorological isomorphism weighting, including the following steps: The future prediction period is divided into multiple consecutive future prediction days, and each future prediction day includes 24 hourly prediction times. Cloud cover is divided into four levels according to preset cloud cover thresholds, and wind speed is divided into four levels according to preset wind speed thresholds. For any given future forecast date, the number of hours each cloud cover level and each wind speed level occurs within 24 hours of that future forecast date is counted and divided by 24 to obtain the time proportion of each level within the future forecast date. The time proportions of the four cloud cover levels and the four wind speed levels are arranged sequentially to form an eight-dimensional forecast weather feature vector for that future forecast date. The same method is used to construct past... N The eight-dimensional measured weather feature vector of each historical day in Tianzhong; Calculate the first The eight-dimensional forecast weather feature vector for the predicted future date and the first The cosine similarity between the eight-dimensional measured weather feature vectors of the historical days is used as the first... The historical day for the first Meteorological isomorphic weights for each future forecast date ; The expression for the seasonal baseline curve is: ; In the formula, N The total number of historical days collected; The predicted date number; Historical date sequence; It is the hour number within a day, and its value is an integer from 0 to 23; For the first One future prediction day The seasonal baseline curve value at that time; For the first A historical day Historical road surface temperature inversion values ​​at that time.

[0014] Preferably, in step S3, the predicted pavement temperature inversion value and the baseline curve for the predicted time period are fused based on a dynamic fusion weight that decays with the forecast lead time. The expression is as follows: ; ; ; In the formula, Forecast lead time, which is the duration calculated from the start of the forecast; For forecast timeliness Dynamic fusion weights at the location; This is the time constant for predicting the decay characteristics; For meteorological stability operators; For forecast timeliness The final predicted road surface temperature at the location; For forecast timeliness The predicted pavement temperature inversion value at the location; To provide a baseline curve for the forecast period in terms of forecast lead time The value at; The above expression applies to the case where the forecast lead time is greater than zero; where, for Add 1 to the integer part of the quotient obtained by dividing by 24. To predict the starting hour and The remainder when the sum is divided by 24 The number of whole hours within the day corresponding to the predicted starting point, with a value ranging from 0 to 23; When the value is zero, the final predicted road surface temperature is taken as the inversion value of the current road surface temperature.

[0015] Preferably, the existence of a source of moisture on the road surface in step S4 means that the prediction time corresponding to the final predicted road surface temperature meets at least one of the following conditions: The meteorological forecast data indicates that precipitation will occur at the predicted time. The virtual water film thickness at the predicted time is greater than zero, and the virtual water film thickness is updated hourly according to the meteorological forecast data during the future prediction period. The dew point temperature calculated based on the air temperature and relative humidity at the predicted time in the meteorological forecast data is higher than the final predicted road surface temperature at the predicted time, and the relative humidity at the predicted time is greater than 90%.

[0016] The advantages of this invention include at least the following: 1. By applying radiation temperature difference correction, wind cooling correction and evaporative cooling correction to the air temperature, the road surface temperature calculation is upgraded from a simple statistical regression to a physical inversion process based on the road surface energy balance; it overcomes the defect of existing empirical formulas that ignore the radiation anomaly in high humidity sunny days and the deviation of radiation behavior in dry cloudy days, so that the model can still maintain physical consistency under atypical meteorological conditions. 2. By updating the virtual water film thickness hourly and mapping it to the wetting coefficient, the road surface moisture status is embedded into the evaporative cooling correction term in real time. Under dry road surface conditions, the wetting coefficient approaches zero to suppress non-existent evaporative cooling. After precipitation or under high humidity conditions, the impact of latent heat of evaporation on the road surface energy balance is taken into account, avoiding false cooling predictions under dry road surface conditions. 3. By utilizing physical boundary constraints such as wet-bulb temperature, dew point or frost point temperature, and snow melting temperature stagnation, the pseudo-true value of road surface temperature is extracted online. This allows the model parameters to be continuously corrected as the road surface ages and the seasons change, even when there are no road surface temperature sensors on the target road section. This overcomes the problem that existing technologies must rely on long-term measured data for offline calibration and quickly fail in road sections where sensors are missing. 4. By integrating the time-by-time inversion results based on physical mechanisms with the seasonal baseline curve based on historical data mapping, the accuracy of the physical response to short-term meteorological fluctuations is preserved, while the drift error of the pure physical model in long-term forecasts is suppressed by long-term statistical rhythm constraints. Attached Figure Description

[0017] Figure 1 This is a flowchart of a method according to an embodiment of the present invention; Figure 2 This is a comparison curve of the predicted road surface temperature and the measured temperature in an embodiment of the present invention. Detailed Implementation

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

[0019] This invention provides a method for predicting long-term road icing based on meteorological parameter inversion. Taking the prediction of a certain highway section over the next 7 days (168 hours) as an example, the technical solution of this invention will be described in detail.

[0020] like Figure 1 As shown, the method in this embodiment includes four main steps, from step S1 to step S4.

[0021] Step S1: Real-time data acquisition and road surface temperature inversion.

[0022] The purpose of this step is to calculate the road surface temperature inversion value using conventional meteorological observation data, based on a first-order linearized approximation model of the surface energy balance principle, in the absence of a road surface temperature sensor.

[0023] First, meteorological observation data for the target road segment at the current time is collected. In this embodiment, the required meteorological observation data includes the following six items: Measured air temperature, in degrees Celsius (°C), is recorded as follows: Wind speed, in m / s, symbol . Relative humidity, expressed as a percentage, with the symbol . Precipitation data is used to determine whether precipitation exists at the current moment; its value is either 0 or 1, denoted as . Cloud cover data, expressed as a percentage of total cloud cover, or coded as sunny, partly cloudy, overcast, etc., is used to determine cloud cover levels; Current time information is used to determine the daytime and nighttime phases and calculate the solar altitude angle, denoted by [symbol missing]. The aforementioned meteorological observation data can be obtained in real time from meteorological stations near the target road section.

[0024] Road surface temperature is primarily influenced by air temperature, but deviations can occur due to radiation and wind. During the day, shortwave solar radiation causes road surface temperature to be higher than air temperature, while at night, longwave radiation loss from the ground causes road surface temperature to be lower than air temperature. Additionally, convective cooling from wind and the absorption of latent heat from road surface moisture evaporation also affect road surface temperature. Based on these physical mechanisms, this embodiment uses the following simplified formula to calculate the road surface temperature inversion value: ; In the above formula, This is the current road surface temperature inversion value; This is the actual measured temperature; This is a correction term for radiation temperature difference; For air-cooled correction; This is a correction item for evaporative cooling. The three correction items are explained below.

[0025] The calculation of the radiation temperature difference correction term requires first distinguishing between daytime and nighttime periods. In this embodiment, the division between daytime and nighttime is not based on fixed times, but rather on sunrise and sunset times calculated using astronomical formulas based on the latitude and longitude information of the target road segment. When the current time falls between sunrise and sunset, it is determined to be the daytime period; when the current time falls between sunset and the sunrise time of the next day, it is determined to be the nighttime period. In practical applications, to simplify calculations, an approximate fixed time period can also be used as the daytime period under specific regional and seasonal conditions.

[0026] The intensity of solar shortwave radiation is affected not only by atmospheric transparency but also by the sun's geometric position. This embodiment calculates the solar altitude angle using the latitude and longitude information of the target road segment and the current time. And introduce solar elevation angle modulation factor It is used to characterize the differences in radiative heating effects under different solar positions during the daytime phase.

[0027] Solar altitude angle Calculate using the following steps: First, determine the year-end based on the current date. n , n This is the day number starting from January 1st of the current year, with values ​​ranging from 1 to 365 or 366.

[0028] Calculate the solar declination angle : .

[0029] Calculating true solar time : ; In the formula, The local standard time for the target road segment is in hours. The longitude of the target road segment is in degrees, with east longitude positive. The longitude of the standard meridian of the time zone, in degrees; The time difference, in minutes, is calculated using the following formula: ; ; In the formula, The angle is the solar angle, and the unit is degrees.

[0030] Calculate the solar hour angle : .

[0031] Calculate the solar altitude angle α: ; In the formula, The latitude of the target road segment is expressed in degrees, with positive values ​​for North latitude.

[0032] Solar elevation angle modulation factor Determine by the following formula: when At that time, the sun is below the horizon, and the contribution of solar shortwave radiation is negligible, making ;when season , which serves as a normalized characterization of solar radiation intensity.

[0033] Therefore, during the daytime phase, the radiation temperature difference correction term is calculated using the following formula: ; During the nighttime phase, there is no solar shortwave radiation input, and the radiation temperature difference correction term is calculated using the following formula: .

[0034] In the above two formulas, This is the radiation correction factor, used to characterize the heating effect of shortwave solar radiation on the road surface; The diurnal radiation coefficient is used to characterize the cooling effect of long-wave radiation loss from the ground surface at night on the road surface. Relative humidity; This is the solar elevation angle modulation factor.

[0035] By introducing solar altitude angle modulation, the model can reflect the continuous changes in solar radiation intensity during sunrise, noon, and sunset, avoiding overestimation at sunrise and sunset and underestimation at noon caused by giving the same radiation correction amount throughout the daytime.

[0036] It should be noted that in road application environments where direct solar radiation observation data is lacking, relative humidity is used to characterize the effect of atmospheric transparency, and solar altitude angle is used to characterize the effect of the sun's geometric position. Together, they constitute a simplified physical description of daytime radiation.

[0037] The above formula uses relative humidity as a dimensionality-reduced proxy variable for cloud cover. Under typical meteorological scenarios of low humidity on sunny days and high humidity on cloudy days, there is a good inverse correlation between relative humidity and solar radiation flux, which can effectively characterize the strength of the radiative heating effect.

[0038] However, in atypical meteorological scenarios such as high-humidity sunny days or low-humidity cloudy days, relative humidity cannot be completely equated to cloud cover, leading to errors in the aforementioned proxy relationship. To address this issue, this embodiment introduces a dynamic compensation term for atypical daytime weather conditions. and atmospheric window transparency correction operator Cross-validation was performed by combining weather phenomenon codes with relative humidity.

[0039] Under typical meteorological conditions, weather phenomenon codes are consistent with sky conditions characterized by relative humidity. At this time, the dynamic compensation item for atypical daytime weather... and atmospheric window transparency correction operator When all values ​​are set to zero, the radiation correction formula degenerates into a basic radiation correction model based on relative humidity and solar altitude angle modulation factors. Under atypical meteorological conditions, the corresponding compensation parameters are triggered according to the combination relationship between weather phenomenon codes and relative humidity.

[0040] This embodiment sets compensation parameters for different atypical meteorological scenarios, including positive compensation parameters for daytime conditions under sunny and humid conditions. Negative compensation parameters during the day under cloudy and low humidity conditions Nighttime negative compensation parameters under clear and dry nighttime conditions And the positive compensation parameters for nighttime under conditions of high humidity or light fog. .

[0041] The aforementioned compensation parameters are predetermined using a global reference database. Specifically, historical data of reference road sections equipped with pavement temperature monitoring devices are extracted, and the reference road sections are classified according to climate zone type, pavement material, and micro-topographic features. Within each category, corresponding atypical meteorological samples are extracted based on weather phenomenon codes and relative humidity thresholds. The base pavement temperature inversion value for each sample is calculated under the condition that the compensation term is zero. And calculate the road surface temperature inversion residual using the following formula: ; In the formula, For the first k Road surface temperature inversion residuals for each sample; For the first k The measured road surface temperature corresponding to each sample; The base pavement temperature inversion value is calculated when the compensation term is zero.

[0042] For the same type of atypical meteorological samples, the median or truncated mean of the residuals of this type of samples is used as the initial value of the corresponding compensation parameter; alternatively, the corresponding compensation parameter can be determined by minimizing the sum of squares of the road surface temperature inversion residuals of this type of samples. The selectable range of compensation parameters is determined by using a preset quantile interval of the residuals of this type of historical samples, for example, using the interval corresponding to the 10th percentile to the 90th percentile, in order to reduce the impact of individual abnormal samples on the parameter determination results.

[0043] When the weather phenomenon is coded as sunny and the relative humidity (RH) is greater than 70%, it indicates that the sky is actually clear but the air humidity is high. Using relative humidity alone as a cloud cover proxy will underestimate the heating effect of solar shortwave radiation. In this case, the daytime positive compensation parameter under sunny and high humidity conditions is triggered. and order ,and >0, to compensate for the underestimation of radiant heating caused by excessively high relative humidity.

[0044] When the weather phenomenon is coded as cloudy and the relative humidity (RH) is less than 30%, it indicates that there is cloud cover but the air is relatively dry. Using relative humidity alone as a proxy for cloud cover will overestimate the heating effect of solar shortwave radiation. This triggers the daytime negative compensation parameter under cloudy, low-humidity conditions. and order ,and <0, to correct the overestimation of radiative heating effect.

[0045] For nighttime scenarios, this embodiment employs a longwave radiation prediction model based on relative humidity and introduces an atmospheric window transparency correction operator. When it is nighttime, the weather phenomenon is coded as clear, and the relative humidity (RH) is less than 30%, it indicates that there are few clouds and the air is dry, the atmospheric window is highly transparent, and the loss of longwave radiation from the surface is enhanced. At this time, the nighttime negative compensation parameter under clear and dry nighttime conditions is triggered. and order ,and This will increase the temperature correction range corresponding to nighttime long-wave radiation cooling and reduce the possibility of underreporting icing risk.

[0046] When it is nighttime and the relative humidity (RH) is greater than 85%, or when the weather phenomenon code indicates light fog, it indicates a high water vapor content in the near-surface layer, enhanced atmospheric back radiation, and possibly the release of latent heat of water vapor condensation. This triggers the nighttime positive compensation parameters under conditions of high humidity or light fog. and order ,and >0, in order to reduce the overestimation of the heat dissipation effect of long-wave radiation at night by the basic model and reduce the possibility of false alarms for icing warnings.

[0047] For target road sections without any pavement temperature sensors, the K-nearest neighbor algorithm is used to select the best-matching reference road section parameter cluster from the global reference database based on the climate zone type, pavement material, and micro-topographic features of the target road section. The various compensation parameters in this parameter cluster are then used as the initial compensation parameters for the target road section. During the long-term operation of the target road section, the corresponding compensation parameters can be adaptively fine-tuned by incorporating natural meteorological calibration events (described later).

[0048] The wind-cooling correction term characterizes the effect of wind convection cooling on road surface temperature and is calculated using the following formula: ; In the formula, This is the correction factor for air cooling; This represents the wind speed at the current moment.

[0049] The evaporative cooling correction term characterizes the cooling effect of latent heat absorption on pavement temperature during moisture evaporation, and is calculated using the following formula: ; In the formula, The road surface wettability coefficient ranges from 0 to 1. Let be the evaporative cooling coefficient. The above equation is a first-order linearized approximation of latent heat cooling: the natural convection evaporation term in the water film budget. This is neglected in the energy budget because the temperature drop corresponding to wind speeds less than 1 m / s is typically less than 0.2℃, which is below the resolution of this first-order linearized model; the change in saturated vapor pressure with air temperature is... The evaporative cooling coefficient is absorbed within the freezing temperature range of 10℃ to 10℃. , Updated online by natural calibration events, it adaptively adjusts with seasonal temperature shifts, thus the evaporative cooling correction term no longer explicitly includes saturated vapor pressure. Road surface wettability coefficient. The value of is not determined by hardware sensors, but is dynamically generated through the following virtual water film dynamics equation.

[0050] This embodiment introduces a road surface hydrological budget balance equation based on Dalton's law of evaporation, by tracking the thickness of the virtual water film. This method controls the on / off state of the evaporative cooling effect, thus solving the problem of spurious calculations of the evaporative cooling term under dry road surface conditions. The specific process is as follows.

[0051] During the precipitation phase, when meteorological data identifies the precipitation status, the virtual water film thickness... Assigned initial value In this embodiment The value is set to 1.0 mm, representing the maximum water film thickness constant of the road surface. At this point, the road surface wetting coefficient... Setting it to 1 indicates that the road surface is completely wet.

[0052] After the precipitation ends, the system enters a dry phase. Based on Dalton's law of evaporation, the system calculates the... t Evaporation rate of road surface water film at any time To clarify the physical meaning and dimensions of each parameter, the saturated vapor pressure corresponding to the current air temperature is first calculated using the Magnus approximation formula: ; In the formula, For the first The saturated vapor pressure corresponding to the temperature at a given time, in units of ; For the first t The temperature at that moment is in °C; 17.67 is a dimensionless empirical coefficient; 243.5 is the Magnus empirical temperature constant corresponding to the Celsius temperature, in °C.

[0053] The calculation is based on saturated water vapor pressure, relative humidity, and wind speed. t Evaporation rate of road surface water film at any given time: ; In the formula, For the first The rate of evaporation of the water film on the road surface at any given time, in units of ; The natural convection evaporation coefficient under windless conditions is given in units of . ; The wind speed mass transfer coefficient is expressed in units of 1000 m / s. ; For the first t Wind speed at any given moment, in units of ; For the first t Relative humidity at any given time, expressed as a percentage.

[0054] in, Units and Similarly, the relative humidity deficit term is a dimensionless quantity, and the product of all terms results in the evaporation rate. The unit is .

[0055] The preferred value range is... to In this embodiment, it is preferred to take ; The preferred value range is... to In this embodiment, it is preferred to take . and The initial value can also be selected from the corresponding parameter cluster in the global reference database based on the climate zone type, road surface material and micro-topographic features of the target road section, or it can be calibrated using the water film thickness change data after the precipitation ends of the reference road section, with the goal of minimizing the error between the calculated and observed water film thickness values.

[0056] The virtual water film thickness at the current moment as time progresses. Update according to the following integral decay equation: ; In the formula, This refers to the elapsed time after the end of the precipitation, in units of 1200 minutes. ; The unit is The unit of the integral variable is Therefore, the unit of the integral result is ,and and The units are consistent.

[0057] Road surface wettability coefficient The mapping relationship with the virtual water film thickness is as follows: ; When integral decay leads to virtual water film thickness When the coefficient of wettability drops to zero, it indicates that the road surface is completely dry, at which point the road surface wettability coefficient is zero. If the value is zero, the system automatically cuts off the calculation of the evaporative cooling correction term. This mechanism ensures that the evaporative cooling effect only participates in the correction of road surface temperature when there is actual moisture on the road surface, avoiding false cooling calculations under dry road surface conditions.

[0058] In the aforementioned pavement temperature retrieval model, the radiation correction coefficient A, diurnal radiation coefficient B, wind cooling correction coefficient C, and evaporative cooling coefficient D are key model parameters. Traditional methods require the deployment of pavement temperature monitoring equipment on the target road section to collect real pavement temperature data over a period of time and perform offline parameter calibration. To eliminate the dependence on measured pavement temperature sensors and resolve the logical contradiction between traditional model parameter calibration and sensorless application scenarios, this embodiment proposes an unsupervised adaptive calibration method based on physical boundary constraints of natural meteorological events.

[0059] This embodiment first obtains initial parameters based on prior knowledge of a multi-source heterogeneous feature space. A global reference database is constructed, extracting road segment data from high-precision traffic and meteorological stations across the country. Climate zone type, road surface material, and micro-topographic features are used as three-dimensional labels. Historical measured data is used to invert and generate basic parameter clusters covering different scenarios. Climate zone type, road surface material, and micro-topographic features are all taken from preset enumerated values. For example, climate zone type is taken as temperate, warm temperate, subtropical, etc., road surface material is taken as asphalt or cement concrete, and micro-topographic features are taken as plains, bridges, tunnel entrances, high embankments, deep excavations, etc. For target warning road segments without sensors, the system uses the inherent static label of the road segment and takes the Hamming distance according to the three-dimensional label. The same value is counted as 0 and different values ​​as 1. The system extracts the K reference road segments with the smallest distance in the global database through the K nearest neighbor algorithm, and uses the average parameter value of these K reference road segments as the initial model parameter set. K is preferably between 3 and 7, and in this embodiment, it is taken as 5. The basic parameter set includes radiation correction factor A, diurnal radiation factor B, air cooling correction factor C, evaporative cooling factor D, and natural convection evaporation factor. Wind speed mass transfer coefficient and each compensation parameter , , , The parameter cluster also records the historical distribution of each parameter under this category, which is used to determine the upper and lower bounds in the parameter calibration process described later.

[0060] After obtaining the initial parameters, the system can operate without human intervention for a long time. By analyzing regular weather forecast data in real time, it automatically captures natural calibration events with the constraints of the laws of physical thermodynamics to extract the pseudo-true value of the road surface temperature, and updates only the parameters dominated by the event each time.

[0061] In this embodiment, the pseudo-true value of road surface temperature is preferentially derived from natural meteorological events with mandatory physical constraints; in the period when there are no available physical constraint events, the median measured road surface temperature of reference road sections with the same climate zone type, the same road surface material, and the same micro-topographic features in the global reference database under the same solar altitude angle and similar meteorological conditions is used as the reference road surface temperature; similar meteorological conditions refer to the temperature difference not exceeding 2℃ and the same cloud cover level and wind speed level.

[0062] In this embodiment, the snow or ice condition on the road surface is estimated by meteorological observation data: when precipitation data indicates precipitation and the temperature is not higher than 0°C, it is determined that there is snow or ice on the road surface and a snow accumulation sign is set; after the snow accumulation sign is set, when the temperature rises from not higher than 0°C to higher than 0°C, it is determined that melting has started; when the cumulative positive accumulated temperature since the start of melting reaches the preset melting accumulated temperature threshold, it is determined that melting has ended and the snow accumulation sign is cleared.

[0063] Natural calibration events include the following five types: Evaporative cooling parameter calibration event: When the system detects the end of precipitation and the virtual water film thickness. Greater than zero, road surface wettability coefficient When the wind speed is greater than or equal to 0.8 and the wind speed is greater than 2 m / s, the road surface is in a nearly completely wet state. Heat exchange on the road surface is dominated by latent heat, and the road surface temperature is forcibly flattened to the wet-bulb temperature of the atmosphere. In this case, the pseudo-true value for the road surface temperature is set to the wet-bulb temperature calculated based on the Stull empirical formula. The main update is to the evaporative cooling coefficient D, while other parameters remain unchanged.

[0064] The empirical formula for Stull is: ; In the formula The value represents the wet-bulb temperature in °C; the arctangent function values ​​are expressed in radians; this formula applies to air temperatures under standard atmospheric pressure. Conditions of 20℃ to 50℃ and relative humidity of 5% to 99%.

[0065] Nighttime radiation parameter calibration event: When the system identifies the nighttime phase, the weather phenomenon is coded as clear, the relative humidity (RH) is less than 60%, there is no precipitation, the virtual water film thickness is zero, and the wind speed is less than 2 m / s, the loss of long-wave radiation from the ground surface is the main correction factor, and the effects of wind cooling and evaporative cooling can be ignored. If the meteorological observation record indicates that condensation or frost has occurred at this time, the road surface temperature has dropped to the dew point or frost point. The pseudo-true value of the road surface temperature is set as follows: when the air temperature is above 0℃, the dew point temperature calculated from the air temperature and relative humidity is used; when the air temperature is not above 0℃, the corresponding frost point temperature is used. Dew point temperature is calculated using the following formula: ; Frost point temperature is calculated using the following formula: ; ; In the formula, This refers to the dew point temperature. These are the frost point temperatures, all in °C. Determined by air temperature and relative humidity, the temperature term uses the water surface coefficient from the Magnus formula, corresponding to the definition of relative humidity using the saturated water vapor pressure at the water surface as a reference. The dew point temperature formula is the inverse solution of the Magnus formula under the water surface coefficient; the frost point temperature formula uses the same actual water vapor pressure as input and is solved inversely according to the ice surface coefficient, therefore both formulas share the same... The only difference is the coefficients of the outer layer of the inverse solution; the two cannot be used interchangeably.

[0066] The two equations are the inverse solutions of the Magnus formula under the water surface coefficient and the ice surface coefficient, respectively.

[0067] If there are no records of condensation or frost, the reference road surface temperature will be used. This event primarily updates the diurnal radiative coefficient B, while other parameters remain unchanged.

[0068] Air-cooled parameter calibration event: When the system identifies no precipitation, zero virtual water film thickness, and the same cloud cover level for 6 consecutive hours, with wind speeds between 2 m / s and 8 m / s and a range of at least 3 m / s within this time window, evaporative cooling has been cut off, radiation conditions are relatively stable, and the variation of road surface temperature inversion residuals with wind speed is mainly contributed by the air-cooled correction term. At this time, the difference between the road surface temperature inversion value at the moment of minimum wind speed within this time window and the reference road surface temperature at the same moment is used as the bias correction. The pseudo-true value of road surface temperature at each moment is the sum of the reference road surface temperature at the corresponding moment and the bias correction. A linear regression is performed on the wind speed using the difference between the pseudo-true value of road surface temperature and the air temperature at each moment, and the air-cooled correction coefficient C is updated with the negative value of the regression slope, while other parameters remain unchanged.

[0069] Daytime Radiation Parameter Calibration Event: When the system identifies a daytime period, the weather phenomenon is coded as sunny, there is no precipitation, the virtual water film thickness is zero, and the wind speed is less than 2 m / s, solar shortwave radiation is the primary correction factor. If snow or ice is detected on the road surface and melting begins during this period, the pseudo-true value for the road surface temperature is set to 0℃ based on the physical law that temperature stagnation occurs during snow and ice melting. Considering that meltwater forms a thin layer of liquid water on the road surface and absorbs latent heat of vaporization, this event only uses samples from the first two hours after melting begins for calibration to limit the latent heat from being mistakenly included in the radiation correction coefficient. If there is no snow or ice melting process, the reference road surface temperature is used. This event primarily updates the radiation correction coefficient A; all other parameters remain unchanged.

[0070] Extreme Freezing Constraint Event: When the system detects snow accumulation on the road surface in winter and the temperature rises above 0°C, it utilizes the physical law that the latent heat of snow melting causes temperature stagnation, forcibly setting the pseudo-true value of the road surface temperature to 0°C and maintaining it until the snow melt ends. This event is used to compensate parameters triggered in the current time period. , , or Update the parameters while keeping the rest unchanged; if this event and the daytime radiation parameter calibration event are both met, this event will be processed first.

[0071] Whenever the system identifies a natural meteorological calibration event, it performs an online calibration for the model parameters dominated by that event. To reduce the impact of a single abnormal meteorological observation on the parameter update results, this embodiment employs a projection gradient descent method based on a sliding time window.

[0072] The sliding time window is from the most recent m This consists of several natural calibration events of the same type. In this embodiment, m The optimal value is 10 to 30. When the system identifies a new calibration event of the same type, it adds the event sample to the sliding time window and deletes the oldest sample, ensuring that the number of samples in the sliding time window does not exceed [a certain value]. m .

[0073] For each sample within the sliding time window, the residual between the road surface temperature inversion value under the existing model parameters and the pseudo-true value of the road surface temperature determined by the corresponding physical law is calculated, and the average of the squared residuals of all samples within the time window is used as the loss function: ; In the formula, For the first The loss function value corresponding to the next iteration; This represents the number of samples in the sliding time window; For the first k The road surface temperature inversion value corresponding to each calibration sample; For the first k The pseudo value of the road surface temperature corresponding to each calibration sample.

[0074] Because the radiation correction factor A, diurnal radiation factor B, wind-cooling correction factor C, and evaporative cooling factor D have different dimensions and numerical ranges, the parameters to be calibrated are first normalized to the range of 0 to 1 before performing gradient updates. For any parameter to be calibrated... Its normalized value x Calculate using the following formula: ; In the formula, The current value of the parameter to be calibrated; This is the lower bound of the parameter; This is the upper bound of the parameter; The normalized parameter value ranges from 0 to 1.

[0075] This embodiment uses an adaptive learning rate that gradually decreases with the number of iterations: ; In the formula, For the first The learning rate used in this iteration; The initial learning rate; This is the iteration number in this calibration process. Initial learning rate. The value is preferably between 0.01 and 0.05, and in this embodiment, it is preferably 0.02.

[0076] The parameter update direction is determined by the gradient of the loss function with respect to the normalized parameters. To prevent abrupt parameter changes caused by a single abnormal calibration event, the normalization change of a single parameter in each iteration must not exceed a preset maximum update magnitude. ,Right now: ; in The preferred value is 0.01 to 0.05, indicating that the magnitude of a single parameter update does not exceed 1% to 5% of the allowable value range of the parameter. In this embodiment, 0.03 is preferred. After the normalized parameter update is completed, it is converted into a parameter value with actual physical units.

[0077] After obtaining candidate parameter values ​​in each iteration, the candidate parameter values ​​are restricted to the corresponding physically reasonable range: ; In the formula, These are the candidate parameter values ​​obtained after gradient update; These are the parameter values ​​after the upper and lower bound constraints are completed.

[0078] The upper and lower bounds of each parameter are primarily determined based on the global reference database: Within parameter clusters that match the climate zone type, pavement material, and micro-topographic features of the target road segment, the 5th percentile of the historical distribution of each parameter is used as the lower bound, and the 95th percentile as the upper bound. When there are insufficient reference samples matching the target road segment, the following safety margins are used: Values ​​for A range from 0℃ to 25℃; values ​​for B range from 0℃ to 15℃; and values ​​for C range from 0 to 1.0. The value of D ranges from 0 to 1.5. Where A and B represent the maximum temperature difference coefficients corresponding to the daytime radiative heating effect and the nighttime radiative cooling effect, respectively, in °C; C and D represent the wind-cooling correction coefficient and the evaporative cooling correction coefficient per unit wind speed, respectively, in units of _____. .

[0079] The aforementioned numerical range is used to limit parameter updates and prevent parameters from going out of bounds; the initial parameters of the target road segment are still determined by the matching parameter cluster in the global reference database. In addition, when the inverted pavement temperature is lower than the wet-bulb temperature calculated by the Stull empirical formula at the same time, the wet-bulb temperature is taken as the pavement temperature inversion value to avoid inversion distortion caused by the linear superposition of wind cooling and evaporative cooling under high wind conditions.

[0080] In this calibration process, the parameters are considered to have converged when the following two conditions are met simultaneously and maintained for three consecutive iterations: First, the relative rate of change of the loss function between two adjacent iterations is less than 0.0001; second, the maximum normalized change of all parameters to be calibrated is less than 0.001. The relative rate of change of the loss function is calculated using the following formula: ; In the formula, For the first The relative rate of change of the loss function in each iteration; For the first The loss function value of the next iteration; This is the loss function value from the previous iteration; In this embodiment, a positive number is set to prevent the denominator from being zero. Take 0.000001.

[0081] The maximum number of iterations for each natural weather calibration event is preferably between 20 and 100, and in this embodiment, it is preferably set to 50. When the number of iterations reaches 50 but the convergence condition has not been met, the calibration is stopped, and the parameter value corresponding to the minimum loss function during this iteration is adopted. If the minimum loss function after calibration is still greater than the loss function before calibration, or if there are obvious missing values, outliers, or sudden changes in weather conditions in the calibration sample, the parameter update is abandoned, and the parameter value before calibration is retained.

[0082] As seasons change and road surfaces age, the system continuously triggers the aforementioned calibration events during long-term operation, enabling unsupervised adaptive iteration of parameters constrained by these events, thus achieving calibration-free deployment; natural convection evaporation coefficient With wind speed mass transfer coefficient It does not participate in online calibration and keeps the initial value unchanged.

[0083] Step S2: Construct the seasonal baseline curve.

[0084] To avoid the accumulation of errors in long-term weather forecasts, this step constructs a road surface temperature background curve that reflects recent climate characteristics, serving as a benchmark reference for long-term forecasts.

[0085] Hourly meteorological observation data for the past N days are collected; in this embodiment, the preferred value for N is 7 days. The hourly meteorological observation data includes measured air temperature, wind speed, relative humidity, and precipitation. Using the pavement temperature inversion model in step S1, the historical pavement temperature inversion values ​​for each hour of each day over the past N days are calculated. .

[0086] Traditional methods using arithmetic averaging can easily lead to the cancellation of heterogeneous weather features. For example, the temperature curves for sunny and rainy days differ significantly in shape, and simple averaging will blur the physical characteristics of the two. To address this issue, this embodiment introduces meteorological isomorphism weights, performing similarity matching by extracting weather feature vectors.

[0087] First, obtain hourly cloud cover forecast data and hourly wind speed forecast data for the target road section during the future forecast period, and then obtain past data. N Hourly measured cloud cover data and hourly measured wind speed data for the entire day.

[0088] This embodiment divides the 168-hour forecast period into 7 consecutive forecast days, each forecast day including 24 hourly forecast times. A forecast weather feature vector is constructed for each forecast day, and this feature vector is then compared with past forecasts. N The measured weather feature vectors for each historical day in Tianzhong are compared.

[0089] This embodiment classifies cloud cover into four levels based on the percentage of total cloud cover: Level 1 cloud cover is defined as total cloud cover greater than or equal to 0% and less than 20%; Level 2 cloud cover is defined as total cloud cover greater than or equal to 20% and less than 50%; Level 3 cloud cover is defined as total cloud cover greater than or equal to 50% and less than 80%; and Level 4 cloud cover is defined as total cloud cover greater than or equal to 80% and less than or equal to 100%. When the acquired cloud cover data is represented as a decimal between 0 and 1, the cloud cover data is multiplied by 100% before being classified according to the above thresholds. When meteorological data only provides weather phenomenon codes such as sunny, partly cloudy, and overcast, a correspondence between weather phenomenon codes and cloud cover levels is pre-established, mapping sunny to Level 1 cloud cover, partly cloudy or partially cloudy to Level 2 cloud cover, partly cloudy to Level 3 cloud cover, and overcast or completely cloud-covered to Level 4 cloud cover.

[0090] This embodiment is based on hourly wind speed. V Wind speed is divided into 4 levels: wind speed V Wind speeds less than 2 m / s are classified as Level 1; wind speeds greater than or equal to 2 m / s and less than 4 m / s are classified as Level 2; wind speeds greater than or equal to 4 m / s and less than 6 m / s are classified as Level 3; and wind speeds greater than or equal to 6 m / s are classified as Level 4.

[0091] For any future forecast date, count the number of hours each cloud cover level and each wind speed level occurs within the 24 hours of that forecast date, and divide each by 24 to obtain the time proportion of each level within that forecast date. Arrange the time proportions of cloud cover levels 1 to 4 and wind speed levels 1 to 4 in order to form the forecast weather feature vector for that future forecast date: Forecast weather feature vector = (Level 1 cloud cover percentage, Level 2 cloud cover percentage, Level 3 cloud cover percentage, Level 4 cloud cover percentage, Level 1 wind speed percentage, Level 2 wind speed percentage, Level 3 wind speed percentage, Level 4 wind speed percentage).

[0092] Therefore, the forecast weather feature vector in this embodiment is an 8-dimensional feature vector, including four cloud cover level proportion components and four wind speed level proportion components. Using the same cloud cover level threshold, wind speed level threshold, and proportion statistics method, an 8-dimensional measured weather feature vector for each historical day in the past N days is constructed.

[0093] Since each component in the weather feature vector is the ratio of the number of hours the corresponding level occurs to 24 hours, and each component is a dimensionless value in the range of 0 to 1, the forecast weather feature vector and the historical weather feature vector have the same physical meaning and numerical range in each dimension. Therefore, there is no need to perform minimum-maximum normalization processing before calculating the cosine similarity.

[0094] Calculate the first one respectively d The weather feature vector for the predicted future date and the past N The cosine similarity between the measured weather feature vectors of each historical day in the sky is used, and this cosine similarity is directly used as the correlation between that historical day and the first historical day. d Meteorological isomorphic weights for each future forecast date . The larger the value, the higher the value. i The cloud cover level distribution and wind speed level distribution of the historical day are compared with the first d The closer a future prediction date is to the historical date, the greater the weight of the historical road surface temperature inversion value for that historical date in the calculation of the baseline curve for that prediction date. Cosine similarity can also be replaced by Gaussian kernel similarity or similarity based on normalized Euclidean distance; the negative exponential fusion weight in step S3 can also be replaced by piecewise linear weights that monotonically decrease with the forecast lead time, and both have the same function in the embodiments of this invention.

[0095] No. d The formula for calculating the seasonal baseline curve for each future forecast date is: ; In the formula, d The predicted date number is 1 to 7. From 1 to ; Hourly serial number; For the first Heavenly Historical road surface temperature inversion values ​​for hours. Using the above weighting method, if a historical day is relative to the [number]th [hour] [day], [the following is a calculation method ... d If extreme heterogeneous weather occurs on a future predicted date, then the weight corresponding to that historical date... Approaching zero, it is automatically filtered out, thus extracting a pure background curve that is highly isomorphic to the weather conditions for each future forecast day.

[0096] The daily variation patterns of each future prediction day are mapped to the corresponding prediction time period according to the hour sequence. In this embodiment, the historical road surface temperature inversion values ​​are organized according to the hour of the calendar day, i.e. Represents the hourly rate within a day, ranging from 0 to 23; Seasonal baseline curve What is depicted is the first One future prediction day Background road surface temperature at the time. Due to the forecast lead time. The forecast period is calculated from the prediction start point, which may not be at midnight. Therefore, the forecast lead time must be converted to hours before mapping: Let the hourly number of the predicted start point within the day be denoted as... If the value is an integer between 0 and 23, then the forecast lead time is... ,Pick , That is, the baseline curve for the forecast period is Thus, the background road surface temperature curve for the next 168 hours is obtained.

[0097] If taken directly If the prediction starting point is not at zero, then the day-night phase of the baseline curve will be shifted relative to the actual time. Hours, up to a maximum of 23 hours, will apply historical day and night background temperatures to future daytime periods, contradicting the diurnal variation patterns that should be preserved by diurnal cycle mapping. (Introduction) After that, regardless of when the forecast starts, the peak and trough positions of the baseline curve always align with the actual clock.

[0098] It should be noted that in this embodiment, the future forecast days are divided into 24-hour periods starting from the forecast start point, while the historical days are divided into calendar days. The window start points for the two can be different. Each component of the weather feature vector is the ratio of the number of hours the corresponding level occurs to 24 hours, which represents the proportion of weather composition within the window. Its physical meaning does not change with the window start point, so it does not affect the calculation of cosine similarity and the determination of meteorological isomorphism weights.

[0099] Compared to using the same 24-hour curve for all forecast days, the above-mentioned forecast-day weighting method allows the background curve to change according to the differences in weather conditions on each future day.

[0100] Step S3: Long-term fusion prediction.

[0101] First, based on the measured meteorological observation data from the weather station at the current moment, the current road surface temperature inversion value is calculated using the road surface temperature inversion model in step S1. This value serves as the final predicted road surface temperature at the zero point of the prediction time.

[0102] Then, obtain the hourly temperature forecast data released by the meteorological department for the next 7 days, along with the corresponding wind speed, humidity, and precipitation forecasts. The virtual water film thickness during the forecast period is initialized with the value recursively derived from the observation data at the forecast start time, and updated hourly according to the hourly precipitation forecast: when the precipitation forecast indicates precipitation at that time, the virtual water film thickness is reassigned. And order M =1; When there is no precipitation, substitute the predicted air temperature, relative humidity, and wind speed at that moment into the aforementioned evaporation rate formula, and decrease the virtual water film thickness hourly according to the integral decay equation to obtain the road surface wetting coefficient at each predicted moment. M Substitute the above forecast data and the pavement wetting coefficient at each forecast time into the pavement temperature inversion formula in step S1 to calculate the predicted pavement temperature inversion value for each future hour. .

[0103] To scientifically characterize the evolution of weather forecast errors over time, this embodiment employs a negative exponential decay model based on atmospheric predictability duration to construct dynamic fusion weights. The formula for calculating the dynamic fusion weights is as follows: ; In the formula, For forecast lead time, it is in hours; To predict the decay characteristic time constant, The value of is greater than zero; The physical meaning is the time required for the forecast confidence level to drop to 1 / e, or approximately 37%. Parameter The forecast duration was determined based on historical error statistics of the numerical models of the meteorological bureau in the area where the target road segment is located. By statistically analyzing the root mean square error (RMSE) between the numerically predicted and observed temperatures in the region over the past year, the forecast duration corresponding to when the RMS error reaches 50% of the historical temperature statistical standard deviation was calculated, and this duration was used as the forecast duration. The value of is set as follows: For example, if the root mean square error of a numerical weather prediction model for a certain region increases to 50% of the historical temperature statistical standard deviation when the forecast lead time reaches 120 hours, then is set as follows. =120. For meteorological stability operators, The value of is not less than 1 and not greater than 2. When the meteorological authority issues a blue or higher cold wave warning signal for the area where the target road section is located, or when meteorological forecast data shows that at least one of the following conditions is met: a daily minimum temperature drop of not less than 8℃ in the next 24 hours, a daily minimum temperature drop of not less than 10℃ in the next 48 hours, or a daily minimum temperature drop of not less than 12℃ in the next 72 hours, and the daily minimum temperature after the temperature drop is not higher than 4℃, it is determined that a large-scale strong cold air system or cold wave process exists in the future. At this time, the atmospheric circulation is controlled by a large-scale system, and the forecast certainty is enhanced. Take a value greater than 1 and no greater than 2, such as 1.2, to extend the trust period for the forecast data; An excessively large value would cause the baseline curve to lose its long-term constraint effect; therefore, 2 is set as the upper limit. Under normal meteorological conditions... Set to 1, and the dynamic fusion weight decays exponentially according to the usual method.

[0104] The final formula for predicting road surface temperature is: ; In the formula, To predict the inversion value of road surface temperature; This serves as the baseline curve for the forecast period.

[0105] The above formula is applicable to predicting lead time. The period of time when the value is greater than zero. For That is, the predicted starting time and the final predicted road surface temperature are directly taken from the inverted value of the current road surface temperature. ,Right now This setting ensures that the starting point of the forecast curve is anchored to the best inversion result based on the current measured meteorological observation data, avoiding the propagation of systematic biases in the meteorological forecast data at the initial moment to subsequent forecasts.

[0106] The physical meaning of the above formula is that, in the early stages of forecasting, dynamic fusion weights... Approaching 1, the final predicted road surface temperature is mainly determined by the forecast inversion value; as the forecast lead time increases, The temperature gradually decreases towards zero, eventually predicting that the road surface temperature will gradually transition to the seasonal baseline curve. This mechanism trusts the baseline pattern during normal periods, and automatically triggers weight adjustments and prioritizes responses to forecast trends during forced weather events such as cold waves.

[0107] Step S4: Determine the icing warning based on the combination of road surface temperature and water film condition.

[0108] This step combines the road surface wetting coefficient calculated in step S1. and virtual water film thickness The icing determination logic is constructed based on the final predicted road surface temperature calculated in step S3.

[0109] When the final predicted road surface temperature If the temperature is below 0°C and there is a source of moisture on the road surface at the predicted time, the target road section is deemed to have a risk of icing at that time. The presence of a source of moisture on the road surface includes meeting at least one of the following conditions: First, the precipitation forecast indicates that precipitation is present at that time; second, the virtual water film thickness calculated using the aforementioned method at that time is greater than zero, meaning that the moisture left on the road surface by previous precipitation has not completely evaporated; third, the dew point temperature calculated based on the predicted air temperature and relative humidity at that time is higher than the final predicted road surface temperature at that time, and the relative humidity at that time is greater than 90%, at which point the road surface microenvironment meets the condensation conditions, and frost or dew may form on the road surface even without precipitation. It should be noted that the virtual water film thickness in this embodiment... Initialized solely by precipitation events, the trace water film formed by water vapor condensation does not participate in the update of the virtual water film thickness, nor in the calculation of the evaporative cooling correction term. However, this condensation condition still serves as an independent criterion for determining the freezing risk.

[0110] The above-mentioned judgment logic no longer relies solely on a single temperature threshold, but instead uses both the predicted road surface temperature and the residual thickness of the virtual water film as judgment conditions, thereby achieving a refined early warning of the risk of road icing and black ice.

[0111] Example 2: Evolution of Virtual Water Film under Solid Precipitation Conditions The difference between this embodiment and Embodiment 1 is that the precipitation falls onto the road surface in solid form. In this case, the virtual water film thickness is still initialized as described above. However, the saturated vapor pressure in the evaporation rate calculation The calculation is now based on the ice surface formula, where the coefficients 17.67 and 243.5 in the Magnus formula are replaced with 22.46 and 272.62 respectively. Simultaneously, since relative humidity is referenced to the saturated vapor pressure at the water surface, the humidity deficit term is no longer simplified. , No. The evaporation rate at any given time is calculated using the following formula: ; In the formula, For the first The saturated vapor pressure on the ice surface corresponding to the temperature at that moment; The two values ​​represent the saturated vapor pressure at the water surface corresponding to the air temperature at the same moment. They are calculated using the Magnus formula based on the ice surface coefficient and water surface coefficient, respectively, and the units are both... When the water film on the road surface is in a liquid state and Equal to the above, the above formula degenerates into the aforementioned form. The remaining steps are the same as in Example 1.

[0112] To verify the effectiveness of the method in this embodiment of the invention, a seasonal baseline curve was constructed using hourly measured meteorological data from February 22 to March 1, 2025. Hourly measured meteorological data from March 1 to March 8 were used instead of the forecast data for the same period as input to predict the road surface temperature for this 168-hour period, and the results were compared and analyzed with the measured road surface temperature data for the same period. The purpose of using measured data instead of forecast data is to isolate the inherent error of numerical forecasting and to separately test the accuracy of the road surface temperature inversion and curve fusion process. In practical applications, the error when inputting forecast data will be higher than the results of this verification.

[0113] This exemplary calculation verification embodiment selects meteorological data from a meteorological station in a certain city from March 1 to March 8, 2025 at 0:00, and obtains the measured road surface temperature data of the corresponding road section during the same period to verify the input data processing, parameter calling, hourly calculation and error evaluation process of the road surface temperature inversion method of the present invention.

[0114] Meteorological data includes air temperature, wind speed, relative humidity, precipitation data, cloud cover data, and time information. Measured road surface temperature data serves as the true reference in this verification embodiment. The meteorological data and the measured road surface temperature data are correlated using the same time base.

[0115] The above data is only used to illustrate the implementation and verification process of the present invention and does not limit the geographical scope of application of the method of the present invention.

[0116] Based on hourly meteorological data, the method of this invention is used to predict road surface temperature for the next 168 hours. The specific process includes: first, calculating the current road surface temperature inversion value using a road surface temperature inversion model based on the meteorological observation data at the current moment as the prediction starting point; then, calculating the predicted road surface temperature inversion values ​​for each future moment using the road surface temperature inversion model based on the forecasted air temperature, wind speed, relative humidity, and precipitation conditions; constructing a seasonal baseline curve based on the aforementioned historical meteorological data; fusing the forecasted inverted road surface temperature with the seasonal baseline curve according to the dynamic fusion weight corresponding to the prediction duration to obtain the final predicted road surface temperature value for the next 168 hours; and determining whether there is a risk of road icing at each moment based on whether the predicted road surface temperature does not exceed 0℃ and in combination with precipitation, relative humidity, and precipitation conditions in the previous period.

[0117] The 168-hour pavement temperature prediction curve obtained in this embodiment of the invention is compared with the measured pavement temperature curve during the same period, as follows: Figure 2 As shown in the figure. The results show that the predicted curve obtained by the method of the present invention has good consistency with the measured road surface temperature curve in terms of overall trend, diurnal temperature variation and identification of low temperature periods, and can accurately reflect the road surface temperature change process in the next 7 days.

[0118] To quantitatively evaluate the prediction accuracy, the Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) are used as evaluation indicators. Calculations show that the method of this embodiment of the invention achieves the following prediction results on the validation sample: MAE = 2.7℃ and RMSE = 3.6℃. These results demonstrate that, without relying on long-term historical pavement temperature samples for training and without requiring the deployment of pavement contact sensors on the target road section, the method of this embodiment can still effectively predict pavement temperature for the next 168 hours.

[0119] Different application areas use the same model structure and calculation steps. The model only matches the corresponding initial parameter set from the global reference database based on the climate zone type, road surface material and micro-topographic features of the target road segment. During operation, the model is continuously adjusted by natural meteorological calibration events, thereby achieving model adaptation to different road environments.

[0120] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described; only preferred embodiments of the present invention are illustrated. The descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the present invention. As long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification.

[0121] It should be noted that those skilled in the art can make various modifications and improvements without departing from the inventive concept, and these all fall within the scope of protection of this invention. Therefore, the scope of protection of this invention should be determined by the appended claims.

Claims

1. A method for predicting long-term road icing based on meteorological parameter inversion, characterized in that, Includes the following steps: Step S1: Collect meteorological observation data of the target road section at the current time, and apply radiation temperature difference correction term, wind cooling correction term and evaporative cooling correction term to the air temperature based on the meteorological observation data to obtain the current road surface temperature inversion value; The virtual water film thickness is initialized based on the precipitation data in the meteorological observation data, and the virtual water film thickness is updated hourly based on the evaporation rate. When the virtual water film thickness drops to zero, the evaporation cooling correction term is cut off. Step S2: Collect hourly meteorological observation data for the past N days, calculate the historical road surface temperature inversion value for each historical day, obtain meteorological forecast data, and generate a seasonal baseline curve by weighting the historical road surface temperature inversion value based on the similarity between the forecast weather characteristics and the weather characteristics of each historical day. Map the seasonal baseline curve to the future prediction period according to the daily cycle to obtain the prediction period baseline curve. Step S3: Calculate the predicted road surface temperature inversion value for each future time based on the meteorological forecast data. Using the current road surface temperature inversion value as the prediction starting point, fuse the predicted road surface temperature inversion value with the baseline curve of the prediction period to obtain the final predicted road surface temperature. Step S4: When the final predicted road surface temperature does not exceed 0℃ and there is a source of moisture on the road surface, it is determined that the target road section has a risk of icing.

2. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 1, characterized in that: The expressions for the current road surface temperature inversion value, the air-cooling correction term, and the evaporative cooling correction term in step S1 are as follows: ； ； ； In the formula, This is the current road surface temperature inversion value; Temperature from meteorological observation data; This is a correction term for radiation temperature difference; For air-cooled correction; This is a correction term for evaporative cooling. This is the correction factor for air cooling; The coefficient of performance is the evaporative cooling rate. Wind speed in meteorological observation data; Relative humidity in meteorological observation data; This is the road surface wettability coefficient.

3. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 2, characterized in that: The radiation temperature difference correction term is calculated separately for the daytime and nighttime periods, including the following steps: The sunrise and sunset times are calculated based on the latitude and longitude information of the target road segment. If the current time is between the sunrise and sunset times, it is determined to be the daytime stage; otherwise, it is determined to be the nighttime stage. Calculate the solar altitude angle based on the latitude and longitude information of the target road segment and the current time. And determine the solar elevation angle modulation factor. ,when hour ,when hour ; During the daytime phase, the expression for the radiation temperature difference correction term is: ； During the nighttime phase, the expression for the radiation temperature difference correction term is: ； In the formula, The solar altitude angle; Solar elevation angle modulation factor; This is the radiation correction factor; The diurnal radiation coefficient; This is a dynamic compensation item for atypical daytime weather conditions; Operator for atmospheric window transparency correction.

4. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 3, characterized in that: The daytime atypical meteorological dynamic compensation term and the atmospheric window transparency correction operator are valued in the following manner: During the daytime phase, the weather phenomenon in the meteorological observation data is coded as sunny and When it is greater than 70%, set Positive compensation parameters for daytime conditions under sunny and high humidity conditions ,and ; During the daytime phase, the weather phenomenon is coded as cloudy and When less than 30%, set Negative compensation parameters for daytime conditions under cloudy and low humidity conditions ,and ; During the nighttime phase, the weather phenomenon is coded as clear and When less than 30%, set The negative compensation parameter for nighttime under clear and dry nighttime conditions. ,and ; During the nighttime phase and When the weather phenomenon is greater than 85% or coded as light fog, set Nighttime positive compensation parameters under conditions of high humidity or light fog. ,and ; If none of the above four conditions are met, and All are set to zero.

5. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 2, characterized in that: The virtual water film thickness and the road surface wettability coefficient are determined as follows: When precipitation data in the meteorological observation data is identified as precipitation, the virtual water film thickness... Take initial value The road surface wettability coefficient Take 1; After the precipitation ends, the Magnus approximation formula is used to calculate the first... The saturated vapor pressure corresponding to the temperature at that moment: ； Calculate the first according to Dalton's law of evaporation. Evaporation rate of road surface water film at any given time: ； Update the virtual water film thickness according to the integral decay equation: ； The mapping relationship between the road surface wettability coefficient and the virtual water film thickness is as follows: ； In the formula, This is the elapsed time since the end of the precipitation. The time corresponding to the current moment value; The virtual water film thickness; This is the constant for the maximum water film thickness on the road surface; For the first The saturated vapor pressure corresponding to the temperature at a given moment; For the first Temperature at any given moment; For the first The rate of evaporation of the water film on the road surface at any given time; The coefficient of evaporation due to natural convection; The mass transfer coefficient is the wind speed. For the first Wind speed at any given moment; For the first Relative humidity at any given time.

6. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 4, characterized in that: The radiation correction coefficient, the diurnal radiation coefficient, the wind-cooling correction coefficient, the evaporative cooling coefficient, the daytime atypical meteorological dynamic compensation term, and the atmospheric window transparency correction operator are subjected to unsupervised adaptive calibration based on physical boundary constraints of natural meteorological events, including the following steps: Identify natural calibration events from the meteorological observation data, where the road surface temperature is constrained by physical laws to a known value; Based on the aforementioned natural calibration events, pseudo-true values ​​of road surface temperature are extracted; The residual is calculated by comparing the current road surface temperature inversion value with the road surface temperature pseudo value. The parameters are updated with the goal of minimizing the residual using the projection gradient descent method based on a sliding time window. Each calibration only updates the parameters corresponding to the natural calibration event. The sliding time window consists of the most recent preset number of the same type of natural calibration events, and the projected gradient descent method projects the candidate parameter values ​​to the interval defined by the lower and upper bounds of the updated parameters after each iteration.

7. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 6, characterized in that: The natural calibration events include the following five types: Evaporative cooling parameter calibration event: When the virtual water film thickness is greater than zero, the road surface wetting coefficient is not less than the preset wetting threshold, and the wind speed is greater than 2m / s after the precipitation ends, the pseudo-true value of the road surface temperature is taken as the wet-bulb temperature calculated according to the Stull empirical formula, and the evaporative cooling coefficient is updated accordingly. Nighttime radiation parameter calibration event: When the nighttime phase is in progress, the weather phenomenon is coded as sunny, the relative humidity is less than the preset humidity threshold, there is no precipitation, the virtual water film thickness is zero and the wind speed is less than 2m / s, and the meteorological observation record indicates the occurrence of dew or frost, the pseudo-true value of the road surface temperature is taken as the dew point temperature when the air temperature is higher than 0℃ and as the frost point temperature when the air temperature is not higher than 0℃, and the diurnal radiation coefficient is updated accordingly. Wind-cooled parameter calibration events: no precipitation, the virtual water film thickness is zero, the cloud cover in the meteorological observation data remains at the same cloud cover level within a preset time period, the cloud cover level is divided according to a preset cloud cover threshold and the wind speed appears in a stage not less than a preset range within a preset interval, the difference between the current road surface temperature inversion value at the moment of minimum wind speed in the stage and the reference road surface temperature at the same moment is used as the bias correction amount, the sum of the reference road surface temperature at each moment in the stage and the bias correction amount is used as the pseudo value of the road surface temperature at the corresponding moment, the difference between the pseudo value of the road surface temperature and the air temperature is used to perform a linear regression on the wind speed, and the negative value of the regression slope is used to update the wind-cooled correction coefficient accordingly, the reference road surface temperature is the median of the measured road surface temperature of the reference road section under the same solar altitude angle and similar meteorological conditions, and the climate zone type, road surface material and micro-topographic features of the reference road section are the same as those of the target road section; Daytime radiation parameter calibration event: When the daytime phase is in progress, the weather phenomenon is coded as sunny, there is no precipitation, the virtual water film thickness is zero, the wind speed is less than 2m / s and the temperature is higher than 0℃, and the road surface snow or ice begins to melt, the pseudo-true value of the road surface temperature is set to 0℃, and the radiation correction coefficient is updated accordingly. Extreme freezing constraint event: When snow is detected on the road surface and the air temperature rises from no higher than 0°C to higher than 0°C, the pseudo-true value of the road surface temperature is set to 0°C and continues until the snow melt ends. The corresponding daytime atypical meteorological dynamic compensation item or atmospheric window transparency correction operator triggered in the current period is updated, and the extreme freezing constraint event is executed before the daytime radiation parameter calibration event.

8. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 1, characterized in that: The similarity mentioned in step S2 is cosine similarity. The historical road surface temperature inversion values ​​are weighted using meteorological isomorphism weighting, and the steps include: The future prediction period is divided into multiple consecutive future prediction days, and each future prediction day includes 24 hourly prediction times. The cloud cover is divided into four cloud cover levels according to the preset cloud cover threshold, and the wind speed is divided into four wind speed levels according to the preset wind speed threshold. For any given future forecast date, the number of hours each cloud cover level and each wind speed level occurs within 24 hours of that future forecast date is counted and divided by 24 to obtain the time proportion of each level within the future forecast date. The time proportions of the four cloud cover levels and the four wind speed levels are arranged sequentially to form an eight-dimensional forecast weather feature vector for that future forecast date. The same method is used to construct past... N The eight-dimensional measured weather feature vector of each historical day in Tianzhong; Calculate the first The eight-dimensional forecast weather feature vector for the predicted future date and the first... The cosine similarity between the eight-dimensional measured weather feature vectors of the historical days is used as the first... The historical day for the first Meteorological isomorphic weights for each future forecast date ; The expression for the seasonal baseline curve is: ； In the formula, N The total number of historical days collected; The predicted date number; Historical date sequence; It is the hour number within a day, and its value is an integer from 0 to 23; For the first One future prediction day The seasonal baseline curve value at that time; For the first A historical day Historical road surface temperature inversion values ​​at that time.

9. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 8, characterized in that: In step S3, the predicted pavement temperature inversion value and the baseline curve for the predicted time period are fused based on a dynamic fusion weight that decays with the forecast lead time. The expression is as follows: ； ； ； In the formula, Forecast lead time, which is the duration calculated from the start of the forecast; For forecast timeliness Dynamic fusion weights at the location; This is the time constant for predicting the decay characteristics; For meteorological stability operators; For forecast timeliness The final predicted road surface temperature at the location; For forecast timeliness The predicted pavement temperature inversion value at the location; To provide a baseline curve for the forecast period in terms of forecast lead time The value at; The above expression applies to the case where the forecast lead time is greater than zero; where, for Add 1 to the integer part of the quotient obtained by dividing by 24. To predict the starting hour and The remainder when the sum is divided by 24 The number of whole hours within the day corresponding to the predicted starting point, with a value ranging from 0 to 23; When the value is zero, the final predicted road surface temperature is taken as the inversion value of the current road surface temperature.

10. The method for predicting long-term road icing based on meteorological parameter inversion according to claim 1, characterized in that: The existence of a source of moisture on the road surface mentioned in step S4 means that the prediction time corresponding to the final predicted road surface temperature meets at least one of the following conditions: The meteorological forecast data indicates that precipitation will occur at the predicted time. The virtual water film thickness at the predicted time is greater than zero, and the virtual water film thickness is updated hourly according to the meteorological forecast data during the future prediction period. The dew point temperature calculated based on the air temperature and relative humidity at the predicted time in the meteorological forecast data is higher than the final predicted road surface temperature at the predicted time, and the relative humidity at the predicted time is greater than 90%.

Images