An Automatic Reprocessing Method for FY-3D MERSI L1B Data

Abstract

Aiming at the MERSIL1B orbital observation data received and preprocessed by the Fengyun-3D (FY-3D) meteorological satellite, this invention proposes an automated reprocessing method for FY-3D MERSIL1B data. It uses a Python development platform and GDAL raster processing technology to intensively manage modules such as data detection and transmission, radiometric calibration, spectral reflectance / brightness temperature calculation, equal-latitude and longitude projection, regional cropping, data format conversion, and archiving, realizing the close connection and automated operation of each processing module. The processing results show that this method has the advantages of high data reading efficiency and less processing time consumption. It can automatically process and archive multi-channel images of equal-latitude and longitude grids with spatial resolutions of 1 km and 250 m. The image quality is basically the same as that of the original processing technology and can be applied to the satellite remote sensing monitoring system, thus effectively supporting the meteorological remote sensing service business in this region.

Description

Technical Field

[0001] The present invention relates to the technical field of FY-3D MERSI L1B data reprocessing techniques, and particularly relates to an automated reprocessing method for FY-3D MERSI L1B data. Background Art

[0002] FengYun-3D (FY-3D) was launched at the Taiyuan Satellite Launch Center in November 2017. It is China's second-generation polar-orbiting meteorological satellite. Together with FengYun-3C (FY-3C) and FengYun-3E (FY-3E), it forms China's polar-orbiting meteorological satellite dawn, morning, and afternoon observation service networks, providing a complete global coverage of remote sensing observation data for numerical weather prediction assimilation services every 6 hours. FY-3D is equipped with 10 sets of remote sensing observation instruments, including the second-generation Medium Resolution Spectral Imager (MERSI-II), the second-generation Microwave Temperature Sounder (MWTS-II), the Microwave Radiation Imager (MWRI), the second-generation Microwave Humidity Sounder (MWHS-II), the Space Environment Monitor (SEM), the Global Navigation Satellite Occultation Sounder (GNOS), the Infrared Hyperspectral Atmospheric Sounder (HIRAS), the Near-Infrared Hyperspectral Greenhouse Gas Monitor (GAS), the Wide-Angle Aurora Imager (WAI-I), and the Ionospheric Photometer (IPM). The observation data are widely used in the fields of space weather, atmospheric radiation, land surface ecology, numerical weather prediction, etc. MERSI-II is widely used for the production of monitoring and evaluation products for regional atmospheric aerosols, clouds, fire monitoring, and land surface ecology. At present, domestic remote sensing application research based on FY-3D MERSI data has been carried out, achieving fire monitoring research based on FY-3D MERSI-II far-infrared data (Zheng Wei, Chen Jie, etc., 2020) and AOD remote sensing inversion (Chen Hui, etc., 2022), studying the cloud detection model of FY-3D MERSI based on deep learning (Qu Jianhua, etc., 2019), and conducting quality assessment on the FY-3D MERSI NDVI product (Wang Yuanyuan, etc., 2022). Provincial meteorological departments mainly use FY3D MERSI-II data products to carry out local meteorological remote sensing real-time monitoring service operations.

[0003] Since the FY-3D MERSI L1B data does not have equal longitude and latitude projection information, the spectral band data is mainly stored as integer radiance gray values (Digital Number, DN). For non-remote sensing technology professionals to analyze and apply this data, radiation calibration, pixel reflectance / brightness temperature calculation, reprojection, geometric correction, regional cropping, and image format conversion need to be performed on the data. After reprocessing and processing, multi-channel raster remote sensing image data that is easy to understand and apply is output. Currently, in local operations, the reprocessing and processing method of FY-3D MERSI L1B data requires a large amount of manual intervention and settings, and no automated processing engineering mode for data processing has been formed. The reprocessing efficiency is not high and cannot meet the requirements of the intelligent development of meteorological remote sensing monitoring services. Summary of the Invention

[0004] In view of the above problems, the present invention provides an automated reprocessing method for FY-3D MERSI L1B data. Based on the Network File System (NFS) transfer protocol, it can be docked with the data transfer of the preprocessing platform of the Fengyun-3 provincial utilization station. Python and GDAL raster technology are used to reprocess the orbital data, realizing the coupling connection and trigger-based asynchronous automatic operation of the main modules such as FY-3D MERSI L1B data file detection, metadata and geographic data parsing, radiation calibration, spectral reflectance / brightness temperature calculation, reprojection, geometric correction, and regional cropping, improving the processing efficiency and automation level of FY-3D MERSI data, and automatically generating and archiving multi-channel raster format remote sensing images with two spatial resolutions of 250m and 1km.

[0005] The technical solution to achieve the object of the present invention is as follows:

[0006] An automated reprocessing method for FY-3D MERSI L1B data, characterized by comprising the following steps:

[0007] Step 1: Receive and preprocess the orbital observation data of the FY-3D MERSI L1B transit area in real time through the Fengyun-3 provincial direct receiving station;

[0008] Step 2: Automatically detect the newly received FY-3D MERSI L1B data, and perform automatic integrity check and reading;

[0009] Step 3: Analyze the attributes of the inspected observation data to obtain the metadata information of the spectral band pixels of the data file; perform feature analysis on the metadata information to obtain the observation angle, correction coefficient, and geographic coordinate information;

[0010] Step 4: Calculate the atmospheric apparent reflectance of the spectral band channels or the brightness temperature of the thermal radiation band based on the obtained observation angle and correction coefficient;

[0011] Step 5: Based on the obtained pixel geographic coordinates, by setting affine parameters, perform equi - longitude - latitude grid projection and geometric correction on the image data;

[0012] Step 6: According to the longitude - latitude range of the business observation area, crop the monitoring image data, and use the GDAL raster data processing technology to generate a regional image product in TIF format.

[0013] Further, the specific operation steps of Step 2 include:

[0014] Step 21: Regularly scan the received FY - 3D MERSI L1B data files in the FY - 3D MERSI L1B orbit - time - series observation data directory, identify them, and generate an identification file;

[0015] Step 22: Determine whether there is a new identification file. If not, return to Step 21 to continue scanning; if so, enter Step 23 for a complete check;

[0016] Step 23: Based on the file size, the number of data files, and the time attribute of the identification file, automatically check the integrity of the data identification file based on the scheduled task service. If the file is complete, data reading can be performed; if the file is incomplete, return to Step 21.

[0017] Further, the specific steps of Step 4 include:

[0018] Step 41: Calculate the atmospheric apparent reflectance

[0019] Step 411: After detecting outliers and missing values for the digital quantization values of each - band pixels, perform radiometric calibration on the reasonable pixel DN values through Equation (1) to obtain the calibrated radiometric quantization dn after radiometric calibration, that is:

[0020] dn = DN * slope+intercept (1)

[0021] Where, slope and intercept are the gain and offset coefficients required for radiometric calibration, which can be obtained from the dataset attributes; DN is the radiometric digital quantization value of the channel pixel, and dn is the radiometric value after radiometric calibration;

[0022] Step 412: Calculate the corrected pixel radiance R of this channel band according to the correction coefficient:

[0023] R = Cal2 * dn 2 +Cal1 * dn+Cal0 (2)

[0024] Where, Cal2, Cal1, and Cal0 are correction coefficients;

[0025] Step 413: Perform solar zenith angle correction on the average spectral radiation at the top of the atmosphere and calculate the atmospheric apparent reflectance v of the spectral band channel through equation (3). toa :

[0026]

[0027] Among them, E m is the average spectral radiation at the top of the atmosphere, d is the heliocentric astronomical unit, and θ is the solar zenith angle; L sensor is the irradiance at the entrance pupil of the satellite sensor;

[0028] Step 414: Calculate L through equation (4). sensor Substitute equation (4) into equation (3) to obtain equation (5), and quickly calculate the atmospheric apparent reflectance through equation (5):

[0029] L sensor =(R * E m ) / π (4)

[0030]

[0031] Step 42: Calculate the brightness temperature

[0032] Step 421: After detecting outliers and missing values in the thermal radiation band, perform radiometric calibration on the radiation L of the infrared channel and calculate the radiance L of the calibrated band λ :

[0033] L λ =L * slope + intercept (6)

[0034] Among them, L is the uncalibrated radiance quantization value of the channel, slope is the calibration gain, intercept is the calibration offset, and L λ is the calibrated radiance;

[0035] Step 422: Calculate the equivalent blackbody brightness temperature T of the thermal infrared band according to the Planck function conversion formula (7). f :

[0036]

[0037] Among them, C1 and C2 are the integrated calculation coefficients, h = 6.62606876e-34 J.s is the Planck constant, c = 2.99792458e+8 m / s is the speed of light, k = 1.3806503e-23 J / K is the Boltzmann constant, υ is the wave number of this channel, and L λ is the radiance after radiometric calibration;

[0038] Step 423: Perform correction on T through equation (8). fRevise it to obtain the brightness temperature T observed in this band b :

[0039] T b = TBB_a * T f + TBB_b (8)

[0040] where TBB_a and TBB_b are brightness temperature correction coefficients, T f is the equivalent blackbody brightness temperature of the thermal infrared channel, T b is the brightness temperature of the band, and the brightness temperature of the band observed by the satellite remote sensing observation instrument.

[0041] Furthermore, the specific steps of step 5 include:

[0042] Step 51: According to the pixel geographic coordinates of different spatial resolutions read, obtain the maximum latitude Latmax, minimum longitude Lonmin, pixel row height Res_Line, and pixel column width Res_pixel of images with different resolutions;

[0043] Step 52: Set the satellite image affine parameters according to Latmax, Lonmin, Res_Line, and Res_pixel, construct an equi - longitude - latitude projection grid, and establish the relationship between the original image pixel geographic coordinates and the projection grid row - column numbers through equations (9) - (10):

[0044] Lat var = Lat max + lines * Res_line (9)

[0045] Lon var = Lon min + cols * Res_pixel (10)

[0046] where Latmax is the maximum latitude of the satellite image, Latvar is the pixel latitude of the satellite image, Res_line is the row height of the projection grid, and lines is the number of rows of the projection grid; Lonmin is the minimum longitude of the satellite image, Lonvar is the pixel longitude of the satellite image, Res_pixel is the column width of the projection grid, and cols is the number of columns of the projection grid;

[0047] Step 53: Based on the mapping relationship, perform geometric correction and reprojection of the satellite image through gdalWarp.

[0048] Furthermore, the specific steps of step 6 include:

[0049] Step 61: Obtain the longitude - latitude coordinate range information of the area of the FY - 3D - MERSI image to be monitored, and crop the FY - 3D MERSI image data;

[0050] Step 62: Within the cropped area, construct a three-dimensional raster dataset based on the GDAL raster data processing technology;

[0051] Step 64: Write the obtained solar and satellite observation angles of pixels, the calculated atmospheric apparent reflectance of spectral band channels, and the pixel brightness temperature of the thermal infrared band into the three-dimensional raster dataset in sequence according to the size of the image data row and column dimensions;

[0052] Step 65: Output multi-channel raster image data with two spatial resolutions of 1 km and 250 m.

[0053] Compared with the prior art, the present method has the following beneficial effects:

[0054] The present invention reprocesses the orbital data, realizes the coupling connection and trigger-based asynchronous automatic operation of the main modules such as the detection of FY-3D MERSI L1B data files, the parsing of metadata and geographic data, radiometric calibration, the calculation of spectral reflectance / brightness temperature, reprojection, geometric correction, and regional cropping, improves the data processing efficiency and automation level of FY-3D MERSI, and automatically generates and archives multi-channel band raster remote sensing images with two spatial resolutions of 250 m and 1 km; experiments show that the data volume of the multi-channel TIF images processed by the method of the present invention is generally smaller than that before processing, which is convenient for users to read and call; and because the present invention uses the network file transfer protocol NFS to establish a real-time transfer link with the FY-3D storage device, it realizes the data automatic detection and fast reading ability, and the time consumption is significantly reduced, thus improving the data reprocessing efficiency to a certain extent. Brief Description of the Drawings

[0055] Figure 1 It is a flowchart of data automatic detection;

[0056] Figure 2 It is the automatic reprocessing method proposed by the present invention;

[0057] Figure 3 It is a trend chart of the change in the time-consuming of automatic reprocessing;

[0058] Figure 4 It is a reflectance image of band13 (0.709 μm); among them Figure 4 (a) is the image obtained by the original processing technology, Figure 4 (b) is the image obtained by the present invention;

[0059] Figure 5 It is an image of the thermal infrared band band23 (8.550 μm); among them Figure 5 (a) is the image obtained by the original processing technology, Figure 5 (b) is the image obtained by the present invention;

[0060] Figure 6 are the MAE and RMSE for the visible and infrared RSB bands of band1 - band19; where Figure 6 (a) is the MAE, Figure 6 (b) is the RMSE;

[0061] Figure 7 are the MAE and RMSE of the brightness temperature for the long - wave thermal infrared radiation TEB bands in Band20 - Band25; where Figure 7 (a) is the MAE, Figure 7 (b) is the RMSE;

[0062] Figure 8 is the result of fitting the original processing method and the processing method of the present invention; where Figure 8 (a) is the linear regression fitting result for band17, and 8(b) is the linear regression fitting result for band25;

[0063] Figure 9 is a remote sensing image with a spatial resolution of 250m; where Figure 9 (a) is band3 (0.650μm), Figure 9 (b) is band24 (10.8μm). Specific implementation manners

[0064] In order to enable those of ordinary skill in the art to better understand the technical solution of the present invention, the technical solution of the present invention will be further described below in conjunction with the accompanying drawings and implementation cases.

[0065] 1. Introduction to FY - 3D MERSI L1B data

[0066] The data to be processed in the present invention is the FY - 3D MERSI L1B data received and pre - processed by the FY - 3 provincial direct - receiving station system. MERSI of FY - 3 is the second - generation medium - resolution spectral imager, and there are a total of 25 spectral channels (see Table 1). Compared with the previous generation MERSI - I, 5 spectral channels such as cirrus, water vapor, and land surface temperature in the atmospheric window are added.

[0067] The FY-3D MERSI L1B orbital overpass observation data received and preprocessed by the Fengyun-3 provincial direct receiving stations mainly include four HDF5 format data files, namely the calibrated observation data set with a spatial resolution of 1 km, the geographic coordinate data with a spatial resolution of 1 km, the calibrated observation data set with a spatial resolution of 250 m, and the geographic coordinate data with a spatial resolution of 250 m (see Table 2). Data descriptions can be queried from the data service website of the National Satellite Meteorological Center. The calibrated observation data set with a spatial resolution of 1 km mainly contains 4 subsets, the data set with a spatial resolution of 250 m contains 6 subsets, and there are geographic coordinate files with different resolutions. Table 1 shows the spectral band information of MERSI-I and MERSI-II, and Table 2 shows the information of FY-3D MERSI L1B data files.

[0068] Table 1 MERSI-I and MERSI-II Spectral Information

[0069]

[0070]

[0071] Table 2 FY-3D MERSI L1B File Information

[0072]

[0073]

[0074] Next, the reprocessing method proposed in the present invention will be described from the aspects of reflectance calculation in the Reflective Solar Bands (RSB), brightness temperature calculation in the Thermal Emissive Bands (TEB), equidistant cylindrical projection, and raster data generation. It is based on the Python development platform and the gdal processing module to achieve intensive management and automated operation of data in the processing module.

[0075] 2. Reflectance Calculation

[0076] The data storage types of the FY-3D MERSI L1B-level data files with different resolutions (1 km, 250 m) are all UInt16 pixel spectral brightness DN values. It is not convenient to directly apply them for remote sensing monitoring services. After radiometric calibration of the DN values of each channel data, the corresponding reflectance or radiance brightness temperature needs to be calculated and then applied to the inversion models of basic observation factors such as vegetation indices and land surface temperature.

[0077] The reflectance refers to the ratio of the radiance reflected by the remotely sensed object received by the satellite sensor to the solar radiance absorbed by the remotely sensed object, which characterizes the absorption and reflection capabilities of the remotely sensed object to solar radiation. The reflectance mentioned in this paper is the atmospheric apparent reflectance, which is the sum of the spectral reflectances of the atmosphere and the remotely sensed ground object received by the satellite sensor. The visible and near-infrared channels of the Reflective Solar Bands (RSBs) of the MERSI instrument on FY-3D have a total of 19 channels, including 4 channels with a spatial resolution of 250m and 15 channels with a spatial resolution of 1km. Calculating the atmospheric apparent reflectance of the spectral band requires a large amount of computing time. For each band, the Digital Number (DN) value of the pixel, after being detected and processed for outliers and missing values, is radiometrically calibrated for the reasonable DN value through Equation (1) to obtain the calibrated radiometric quantization dn, as follows:

[0078] dn = DN * slope + intercept (1)

[0079] Among them, slope and intercept are the gain and offset coefficients required for radiometric calibration, which can be obtained from the dataset attributes; DN is the radiometric digital quantization value of the channel pixel, and dn is the radiometric quantization after radiometric calibration;

[0080] After radiometric calibration, read the file correction coefficient and calculate the corrected pixel radiance R of the channel band, as follows:

[0081] R = Cal2 * dn 2 + Cal1 * dn + Cal0 (2)

[0082] Among them, Cal2, Cal1, and Cal0 are correction coefficients, which can be obtained from the L1 data file;

[0083] After calculating the corrected radiation R of the band channel, perform solar zenith angle correction on the average spectral radiation E at the top of the atmosphere through Equation (3) m and calculate the atmospheric apparent reflectance ρ of the spectral band channel toa :

[0084]

[0085] Among them, E m is the average spectral radiation at the top of the atmosphere, d is the heliocentric astronomical unit, θ is the solar zenith angle, and these variables can be obtained from the dataset attributes of the L1B file dataset; L sensor is the irradiance at the entrance pupil of the satellite sensor,

[0086] L sensor and the corrected radiation R, the average spectral radiation E at the top of the atmosphere mThere is a relationship as shown in Equation (4), and L can be calculated through this relationship sensor :

[0087] L sensor = (R * E m ) / π (4)

[0088] Substitute Equation (4) into Equation (3) and simplify to obtain the calculation formula for the apparent atmospheric reflectance of the channel (5). Convert the relationship between the apparent atmospheric reflectance ρ toa and L sensor into the relationship with the radiance R for channel calibration and correction. The apparent atmospheric reflectance ρ toa can be quickly calculated through Equation (5):

[0089]

[0090] 3. Brightness temperature calculation

[0091] When the spectral radiance of an object is equal to the spectral radiance of a certain blackbody, the physical temperature of the blackbody can be used to characterize the brightness temperature (i.e., the bright temperature) of the object. Although the bright temperature has the dimension and unit of temperature, it does not have the physical meaning of temperature and is a synonym for the radiance of the object. There are a total of 6 thermal infrared radiation spectral TEB bands in FY-3D MERSI-II, among which there are 2 thermal infrared bands with a spatial resolution of 250 m. The data of these thermal infrared bands are stored as uncalibrated radiance integer quantization values L. After completing the detection and processing of outliers and missing values in the thermal radiation bands, perform radiometric calibration on the infrared channel radiance L λ , and the solution process is as follows:

[0092] L λ = L * slope + intercept (6)

[0093] where L is the uncalibrated radiance quantization value of the channel, slope is the gain of calibration, intercept is the offset of calibration, and L λ is the calibrated radiance;

[0094] The unit of the radiance data of the thermal infrared bands here is mW / (m 2 *cm -1 *sr), but in the atmospheric correction or land surface temperature inversion algorithm models, the unit of radiance is sometimes W / (m 2 *μm*sr) or μW / (cm 2 *sr*nm). According to the unit of the specific application scenario of the radiance, it is necessary to perform the conversion of the wave number υ and the wavelength λ and the unit transformation when necessary. This unit transformation adopts the publicly disclosed method in the existing technology and will not be elaborated here.

[0095] According to the Planck function conversion formula (7), calculate the equivalent blackbody brightness temperature T in the thermal infrared band f :

[0096]

[0097] where C1 and C2 are the integrated calculation coefficients, h = 6.62606876e-34 J·s is the Planck constant, c = 2.99792458e+8 m / s is the speed of light, k = 1.3806503e-23 J / K is the Boltzmann constant, is the wave number of this channel, and L λ is the radiance after radiometric calibration.

[0098] After completing the calculation of the equivalent blackbody brightness temperature T of the thermal infrared channel through Equation (7) f and then, correct T through Equation (8) f to obtain the brightness temperature T observed in this band b :

[0099] T b = TBB_a * T f + TBB_b (8)

[0100] where TBB_a and TBB_b are the channel brightness temperature correction coefficients, T f is the equivalent blackbody brightness temperature of the thermal infrared channel, T b is the band brightness temperature, and is the band brightness temperature observed by the satellite remote sensing observation instrument.

[0101] The main codes for calculating the reflectance in the solar reflection band and the brightness temperature in the thermal radiation band are as follows:

[0102]

[0103]

[0104]

[0105] 4. Equal latitude-longitude reprojection

[0106] Since the FY-3D MERSI L1B data image does not have equal latitude-longitude projection information, after the reflectance / radiance temperature calculation and processing of the spectral band data, it is necessary to reproject the image data. Based on the latitude-longitude information of the original orbit data, perform equal latitude-longitude projection on the image. Read the geographic coordinates of the FY-3D L1B 1 km and 250 m resolution image data respectively, and obtain the maximum latitude Lat max and the minimum longitude Lon min and the row height Res_Line and column width Res_pixel of the image pixels. According to these parameters, set the affine parameters of the satellite image and construct an equal latitude-longitude projection grid as follows.

[0107] Lat var = Lat max + lines * Res_line(9)

[0108] Lon var = Lon min + cols * Res_pixel(10)

[0109] Among them, Lat max is the maximum latitude of the satellite image, Lat var is the pixel latitude of the satellite image, Res_line is the row height of the projection grid, and lines is the number of rows of the projection grid; among them, Lon min is the minimum longitude of the satellite image, Lon var is the pixel longitude of the satellite image, Res_pixel is the column width of the projection grid, and cols is the number of columns of the projection grid;

[0110] The mapping relationship between the geographic coordinates of the original satellite image and the row and column numbers of the new projected grid image is established through equations (9) and (10), and then the geometric correction and reprojection of the satellite image can be completed using gdalWarp.

[0111] 5. Generation of multi-channel raster images

[0112] The range of the FY-3D MERSI image data cropping area set by the present invention is 104.0E to 112.0E and 26.0N to 34.0N, covering parts of the upper reaches of the Yangtze River in Chongqing and some areas of neighboring provinces. Through the GDAL raster data processing technology, auxiliary data sets such as the observation angles of the sun and satellite, the reflectance in the visible and near-infrared bands, and the radiance temperature data set in the thermal infrared band are recombined together to construct three-dimensional data sets with a resolution of 1 km and a size of 29 * 800 * 800 and a resolution of 250 m and a size of 10 * 3200 * 3200, and multi-channel raster image data with two spatial resolutions of 1 km and 250 m are processed and output. Table 3 shows the information of the processed output image files, where PRJ indicates that reprojection processing has been performed, L2 is the data processing level, YYYYMMDDHHmm is the satellite orbit data reception date, GLL is the equidistant longitude and latitude projection, and information such as the satellite, payload, and spatial resolution can also be found in the file.

[0113] Table 3 TIF image product information

[0114]

[0115] 6. Data automatic processing technology

[0116] Regarding the automatic reprocessing technology of FY-3D MERSI L1 data, it is first necessary to realize the automatic detection of L1 data and data integrity check. This step is mainly based on the Linux system network file transfer protocol (NFS) and scheduled task (Crond) service, so as to realize the multi-frequency automatic detection of FY-3D MERSI L1B data. The automatic data detection process is shown in the appendix Figure 1 . When there is an update of L1B data in the direct receiving station, the program deployed on the platform side of the direct receiving station will detect it regularly and generate a data identification file. When the automatic detection module scans the new data identification file, it will check the integrity of the data file through logical conditions such as file size and number of files, and automatically read the data information of the data file.

[0117] The automatic reprocessing technology of FY3D MERSI L1 data can be closely connected with the provincial direct receiving station receiving and preprocessing platform. Using the Python development platform, combined with the HDF5 file and GDAL raster data processing modules, the data detection and reading, radiometric calibration, reflectance / brightness temperature calculation, equal latitude and longitude projection, geometric correction, regional cropping and raster image product generation and other processing modules are coupled and intensively managed, realizing an asynchronous condition triggering operation mode between modules, and improving the intelligent level of the reprocessing of FY-3D MERSIL1B-level data. The data automatic reprocessing technology process is shown in Figure 2 , the input end is the FY-3D MERSIL1B transit orbit data, and the output is the reprocessed FY-3D MERSI multi-channel raster image.

[0118] Embodiment

[0119] In order to further verify the effectiveness of the automatic processing method proposed by the present invention, the FY-3D MERSI L1B data at 13:33 on August 14, 2020 with complete reception was selected. At this time point, the central and western parts of Chongqing were mainly clear skies, with less cloud cover, and the image data in the visible light band and thermal infrared band were relatively clear, which was convenient for accurately analyzing the image quality. The image quality of the method of the present invention was compared and analyzed.

[0120] 1. Data volume before and after processing

[0121] The average data sizes before and after being processed by the method of the present invention are shown in Table 4. The data volume of the processed multi-channel TIF image as a whole is smaller than that before processing. The data volume of the 1km spatial resolution data decreased by about 117M after processing, and the data volume of the 250m spatial resolution decreased by about 206M after processing. In order to facilitate technicians to quickly read, the present invention did not convert the data to integer type for binary storage like the existing original processing method. The spectral band data is stored in real floating-point type. Generally, the data volume after processing does not decrease significantly compared with the original processing method, but it is convenient for users to read and call.

[0122] Table 4 Data size before and after processing (unit: MB)

[0123] Table 1 Data size before and after processing(unit:MB)

[0124]

[0125] 2. Data processing takes time

[0126] The time consumption and time efficiency of FY-3D MERSI L1 data automatic reprocessing are shown in the figure below. Figure 3 It can be seen that the processing time of the present invention ranges from 79 to 128 seconds, with an average running time of 103.14 seconds. Due to the actual performance of the server and network fluctuation factors, the data processing time shows certain fluctuations.

[0127] In order to analyze the time-consuming of the reprocessing process in depth, the time-consuming of 5 data processing modules is compared (Table 5). In the original processing mode, the satellite data files need to be manually detected and manually transmitted, which is time-consuming. The present invention uses the network file transfer protocol NFS to establish a real-time transmission link with the FY-3D storage device, realizes data automatic detection and fast reading capabilities, and the time-consuming is significantly reduced. And the time-consuming of other processing nodes of the present invention is also substantially better than the original processing mode. On the whole, the time-consuming of the present invention is reduced by about 111 seconds compared with the original processing mode, which improves the data reprocessing efficiency to a certain extent.

[0128] Table 5 Average time consumption of main processing modules of original processing method and new processing technology (unit: seconds)

[0129] Table 7 the main models average runtime between old method and new technology (unit: seconds)

[0130]

[0131] 3. Image quality comparison

[0132] The image quality comparison and analysis was carried out using the FY-3D MERSI L1B data at 13:33 on August 14, 2020 as an example. At that time, the central and western parts of Chongqing were mainly sunny with few clouds overhead, and the image data of the visible light channel and thermal infrared channel were relatively clear, making it easier to compare image quality.

[0133] Band13 (0.709μm) is the visible light band, which can be used for the calculation of the vegetation index NDVI. Band23 (8.550μm) is the thermal infrared radiation band, which can be used for the retrieval of atmospheric water vapor. In this invention, these two bands are selected to conduct an intuitive comparison of the remote sensing image quality. From the spatial distribution of the two images of Band13 (0.709μm) of the original processing method and the technology of this invention (as Figure 4 shown), the distribution of the observed clouds is concentrated in the eastern region of Chongqing. The location of the main urban area of Chongqing is clearly visible, and the geographical location of the main urban area is accurately specified. After the correction of the solar altitude angle by the new processing method, the image is brighter and the overall image is clearer. From the comparison of the thermal infrared band Band23 (8.550μm) images of the original processing method and the technology of this invention (as Figure 5 shown), affected by the mountainous terrain in the southwestern region, it is found that the variation range of the radiance brightness temperature of the image channels is obvious. The variation of the brightness temperature of the remote sensing image in the central and western regions of Chongqing is basically the same. In the central and western regions of Chongqing, the radiance brightness temperature of the infrared channel is basically the brightness temperature of the ground object, and the brightness temperature value is on the high side.

[0134] By calculating the mean absolute error MAE and root mean square error RMSE of the pixels of the spectral reflectance solar light band RSB and the thermal radiation band TEB, the image data errors of the two processing methods are further analyzed. From the MAE and RMSE of the RSB band of visible light and infrared from Band1 to Band19 (as Figure 6 shown), it is found that the MAE and RMSE of Band5 are the smallest, and the spectral band errors of Band8 - Band15 are relatively large. Generally, the reflectance MAE of the RSB band does not exceed 0.07, and the RMSE is within 0.08. From the MAE and RMSE of the brightness temperature of the long-wave thermal infrared radiation TEB band from Band20 to Band25 (as Figure 7 shown), it is found that the radiation brightness temperature error is generally small, its MAE does not exceed 0.53, and the RMSE does not exceed 0.8. Among them, the error of Band22 is the smallest, below 0.01. Due to the loss of data precision and the difference in interpolation algorithms during the calculation and processing, there are certain errors in the image data between the new processing technology and the original processing method, but the errors of each spectral band are controlled within an acceptable range. Here, two bands, Band17 and Band25, are selected to fit the results of the two methods of the original processing method and the new processing technology, and it is found that there is a linear correlation, and the correlation coefficients are all around 0.98 (as Figure 8 shown).

[0135] This invention can also output band raster remote sensing images with a spatial resolution of 250m. Taking the image data of the red light of Band3 (0.650μm) and the thermal infrared radiation band of Band24 (10.8μm) as an example, the quality of the image data with a spatial resolution of 250m is analyzed. Figure 9It is the reflectance image of visible light with a spatial resolution of 250m and the radiation brightness temperature image of the thermal infrared band. Generally speaking, the texture of the visible light image is finer, the accuracy of cloud amount and cloud spatial distribution is higher, and the geographical positions of mountains and rivers in the main urban area of Chongqing are accurately pointed. Through the variance and standard deviation of the two bands (Table 6), the standard deviation of band3 is 0.256, and the standard deviation of band24 is 14.134. Due to the ineffective removal of the influence of the atmosphere and clouds, the difference between the maximum and minimum values of the data is large. However, the standard deviation distribution of the radiation brightness temperature in the infrared band can basically reflect the characteristics of large topographic changes in Chongqing and its surrounding mountains.

[0136] Table 6 Data information of band3 and band24

[0137]

[0138] The content not described in detail in this specification belongs to the prior art well-known to those skilled in the art. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. An automatic reprocessing method for FY-3D MERSI L1B data, characterized in that, It includes the following steps: Step 1: Receive and preprocess the FY-3D MERSI L1B orbital observation data in the transit area in real time through the Fengyun-3 provincial direct receiving stations; Step 2: Automatically detect the newly received FY-3D MERSI L1B data, and perform automatic integrity check and reading; Step 3: Analyze the attributes of the inspected observation data to obtain the metadata information of the spectral band pixels of the data file; perform feature analysis on the metadata information to obtain the observation angle, correction coefficient, and geographic coordinate information; Step 4: Calculate the atmospheric apparent reflectance of the spectral band channels and the brightness temperature of the thermal radiation band based on the obtained observation angle and correction coefficient; Step 5: Based on the obtained pixel geographic coordinates, set the affine parameters to perform equi-rectangular grid projection and geometric correction on the image data; Step 6: Crop the monitoring image data according to the longitude and latitude range of the business observation area, and use the GDAL raster data processing technology to generate a TIF format regional image product; The specific steps of Step 6 include: Step 61: Obtain the longitude and latitude coordinate range information of the FY-3D-MERSI image area to be monitored, and crop the FY-3D MERSI image data; Step 62: In the cropped area, construct a three-dimensional raster dataset based on the GDAL raster data processing technology; Step 64: Write the obtained solar and satellite observation angles of pixels, the calculated atmospheric apparent reflectance of the spectral band channels, and the pixel brightness temperature of the thermal infrared band into the three-dimensional raster dataset in sequence according to the row and column dimensions of the image data; Step 65: Output multi-channel raster image data with two spatial resolutions of 1 km and 250 m.

2. The automatic reprocessing method for FY-3D MERSI L1B data according to claim 1, wherein The specific operation steps of Step 2 include: Step 21: Regularly scan the received FY-3D MERSI L1B data files in the FY-3D MERSI L1B orbital time series observation data directory, identify them, and generate identification files; Step 22: Determine whether there is a new identification file. If not, return to Step 21 to continue scanning; if so, enter Step 23 for integrity check; Step 23: Automatically check the integrity of the data identification file based on the file size, the number of data files, and the time attribute of the identification file by using the scheduled task service. If the file is complete, the data can be read; if the file is incomplete, return to Step 21.

3. The automatic reprocessing method for FY-3D MERSI L1B data according to claim 2, characterized in that, The specific steps of Step 4 include: Step 41: Calculate the atmospheric apparent reflectance Step 411: After detecting outliers and missing values for the digital quantization values of each band pixel, perform radiometric calibration on the reasonable pixel DN values through Equation (1) to obtain the calibrated radiometric quantization after radiometric calibration dn , that is: (1) where slope and intercept are the gain and offset coefficients required for radiometric calibration, which can be obtained from the dataset attributes; DN is the radiometric digital quantization value of the channel pixel, and dn is the radiometric value after radiometric calibration; Step 412: Calculate the corrected pixel radiation R of this channel band according to the correction coefficient: (2) where Cal2, Cal1, and Cal0 are the correction coefficients; Step 413: Carry out solar altitude angle correction on the average spectral radiation at the top of the atmosphere and calculate the atmospheric apparent reflectance of the spectral band channels through equation (3) : (3) where E m is the average spectral irradiance at the top of the atmosphere, d is the astronomical unit, and θ is the solar zenith angle; L sensor is the radiance at the entrance pupil of the satellite sensor; Step 414: Calculate L using Equation (4) sensor , substitute Equation (4) into Equation (3) to obtain Equation (5), and quickly calculate the atmospheric apparent reflectance using Equation (5): (4) （5）； Step 42: Calculate the brightness temperature Step 421: After detecting outliers and missing values in the thermal radiation band, perform radiometric calibration on the infrared channel radiation L, and calculate the radiance L of the band after radiometric calibration λ : （6） Among them, L is the radiance quantization value of the uncalibrated channel, slope is the calibration gain, intercept is the calibration offset, and L λ is the calibrated radiance; Step 422: Calculate the equivalent blackbody brightness temperature T in the thermal infrared band according to the Planck function conversion formula (7) f : (7) Among them, C1 and C2 are the integrated calculation coefficients, h = 6.62606876e-34 J.s is the Planck constant, c = 2.99792458e+8 m / s is the speed of light, k = 1.3806503e-23 J / K is the Boltzmann constant, is the wave number of this channel, L λ is the radiance after radiometric calibration; Step 423: Correct T through equation (8) f to obtain the brightness temperature T b observed in this band: (8) Among them, TBB_a and TBB_b are the brightness temperature correction coefficients, T f is the equivalent blackbody brightness temperature of the thermal infrared channel, T b is the band brightness temperature, which is the band brightness temperature observed by the satellite remote sensing observation instrument.

4. The automatic reprocessing method for FY-3D MERSI L1B data according to claim 3, wherein The specific steps of Step 5 include: Step 51: Obtain the maximum latitude Lat, minimum longitude Lon max of the images with different resolutions, as well as the row height Res_Line and column width Res_pixel of the image pixels according to the read geographical coordinates of the pixels with different spatial resolutions. max , minimum longitude Lon min and the row height Res_Line and column width Res_pixel of the image pixels; Step 52: Set the satellite image affine parameters according to Lat max , Lon min , Res_Line, and Res_pixel, construct an equatorial and meridional projection grid, and establish the relationship between the geographical coordinates of the original image pixels and the row and column numbers of the projection grid through equations (9)-(10): (9) (10) Among them, Lat max is the maximum latitude of the satellite image, Lat var is the pixel latitude of the satellite image, Res_line is the row height of the projection grid, and lines is the number of rows of the projection grid; Among them, Lon min is the minimum longitude of the satellite image, Lon var is the pixel longitude of the satellite image, Res_pixel is the column width of the projection grid, and cols is the number of columns of the projection grid; Step 53: Perform geometric correction and reprojection of the satellite image through gdalWarp based on the relationship between the original image pixel geographic coordinates and the projection grid row and column numbers constructed.

Images