A method for early warning of dengue fever import risk in border areas

Abstract

The application discloses a dengue fever import risk early warning method for border areas, relates to infectious disease monitoring and early warning and remote sensing data processing, and comprises the following steps: acquiring monthly average temperature and relative humidity, and calculating local propagation suitability corrected by humidity; extracting annual night light effective pixels in a buffer area on both sides of a national boundary line, and calculating spatial exposure weight according to average brightness on both sides; normalizing monthly cases of neighboring countries, applying a monthly lag, multiplying the spatial exposure weight and the local propagation suitability, and obtaining an import risk score. The early warning level is determined in combination with a local effective reproduction number, and the upgrade takes effect in the current month, and the downgrade falls back gradually by month. The method combines annual night light spatial information with monthly epidemic situation and meteorological changes, and is used for risk early warning of land border areas lacking detailed cross-border population flow data.

Description

Technical Field

[0001] This invention relates to the fields of infectious disease monitoring and early warning and remote sensing data processing technology, specifically to a method for quantifying and classifying the risk of dengue fever importation into border areas using nighttime light remote sensing data, epidemic data from neighboring countries, and local meteorological data. Background Technology

[0002] Dengue fever control in border areas requires a comprehensive consideration of the epidemic situation in neighboring countries, cross-border exposure conditions, and local mosquito-borne transmission environments. For land border areas, detailed border inspection and customs clearance data and cross-border population flow survey data may be difficult to obtain continuously, thus limiting input risk assessment methods that rely on such data. Using only the number of cases in neighboring countries or local meteorological conditions is also insufficient to simultaneously describe the pressure of the epidemic from abroad, the spatial distribution of human activities along the border, and local transmission conditions.

[0003] Annual nighttime light remote sensing imagery can provide information on the spatial distribution of human activity on both sides of the border. This information has a different temporal resolution than monthly epidemic data, necessitating a clear definition of the respective roles of spatial and monthly dynamic information, and the establishment of data processing relationships corresponding to local transmission conditions. This is crucial for developing monthly calculable risk indicators in the absence of detailed cross-border population movement data. Summary of the Invention

[0004] The technical problem to be solved by this invention is how to combine annual nighttime light remote sensing data with monthly epidemic and meteorological data in land border areas where detailed cross-border population flow data is lacking, quantify the monthly changes in spatial exposure differences and input risks, and generate graded early warning results in combination with local transmission status.

[0005] To address the aforementioned technical problems, this invention provides a method for early warning of dengue fever importation risks in border areas, comprising the following steps: S1. Obtain the monthly average temperature and monthly average relative humidity of each spatial unit in the study area. Use the Briere temperature function to determine the mosquito bite rate and external incubation period progress rate. Apply the humidity correction factor to the temperature-dependent mosquito mortality rate. Calculate and normalize the local transmission suitability index of each spatial unit. S2. Acquire annual nighttime light remote sensing images. Delineate strip buffer zones along the national boundary lines corresponding to each spatial unit on the Chinese side and the neighboring country side respectively. Extract effective pixels with radiance greater than zero within each buffer zone and calculate the average brightness on the Chinese side and the average brightness on the neighboring country side respectively. When the study area includes multiple spatial units, calculate the geometric mean of the average brightness on both sides of each spatial unit and perform range normalization among the multiple spatial units to obtain the spatial exposure weight of each spatial unit. When the study area is a single county-level spatial unit, divide the average brightness on the neighboring country side by the sum of the average brightness on both sides to obtain the spatial exposure weight of that spatial unit. S3, obtain the monthly dengue fever case count of neighboring countries, and normalize the monthly dengue fever case count into the epidemic pressure index of neighboring countries; S4. For each spatial unit, multiply the neighboring country's epidemic pressure index (lagging by a preset number of months), the spatial exposure weight of the spatial unit, and the local transmission suitability index for the current month to obtain the input risk score for the current month; the spatial exposure weight remains unchanged within the nighttime light data update interval, and the neighboring country's epidemic pressure index and the local transmission suitability index are updated monthly. S5. Estimate the local effective reproduction number based on the local dengue fever case time series, and jointly determine the input risk score and the local effective reproduction number according to a preset threshold to generate an early warning level.

[0006] Further, in step S1, the local propagation suitability index is obtained by normalizing the original propagation suitability value as shown in the following formula:

[0007] in, The initial values ​​for propagation suitability are determined by temperature and humidity. The average monthly temperature is expressed in °C. The monthly average relative humidity is expressed as a percentage. Mosquito bite rate, Mosquito mortality rate corrected for humidity. The external latency progression rate, and These represent the transmission efficiencies from mosquitoes to humans and from humans to mosquitoes, respectively. The natural death rate of the population. The recovery rate of infected individuals is given; the mosquito bite rate and the external incubation period progression rate are respectively expressed using the Briere function, and the temperature-dependent mosquito mortality rate is the reciprocal of the average adult lifespan in quadratic form.

[0008] In the above formulas, the units for mosquito bite rate, humidity-corrected mosquito mortality rate, external incubation period progression rate, natural mortality rate, and recovery rate of infected individuals are all daily. -1 Both transmission efficiencies and transmission suitability raw values ​​are dimensionless. Subscripts m and v are used to identify mosquito-borne parameters, subscript h to identify human-borne parameters, and vh and hv to identify the transmission directions from mosquito to human and from human to mosquito, respectively. All rates use a consistent time scale, and the mosquito mortality rate is a positive value.

[0009] Furthermore, the humidity-corrected mosquito mortality rate is:

[0010] in, Temperature-dependent mosquito mortality rate, For reference relative humidity, the value is 80%; The humidity sensitivity coefficient ranges from 0.005 to 0.01, and the relative humidity is expressed as a percentage. Mosquito mortality rate corrected for humidity. and These are the monthly average temperature and the monthly average relative humidity, respectively. This represents an exponential function with the natural constant as its base.

[0011] With a relative humidity of 80%, substitute the value 80 into the above formula. The humidity sensitivity coefficient is calculated as a percentage point. -1 The calculation makes the exponential term dimensionless; the mosquito mortality rate before and after humidity correction is in daily units. -1 .

[0012] Furthermore, before normalization, the original value of the propagation suitability is multiplied by a mosquito density correction multiplier; the mosquito density correction multiplier is 1 when there is no mosquito monitoring data, and is the ratio of the current month's Breteau index to the Breteau index reference value when Breteau index monitoring data is available, or is the ratio of the current month's mosquito attractant index to the mosquito attractant index reference value when mosquito attractant oviposition index monitoring data is available.

[0013] The Breteau index reference value can be taken as 5, and the recommended reference value for the mosquito attractant and oviposition index is also 5. The mosquito attractant and oviposition index is substituted as a percentage value, for example, 5% is substituted with the value 5. When calculating the mosquito vector density correction multiplier, the monitored value and the corresponding reference value use the same numerical scale, and the resulting multiplier is dimensionless.

[0014] Furthermore, in step S2, the width of each of the two side strip buffer zones is 50 km; the spatial resolution of the annual nighttime light remote sensing image is 500 m, and the spatial exposure weight is updated every 1 to 3 years based on the annual nighttime light remote sensing image.

[0015] Furthermore, in step S3, the number of dengue fever cases in the neighboring country in the current month is divided by the maximum number of dengue fever cases in the past three years to obtain the epidemic pressure index of the neighboring country; for the missing data in the number of dengue fever cases in the month, linear interpolation is used to fill in the missing data, and the data quality is marked in the early warning output.

[0016] Furthermore, in step S4, the preset number of months is 1 to 3 months.

[0017] Furthermore, the method also includes: setting a rolling window with a length of 4 to 6 months; calculating the Pearson correlation coefficient between the input risk score and the number of local monthly dengue fever cases within the rolling window at each integer month lag from zero to a preset maximum lag order; taking the largest correlation coefficient as the synchronicity indicator; and obtaining the corresponding optimal lag month; wherein the maximum lag order is 0 to 3 months; and simultaneously providing the synchronicity indicator and the optimal lag month in the early warning output results to characterize the temporal correlation between input risk and local cases.

[0018] Furthermore, in step S5, the warning levels include Level I, Level II, Level III, and Level IV, with the risk level decreasing sequentially. The specific determination rules are as follows: When the input risk score is greater than 0.15 and the number of effective local regenerations is greater than 1.5, it is judged as Level I.

[0019] When the input risk score is greater than 0.15 and the local effective regeneration count is greater than 1.0 and less than or equal to 1.5, it is judged as Level II.

[0020] When the input risk score is greater than 0.08 and less than or equal to 0.15, and the local effective regeneration count is greater than 1.0, it is judged as Level III.

[0021] When the input risk score is less than or equal to 0.08, or the local effective reproduction number is less than or equal to 1.0, it is judged as Level IV.

[0022] Furthermore, the level determined according to the preset threshold is used as the condition level for the current month, and the initial warning state is set to level IV; when the risk level of the current month's condition level is higher than the current release level, the level is upgraded to the current month's condition level; when the risk level of the current month's condition level is lower than the current release level, it is reduced by one level each month in the order of level I, level II, level III, and level IV, until the current month's condition level is reached; when the two are the same, the current release level is maintained.

[0023] When the local effective reproduction number is missing, the upgrade will not be triggered for that month, and the gradual reduction will be performed with Level IV as the condition level for that month; when the 95% range of the local effective reproduction number contains 1, the propagation trend will be marked as uncertain in the output.

[0024] The present invention has the following beneficial effects: This invention extracts effective nighttime light pixels within the buffer zones on both the Chinese and neighboring country sides along the national border. The average brightness on both sides is converted into spatial exposure weights, which are then multiplied by a lagged neighboring country epidemic pressure index and a monthly local transmission suitability index. This ensures that the risk score simultaneously incorporates spatial information on human activity on both sides of the border, changes in the overseas epidemic situation, and local transmission conditions. Even in the absence of detailed border inspection data or cross-border population flow survey data, monthly risk calculations can still be completed using remote sensing, case data, and meteorological data, reducing reliance on population flow data that is difficult to obtain continuously.

[0025] This invention sets up weight calculation methods for multiple spatial units and a single county-level spatial unit: for multiple spatial units, the geometric mean of the average brightness on both sides is used and the spatial range is normalized to form a comparable relative exposure weight within the study area; for a single county-level spatial unit, the ratio of the average brightness on the neighboring country side to the sum of the average brightness on both sides is used, without the need for spatial range normalization on a single spatial unit, making the nighttime light weight calculation applicable to the above two study scales.

[0026] This invention uses the Briere temperature function to describe mosquito bite rate and external incubation period progression rate, and adjusts the temperature-dependent mosquito mortality rate using a humidity correction factor, so that transmission suitability can reflect the combined effects of temperature and humidity. When mosquito vector monitoring data is available, the ratio of the Breteau index or the mosquito attractant / oviposition index to the reference value is further used for correction, incorporating local mosquito density changes into the calculation; when mosquito vector monitoring data is unavailable, the correction multiplier is set to 1, and transmission suitability can still be calculated based on meteorological data.

[0027] This invention keeps the spatial exposure weight constant within the nighttime light data update interval and updates the neighboring country epidemic pressure index and local transmission suitability index monthly. This allows annual remote sensing data and monthly dynamic data to participate in risk calculation according to their respective time scales, avoiding the repeated weight extraction of the same annual imagery each month. By setting a 1-3 month lag for the neighboring country epidemic pressure, a clear temporal correspondence can be established between changes in the overseas epidemic and the local transmission conditions in subsequent months.

[0028] This invention combines input risk scores with local effective reproduction numbers for early warning level determination, enabling the output to simultaneously reflect both the level of cross-border input risk and the local transmission status. By employing a release rule that upgrades the warning level monthly when risk increases and gradually reduces it by one level each month when risk decreases, it balances timely response to rising risks with the continuity of the early warning status, limiting the risk of a sudden drop in warning levels due to a single month's decline in indicators. Clear processing and labeling rules are established for cases where effective reproduction numbers are missing or where the 95% interval contains 1, ensuring that the early warning output reflects both the availability of estimated information and the uncertainty of the transmission trend.

[0029] This invention utilizes lag-based rolling window correlation analysis to synchronously output synchronicity indicators and the optimal lag month number according to the warning level. This allows the warning results to simultaneously provide the linear correlation between the input risk and local cases within the selected window, along with their corresponding monthly offsets. This provides a basis for analyzing the temporal correspondence of warning results and calibrating lag parameters. Missing data on neighboring countries' cases are imputed using linear interpolation, and the data quality is indicated in the warning output, enabling users to identify whether the data upon which the results are based has been imputed. Attached Figure Description

[0030] Figure 1 This is a diagram showing the correspondence between the months of the early warning system and the time distribution of the local epidemic situation.

[0031] Figure 2 This is a schematic diagram of the spatial distribution of nighttime lights in the 50 km strip buffer zone along the China-Laos border, as presented in this invention.

[0032] Figure 3 This is a schematic diagram of the data processing and early warning process of the present invention.

[0033] Figure 4 This is a schematic diagram illustrating the nighttime light data processing on both sides of the border and the two spatial weight calculation methods of the present invention. Detailed Implementation

[0034] Example 1: Early warning method for the risk of dengue fever importation in border areas: This embodiment combines Figure 3 and Figure 4 This section describes the data processing procedure of the method. The study area can be a single county-level spatial unit or composed of multiple spatial units distributed along the border. The warning results are calculated and output according to the following steps.

[0035] like Figure 3 As shown, monthly average temperature and relative humidity, annual nighttime light imagery, and monthly cases in neighboring countries are used to calculate local transmission suitability, spatial exposure weight, and lagged epidemic pressure, respectively. These three factors are multiplied to obtain a monthly input risk score. This score is then jointly evaluated with the local effective reproduction number, and escalation and gradual decline rules are applied to output the warning level. Lagged rolling correlation analysis serves as an auxiliary branch, calculating synchronicity indicators and the optimal lag month, which are output synchronously with the warning level for time-series correlation interpretation and lag parameter calibration.

[0036] The study area consists of one or more county-level spatial units, divided into equidistant segments along the border. The time unit is the month. To indicate the month, use This represents a spatial unit. Meteorological data, local case data, and corresponding national border segments are correlated within each spatial unit; the following formulas for transmissibility and synchronicity are calculated for the same spatial unit, omitting spatial subscripts. Parentheses in the formulas denote the independent variables of the function, and subscripts are identifiers rather than multiplication terms. , and These respectively represent the Chinese side, the neighboring country side, and the local side. Indicates a reference value. and These represent the minimum and maximum values, respectively. This indicates that the data has been corrected for mosquito density. This represents nighttime lighting. Calculations involving different spatial units use their respective corresponding data; data involving the same product or ratio use a consistent time and numerical scale.

[0037] Input data includes: monthly average temperature Units are in °C; monthly average relative humidity. Recorded as a percentage; annual nighttime light remote sensing imagery; monthly dengue fever cases in neighboring countries. ; and local case time series for effective reproduction number estimation and synchronicity analysis. For example, 80% relative humidity is calculated by substituting 80 into the humidity correction formula.

[0038] Annual nighttime light imagery can be obtained using NPP-VIIRS-like annual composite imagery, with a spatial resolution of 500 m, a coordinate system of WGS84, and radiance units of nW·cm. -2 ·sr -1 COVID-19 data from neighboring countries can be obtained from the OpenDengue database, reports from the health ministries or disease control agencies of neighboring countries, and WHO regional epidemic reports. The selected data should correspond to monthly case counts.

[0039] Local transmission suitability calculation: The Briere temperature function takes the following form between the lower and upper temperature limits:

[0040] in, For temperature The rate of the corresponding biological process, in days. -1 ; The average monthly temperature is expressed in °C. The temperature response coefficient of the corresponding biological process is given in Celsius (days) on this numerical temperature scale. -1 ·℃ -5 / 2 ; , These represent the lower and upper temperature limits for the biological process, both in °C. Outside this temperature range, the Briere function is set to 0. The mosquito bite rate and the external incubation period progression rate are respectively:

[0041]

[0042] in, For temperature The bite rate of a single mosquito per order. The incubation period progression rate of the virus within mosquitoes is given in days. -1 ; and Temperature response coefficients for the bite rate function and the external latency progression rate function, respectively, in units of... same; and These represent the lower and upper temperature limits of the biting process, respectively. and These represent the lower and upper temperature limits of the external incubation period, respectively, both in °C. (Subscript) Identify the bite process. The incubation period progression outside the identification stage. The bite rate function is taken as follows. , ℃ ℃; External latency progression rate function is taken as , ℃ ℃. Both rates are in daily units.

[0043] The average lifespan of temperature-dependent adults is calculated using a quadratic function, and the mortality rate is calculated within the applicable temperature range where the lifespan function is positive:

[0044]

[0045] in, For temperature The average lifespan of adult mosquitoes, expressed in days; The mosquito mortality rate at this temperature is not corrected for humidity, in days. -1 ; The temperature response coefficient is a function of lifetime, expressed in days·℃. -2 ; and These are the lower and upper temperature limits of the lifetime function, respectively, in °C, with subscripts... Identify the lifespan process. Parameter selection. , ℃, At 25°C, the average lifespan calculated based on the above parameters is approximately 29.84 days, and the temperature-dependent daily mortality rate is approximately 0.0335.

[0046] Relative humidity corrects for mortality rates using the following multipliers:

[0047]

[0048] in, For dimensionless humidity correction multipliers, Represents an exponential function with the natural constant as its base; The average relative humidity for the month. For reference relative humidity, both values ​​are substituted as percentages; The humidity sensitivity coefficient is expressed as a percentage point on this substitution scale. -1 Calculate, make the index It is dimensionless; For temperature and relative humidity Mosquito mortality rate under combined effects, in units of daily -1 Reference relative humidity Take 80%, humidity sensitivity coefficient The value is set to 0.005 to 0.01, with a default value of 0.008. When the relative humidity is equal to the reference value, the multiplier is 1; when it is lower than the reference value, the multiplier is greater than 1; and when it is higher than the reference value, the multiplier is less than 1.

[0049] The original value of propagation suitability is calculated using the following formula:

[0050] in, The values ​​are raw values ​​of transmissibility without mosquito density correction and normalization, and are uniformly adopted at all rates. -1 Furthermore, when the propagation efficiency is dimensionless, this value is dimensionless; , and These are the biting rate, the external incubation period progression rate, and the humidity-corrected mosquito mortality rate, respectively. The efficiency of mosquito transmission to humans The efficiency of human-to-mosquito transmission, both are dimensionless and range from 0.5 to 0.75; subscript and They represent mosquito vectors and humans, respectively. and These indicate the two directions of transmission. Natural mortality rate. Approximately 0.000016 days -1 Recovery rate of infected individuals The value range is 0.14 to 0.25 days. -1 All rates use the same diurnal scale. The squared term of the bite rate corresponds to the two-stage transmission chain of human-to-mosquito and mosquito-to-human bites.

[0051] Table 1. Meaning of each parameter in the P-index calculation formula.

[0052] Wherein, "normalized 0-1" corresponds to this embodiment. The unnormalized raw value is The typical bite rates listed in the table differ from the results calculated using the function coefficients above. Specific calculations should use the functions and coefficients explicitly listed above; for example, the bite rate at 22–28℃ is approximately 0.1635–0.2880 days. -1 .

[0053] When mosquito vector monitoring data is available, the ratio of the Breteau index for the current month to its reference value, or the ratio of the mosquito attractant and oviposition index for the current month to its reference value, can be used as the mosquito vector density correction multiplier. The Breteau index can be taken as a reference value of 5; the mosquito attractant and oviposition index reference value is... A value of 5 is recommended. The ovitrap index, also known as the mosquito attractant index, is recorded and substituted as a percentage, for example, 5% corresponds to the value 5. This recommended value refers to the transmission control threshold of an ovitrap index below 5 in the "Guidelines for Monitoring Aedes Mosquitoes, the Vector of Dengue Fever" issued by the Chinese Center for Disease Control and Prevention, and serves as a benchmark for density correction. When using the ovitrap index, The numerator and denominator use the same numerical scale. When there is no mosquito-borne disease monitoring data, .in, For the first The mosquito ovipositor index for the month This is a fixed reference value for the indicator. For the first The dimensionless mosquito density correction multiplier for the month, subscript Mosquito density correction was applied. The original value for transmission suitability after correction is:

[0054] in, For the first Original value of transmission suitability after correction for mosquito density during menstruation; and The first The average monthly temperature and average monthly relative humidity of the month. The original values ​​for the current month are obtained by substituting the two into the propagation suitability formula. It is dimensionless. Range normalization is applied to the corrected sequence to obtain... :

[0055] in, For the first The dimensionless local transmission suitability index for the month; and The minimum and maximum values ​​of the original corrected propagation suitability values ​​within the reference time series used for normalization are given; the calculation requires the maximum value to be greater than the minimum value. Within this reference series, The value of is between 0 and 1.

[0056] like Figure 4 As shown, using the national border line corresponding to the spatial unit as a reference line, a 50 km wide strip buffer zone is delineated on both the Chinese and neighboring country sides. For annual imagery... , will satisfy the buffer pixels As valid pixels, the sets of valid pixels on the Chinese side and the neighboring country side are denoted as follows: and The average brightness on both sides is:

[0057]

[0058] in, For the annual nighttime light images of the selected year, Identify the year of the image. For cell index, This represents the radiance of that pixel; and These are the sets of valid pixels within the buffer zones on the Chinese side and the neighboring country side, respectively, with set symbols. With propagation suitability function Distinguish by the formula it appears in; and The vertical lines represent the number of valid pixels on both sides, and the vertical lines represent the number of elements in the set. and These represent summation over all valid cells in the corresponding set; and These are the average radiance of the effective pixels on both sides, and their units are... The values ​​are the same, both being nW·cm -2 ·sr -1 The above average values ​​are calculated with the number of effective pixels as the denominator, and the calculation requires that both sets of effective pixels are not empty. For multiple spatial cells, the geometric mean of the average brightness on both sides is calculated separately:

[0059] in, For spatial unit indexing, and The first Average effective pixel radiance of each spatial unit on the Chinese side and the neighboring country side; The value represents the geometric mean of the average brightness on both sides of the unit, with the subscript... The geometric mean is indicated, and its unit remains nW·cm. -2 ·sr -1 Range normalization is performed between the geometric means of all spatial units within the study area:

[0060] in, For the first Dimensionless spatial exposure weight of nighttime light in each spatial unit; and These are the minimum and maximum geometric mean values ​​for all spatial units, respectively, with the maximum value required to be greater than the minimum. The resulting weights reflect the relative position of each spatial unit within the study area.

[0061] When there is only one county-level spatial unit, the spatial exposure weight of that unit is calculated using the following formula:

[0062] in, The dimensionless spatial exposure weight for this single county-level spatial unit. and These are the average effective pixel brightness on the Chinese side and the neighboring country side of the unit, respectively, and their sum should be greater than zero. This weight describes the proportion of the average brightness on the neighboring country side to the sum of the average brightness on both sides. The spatial exposure weight is updated every 1 to 3 years based on annual night light data, and is used as a fixed weight for monthly risk calculations between two updates.

[0063] The neighboring country Monthly case count Divide by the highest monthly number of cases in the past three years Calculate the COVID-19 stress index of neighboring countries:

[0064] in, For the first The dimensionless neighboring country's epidemic pressure index for the month This represents the number of dengue fever cases in neighboring countries that month. The maximum monthly case count within the three-year reference period is used, with case counts as an example. The maximum monthly case count should be a positive value. Linear interpolation is used for missing monthly case data from neighboring countries, and data quality is indicated in the warning results. A lag parameter is set. Its value ranges from 1 to 3 months, with a default value of 1 month.

[0065]

[0066] in, For the first The spatial unit in the first Monthly dimensionless input risk score, subscript Identify input risks; For relative to the first Monthly lag The COVID-19 stress index of neighboring countries over the past month. The lag parameter is a non-negative integer month; in this implementation, it is taken as 1 to 3 months. Assign a nighttime light exposure weight to this spatial unit; For this spatial unit The local transmission suitability index for the month is shown here, with spatial subscripts omitted. When all three factors are in the range of 0 to 1, their product is also in the range of 0 to 1. This product is used to represent the relative score under the combined effects of overseas epidemic pressure, nighttime light exposure weight, and local transmission suitability; when any factor is 0, the score is 0.

[0067] Within the same spatial unit, The length of the scrolling window is one month. Choose 4 to 6 months, with 6 months as the default; The maximum lag order is set to 0–3 months. For each integer month, the lag is... Calculation window and The Pearson correlation coefficient. For simplicity, let's denote it as... , , , Let represent the mean values ​​of the data within these two windows, respectively.

[0068] in, As of the date Months, lagging behind integer months The Pearson correlation coefficient calculated below is dimensionless, and its value ranges from -1 to 1 when the definition holds true; The number of months included in the scrolling window. For the index of the sum of months within the window, from to , with input risk lag parameter Different meanings; Take 0 to 1 in sequence integers, The maximum number of months of lag for comparison. For this spatial unit Monthly input risk score, For the first The number of local cases in a month, for example; and These two groups each contain The arithmetic mean of the data varies with the selected window. Also Changes; Summation symbol This represents the cumulative sum of data for each month within the window. The calculation requires that data for the desired months be available, and that the sum of the squared deviations of both sets of data from the mean within the window be greater than zero, so that the denominator is non-zero.

[0069] The maximum correlation coefficient under each lag is taken as the synchronicity indicator. The lag order that reaches this maximum value is taken as the optimal lag month. The Pearson correlation coefficient reflects the linear co-variance of two sets of data around their respective means, and is not used to determine whether the absolute levels of the number of cases and the risk score are the same. In the early warning output, this synchronicity indicator is included. and optimal lag months Provided in tandem with the corresponding month's alert level to characterize the temporal correlation between imported risk and local cases, and to... The rate determination provides auxiliary information. Among them, As of the date The maximum value among the lagged correlation coefficients compared in the month is dimensionless. To reach this maximum value, an integer month lag is required, in months, with an asterisk indicating the optimal value rather than multiplication. and All are relative to what is already included. Additional offset of the lagged input risk sequence. The synchronicity index characterizes the degree of linear correlation under the selected lag, and the optimal lag month number characterizes the monthly offset corresponding to this correlation.

[0070] The effective reproduction number was estimated using the EpiEstim method based on local case time series. The interval distribution used was a Gamma distribution with a mean of approximately 15 days and a standard deviation of approximately 5.8 days. The estimation window was 30 days, and the prior was Gamma(1,1). The effective reproduction number is presented as the corresponding monthly value. Participation level determination. Among them, For the corresponding number The monthly local effective reproduction number represents the average number of secondary infections caused by each infected person under the then-current transmission conditions. It is dimensionless and has a subscript. The number of valid reproductions is identified; it is related to the input risk score. For different indicators. Representing the gamma distribution, the mean and standard deviation of the interval distribution are in daily terms; prior... The two values ​​in the table represent the prior parameters. Under common shape-scale or shape-rate conventions, these values ​​all correspond to the same distribution with a shape parameter of 1 and a mean of 1. The threshold table below is simplified and omitted. Spatial and monthly independent variables and The independent variable is the month, and all comparisons are made for the same spatial unit and the same month.

[0071] The risk score and the number of effective recurrences will be combined to determine the condition level for the current month, as follows: Table 2. Level IV Warning Table

[0072] Levels I, II, and III correspond to extremely high, high, and moderate alert states, respectively, while Level IV corresponds to routine monitoring. All conditions for Levels I to III in Table 2 must be met simultaneously; for Level IV, only one of the two conditions needs to be met. The risk descriptions and recommended response measures in the table correspond to each alert level.

[0073] The initial warning status is Level IV. When the monthly condition level is higher than the current release level, it will be immediately upgraded to the current month's condition level, allowing for upgrades across levels; when the monthly condition level is lower than the current release level, it will drop by a maximum of one level per month in the order of Level I → Level II → Level III → Level IV, and will remain at the current month's condition level; when the two are the same, the current release level will be maintained.

[0074] When the effective reproduction number cannot be estimated due to reasons such as consecutive zero reports, it is recorded as missing. No upgrade is triggered for that month, and the condition level is set to Level IV for a gradual decline. Therefore, months with missing effective reproduction numbers may still retain a transitional release level due to the previous month's status. If the 95% range of the effective reproduction number contains 1, an uncertain propagation trend is displayed.

[0075] Example 2: Application in Mengla County, Yunnan Province, on the border with Laos: This embodiment uses the method of Example 1, taking Mengla County as a single county-level research unit, to conduct historical calculations from January 2023 to December 2025. Mengla County is located in Xishuangbanna Dai Autonomous Prefecture, Yunnan Province, bordering Laos. During the study period, there were a total of 3466 local cases, including 2368 cases in 2023, 1059 cases in 2024, and 39 cases in 2025; during the same period, there were a total of 63815 cases in Laos. Data content and sources are shown in Table 3.

[0076] Table 3 Data Sources and Explanations for the Examples

[0077] Using Mengla County, Yunnan Province as a single county-level spatial unit, historical calculations were performed for 36 months from January 2023 to December 2025. Cases from neighboring countries were represented by monthly national-level cases from Laos, while local cases were represented by monthly reported cases from Mengla County. Monthly average temperature and monthly average relative humidity were calculated using meteorological data; the mosquito vector density correction multiplier was set to 1.

[0078] like Figure 2 As shown, a 50 km buffer zone was delineated on both sides of the China-Laos border along Mengla County, and the average brightness of effective pixels was extracted from the annual imagery. In the figure, B_cn and B_nb represent the buffer zones on the Chinese and Lao sides, respectively, and nighttime light intensity represents the radiance of the annual remote sensing imagery. The location of the border crossing is used as regional background information and is not an independent input variable for risk scoring. The pixel extraction and weight calculation process within the buffer zones is as follows: Figure 4 As shown. The annual average brightness and weight are as follows; both brightness and weight are values ​​under display resolution.

[0079] Table 4. Nighttime light intensity and spatial exposure weights in the 50 km buffer zone along the China-Laos border (annual)

[0080] In this embodiment, the highest monthly case count of 9988 cases from the three-year data was used as the denominator for the normalization of the epidemic pressure in neighboring countries, with a lag parameter of one month. The original value of transmission suitability was then corrected for mosquito density and normalized for range. The rolling window length for the synchronicity analysis is 6 months, and the maximum lag order is 3 months.

[0081] The complete calculation results for the 36 months from 2023 to 2025 are shown in Table 5. As a local transmission suitability index, The COVID-19 stress index of neighboring countries for that month. To input risk scores, This represents the number of effective local regenerations. (Calculation) Use a 1-month lag Values ​​cannot include the current month's values ​​from the same row in the table. The value is used directly for the monthly rating.

[0082] Table 5 Monthly Calculation Results of the Example

[0083] The "-" in the column indicates that the effective reproduction number for that month is missing, and the "-" in the warning level column indicates that there are no warnings at levels I to III, corresponding to level IV normal monitoring status; the two columns have different meanings. The warning level is the release level after the rule of immediate effect for upgrades and gradual decline for downgrades is implemented.

[0084] The above case numbers are counted per case. The asterisked 10.00 is the displayed value after being truncated to the 0.1–10 range. Adding the listed monthly cases together, the number of cases in Mengla in 2023, 2024, and 2025 were 2368, 1059, and 39 respectively, while the number of cases in Laos were 32084, 19486, and 12245 respectively.

[0085] A Level III warning was issued in June 2023, Level I in July and August, Level II in September, Level III in October, and no Level I-III warnings from November onwards. In August 2024, it was Level III, Level I in September, Level II in October, Level III in November, and no Level I-III warnings from December onwards. No Level I-III warnings were issued throughout 2025. The above time distribution corresponds to Table 5 and... Figure 1 The historical results are shown.

[0086] For example, in August 2023, the lagged epidemic pressure was calculated using 9,988 cases from neighboring countries in July 2023. Substituting a weight of 0.4217 and a monthly transmissibility of 0.95, the input risk score was approximately 0.4006. Based on the effective reproduction number of 1.54 in the table, the condition level for that month was determined to be Level I. In September 2023, the effective reproduction number was 0.71, corresponding to a condition level of Level IV. However, the level dropped one level from Level I in August, and the level issued in September was Level II.

[0087] In September 2024, the lagged epidemic pressure was calculated using 4,487 cases in Laos in August. Substituting the annual weight of 0.3839 and the appropriateness of transmission for the month of 0.93, the input risk score was approximately 0.1604. Based on the effective reproduction number of 1.84 in the table, the condition level for that month was Level I. The above example illustrates the calculation relationship between the risk score and the state transition rules.

[0088] Figure 1 The monthly input risk score, local case count, and warning month of the above embodiment are displayed on the same time axis. The bars represent the number of reported cases in Mengla County in the current month, corresponding to the left axis; the line represents the input risk score, corresponding to the right axis; the background color band represents the month in which a Level I, Level II, or Level III warning was actually issued. The 0.15 and 0.08 horizontal lines in the figure are the risk score thresholds. The issuance level also needs to be determined in combination with the local effective reproduction number and the step-down rule, and cannot be determined solely by whether the line exceeds the horizontal line.

[0089] Historical data was used to calibrate the humidity sensitivity coefficient, the lag month for epidemic pressure, the rolling window length, and the risk classification threshold. The candidate range for the humidity sensitivity coefficient was 0.005–0.01, the lag month was 1–3 months, and the rolling window length was 4–6 months.

[0090] The validation criteria included the time deviation between the warning month and the actual outbreak month, the Spearman rank correlation between the peak risk score and the peak incidence, classification sensitivity and specificity, the area under the receiver operating characteristic curve, and the number of months of advance warning. An outbreak was defined as a monthly incidence rate exceeding the 90th percentile of the past three years, and a one-year retention cross-validation method was used: parameters were calibrated using the remaining years each time, and validation was performed using the retained year, rotating year by year.

Claims

1. A method for early warning of dengue fever importation risk in border areas, characterized in that, Includes the following steps: S1. Obtain the monthly average temperature and monthly average relative humidity of each spatial unit in the study area. Use the Briere temperature function to determine the mosquito bite rate and external incubation period progress rate. Apply the humidity correction factor to the temperature-dependent mosquito mortality rate. Calculate and normalize the local transmission suitability index of each spatial unit. S2. Acquire annual nighttime light remote sensing images. Delineate strip buffer zones along the national boundary lines corresponding to each spatial unit on the Chinese side and the neighboring country side respectively. Extract effective pixels with radiance greater than zero in each buffer zone and calculate the average brightness on the Chinese side and the average brightness on the neighboring country side respectively. When the study area includes multiple spatial units, calculate the geometric mean of the average brightness on both sides of each spatial unit and perform range normalization among the multiple spatial units to obtain the spatial exposure weight of each spatial unit. When the study area is a single county-level spatial unit, the spatial exposure weight of the spatial unit is obtained by dividing the average brightness of the neighboring country side by the sum of the average brightness of both sides. S3, obtain the monthly dengue fever case count of neighboring countries, and normalize the monthly dengue fever case count into the epidemic pressure index of neighboring countries; S4. For each spatial unit, multiply the neighboring country's epidemic pressure index (lagging by a preset number of months), the spatial exposure weight of the spatial unit, and the local transmission suitability index for the current month to obtain the input risk score for the current month; the spatial exposure weight remains unchanged within the nighttime light data update interval, and the neighboring country's epidemic pressure index and the local transmission suitability index are updated monthly. S5. Estimate the local effective reproduction number based on the local dengue fever case time series, and jointly determine the input risk score and the local effective reproduction number according to a preset threshold to generate an early warning level.

2. The method according to claim 1, characterized in that, In step S1, the local propagation suitability index is obtained by normalizing the original propagation suitability value as shown in the following formula: in, The initial values ​​for propagation suitability are determined by temperature and humidity. The average monthly temperature is expressed in °C. The monthly average relative humidity is expressed as a percentage. Mosquito bite rate Mosquito mortality rate corrected for humidity. The external latency progression rate, and These represent the transmission efficiencies from mosquitoes to humans and from humans to mosquitoes, respectively. The natural death rate of the population. The recovery rate of infected individuals is given; the mosquito bite rate and the external incubation period progression rate are respectively expressed using the Briere function, and the temperature-dependent mosquito mortality rate is the reciprocal of the average adult lifespan in quadratic form.

3. The method according to claim 2, characterized in that, The humidity-corrected mosquito mortality rate is: in, Temperature-dependent mosquito mortality rate For reference relative humidity, the value is 80%; The humidity sensitivity coefficient ranges from 0.005 to 0.01, and the relative humidity is expressed as a percentage. Mosquito mortality rate corrected for humidity. and These are the monthly average temperature and the monthly average relative humidity, respectively. This represents an exponential function with the natural constant as its base.

4. The method according to claim 2, characterized in that, Before normalization, the original value of the propagation suitability is multiplied by a mosquito density correction multiplier; the mosquito density correction multiplier is 1 when there is no mosquito monitoring data, and is the ratio of the current month's Breteau index to the Breteau index reference value when Breteau index monitoring data is available, or the ratio of the current month's mosquito attractant index to the mosquito attractant index reference value when mosquito attractant oviposition index monitoring data is available.

5. The method according to claim 1, characterized in that, In step S2, the width of the two side strip buffer zones is 50 km each; the spatial resolution of the annual nighttime light remote sensing image is 500 m, and the spatial exposure weight is updated every 1 to 3 years based on the annual nighttime light remote sensing image.

6. The method according to claim 1, characterized in that, In step S3, the number of dengue fever cases in the neighboring country in the current month is divided by the maximum number of monthly dengue fever cases in the past three years to obtain the epidemic pressure index of the neighboring country; for the missing data in the number of monthly dengue fever cases, linear interpolation is used to fill in the missing data, and the data quality is marked in the early warning output.

7. The method according to claim 1, characterized in that, In step S4, the preset number of months is 1 to 3 months.

8. The method according to claim 1, characterized in that, Also includes: A rolling window with a length of 4 to 6 months is set. At each integer month lag from zero to a preset maximum lag order, the Pearson correlation coefficient between the lagging input risk score and the local monthly dengue fever case count within the rolling window is calculated. The largest correlation coefficient is taken as the synchronicity indicator, and the corresponding optimal lag month is obtained to analyze the temporal correlation between input risk and local cases. The maximum lag order is 0 to 3 months. The synchronicity indicator and the optimal lag month are simultaneously provided in the early warning output results to characterize the temporal correlation between input risk and local cases.

9. The method according to claim 1, characterized in that, In step S5, the warning levels include Level I, Level II, Level III, and Level IV, with the risk level decreasing sequentially. When the input risk score is greater than 0.15 and the local effective recurrence count is greater than 1.5, it is determined to be Level I; when the input risk score is greater than 0.15 and the local effective recurrence count is greater than 1.0 and less than or equal to 1.5, it is determined to be Level II; when the input risk score is greater than 0.08 and less than or equal to 0.15, and the local effective recurrence count is greater than 1.0, it is determined to be Level III; when the input risk score is less than or equal to 0.08, or the local effective recurrence count is less than or equal to 1.0, it is determined to be Level IV.

10. The method according to claim 9, characterized in that, The level determined according to the preset threshold will be used as the condition level for the current month, and the initial warning status will be set to Level IV; when the risk level of the current month's condition level is higher than the currently issued level, the current month's condition level will be upgraded to the current month's condition level. When the risk level of the current month's condition level is lower than the current release level, it will be reduced by one level each month in the order of Level I, Level II, Level III, and Level IV, until the current month's condition level is reached; When both are the same, maintain the current release level; when the local effective reproduction number is missing, the upgrade will not be triggered for the month, and the gradual reduction will be performed with Level IV as the condition level for the month; when the 95% range of the local effective reproduction number contains 1, the propagation trend will be marked as uncertain in the output.

Images