Border fire smoke extraction method based on stationary satellite multispectral index coupling

Abstract

The application discloses a border fire smoke extraction method based on stationary satellite multispectral index coupling and belongs to the technical field of satellite remote sensing fire monitoring. The method comprises the following steps: acquiring Kuifah satellite multi-temporal images, external fire point space-time information and meteorological wind field data, performing radiation calibration, geometric correction, atmospheric correction and interference mask generation on the acquired data, and then outputting normalized multi-band data; defining a rectangular buffer area with the fire point position as the center to perform space-time screening; calculating seven kinds of spectral indexes, i.e., NDVI, SAVI, CI, BAI, NDB N , NDR N and NDR SW , identifying smoke pixels through hierarchical collaborative discrimination and adaptive window filtering, and outputting a binary mask; calculating the centroid offset and area overlap rate based on the smoke masks of consecutive time periods to realize time sequence tracking, and combining a ground-atmosphere wind speed conversion model and a self-calibration mechanism to predict the smoke spreading trend; and finally outputting a standardized raster product in the GeoTIFF format. The application can effectively meet the near-real-time needs of fire emergency response and atmospheric environment monitoring and other businesses.

Description

Technical Field

[0001] This invention relates to the field of satellite remote sensing fire monitoring technology, and more specifically, to a method for extracting smoke and dust from border fires based on geostationary satellite multispectral index coupling. Background Technology

[0002] Fire smoke not only destroys ecosystems and threatens community safety, but also emits large amounts of fine particulate matter and harmful gases, severely impacting public health and the atmospheric environment. This is especially true in border regions, which are sparsely populated and have weak transportation and ground-based monitoring infrastructure, often resulting in delayed fire detection and response. Furthermore, border fires easily trigger cross-border smoke transport, threatening the atmospheric quality and ecological security of neighboring countries and posing a serious challenge to cross-border environmental diplomacy and emergency coordination. Therefore, rapid and accurate detection and dynamic tracking of smoke from border fires are of significant and unique practical importance for fire emergency response, cross-border joint prevention and control, and atmospheric environmental assessment.

[0003] Satellite remote sensing technology, due to its wide coverage, high timeliness, and multispectral observation capabilities, has become an important means of monitoring fire smoke and dust. The new generation of geostationary meteorological satellites, represented by Himawari-8 / 9, possess the capability to conduct comprehensive observations every 10 minutes and cover 16 bands from visible light to far-infrared, providing an ideal data source for high temporal resolution dynamic monitoring of fires and exhibiting significant advantages in monitoring the rapid evolution of smoke and dust.

[0004] In existing technologies, fire detection methods are relatively mature, and smoke identification methods have also been studied, but the following shortcomings still exist: Smoke detection and fire detection are disconnected. Existing smoke recognition algorithms often operate independently of fire information, lacking a spatiotemporal coordination mechanism. Performing full-map detection without fire information is not only computationally inefficient but also prone to misclassifying non-fire smoke such as industrial emissions and dust as fire smoke. This misclassification problem is particularly prominent in border areas due to the lack of real-world ground verification.

[0005] The problem of spectral confusion between smoke and dust, clouds, and ground background is prominent. Wildfire smoke and different types of clouds have similar spectral responses in remote sensing images, and existing threshold discrimination methods based on single or a few spectral indices are difficult to distinguish effectively, especially in cloud edges, fragmented clouds, and thin cloud areas where the misclassification rate is high. At the same time, there is a lack of systematic suppression methods for background interference such as vegetation and bare soil.

[0006] The ability to predict the spread of smoke and dust is insufficient. Existing methods mostly stop at static distribution identification and lack dynamic prediction capabilities. The few methods that do involve prediction rely on complex fluid dynamics models, which have high computational costs and slow response speeds, making it difficult to meet near real-time emergency needs. At the same time, existing predictions directly use ground wind speed to extrapolate displacement, ignoring the difference between the wind speed at height and the ground wind speed after the smoke and dust are lifted, resulting in low prediction accuracy.

[0007] The existing algorithms lack adaptability for operational deployment. Most are experimental studies in specific regions, with fixed threshold parameters and limited generalization ability, making them difficult to adapt to changes in spectral characteristics across different seasons and geographical areas. Furthermore, Himawari-8 / 9 full-panel imagery, with individual scenes ranging from 500 to 800 MB, has a high computational load for pixel-by-pixel processing, making it difficult to meet near real-time response requirements, especially in border areas where rapidly changing fire situations demand stringent timeliness.

[0008] In summary, there is an urgent need for a smoke extraction and tracking analysis method that can work in conjunction with fire point information, possess multispectral index joint discrimination capabilities, enable dynamic prediction of smoke spread, and adapt to operational deployment. Summary of the Invention

[0009] The purpose of this invention is to provide a method for extracting smoke and dust from border fires based on geostationary satellite multispectral exponential coupling, so as to solve the problems mentioned in the background art.

[0010] To achieve the above objectives, the present invention provides the following technical solution: A method for extracting smoke from border fires based on geostationary satellite multispectral index coupling includes the following steps: Step 1: Data acquisition and preprocessing. Acquire multi-temporal satellite remote sensing image data of the target area, spatiotemporal information of external fire points, and meteorological wind field data. Perform radiometric calibration, geometric correction, atmospheric correction, and interference mask generation on the satellite remote sensing image data, and output standardized multi-band reflectivity and brightness temperature data. Step 2: Data filtering. A spatial buffer area is delineated centered on the fire location in the external fire spatiotemporal information, and the matching satellite observation times are filtered based on the fire occurrence time. Effective image data within the target detection range is extracted from the preprocessed multi-band data. Step 3: Spectral index calculation. Based on the effective image data, calculate at least a variety of spectral indices, including the normalized vegetation index, soil-adjusted vegetation index, cloud detection index, fire area index, normalized "blue-near-infrared" index, normalized "red-near-infrared" index, and normalized "red-shortwave infrared" index. Step 4: Smoke and dust identification. Multiple spectral indices are collaboratively discriminated. Pixels that meet the vegetation index threshold condition are removed, and pixels that meet the cloud index threshold condition are marked as cloud masks. Pixels that simultaneously meet multiple smoke and dust enhancement index threshold conditions are marked as smoke and dust candidate pixels. Adaptive window filtering is performed on the smoke and dust candidate pixels to remove isolated noise points and fill small holes, and outputs a binary mask of fire smoke and dust. Step 5: Combustion trajectory tracking. Based on the binary smoke mask of multiple consecutive time periods, calculate the centroid offset and area overlap rate of the smoke mask in consecutive time periods, track the spatial displacement and morphological evolution of the smoke area between consecutive observation time periods, and combine meteorological wind field data and "ground-atmosphere" wind speed conversion model to predict the smoke spread trend and generate fire combustion trajectory tracking results. Step 6: Output data results. Output the binary mask of fire smoke and dust and the fire combustion trajectory tracking results as a standardized geospatial raster product.

[0011] Preferably, the satellite remote sensing image data is HSD or NetCDF format full-view observation data from the Himawari 8 / 9 geostationary meteorological satellite, with multiple time phases including at least two consecutive observation times, each time phase being 10 minutes apart.

[0012] Preferably, atmospheric correction is performed by inverting the true surface reflectance pixel by pixel by calling the pre-calculated lookup table of the 6S radiative transfer model; interference mask generation includes: using the cloud detection product built into the Sunflower satellite or generating a cloud mask based on band combination logic, using the water vapor absorption band to identify high water vapor areas, and excluding bright deserts, urban building areas, and snow and ice non-smog surface types.

[0013] Preferably, the spatial buffer zone is a rectangular buffer zone with a side length of 100 kilometers centered on the fire point.

[0014] Preferably, the calculation methods for various spectral indices are as follows: ; ; ; ; ; ; ; Where NIR is the near-infrared band reflectance, R is the red band reflectance, G is the green band reflectance, B is the blue band reflectance, and SWIR1 and SWIR2 are the reflectances of the first and second short-wave infrared bands, respectively.

[0015] Preferably, in the collaborative discrimination: when the NDVI or SAVI value exceeds a preset threshold, the corresponding pixel is determined to be a vegetation-dominant area and is masked; an initial cloud mask is generated using the CI index, and the BAI index is introduced to refine the distinction between clouds and dust through combined criteria; candidate pixels simultaneously satisfy NDB N NDR N and NDR SW Multiple index threshold conditions are used to mark as candidate pixels for smoke and dust.

[0016] Preferably, the adaptive window filtering includes: taking each candidate smoke pixel as the center, calculating the proportion of smoke pixels in the initial window; if the proportion is lower than a preset threshold, gradually expanding the window size until the condition is met or the maximum window limit is reached; if the proportion is still lower than the threshold within the maximum window range, it is determined to be an isolated noise point and removed; performing confidence-weighted smoothing on the retained smoke pixels, and using morphological closing operation for filling and smoothing processing.

[0017] Preferably, the "ground-atmosphere" wind speed conversion model is as follows: ;where V sky For atmospheric wind speed, V land H represents the ground wind speed. sky H is the atmospheric wind height. land The ground wind height is taken as 10 meters, which is the standard meteorological observation height, and k is the conversion coefficient. The prediction of the smoke and dust spread trend uses the linear extrapolation method to predict the spatial distribution of smoke and dust at the next observation time. After the actual smoke and dust detection results of the next time are obtained, the Euclidean distance between the predicted smoke and dust centroid and the actual smoke and dust centroid is calculated as the position offset error. Calibration is triggered only when the position offset error is not less than 0.5 km and the ratio of the predicted area to the actual area is between 0.5 and 2.0. The conversion coefficient k is updated online using the moving average method.

[0018] Preferably, the standardized geospatial raster product is in GeoTIFF format and adopts a unified naming convention of <satellite type>_<date>_<time>_<product type>_<file extension>; the fire trajectory tracking results are output in the form of hierarchical coded raster images, where different pixel values ​​represent non-smoke areas, the current smoke location, and the predicted smoke spread location, respectively.

[0019] A border fire smoke extraction system based on geostationary satellite multispectral index coupling includes a module for performing a border fire smoke extraction method based on geostationary satellite multispectral index coupling.

[0020] Compared with the prior art, the beneficial effects of the present invention are as follows: (1) This invention introduces external fire point spatiotemporal information for pre-data filtering, delineates a rectangular buffer area with a side length of 100 kilometers centered on the fire point location, and only filters satellite observation times matching the time before and after the fire point occurrence, compressing the effective data volume involved in the calculation to less than 1% of a single full-panel image, reducing the data volume processed per batch from 500-800MB to less than 5MB, and shortening the processing time from minutes to seconds, significantly reducing the computational load and meeting the near real-time business needs in fire emergency response scenarios. At the same time, the GeoTIFF format raster product output by this invention completely preserves the projection information of the original image, supports direct reading and visualization by mainstream GIS software, and the unified product naming standard facilitates batch management and time-series retrieval. The hierarchical coding design allows users to intuitively compare the spatial relationship between source smoke and predicted smoke in the same view, and has good adaptability for business deployment.

[0021] (2) This invention constructs a normalized vegetation index (NDVI), soil-adjusted vegetation index (SAVI), cloud detection index (CI), burn area index (BAI), and normalized blue-near-infrared index (NDB). N ), Normalized Red-Near Infrared Index (NDR) N ) and Normalized Red-Shortwave Infrared Index (NDR) SW The collaborative discrimination system of seven spectral indices, including background suppression, cloud refinement, and smoke enhancement, is implemented through a hierarchical filtering architecture. This effectively solves the prominent problem of spectral confusion between smoke, clouds, and ground background in existing technologies. On this basis, an adaptive window filtering method is introduced to dynamically adjust the size of the filtering window according to the distribution density of smoke pixels in the local neighborhood, effectively eliminating isolated noise points. At the same time, confidence-weighted smoothing and morphological closing operations are used to maintain the edge clarity and spatial integrity of the smoke area.

[0022] (3) Based on the binary mask of fire smoke in multiple consecutive time periods, this invention realizes the quantitative tracking of the spatial displacement and morphological evolution of the smoke area between consecutive observation time periods by calculating the centroid offset and area overlap rate, filling the technical gap of the existing technology that stops at static spatial distribution identification in smoke detection; on this basis, for the first time, a "ground-atmosphere" wind speed conversion model is introduced to convert the ground wind speed into the atmospheric wind speed at the height of the smoke, and the linear extrapolation method is used to realize the rapid prediction of the spatial distribution of smoke at the next observation time (10 minutes later), and a self-calibration mechanism is designed - after the real smoke detection results of the next time period are obtained, the Euclidean distance between the predicted smoke centroid and the real smoke centroid is calculated as the position offset error, and the conversion coefficient k is updated online by using the moving average method, realizing the continuous optimization from the initial empirical parameters to regional adaptation. Attached Figure Description

[0023] Figure 1 This is a flowchart of the smoke and dust extraction algorithm of the present invention; Figure 2 It is a true-color composite satellite image; Figure 3 This is a false-color composite satellite image; Figure 4 This is a graph showing the results of fire smoke detection (smoke - white; non-smoke - black). Figure 5 The results of fire smoke spread prediction (initial location - purple; predicted spread location - pink; non-smoke - black). Detailed Implementation

[0024] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0025] Example: Please see Figure 1 A method for extracting smoke from border fires based on geostationary satellite multispectral index coupling includes the following steps: Step 1: Data acquisition and preprocessing. Acquire multi-temporal satellite remote sensing image data of the target area, spatiotemporal information of external fire points, and meteorological wind field data. Perform radiometric calibration, geometric correction, atmospheric correction, and interference mask generation on the satellite remote sensing image data, and output standardized multi-band reflectivity and brightness temperature data.

[0026] 1. Data Acquisition: Satellite remote sensing imagery data is full-panel observation data in HSD or NetCDF format from the Himawari 8 / 9 geostationary meteorological satellite. Multi-temporal phases include at least two consecutive observation times, with each time interval being 10 minutes. Raw data is downloaded from JAXA's (Japan Aerospace Exploration Agency) P-Tree system or JMA's data distribution server via FTP / HTTPS protocol. Downloaded data is stored in raw HSD format or converted to NetCDF format. Each complete full-panel image contains observation data for 16 bands, with a file size of approximately 500–800 MB.

[0027] Acquisition of external fire point spatiotemporal information External fire point spatiotemporal information can be derived from real-time output results of satellite fire point detection systems (such as Himawari-8 / 9 fire point detection products), official fire point products from satellites such as MODIS / VIIRS, or fire point spatiotemporal information manually calibrated by users according to actual needs. For border areas, external fire point information can also be supplemented by fire monitoring data shared by neighboring countries or fire point coordinates provided under bilateral cooperation mechanisms. External fire point information supports input in three formats: JSON, CSV, and Shapefile. Each record must contain at least three fields: longitude (unit: degrees, value range: -180 to 180), latitude (unit: degrees, value range: -90 to 90), and observation time (UTC, format: YYYY-MM-DDTHH:MM:SS). For fire point data that supports confidence level input, subsequent processing is triggered only for fire points with a confidence level of "medium" or "high".

[0028] Acquisition of meteorological wind field data Meteorological wind field data were obtained from ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) or CMA-GD regional numerical weather prediction products from the China Meteorological Administration. Wind speed V at a height of 10 meters above the ground was acquired. land (Unit: m / s) and wind direction θ (unit: degrees, with true north as 0° and clockwise as positive), as well as the wind speed at the 850 hPa isobaric surface, are used as reference benchmarks for atmospheric height wind speed. The temporal resolution of the wind field data is no less than 1 hour, and the spatial resolution is no less than 0.25° × 0.25°. The data is resampled to the same spatial grid as the Sunflower imagery through bilinear interpolation.

[0029] 2. Radiometric calibration: Read the digital count values ​​(DN values) for each band from the original HSD or NetCDF format file. In the standard data products of the Himawari-8 / 9 satellite, the DN value of each band pixel is stored as an unsigned 16-bit integer (UInt16).

[0030] For the visible and near-infrared bands (Band 1 to Band 6), the DN values ​​are converted to TOA reflectance using the calibration lookup table published by JMA. The conversion formula is: ; where OFFSET and GAIN are calibration coefficients specific to each band, read from the data header file; θ s The zenith angle is provided by the observation geometry parameters in the header file. When the DN value is equal to the padding value (e.g., 65535), the cell is marked as invalid and does not participate in any subsequent calculations.

[0031] For the mid-infrared and thermal infrared bands (Band 7–Band 16), the DN values ​​are converted to brightness temperature. First, the DN values ​​are converted to radiance L (unit: W·m) using a scaling factor.- ²·sr - ¹·μm - ¹), and then use the inverse Planck function to calculate the brightness temperature T. B : ; Where c1 and c2 are Planck's constants, and ν is the center wavenumber (cm) of each band. - ¹), the ν values ​​for each band are provided by the band parameter table published by JMA.

[0032] The calibration process described above employs a band-by-band block parallel strategy: the entire image is divided into several data blocks of 512×512 pixels. A Python multiprocessing pool (multiprocessing.Pool) is used to independently perform calibration calculations on each data block, ensuring no interference between processes. The calculation results are then merged back into the complete image according to the original row and column numbers. After radiometric calibration, standardized atmospheric top reflectivity and brightness temperature data are output.

[0033] 3. The purpose of geometric correction is to eliminate geometric deformation caused by changes in the Earth's curvature, rotation, and satellite scanning attitude, and to ensure spatial consistency among multi-band and multi-temporal data.

[0034] First, key parameters such as the projection center latitude and longitude, column offset, row offset, and projection scale factor are read from the data header file. Himawari-8 / 9 uses the GCL projection (grid projection), and the latitude and longitude coordinates of each pixel are determined by the following projection formula: , ; Where col and row are the column and row numbers of the pixel, COFF and LOFF are the column and row offsets, and CFAC and LFAC are the column and row scale factors. All of these parameters are extracted from the data header file. After converting the image coordinate system to the Earth fixed grid coordinate system using the above formula, each pixel is then mapped onto the standard geographic coordinate grid through inverse projection.

[0035] In the resampling stage, a bicubic convolution interpolation method is employed, which effectively suppresses aliasing while maintaining edge sharpness. For registration between multi-temporal images, mutual information maximization is used as the registration quality evaluation index. When the registration error exceeds 0.5 pixels, the image at that time is marked and discarded in subsequent processing or subjected to secondary fine registration. Geometric correction ensures accurate spatial alignment of each band and multi-temporal images, guaranteeing the reliability of subsequent band calculations and temporal analysis.

[0036] 4. Atmospheric correction aims to eliminate the attenuation effects of aerosol scattering, water vapor absorption, and molecular Rayleigh scattering on surface reflectance, converting TOA reflectance into true surface reflectance.

[0037] This implementation employs a pixel-by-pixel inversion method based on a pre-computed lookup table (LUT) of the 6S radiative transfer model. The LUT is constructed using five input dimensions as indices: aerosol optical thickness (AOD, range 0–2.0, step size 0.05), solar zenith angle (range 0°–80°, step size 5°), observed zenith angle (range 0°–80°, step size 5°), relative azimuth angle (range 0°–360°, step size 20°), and atmospheric model (six standard atmospheric models: tropical, mid-latitude summer, mid-latitude winter, subpolar summer, subpolar winter, and US62). A multidimensional LUT data structure is constructed using Python's xarray and scipy.interpolate libraries and stored in NetCDF format for subsequent real-time retrieval.

[0038] For each pixel, the corresponding atmospheric model is determined based on its solar zenith angle, observed zenith angle, relative azimuth angle, and the season in which the image was acquired. The eight nearest neighbors in the four-dimensional parameter space are retrieved from the LUT, and the correction coefficients for that pixel are calculated using a three-dimensional linear interpolation method. The TOA reflectance is divided by the correction coefficients to obtain the true surface reflectance. For pixels with physical anomalies (such as reflectance less than 0 or greater than 1) during the inversion process, their surface reflectance is truncated to 0.01 and 0.99 respectively, and then marked and recorded.

[0039] 5. To reduce the interference of clouds, high water vapor areas, and bright ground surfaces on smoke and dust detection, multiple types of mask layers need to be generated.

[0040] The Cloud Mask (CLM) product from the Himawari satellite standard product is used preferentially. After reading the CLM data, pixels marked as "cloud", "cloud shadow" and "probably cloud" are directly included in the cloud mask.

[0041] For areas not covered by the CLM product, or situations where the CLM indicates clear skies but thin clouds are present, band combination logic is used for supplementary judgment. The specific judgment rules are as follows: (1) Cloud identification: When the ratio of the reflectance of Band 3 to that of Band 4 is less than 0.85 and the brightness temperature of Band 15 is less than 240K, it is identified as a cloud pixel.

[0042] (2) High-brightness surface identification: When the reflectivity of Band 4 is greater than 0.9 and the difference between the brightness temperature of Band 7 and Band 14 is less than 5K, it is identified as a high-brightness desert or snow pixel.

[0043] (3) Water vapor region identification: High water vapor regions are identified using water vapor absorption bands (Band 13, center wavelength 6.2 μm; Band 14, center wavelength 7.0 μm). When the brightness temperature difference between Band 13 and Band 15 is less than 3 K, it is identified as a high water vapor region.

[0044] All the above masking results are merged into a single mask layer, where a cell with a mask value of 1 represents an interference area that needs to be excluded and is not included in any calculations in subsequent smoke and dust detection.

[0045] After the above preprocessing steps, the raw sunflower data is converted into a high-quality input product with a unified spatial reference, true surface reflectance / brightness temperature, and effective detection range.

[0046] Step 2: Data filtering. A spatial buffer area is delineated centered on the fire location in the external fire spatiotemporal information, and matching satellite observation times are filtered based on the fire occurrence time. Effective image data within the target detection range is extracted from the preprocessed multi-band data. The spatial buffer area is a rectangular buffer area with a side length of 100 kilometers centered on the fire location.

[0047] 1. Analysis of external fire point information Read the spatiotemporal information of fire points from external input. When multiple fire points exist (e.g., multiple combustion points in the same area), process each fire point separately. For fire point data with confidence level information, only trigger the smoke detection process for fire points with a confidence level of "medium" or "high", and do not process fire points with low confidence levels.

[0048] When there are multiple fire points within the spatial buffer area, a geometric union operation is performed on the buffer areas of each fire point. The overlapping areas are then merged and used as a unified target detection range to avoid redundant calculations of the overlapping areas.

[0049] 2. Specific Implementation Methods of Spatial Filtering A rectangular buffer zone with a side length of 100 kilometers is constructed, centered on the input fire location. In the pixel coordinate system, this buffer zone is constructed as follows: first, the latitude and longitude coordinates of the fire location are converted into image row and column numbers (col) using the GCL projection formula. c row c Then, the cell span of the buffer in the row and column directions is calculated separately. The spatial resolution in the row and column directions is read from the data header file (usually 1 km / pixel, i.e., each pixel corresponds to approximately 1 km × 1 km of ground). From this, the row and column ranges of the rectangular buffer can be determined as: [row c -50, row c +50] and [col] c -50, colc +50] (in pixels). The boundaries of the buffer are cropped to ensure that the effective data range of the original image is not exceeded (0 ≤ row ≤ total number of rows in the image, 0 ≤ col ≤ total number of columns in the image).

[0050] 3. Specific Implementation Methods of Time-Based Filtering Using the observation time of the fire point as the baseline time T0, image data from seven satellite observation times before and after time T0 are selected, namely T0-2, T0-1, T0, T0+1, T0+2, T0+3, and T0+4 (corresponding to a 70-minute time window).

[0051] In each candidate time period, the cloud coverage ratio within the rectangular buffer after spatial filtering is calculated. If the cloud coverage ratio of a certain time period is less than 50%, that time period is selected first; if the cloud coverage ratio of all time periods is greater than 50%, the time period with the lowest cloud coverage ratio is selected for subsequent processing.

[0052] After the above data screening, effective pixel data within the target detection range is extracted from the preprocessed multi-band data. The number of effective pixels participating in subsequent calculations usually does not exceed 1% of the total number of pixels in the entire image (based on the total image of approximately 27.5 million pixels, the effective pixels after screening are approximately 200,000 pixels), which significantly reduces the computational load. This is particularly crucial for the operational use of monitoring stations in border areas where communication and computing resources are relatively limited.

[0053] Step 3: Spectral index calculation. Based on the effective image data, calculate at least a variety of spectral indices, including the normalized vegetation index, soil-adjusted vegetation index, cloud detection index, fire area index, normalized "blue-near-infrared" index, normalized "red-near-infrared" index, and normalized "red-shortwave infrared" index. The correspondence between the bands involved in the calculation of each index and the Himawari-8 / 9 satellite band numbers is shown in Table 1.

[0054] Table 1. Band correspondences involved in spectral index calculation The calculation methods for various spectral indices are as follows: (1) Normalized Difference Vegetation Index (NDVI): ; (2) Soil-Adjusted Vegetation Index (SAVI): ; (3) Cloud Detection Index (CI): ; (4) Burning Area Index (BAI): ; (5) Normalized Blue-Near Infrared Index (NDB) N ): ; (6) Normalized Red-Near Infrared Index (NDR) N ): ; (7) Normalized Red-Shortwave Infrared Index (NDR) SW ): ; Where NIR is the near-infrared band reflectance, R is the red band reflectance, G is the green band reflectance, B is the blue band reflectance, and SWIR1 and SWIR2 are the reflectances of the first and second short-wave infrared bands, respectively.

[0055] In practical calculations, numerical stability processing is required: when the absolute value of the denominator in any exponent calculation is less than 1 × 10⁻⁶. -6 If the pixel is not found, it is marked as invalid (NaN) and will not be included in subsequent thresholding. Simultaneously, the calculation results are dynamically truncated; for example, the NDVI calculation results are truncated to the [-1, 1] interval, and any values ​​exceeding this interval are directly assigned boundary values ​​(-1 or 1). All spectral index calculation results are saved as temporary GeoTIFF files using Float32 data type for use in subsequent steps.

[0056] All the above index calculations adopt a block-parallel strategy: the effective data area is divided into data blocks of 512×512 pixels, all seven indices are calculated independently for each data block, and parallel processing between data blocks is achieved using a Python multi-process pool (the number of processes is the number of CPU cores minus 1 by default). Finally, the calculation results of each block are merged into a complete index raster according to the original position.

[0057] Step 4: Smoke and dust identification. Multiple spectral indices are collaboratively discriminated. Pixels that meet the vegetation index threshold condition are removed, and pixels that meet the cloud index threshold condition are marked as cloud masks. Pixels that simultaneously meet multiple smoke and dust enhancement index threshold conditions are marked as smoke and dust candidate pixels. Adaptive window filtering is performed on the smoke and dust candidate pixels to remove isolated noise points and fill small holes, and outputs a binary mask of fire smoke and dust.

[0058] (1) Background suppression Pixels with vegetation and bare ground are removed using NDVI and SAVI indices. In this embodiment, the NDVI threshold is set to 0.4: pixels with an NDVI value greater than 0.4 are identified as high-coverage vegetation pixels and are masked; the SAVI threshold is set to 0.3: pixels with a SAVI value greater than 0.3 are identified as sparse vegetation or shrub pixels and are masked. These thresholds are default values ​​and can be adjusted via external configuration files in actual operation. For example, in arid and semi-arid regions, the NDVI threshold can be appropriately lowered to 0.3 to avoid excessive removal, while in humid forest regions it can be increased to 0.5.

[0059] (2) Cloud exclusion First, an initial cloud mask is generated using the CI index. When the CI value is greater than 1.5, the corresponding pixel is marked as a cloud candidate pixel. Based on this, the BAI index is introduced to perform a second discrimination of the cloud candidate pixels: if the BAI value is greater than 500 and the CI value is greater than 1.5, it is determined to be a cloud pixel and is excluded; if the CI value is greater than 1.5 but the BAI value is less than 500, it is determined to be a high-reflectivity surface (non-cloud) and is retained to enter the next discrimination stage.

[0060] (3) Smoke and dust enhancement discrimination Simultaneously utilize NDB N NDR N and NDR SW Three smoke enhancement indices are used for joint discrimination. The default threshold ranges for the three indices are: NDB N ∈ (0.1, 0.6); NDR N ∈ (-0.3, 0.3); NDR SW ∈ (0.0, 0.5); Candidate pixels must simultaneously meet the above three threshold conditions to be marked as smoke / dust candidate pixels. The logical AND condition of the three indices effectively reduces the risk of misclassification by a single index. Taking desert areas as an example, desert pixels in NDB... N There may be overlap with smoke and dust, but NDR SW Values ​​greater than 0.6 are usually excluded by this condition.

[0061] All the above thresholds are provided in an external JSON configuration file, which organizes the parameter set by a two-level index of "season-geographic region". For example, separate parameter files are maintained for spring (dry and dusty) in North China and summer (humid and vegetated) in South China. The program automatically reads the parameter set corresponding to the current processing region and season when it starts.

[0062] Detailed implementation of adaptive window filtering Adaptive window filtering is used to eliminate isolated noise caused by fragmented clouds, cloud edges, cloud shadows, and complex terrain boundaries, while maintaining the continuity and edge clarity of the smoke and dust area.

[0063] (1) Initial statistics For each candidate smoke pixel, an initial window of 3×3 pixels (corresponding to approximately 9 km² of ground) is taken. The number of pixels marked as smoke candidate pixels within the window is counted, and the proportion p of the candidate pixels to the total number of pixels in the window is calculated.

[0064] (2) Adaptive window expansion If p ≥ 30%, the current pixel is considered to be within a continuous smoke and dust area and is retained. If p < 30%, the window size is successively expanded to 5×5, 7×7, 9×9, 11×11, 13×13, and 15×15, with the percentage p recalculated after each expansion. Once p reaches 30% or more, expansion stops, and the current pixel is retained. If p is still less than 30% after the window is expanded to 15×15, the current pixel is determined to be an isolated noise point and is discarded.

[0065] (3) Confidence-weighted smoothing For the retained smoke and dust pixels, a Gaussian kernel (σ=1.5) is used to calculate the neighborhood weights using a 3×3 neighborhood as the window. This weighted smoothing is then performed by combining the comprehensive confidence scores of each neighboring pixel in the spectral index discrimination stage (the average of the normalized distances of the three smoke enhancement indices from their respective thresholds). The smoothed pixel values ​​reflect the smoke and dust probability in the local neighborhood, resulting in a more natural transition at the boundaries of the smoke and dust region.

[0066] (4) Morphological post-processing A morphological closing operation (dilation followed by erosion) is performed on the filtered result using a 5×5 square structuring element, iterated once. The closing operation can fill the small pores inside the smoke region and connect broken smoke plume fragments without changing the overall shape of the smoke region. After the closing operation, the result is re-binarized (threshold 0.5), and the final fire smoke binary mask is output, where a pixel value of 1 represents smoke and a pixel value of 0 represents non-smoke.

[0067] Step 5: Combustion trajectory tracking. Based on the binary smoke mask of multiple consecutive time periods, calculate the centroid offset and area overlap rate of the smoke mask in consecutive time periods, track the spatial displacement and morphological evolution of the smoke area between consecutive observation time periods, and combine meteorological wind field data and "ground-atmosphere" wind speed conversion model to predict the smoke spread trend and generate fire combustion trajectory tracking results.

[0068] 1. Specific implementation methods of time series correlation analysis Based on the binary smoke mask of multiple consecutive time periods (at least two consecutive observation time periods), the centroid offset and area overlap rate of the smoke mask in consecutive time periods are calculated to track the spatial displacement and morphological evolution of the smoke region between consecutive observation time periods.

[0069] (1) Centroid calculation For each time step of the smoke mask, calculate its centroid coordinates (weighted average of row and column numbers): ; ; where, value i For the i-th pixel, row is a binary value (0 or 1). i and col i These are the row and column numbers of the pixel.

[0070] (2) Calculation of centroid offset Calculate the Euclidean distance between the centroids of two adjacent time-series dust particles as the centroid offset (unit: pixel), then multiply it by the spatial resolution to convert it into physical distance (unit: km): .

[0071] (3) Calculation of area overlap rate Calculate the spatial overlap of smoke and dust regions between adjacent time intervals. For two binary masks A and B, the area overlap ratio is defined as the intersection-over-union ratio (IoU): ;in, Let be the number of pixels at the intersection of the smoke and dust masks from the two time intervals. The IoU value represents the number of pixels in the union. A higher IoU value indicates greater spatial stability and less morphological change in the smoke and dust region.

[0072] Specific Implementation Methods for Smoke and Dust Spread Prediction By combining meteorological wind field data with a "ground-atmosphere" wind speed conversion model, the spread trend of smoke and dust is predicted, and the fire combustion trajectory tracking results are generated.

[0073] (1) Wind field data preprocessing The acquired meteorological wind field data was resampled to the same spatial grid as the sunflower image using bilinear interpolation to ensure that each smoke pixel location corresponds to an independent wind vector. For locations far from the fire point (beyond the buffer boundary), nearest neighbor extrapolation of the wind field data at the buffer boundary was used for padding.

[0074] (2) Ground-atmosphere wind speed conversion The following empirical model is used to convert the ground wind speed into the atmospheric wind speed at the height of the smoke and dust: ;where V sky For atmospheric wind speed, Vland H represents the ground wind speed. sky H is the atmospheric wind height. land H represents the ground wind height, k is the conversion factor; land Take 10 meters (standard meteorological observation altitude), H sky Based on the vegetation type of the fire area, the following settings are applied: 500 meters for grassland / shrub fires, 800 meters for forest fires, and 650 meters for other types (default value). The initial value of the conversion factor k is set to 0.143, which is derived from the logarithmic law approximation of the near-surface wind speed profile. The wind direction is considered constant in the vertical direction (i.e., the smoke drift direction is consistent with the surface wind direction), and only the wind speed magnitude is corrected for height.

[0075] (3) Linear extrapolation method for prediction Using the current time (T0) binary smoke and dust mask as the initial field, predict the spatial distribution of smoke and dust at time T0+10 (Δt=600 seconds). For each smoke and dust source pixel with a value of 1 in the mask, assign the corresponding wind speed V to its location. sky Converting wind direction θ to a displacement vector in Cartesian coordinates: , ; Displacement is measured in pixels (metric units removed, floating-point values ​​are used based on pixel spatial resolution). To avoid jagged displacement artifacts caused by integer rounding, a bicubic interpolation kernel is used to distribute the "dust mass" of source pixels to neighboring pixels in the target grid. The dust mass (value 1) of each source pixel is distributed to a 4×4 neighborhood around the target coordinates using a bicubic weighting, with the sum of the weights normalized to 1. After energy distribution for all source pixels, the cumulative dust mass in the target grid is normalized to ensure that the total dust mass of the predicted dust pixels is equal to the total dust mass of the source pixels (allowing for ±2% fluctuation due to diffusion).

[0076] 3. Specific implementation methods of the self-calibration mechanism Once the next real dust detection result is obtained, the self-calibration process is triggered: (1) Calculate the centroid coordinates (weighted average of cell row and column numbers) of the predicted smoke and the actual smoke respectively.

[0077] (2) Calculate the Euclidean distance between the predicted dust centroid and the actual dust centroid as the position offset error E (unit: km).

[0078] (3) Calibration is triggered only when E ≥ 0.5km and the ratio of predicted area to actual area is between 0.5 and 2.0 (excessive area difference indicates a sudden change in wind field or fire situation, and is not suitable for calibration).

[0079] (4) Derive the optimal conversion coefficient k from the current error E. optimalAssuming the error is mainly caused by wind speed deviation, under the condition that the wind direction remains unchanged, the wind speed correction is proportional to the displacement correction, k optimal Calculated based on the ratio of the inverse wind speed to the actual wind speed.

[0080] (5) Update the value of k using the moving average method: k new = 0.7 × k current + 0.3 × k optimal .

[0081] Each independent fire event maintains a unique k value. After the fire event ends (the event is considered to have ended after three consecutive smoke and dust detection results), the final k value is written into the regional historical parameter database for subsequent fire events in the same region to inherit and use.

[0082] Step 6: Data Output. The results of the fire smoke binary mask and fire trajectory tracking are output as standardized geospatial raster products. All output products are in GeoTIFF format. Georeferenced information (GeoTransform six parameters and ProjectionRef projection string) is copied from the input sunflower image to ensure that the output products can be directly overlaid and analyzed with other geospatial data in GIS software. The binary mask products are stored in UInt8 data type, and the hierarchical coded products are also stored in UInt8 data type. LZW lossless compression is used to reduce file size. The standardized geospatial raster products are in GeoTIFF format and use the unified naming convention <satellite type>_<date>_<time>_<product type>_<file extension>. The fire trajectory tracking results are output in the form of hierarchical coded raster images, where different pixel values ​​represent non-smoke areas, the current smoke location, and the predicted smoke spread location, respectively.

[0083] Specific examples are as follows: Smoke and dust detection product: H9_20240514_0430_smoke.tif; Spread prediction product: H9_20240514_0430_fcst.tif; Smoke and dust confidence level product: H9_20240514_0430_prob.tif; Multi-band spectral index product: H9_20240514_0430_indices.tif; The date format is YYYYMMDD, and the time format is HHMM (UTC time, 24-hour format).

[0084] Hierarchical coding rules The fire trajectory tracking results are output as a hierarchical coded raster image, where different pixel values ​​represent non-smoke areas, the current smoke location, and the predicted smoke spread location, respectively. The pixel values ​​are defined as follows: 0: Non-smoke / dust-free area; 1: Initial location of smoke and dust (smoke and dust detected in the current time); 2: Predicted location of smoke and dust spread (the area predicted to be covered at the next observation time); This hierarchical coding design allows users to intuitively compare the spatial relationship between source dust and predicted dust in the same view, and is especially convenient for emergency command centers on both sides of the border to obtain consistent dust dynamic information simultaneously.

[0085] Output of auxiliary products In addition to the core products mentioned above, the following auxiliary products may also be selectively provided: (1) Spectral index feature file: All seven spectral indices (NDVI, SAVI, CI, BAI, NDB) N NDR N NDR SW Each band is stored as a GeoTIFF multi-band file, with the band order as shown above, and each band is stored in Float32 format.

[0086] (2) Smoke confidence map: The continuous probability values ​​(0-1) output by the spectral index collaborative discrimination link are generated after adaptive filtering, reflecting the reliability of smoke identification, and stored in Float32 format.

[0087] (3) Prediction deviation statistics table: CSV format text file, which records the predicted area (km²), the actual area (km²), the centroid offset error (km), the area error percentage, the currently used k value and the updated k value for each hour, and is used for algorithm accuracy evaluation and continuous optimization.

[0088] A fire combustion trajectory tracking and analysis system based on multi-source satellite time-series data includes a module for performing a border fire smoke extraction method based on geostationary satellite multispectral index coupling.

[0089] Taking a typical fire incident as an example, the relevant data results of the fire smoke detection task are presented. The algorithm is based on continuous multi-temporal observation data from the Himawari-9 satellite, and generates a series of process documents and final products according to the aforementioned preprocessing and algorithm flow.

[0090] Figure 2-3 The algorithm outputs local area images during the data filtering stage. After filtering the original full-disk image using fire point spatiotemporal information, only the image area within a rectangular range with a radius of 50 kilometers centered on the fire point is retained. Figure 2The satellite data image is a true-color composite (Band 3, 2, 1 combination). In this color scheme, clouds and smoke are not distinguishable and appear as grayish-white or bluish-white. Figure 3 The satellite data images are composited in false color (Band 3, 4, 6). In this color scheme, the clouds are still grayish-white or bluish-white, while the smoke and dust show a highly distinguishable yellow-orange color, which can provide a basis for the original distribution of the smoke and dust.

[0091] Figure 4 This is the binary smoke mask output by the algorithm. Smoke areas are represented by white pixels, while non-smoke areas remain as black pixels. As can be seen from the image, the detected smoke area largely matches the visually interpreted plume outline, and the overall texture features of the smoke are clearly visible. The algorithm successfully removed... Figure 3 The cloud regions mixed with smoke and dust retained the smoke and dust signals relatively intact.

[0092] Figure 5 This is a multi-level encoded raster image, showing the predicted smoke coverage area for the next time interval (T0+10 min) based on the smoke detection results at the current time interval (T0). In the legend, the purple area represents the initial location of the smoke detected at the current time interval, and the pink area represents the predicted smoke spread location. The prediction results use ground wind speed (tested wind speed value 5.2 m / s, wind direction test value 60°, ground height 2 meters) as input conditions, combined with "ground..." The average spread speed of smoke and dust in the atmosphere, derived from the atmospheric wind speed conversion model (at an atmospheric wind height of 800 meters), achieved the task of predicting the spread of smoke and dust. Spatially, the predicted smoke and dust drifted approximately 8 kilometers in a direction 60 degrees south of west, exhibiting a certain degree of objective lateral diffusion.

[0093] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for extracting smoke and dust from border fires based on geostationary satellite multispectral index coupling, characterized in that, Includes the following steps: Step 1: Data acquisition and preprocessing. Acquire multi-temporal satellite remote sensing image data of the target area, spatiotemporal information of external fire points, and meteorological wind field data. Perform radiometric calibration, geometric correction, atmospheric correction, and interference mask generation on the satellite remote sensing image data, and output standardized multi-band reflectivity and brightness temperature data. Step 2: Data filtering. A spatial buffer area is delineated centered on the fire point location in the external fire point spatiotemporal information, and the matching satellite observation times are filtered based on the fire point occurrence time. Effective image data within the target detection range is extracted from the preprocessed multi-band data. Step 3: Spectral index calculation. Based on the effective image data, calculate at least a variety of spectral indices, including the normalized vegetation index, soil-adjusted vegetation index, cloud detection index, fire area index, normalized "blue-near-infrared" index, normalized "red-near-infrared" index, and normalized "red-shortwave infrared" index. Step 4: Smoke and dust identification. The multiple spectral indices are collaboratively discriminated. Pixels that meet the vegetation index threshold condition are removed, and pixels that meet the cloud index threshold condition are marked as cloud masks. Pixels that simultaneously meet multiple smoke and dust enhancement index threshold conditions are marked as smoke and dust candidate pixels. Adaptive window filtering is performed on the smoke and dust candidate pixels to remove isolated noise points and fill small holes, and output a fire smoke and dust binary mask. Step 5: Combustion trajectory tracking. Based on the binary smoke mask of the fire at multiple consecutive time intervals, calculate the centroid offset and area overlap rate of the smoke mask at consecutive time intervals, track the spatial displacement and morphological evolution of the smoke area at consecutive observation time intervals, and combine the meteorological wind field data and the "ground-atmosphere" wind speed conversion model to predict the smoke spread trend and generate the fire combustion trajectory tracking results. Step 6: Output data results. Output the binary mask of fire smoke and dust and the fire combustion trajectory tracking results as a standardized geospatial raster product.

2. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The satellite remote sensing image data is HSD or NetCDF format full-view observation data from the Himawari 8 / 9 geostationary meteorological satellite. The multi-temporal phase includes at least two consecutive observation time periods, with each time period spaced 10 minutes apart.

3. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The atmospheric correction retrieves the true surface reflectance pixel by pixel by calling the pre-calculated lookup table of the 6S radiative transfer model. The interference mask generation includes: using the cloud detection product built into the Sunflower satellite or generating a cloud mask based on band combination logic, using the water vapor absorption band to identify high water vapor areas, and excluding bright deserts, urban building areas, and snow and ice non-smog surface types.

4. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The spatial buffer zone is a rectangular buffer zone with a side length of 100 kilometers centered on the fire point.

5. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The calculation methods for the various spectral indices are as follows: ； ； ； ； ； ； ； Where NIR is the near-infrared band reflectance, R is the red band reflectance, G is the green band reflectance, B is the blue band reflectance, and SWIR1 and SWIR2 are the reflectances of the first and second short-wave infrared bands, respectively.

6. The method for extracting border fire smoke based on geostationary satellite multispectral exponential coupling according to claim 1, characterized in that, In the collaborative discrimination: when the NDVI or SAVI value exceeds a preset threshold, the corresponding pixel is determined to be a vegetation-dominant area and is masked; an initial cloud mask is generated using the CI index, and the BAI index is introduced to finely distinguish between clouds and dust through combined criteria; candidate pixels simultaneously satisfy NDB N NDR N and NDR SW Multiple index threshold conditions are used to mark as candidate pixels for smoke and dust.

7. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The adaptive window filtering includes: taking each candidate smoke pixel as the center, calculating the proportion of smoke pixels in the initial window; if the proportion is lower than a preset threshold, gradually expanding the window size until the condition is met or the maximum window limit is reached; if the proportion is still lower than the threshold within the maximum window range, it is determined to be an isolated noise point and removed; performing confidence-weighted smoothing on the retained smoke pixels, and using morphological closing operation for filling and smoothing processing.

8. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The "ground-atmosphere" wind speed conversion model is as follows: ;where V sky For atmospheric wind speed, V land H represents the ground wind speed. sky H is the atmospheric wind height. land The ground wind height is taken as 10 meters, which is the standard meteorological observation height; k is the conversion coefficient; the prediction of the smoke and dust spread trend uses the linear extrapolation method to predict the spatial distribution position of the smoke and dust at the next observation time; after the actual smoke and dust detection results of the next time are obtained, the Euclidean distance between the predicted smoke and dust centroid and the actual smoke and dust centroid is calculated as the position offset error. Calibration is triggered only when the position offset error is not less than 0.5 km and the ratio of the predicted area to the actual area is between 0.5 and 2.

0. The conversion coefficient k is updated online using the moving average method.

9. The method for extracting border fire smoke based on geostationary satellite multispectral index coupling according to claim 1, characterized in that, The standardized geospatial raster product is in GeoTIFF format and adopts a unified naming convention of <satellite type>_<date>_<time>_<product type>_<file extension>; the fire trajectory tracking results are output in the form of hierarchical coded raster images, where different pixel values ​​represent non-smoke areas, the current smoke location, and the predicted smoke spread location, respectively.

10. A border fire smoke extraction system based on geostationary satellite multispectral index coupling, characterized in that, Includes a module for performing the method according to any one of claims 1 to 9.

Images