Lake surface temperature reconstruction and freeze thawing phenology extraction method based on filtering-automatic identification

Abstract

The invention discloses a lake surface temperature reconstruction and freeze-thaw phenology extraction method based on filtering-automatic identification. The method comprises the steps of obtaining thermal infrared remote sensing data of a target lake in historical time; dividing a freeze-thaw period and a non-freeze-thaw period of the target lake; for a non-freeze-thaw period, acquiring a lake surface temperature time sequence of the non-freeze-thaw period; aiming at the freezing and thawing period, performing nonlinear curve fitting, and determining a freezing starting day, a freezing ending day, an ablation starting day and an ablation ending day of the target lake; according to the method, the continuity and precision of the lake surface temperature time sequence in the non-freezing-thawing period are improved, so that the error of traditional data caused by cold deviation is effectively reduced; and four accurate phenological key dates, namely freezing start and end and ablation start and end, can be scientifically, automatically and quantitatively extracted according to the freezing and thawing period, so that the deviation and low efficiency of subjective visual interpretation are avoided.

Description

Technical Field

[0001] This invention relates to the field of automated inversion technology of lake environmental parameters, specifically to a method for reconstructing lake surface temperature and extracting freeze-thaw phenology based on filtering and automatic identification. Background Technology

[0002] With the rapid development of remote sensing satellites, there are now abundant remote sensing data sources available for monitoring lake temperature and ice condition changes. However, the application of different data types has certain limitations. Some satellite remote sensing data, such as the Landsat series, have high spatial resolution, but their revisit periods are long (e.g., Landsat / TM has a revisit period of 16 days) and are easily affected by atmospheric conditions. Some high temporal resolution sensors, such as NOAA Advanced Very High Resolution (AVHRR), have shorter revisit periods, but lower spatial resolution, which cannot meet the accuracy requirements for monitoring lake changes in smaller lakes or lakes located in mountainous areas. Some satellite products, such as MODIS (Moderate-resolution Imaging Spectroradiometer), can achieve a temporal resolution of four times per day and a spatial resolution of 500 / 1000 meters, which can meet the temporal and spatial resolution requirements for lake surface temperature and ice condition monitoring, and are suitable for retrieving lake temperature changes and ice condition dynamics.

[0003] However, due to the difference between the depth observed by satellite and the depth measured in reality, there is a "cold bias" in the satellite-observed lake surface temperature values. Specifically, the satellite-observed lake surface temperature is considered to be the temperature of a very thin layer of water at a depth of 0-0.2m, while the measured lake surface temperature is generally the average temperature of the water layer at a depth of 0-0.5m. Therefore, there is a negative systematic error between the satellite-observed and measured water temperature values, called the "cold bias value." In addition, the cold bias may also be caused by cloud contamination or undetected clouds in the satellite imagery, which may lead to subsequent algorithms mistakenly identifying the temperature of unidentified cloud tops as the lake surface temperature. The error caused by the "cold bias" can generally reach 0.6-1˚C.

[0004] Lake surface temperature series are currently the primary data source for extracting freeze-thaw phenological information about lakes. Current research often employs the empirical threshold method, which uses remote sensing data to derive the lake surface temperature time series and fixes a temperature threshold (e.g., 0°C). The freezing and thawing states of the lake are determined by whether the daily lake surface temperature is above or below this threshold. When the temperature remains below the threshold, the lake is considered to be in a freezing period; when the temperature rises above the threshold, the lake is considered to be in a thawing period. However, due to differences in the physicochemical properties of lakes and varying environmental conditions in the study area, the threshold has certain specificities. Therefore, if the selected threshold is inappropriate, it can lead to errors in determining the time points of lake ice thawing and freezing periods. Moreover, when remote sensing products are affected by cloud cover, the availability of data during the freezing or thawing periods is limited, making it impossible to accurately determine the time points using the threshold method.

[0005] Therefore, there is currently a lack of lake surface temperature sequences that can reduce "cold bias" and have high accuracy, as well as robust methods for automatically extracting freeze-thaw phenological parameters. Summary of the Invention

[0006] To address the aforementioned issues, this invention provides a method for reconstructing lake surface temperature and extracting freeze-thaw phenology based on filtering and automatic identification.

[0007] A method for reconstructing lake surface temperature and extracting freeze-thaw phenology based on filtering and automatic identification includes the following steps: Acquire historical thermal infrared remote sensing data of the target lake; the thermal infrared remote sensing data includes surface temperature observations of the target lake at least twice a day; the daily surface temperature observations are synthesized to obtain daily observation values ​​and daily synthesized temperature images of the target lake; The freeze-thaw period and non-freeze-thaw period of the target lake are defined; For the non-freeze-thaw period, the original temperature sequence is first calculated based on the daily observations and the daily synthetic temperature image of the target lake. Then, the original temperature sequence is corrected by a filtering algorithm to obtain the lake surface temperature time series during the non-freeze-thaw period. For the freeze-thaw period, a temperature threshold is set, the state of water or ice is determined on the daily synthetic temperature image, and the percentage of lake ice area is calculated daily to form a time series of lake ice area percentage. Nonlinear curve fitting is performed on the freezing and melting stages in the time series of lake ice area percentage. Based on the fitting results, the freezing start date, freezing end date, melting start date, and melting end date of the target lake are determined.

[0008] Note: The above method significantly improves the continuity and accuracy of lake surface temperature time series during the non-freeze-thaw period by using daily multi-temporal thermal infrared data and applying filtering correction, effectively mitigating the errors caused by "cold bias" in traditional data. By using temperature threshold identification and area percentage curve fitting during the freeze-thaw period, the four precise phenological key dates of the start and end of freezing, the start and end of thawing, and the end of thawing can be scientifically and automatically quantified and extracted, avoiding the bias and inefficiency of subjective visual interpretation.

[0009] Furthermore, the thermal infrared remote sensing data includes the surface temperature of the target lake, obtained four times daily via MODIS.

[0010] Note: Utilizing MODIS's intensive thermal infrared observations four times a day can improve the temporal resolution and redundancy of the data; it can also capture the diurnal variation characteristics of lake surface temperature more precisely, providing rich and reliable raw input for subsequent temperature series reconstruction.

[0011] Furthermore, the method for the synthesis process is as follows: Daily surface temperature observations of the target lake include Four; First, the temperature value of the effective pixel in the maximum composite area is calculated using the following formula (1), and then the daily observation value is calculated using formula (2). (1) (2) In equation (1), This represents the temperature value of the effective pixels within the maximum daily composite area. Let be the discriminant function for the effective temperature pixels of the target lake, when hour Take 1, when Invalid Set to 0; For surface temperature observation The total number of effective pixels; in equation (2), For the first Observed values ​​for the day.

[0012] Note: The above synthetic processing method significantly improves the data quality and reliability of daily observations by intelligently identifying and aggregating effective temperature pixels across multiple time periods and setting an effective data volume threshold; and ensures reliable data output every day by maximizing the complementarity of multiple time periods.

[0013] Furthermore, the freeze-thaw period of the target lake is from November to April of the following year; the non-freeze-thaw period is from May to October of the year.

[0014] Note: The above settings clearly define the freeze-thaw period from November to April of the following year, and the non-freeze-thaw period from May to October. This provides a clear and stable time frame for subsequent processing procedures. These settings are reference values ​​for the lake area studied, and the timeframes may differ in other lakes.

[0015] Further, based on daily observations and daily synthetic temperature images of the target lake, the original temperature sequence is calculated, and a filtering algorithm is used to correct the original temperature sequence to obtain a time series of lake surface temperature during the non-freeze-thaw period; including: Based on the daily synthetic temperature images, the average temperature within the target lake area is calculated and arranged in chronological order to form an original temperature sequence; The original temperature sequence was corrected using an upper profile filtering algorithm to obtain a time series of lake surface temperature during the non-freeze-thaw period.

[0016] Furthermore, the upper profile filtering algorithm is shown in the following formula (3): (3) In equation (3), For the revised first Lake surface temperature during non-freeze-thaw periods For the first Observed values ​​for the day.

[0017] Explanation: The above steps calculate the spatial average temperature of the lake to form the original sequence, and then apply an upper profile filtering algorithm for correction. This effectively identifies and suppresses abnormally low-value noise caused by residual cloud pollution or atmospheric disturbances, while accurately preserving the true temperature change trend and peak characteristics of the lake. Compared with traditional smoothing algorithms, this method is specifically optimized for the error characteristic of "negative bias" which is more prevalent in non-freeze-thaw temperature observations, thereby generating a more reliable and high-quality continuous lake surface temperature time series, laying a solid data foundation for subsequent phenological analysis and climate change research.

[0018] Furthermore, the temperature threshold is -0.75 to 0.75. When the daily observed temperature is below -0.75℃, the state is determined to be ice; when the daily observed temperature is above 0.75℃, the state is determined to be water.

[0019] Note: In certain embodiments of this invention, a specific temperature threshold range of -0.75℃ to 0.75℃ is used for automatic water ice state identification. This effectively avoids misjudgments caused by sensor noise or daytime temperature fluctuations near the freezing point (such as misjudging a brief period of low temperature as freezing). It also ensures that the daily water ice identification results are stable and repeatable, providing a reliable and physically meaningful basis for subsequent accurate calculation of lake ice area percentage and subsequent fitting and extraction of key phenological dates. In other embodiments with different regions and environments, this threshold may change, and this invention does not limit it.

[0020] Furthermore, the temperature threshold is determined based on a time series of lake surface temperatures during the non-freeze-thaw period; the determination method includes: The time series of lake surface temperatures during the non-freeze-thaw period is divided into periods of autumn cooling trend and spring warming trend. By fitting, a trend line of temperature change over time is obtained during the autumn cooling trend period; combined with the trend line, the cooling reference date is determined; and the lower quantile value of the lake surface temperature within 5 to 15 days before the cooling reference date is calculated, and the lower quantile value is used as the lower limit of the temperature threshold; wherein, the cooling reference date is the date when the temperature enters a continuous low temperature state after the theoretical freezing point. Determine the turning point date when the temperature begins to rise continuously during the spring warming trend period, calculate the high quantile value of the lake surface temperature within 5 to 15 days after the turning point date, and use the high quantile value as the upper limit of the temperature threshold.

[0021] Note: The above method abandons fixed empirical values ​​and instead performs dynamic and objective calculations based on the target lake's own non-freeze-thaw temperature sequence. By identifying the turning points in autumn and spring temperature trends and extracting the low and high quantiles of the temperature distribution near these turning points as upper and lower limits for thresholds, the final determined thresholds closely reflect the lake's actual phase transition temperature characteristics and interannual climate fluctuations. This data-driven approach significantly improves the specificity and accuracy of water ice state identification, effectively reduces systematic errors that may arise from using uniform thresholds, and ensures the comparability and reliability of phenological extraction results across different years.

[0022] Furthermore, the fitting formula for the nonlinear curve fitting of the freezing stage in the lake ice area percentage time series is shown in equation (4), and the fitting formula for the nonlinear curve fitting of the melting stage is shown in equation (5). (4) In equation (4), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. t represents the percentage of lake ice area in the fitted curve, and t is time. (5).

[0023] In equation (5), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. t represents the percentage of lake ice area in the fitted curve, and t represents time.

[0024] Note: The nonlinear curve fitting of the lake ice area percentage time series described above can accurately reflect the phenological evolution during the freezing and melting of lake ice. This effectively overcomes the uncertainty of subjective interpretation based solely on raw discrete data points, smooths out observation noise, and ultimately ensures that the extracted phenological indicators are accurate, reliable, and have good time series comparability.

[0025] Furthermore, the calculation formulas for the freezing start date, freezing end date, thawing start date, and thawing end date of the target lake are shown in the following equations (6), (7), (8), and (9): Freeze start date: (6) Freeze End Date: (7) Ablation start date: (8) Ablation end date: (9) In equations (6) and (7), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. The date of the freeze start date, The freeze ends on the due date; In equations (8) and (9), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. The end date of ablation The date of completion of the ablation process.

[0026] Note: The above formula provides a unified and quantitative calculation standard, eliminating the subjective arbitrariness of manual interpretation. This ensures a high degree of repeatability and interannual comparability of the extracted results, and realizes a fully automated process from remote sensing data to final phenological indicators, greatly improving the efficiency and scientific rigor of long-term, large-scale lake freeze-thaw monitoring.

[0027] The beneficial effects of this invention are: This invention significantly improves the continuity and accuracy of lake surface temperature time series during non-freeze-thaw periods by filtering and correcting multi-temporal thermal infrared data, effectively mitigating errors caused by "cold bias" in traditional data. By using temperature threshold identification and area percentage curve fitting during the freeze-thaw period, the four precise phenological key dates of freezing start, end, thawing start, and end can be scientifically and automatically quantified and extracted, avoiding the bias and inefficiency of subjective visual interpretation. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the method flow of Embodiment 1 of the present invention; Figure 2 This is a schematic diagram of the method for determining the lower limit of the temperature threshold in Embodiment 2 of the present invention; Figure 3 This is a schematic diagram of the method for determining the upper limit of the temperature threshold in Embodiment 2 of the present invention; Figure 4 This is a daily lake surface temperature time series data plot reconstructed based on MODIS data in Embodiment 1 of the present invention; Figure 5 This is a comparative analysis of MODIS lake surface temperature data and measured data in Embodiment 1 of the present invention; Figure 6 This is a curve fitting diagram of the freezing and thawing periods of a lake in some years in Embodiment 1 of the present invention; Figure 7 This is the verification result of the lake ice / lake water area percentage during the freezing and melting periods of a certain lake in Embodiment 1 of the present invention. Detailed Implementation

[0029] To further illustrate the methods and effects of this invention, the technical solution of this invention will be clearly and completely described below in conjunction with experiments.

[0030] In light of the background technology, it should be understood that remote sensing surface temperature data refers to the spatial distribution data reflecting the surface temperature of land or water bodies obtained by processing the surface thermal radiation information received by satellite sensors through inversion algorithms. It mainly includes daytime and nighttime temperature observations and their corresponding time, location, and quality control information.

[0031] This invention's research revealed that during the non-freeze-thaw period, lake surfaces are open bodies of water with continuous temperature variations. The main problem with data from this period is insufficient accuracy: satellite-retrieved temperatures exhibit a systematic "cold bias," and residual clouds cause noise that leads to sudden temperature drops. Therefore, the processing objective must be to correct these biases to obtain accurate absolute temperature values; the upper profile filtering method is precisely designed to address this goal.

[0032] During the freeze-thaw cycle, the lake surface undergoes a phase transition from ice to water, with temperatures hovering around the freezing point and fluctuating dramatically. The main problems with data from this period are poor availability and misjudgment of state: cloud cover often leads to missing data for key dates, and the spectral and temperature characteristics of thin and broken ice on the lake surface are blurred. Therefore, the processing objective is no longer to pursue absolute accuracy for individual temperature values, but rather to reliably capture the entire lake surface transition process from water to ice or from ice to water and its timing. The method of using thresholding for initial judgment followed by curve fitting aims to robustly extract these transition moments from incomplete and noisy data.

[0033] Example 1: Based on the above, this invention proposes a method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification, including the following steps: like Figure 1 As shown, S101, acquire the thermal infrared remote sensing data of the target lake over historical time; the thermal infrared remote sensing data includes the surface temperature of the target lake observed at least twice a day; the surface temperature observed each day is synthesized to obtain the daily observation value and the daily synthesized temperature image of the target lake. The thermal infrared remote sensing data includes the surface temperature of the target lake, obtained four times a day via MODIS.

[0034] The method for the synthesis process is as follows: Daily surface temperature observations of the target lake include Four; First, the temperature value of the effective pixel in the maximum composite area is calculated using the following formula (1), and then the daily observation value is calculated using formula (2). (1) (2) In equation (1), This represents the temperature value of the effective pixels within the maximum daily composite area. Let be the discriminant function for the effective temperature pixels of the target lake, when hour Take 1, when Invalid Set to 0; For surface temperature observation The total number of effective pixels; in equation (2), For the first Observed values ​​for the day.

[0035] For example, this embodiment first downloads the MODIS LST products from the MODIS website, specifically the Terra and Aqua satellites, namely MOD11A1 and MYD11A1. The spatial resolution is 1 km, the temporal resolution is four times daily, and the version is V006. The MOD11A1 data comes from the Terra satellite, providing observation data twice daily, approximately at 10:30 AM and 10:30 PM. The MYD11A1 data comes from the Aqua satellite, providing daytime and nighttime observation data, approximately at 1:30 PM and 1:30 AM respectively. The remote sensing image is then flooded using the target lake's vector, and after cropping the lake's shape, a 33-pixel buffer zone is selected to remove mixed land and lake pixels. Then, based on the image's quality control file, pixels marked as "good quality" for each observation are selected. Based on the selected effective temperature data, the number of effective temperature pixels, the number of temperature pixels within the mask, the percentage of effective temperature pixels, and the percentage of effective pixels are calculated. According to the statistical results, the temperature value where the percentage of effective temperature pixels is less than 50% is set as a NaN value. Furthermore, by calculating the average of multiple effective observations at the same pixel location, pixel-by-pixel fusion processing is performed to ultimately generate a daily composite temperature image of the lake area with minimized cloud pollution and maximized spatial coverage. Based on this image, a single temperature observation value representing the daily thermal state of the entire lake is calculated. Determining the pixel state of lake ice and water using daily maximum area composite images: During the freezing period, pixels with temperatures below -0.75°C are labeled as lake ice; during the melting period, pixels with temperatures above 0.75°C are labeled as lake water, and otherwise are labeled as lake water pixels.

[0036] Calculate the percentage of effective pixels and percentage of lake ice area The calculation formula is: When the percentage of valid pixels is below 30% and the percentage of pixels marked as lake water within that day's range is 0, it is marked as an outlier and removed. Next, the values ​​of valid pixels from the three days before and after the outlier, as well as from fixed dates throughout the year, are used to imput missing data within the existing lake ice percentage sequence. Furthermore, the freezing and melting phases of lake ice are reassessed using dates with autumn and spring temperatures of 0˚C: the percentage of lake ice area from September in autumn to the date with autumn temperature of 0˚C is set to 0%, while the percentage of lake ice area from February in late winter to the date with spring temperature of 0˚C is set to 100%.

[0037] The results are shown in Table 1; Table 1 compares the effective pixels after synthesis with those from four images acquired by Terra and Aqua. .

[0039] S102. Determine the freeze-thaw period and non-freeze-thaw period of the target lake; The freeze-thaw period of the target lake is from November to April of the following year; the non-freeze-thaw period is from May to October of the year.

[0040] S103. For the non-freeze-thaw period, firstly, based on the daily observations and the daily synthetic temperature image of the target lake, the original temperature sequence is calculated, and the original temperature sequence is corrected by a filtering algorithm to obtain the lake surface temperature time series during the non-freeze-thaw period. The process involves calculating the original temperature sequence based on daily observations and daily composite temperature images of the target lake, and then applying a filtering algorithm to correct the original temperature sequence to obtain a time series of lake surface temperature during the non-freeze-thaw period; including: Based on the daily synthetic temperature images, the average temperature within the target lake area is calculated and arranged in chronological order to form an original temperature sequence; like Figure 4 , Figure 5 As shown, the original temperature sequence is corrected using an upper profile filtering algorithm to obtain a non-freeze-thaw lake surface temperature time series; in addition, in some other embodiments, other filtering methods can be used for processing, and the present invention does not impose strict limitations on this. The upper profile filtering algorithm is shown in the following formula (3): (3) In equation (3), For the revised first Lake surface temperature during non-freeze-thaw periods For the first Observed values ​​for the day.

[0041] For example, taking the data processing of a certain lake from 2002 to 2016 as an example, after identifying its non-freeze-thaw period, the average temperature of all lake pixels is calculated daily based on the synthetic temperature images generated daily during this period, and arranged in chronological order to construct the original daily average lake surface temperature sequence. Subsequently, the original sequence is automatically corrected using an upper profile filtering algorithm. This algorithm compares the daily temperature with the weighted average of the temperatures of the days before and after it, and selects the larger value as the output, thereby effectively suppressing the abnormal low temperature noise caused by residual thin clouds, and finally obtaining a continuous and high-quality non-freeze-thaw period lake surface temperature time series.

[0042] S104. For the freeze-thaw period, a temperature threshold is set, the state of water or ice is determined on the daily synthetic temperature image, and the percentage of lake ice area is calculated daily to form a time series of lake ice area percentage. Nonlinear curve fitting is performed on the freezing and melting stages in the time series of lake ice area percentage. Based on the fitting results, the freezing start date, freezing end date, melting start date, and melting end date of the target lake are determined.

[0043] In this embodiment of the invention, the temperature threshold is set to -0.75~0.75 based on existing experience. That is, when the temperature of the daily observation is below -0.75℃, the state is determined to be ice; when the temperature of the daily observation is above 0.75℃, the state is determined to be water.

[0044] like Figure 6 As shown, the fitting formula for the nonlinear curve fitting of the freezing stage in the lake ice area percentage time series is shown in equation (4), and the fitting formula for the nonlinear curve fitting of the melting stage is shown in equation (5). (4) In equation (4), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. t represents the percentage of lake ice area in the fitted curve, and t is time. (5) In equation (5), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. t represents the percentage of lake ice area in the fitted curve, and t represents time.

[0045] The calculation formulas for the freezing start date, freezing end date, thawing start date, and thawing end date of the target lake are shown in equations (6), (7), (8), and (9) below: Freeze start date: (6) Freeze End Date: (7) Ablation start date: (8) Ablation end date: (9); In equations (6) and (7), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. The date of the freeze start date, The freeze ends on the due date; In equations (8) and (9), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. The end date of ablation The date of completion of ablation; like Figure 7 As shown, by verifying the percentage of lake ice / lake water area during the freezing and melting periods of the lake, it can be proven that this embodiment has good effects.

[0046] Example 2: This example is largely the same as Example 1, except that the temperature threshold is set differently; like Figure 2 and Figure 3 As shown, the temperature threshold is determined based on the time series of lake surface temperature during the non-freeze-thaw period; the determination method includes: The time series of lake surface temperatures during the non-freeze-thaw period is divided into periods of autumn cooling trend and spring warming trend. By fitting, a trend line of temperature change over time is obtained during the autumn cooling trend period; combined with the trend line, the cooling reference date is determined; and the lower quantile value of the lake surface temperature within 5 to 15 days before the cooling reference date is calculated, and the lower quantile value is used as the lower limit of the temperature threshold; wherein, the cooling reference date is the date when the temperature enters a continuous low temperature state after the theoretical freezing point. Determine the turning point date when the temperature starts to rise continuously during the spring warming trend period, calculate the high quantile value of the lake surface temperature within 5 to 15 days after the turning point date, and use the high quantile value as the upper limit of the temperature threshold. The process of determining the lower limit of the temperature threshold specifically includes: First, extract the autumn window (e.g., from September 1 to November 15) from the reconstructed 2022 non-freeze-thaw temperature series; perform linear regression fitting on the temperature data points (date, temperature) within this window to obtain the cooling trend line; Second, determine the theoretical freezing point: According to historical observations of a certain lake, the physical freezing temperature at which supercooled water and ice crystals begin to form in its freshwater portion is about -0.3°C, rather than the ideal 0°C. In other lakes, this temperature may fluctuate. This example is not limited to comparison. Third, calculate the intersection of the cooling trend line and the -0.3°C physical freezing temperature line. Assuming the calculated date is December 3, 2022, take the theoretical freezing point, i.e., December 3, as the starting point and examine the actual temperature sequence. We look for the first day when the daily average temperature remains consistently below -0.3°C for five consecutive days without a significant rebound (e.g., a single day's temperature increase does not exceed 1°C); this starting day is taken as the cooling baseline date. Fourth, backtrack a short window, such as 10 days, from the cooling baseline date, and extract the actual observed lake surface temperature values ​​during these 10 days to form a dataset containing 10 data points; calculate the 10th percentile of this dataset. Assume the temperatures for these 10 days are: "2.1℃, 1.5℃, 0.8℃, 0.3℃, -0.2℃, -0.5℃, -0.7℃, 1.0℃, -0.1℃, 0.5℃".

[0047] The 10th percentile is -0.7°C, which is the lower limit of the temperature threshold. Specifically, it means that the lake's temperature is below -0.7°C for 10% of the time before freezing. During the freeze-thaw cycle processing of that year, any pixel with a temperature below -0.7°C will be classified as ice. It should be understood that the above lower percentile value uses a calculation method known in existing technology, and the calculation process will not be elaborated here. The process of determining the upper limit of the temperature threshold specifically includes: First, take the temperature series for spring (e.g., March 1st to May 15th); apply a change point detection algorithm (such as Bayesian change point detection, Pettitt test, etc.) to this series. The algorithm scans the series to find dates where statistical characteristics change significantly, and after these dates (the turning point date when the temperature begins a sustained rise), the mean or upward trend of the temperature series changes significantly, shifting from a "fluctuating state" to a "sustainable rising state"; assuming the algorithm detects the turning point date as April 10th, 2023; Secondly, taking the inflection point date as the baseline, a short window is selected, such as 10 days, i.e., April 11 to April 20; the actual observed lake surface temperature values ​​during these 10 days are extracted to form a dataset, and the 90th percentile of this dataset is calculated; assuming the temperatures for these 10 days are: 0.5℃, 1.0℃, 1.8℃, 2.5℃, 3.0℃, 1.2℃, 4.0℃, 2.2℃, 3.5℃, 5.0℃.

[0048] Arranging these temperatures from smallest to largest, the 90th percentile is 0.9 * (10 + 1) = 9.9. Interpolating between the ninth and tenth data points, we get: 4.0 + 0.9(5.0 - 4.0) = 4.9℃. This means the upper limit of the temperature threshold is 4.9℃; specifically, during the initial stage of stable warming, the lake's temperature is below 4.9℃ 90% of the time. During the freeze-thaw cycle processing that year, any pixel with a temperature above 4.9°C will be classified as water.

Claims

1. A method for reconstructing lake surface temperature and extracting freeze-thaw phenology based on filtering and automatic identification, characterized in that, Includes the following steps: Acquire historical thermal infrared remote sensing data of the target lake; the thermal infrared remote sensing data includes surface temperature observations of the target lake at least twice a day; the daily surface temperature observations are synthesized to obtain daily observation values ​​and daily synthesized temperature images of the target lake; The freeze-thaw period and non-freeze-thaw period of the target lake are defined; For the non-freeze-thaw period, the original temperature sequence is first calculated based on the daily observations and the daily synthetic temperature image of the target lake. Then, the original temperature sequence is corrected by a filtering algorithm to obtain the lake surface temperature time series during the non-freeze-thaw period. For the freeze-thaw period, a temperature threshold is set, and the state of water or ice is determined by the daily synthetic temperature image based on the temperature threshold. Then, the percentage of lake ice area is calculated daily to form a time series of lake ice area percentage. Nonlinear curve fitting is performed on the freezing and melting stages in the time series of lake ice area percentage. Based on the fitting results, the freezing start date, freezing end date, melting start date, and melting end date of the target lake are determined.

2. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The thermal infrared remote sensing data includes the surface temperature of the target lake, obtained four times a day via MODIS.

3. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The method for the synthesis process is as follows: Daily surface temperature observations of the target lake include Four; First, the temperature value of the effective pixel in the maximum composite area is calculated using the following formula (1), and then the daily observation value is calculated using formula (2). （1） （2） In equation (1), This represents the temperature value of the effective pixels within the maximum daily composite area. Let be the discriminant function for the effective temperature pixels of the target lake, when hour Take 1, when Invalid Set to 0; For surface temperature observation The total number of valid pixels; In equation (2), For the first Observed values ​​for the day.

4. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The freeze-thaw period of the target lake is from November to April of the following year; the non-freeze-thaw period is from May to October of the year.

5. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The process involves calculating the original temperature sequence based on daily observations and daily composite temperature images of the target lake, and then applying a filtering algorithm to correct the original temperature sequence to obtain a time series of lake surface temperature during the non-freeze-thaw period; including: Based on the daily synthetic temperature images, the average temperature within the target lake area is calculated and arranged in chronological order to form an original temperature sequence; The original temperature sequence was corrected using an upper profile filtering algorithm to obtain a time series of lake surface temperature during the non-freeze-thaw period.

6. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The upper profile filtering algorithm is shown in the following formula (3): (3) In equation (3), For the revised first Lake surface temperature during non-freeze-thaw periods For the first Observed values ​​for the day.

7. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The temperature threshold is -0.75 to 0.

75. When the daily observed temperature is below -0.75℃, the state is determined to be ice; when the daily observed temperature is above 0.75℃, the state is determined to be water.

8. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 7, characterized in that, The temperature threshold is determined based on the time series of lake surface temperature during the non-freeze-thaw period; the determination method includes: The time series of lake surface temperatures during the non-freeze-thaw period is divided into periods of autumn cooling trend and spring warming trend. By fitting, a trend line of temperature change over time is obtained during the autumn cooling trend period; combined with the trend line, the cooling reference date is determined; and the lower quantile value of the lake surface temperature within 5 to 15 days before the cooling reference date is calculated, and the lower quantile value is used as the lower limit of the temperature threshold; wherein, the cooling reference date is the date when the temperature enters a continuous low temperature state after the theoretical freezing point. Determine the turning point date when the temperature begins to rise continuously during the spring warming trend period, calculate the high quantile value of the lake surface temperature within 5 to 15 days after the turning point date, and use the high quantile value as the upper limit of the temperature threshold.

9. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 1, characterized in that, The fitting formula for the nonlinear curve fitting of the freezing stage in the lake ice area percentage time series is shown in Equation (4), and the fitting formula for the nonlinear curve fitting of the melting stage is shown in Equation (5). （4） In equation (4), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. t represents the percentage of lake ice area in the fitted curve, and t is time. （5） In equation (5), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. t represents the percentage of lake ice area in the fitted curve, and t represents time.

10. The method for lake surface temperature reconstruction and freeze-thaw phenology extraction based on filtering and automatic identification as described in claim 9, characterized in that, The calculation formulas for the freezing start date, freezing end date, thawing start date, and thawing end date of the target lake are shown in equations (6), (7), (8), and (9) below: Freeze start date: （6） Freeze End Date: （7） Ablation start date: （8） Ablation end date: （9） In equations (6) and (7), The date when 50% ice coverage occurs during the freezing phase. The fitting rate parameter is used during the freezing phase. The date of the freeze start date, The freeze ends on the due date; In equations (8) and (9), This refers to the date when 50% of the ice covers the area during the ablation phase. The fitting rate parameter is used during the ablation phase. The end date of ablation The date of completion of the ablation process.

Images