A method for inversion of polynyas based on multi-channel satellite data

Abstract

The present invention discloses a method for inverting polynyas based on multi-channel satellite data, belonging to the field of satellite image processing technology. The method comprises the following steps: obtaining reflectivity anomalies and initial values ​​of reflectivity anomaly thresholds for visible light radiation brightness values, and obtaining brightness temperature anomaly values ​​(BTA) and initial values ​​of brightness temperature anomaly thresholds for thermal infrared radiation brightness values; secondly, obtaining a pixel credibility matrix; and thirdly, using an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification to adjust the brightness temperature anomaly thresholds and reflectivity anomaly thresholds, calculate the brightness temperature anomaly membership and reflectivity anomaly membership, and determine the comprehensive probability of a pixel belonging to a polynya based on the memberships and the pixel credibility matrix. The method of the present invention utilizes the advantages of visible light data in identifying object boundaries and textures and the advantages of thermal infrared data in distinguishing objects of different temperatures. It can better adapt to the dynamic complexity of polar environments and improve the accuracy of polynya inversion.

Description

Technical Field

[0001] The present invention belongs to the technical field of image processing, and in particular relates to a polynya waterway inversion method based on multi-channel satellite data. Background Art

[0002] Polar (Antarctic and Arctic) polynyas are key elements of the polar sea ice cover, playing a vital role in the polar climate system, ecological environment, and human activities. With the ongoing changes in polar sea ice, the importance of studying polynyas has become increasingly prominent. Polynyas are linear or narrow cracks in the sea ice cover. They form primarily due to external forces such as wind and ocean currents, which cause uneven stress on the sea ice and cause it to break. These external forces cause the sea ice to move and deform. When the stress exceeds the ice's inherent strength, polynyas form. In winter, the large temperature difference between the ocean and atmosphere causes the newly formed polynyas to rapidly lose heat, quickly forming thin ice, which remains distinct from the surrounding thick ice. This difference is reflected not only in its physical structure, such as the thickness and density of thin ice, but also in its thermodynamic properties, such as its temperature conductivity and heat capacity, which significantly influence energy exchange in the polar regions.

[0003] Since the late 1970s, remote sensing images acquired by satellite sensors have been used to detect polynyas in the polar regions. Willmes et al. (2015) used MODIS thermal infrared imagery to detect polynyas using a background threshold difference method. Hoffman et al. (2019) utilized the MODIS thermal infrared band to design a multi-step method for effectively detecting and characterizing polynyas. Qu et al. (2021) proposed an improved algorithm for inverting the Beaufort Polynya using Terra / MODIS thermal infrared temperature anomaly images. Hoffman et al. (2021) combined thermal infrared data from MODIS and VIIRS satellites and used artificial intelligence (AI) technology to explore the potential and challenges of AI in sea ice research and proposed solutions. Qiu et al. (2023) used high-resolution infrared imagery provided by the Thermal Infrared Spectrometer (TIS) on the Sustainable Development Science Satellite-1 (SDGSAT-1) to detect polynyas in the Arctic. This detection method, based on a single near-infrared channel, has certain limitations. Newly formed thin ice channels (less than 5 cm thick) have a low near-infrared emissivity compared to channels not covered by ice, resulting in poor feature separability. Furthermore, the anisotropic reflectivity caused by snow cover obscures channel edge features, increasing the missed detection rate for channels less than 100 m wide. These issues make a single near-infrared channel incapable of meeting full-scene detection requirements.

[0004] Key et al. (1993) used thermal infrared and visible light imagery from the Advanced Very High Resolution Radiometer (AVHRR) to analyze the presence of polynyas, attempting to determine their widths and providing basic data and methodological references for subsequent research. Their method, which employed a fixed threshold for single-parameter segmentation, was limited by the dynamic complexity of the polar environment. This method suffered from a high misjudgment rate in the ice-water transition zone, making it incapable of adapting to the dynamic polar environment.

[0005] In addition, existing technologies only simulate atmospheric effects through the radiation transfer model (LOWTRAN), which has cloud interference problems and is unable to suppress cloud interference in real time. Summary of the Invention

[0006] In order to solve the technical problem that the existing method of single-parameter segmentation by setting a fixed threshold is used in the inversion of ice channels in polar regions, which is limited by the dynamic complexity of the polar environment and has a high misjudgment rate, the present invention proposes an ice channel inversion method based on multi-channel satellite data to solve the above problem.

[0007] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0008] A method for inverting polynyas based on multi-channel satellite data, comprising:

[0009] Step 1: Data preprocessing step, including: obtaining the reflectance anomaly value of the visible light radiation brightness value , and obtain the brightness temperature anomaly value of the thermal infrared radiation brightness value , obtain the initial values ​​of the reflectivity anomaly threshold and the brightness temperature anomaly threshold of the polynya in the satellite data;

[0010] Step 2: Obtain pixel credibility step, calculate the credibility of each pixel in the satellite data, and form a pixel credibility matrix ;

[0011] Step three, using an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification, adjust the brightness temperature anomaly threshold and the reflectivity anomaly threshold, the thresholds include an upper threshold and a lower threshold, calculate the brightness temperature anomaly membership and the reflectivity anomaly membership, the brightness temperature anomaly membership is the membership of the pixel corresponding to the brightness temperature anomaly to the polynya, the reflectivity anomaly membership is the membership of the pixel corresponding to the reflectivity anomaly to the polynya, and judge the comprehensive probability that the pixel belongs to the polynya based on each membership and the pixel credibility matrix. The iterative dual-threshold adaptive optimization algorithm based on fuzzy classification includes:

[0012] Calculate the brightness temperature anomaly center of the current target class and reflectivity anomaly center The target class is polynyas, the brightness temperature anomaly center is the center value of the brightness temperature anomaly threshold, and the reflectivity anomaly center is the center value of the reflectivity anomaly threshold. According to BTA and Calculate the membership of brightness temperature anomalies , according to RA and Calculate reflectivity outlier membership , k represents the number of iterations, j represents the thermal infrared channel number, and i represents the visible light channel number;

[0013] Classify the pixels according to the brightness temperature anomaly upper threshold and the brightness temperature anomaly lower threshold, calculate the inter-class mean of each class, which is the brightness temperature inter-class mean, and update the brightness temperature anomaly upper threshold and the brightness temperature anomaly lower threshold according to the brightness temperature inter-class mean; classify the pixels according to the reflectivity anomaly upper threshold and the reflectivity anomaly lower threshold, calculate the inter-class mean of each class, which is the reflectivity inter-class mean, and update the reflectivity anomaly upper threshold and the reflectivity anomaly lower threshold according to the reflectivity inter-class mean;

[0014] The iteration is terminated when the brightness temperature threshold change and reflectivity threshold change of two adjacent iterations simultaneously meet the following conditions:

[0015] ;

[0016] ;

[0017] Based on the membership degree obtained in the last iteration and the pixel credibility matrix, the comprehensive probability that the pixel belongs to the polynya is calculated. .

[0018] Compared with the existing technology, the advantages and positive effects of the present invention are: the polynya waterway inversion method based on multi-channel satellite data of the present invention, since the thermal infrared channel of the satellite data can better reflect the thermal radiation characteristics of the ground objects, and the visible light channel of the satellite data can better reflect the detailed information of the polynya waterway, utilizes the characteristics of high brightness temperature and low reflectivity of polynya waterways in the polar region, and searches for the brightness temperature anomaly data of the thermal infrared channel and the reflectivity anomaly values ​​of the visible light channel as the processing objects for inverting the polynya waterway, combines the two, utilizes the advantages of visible light data in identifying the boundaries and textures of ground objects, and utilizes the advantages of thermal infrared data in distinguishing ground objects of different temperatures, thereby achieving complementary advantages.

[0019] When judging brightness temperature anomaly data and reflectivity anomaly data, the pixels are classified according to the current brightness temperature anomaly threshold and reflectivity anomaly threshold, and the inter-class mean of each class is calculated. The inter-class mean can reflect the overall distribution of each class of anomalies. The reflectivity anomaly threshold and brightness temperature anomaly threshold are recalculated based on the inter-class mean, so that this scheme can better adapt to the dynamic complexity of the polar environment, reduce the misjudgment rate, and improve the inversion accuracy of the interglacial waterway.

[0020] Other features and advantages of the present invention will become more apparent after reading the detailed description of the embodiments of the present invention in conjunction with the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 is a flow chart of an embodiment of a polynya waterway inversion method based on multi-channel satellite data proposed by the present invention;

[0022] Figure 2a It is a 10.8 μm visible light channel data image in one embodiment of the polynya waterway inversion method based on multi-channel satellite data proposed by the present invention;

[0023] Figure 2b This is an embodiment of the method for inverting polynya waterways based on multi-channel satellite data proposed by the present invention. Visible light channel data image;

[0024] Figure 2c It is a visible light channel data image in one embodiment of the polynya waterway inversion method based on multi-channel satellite data proposed by the present invention;

[0025] Figure 3a An inversion image of a polynya using a small-scale sliding window in an embodiment of the polynya inversion method based on multi-channel satellite data proposed by the present invention;

[0026] Figure 3b An inversion image of a polynya using a mesoscale sliding window in an embodiment of the polynya inversion method based on multi-channel satellite data proposed by the present invention;

[0027] Figure 3c An inversion image of a polynya using a large-scale sliding window in an embodiment of the polynya inversion method based on multi-channel satellite data proposed by the present invention;

[0028] Figure 4 In one embodiment of the method for inverting polynya waterways based on multi-channel satellite data proposed by the present invention, Figure 3a-3c The fused inversion image of the polynya;

[0029] Figure 5The present invention provides a sea ice density image in one embodiment of a polynya inversion method based on multi-channel satellite data;

[0030] Figure 6 This is an image of an Arctic polynya inverted in an embodiment of the polynya inversion method based on multi-channel satellite data proposed by the present invention. DETAILED DESCRIPTION

[0031] The specific embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0032] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0033] It should be noted that, in the description of the present invention, the terms "upper", "lower", "left", "right", "vertical", "horizontal", "inside", "outside" and the like indicating directions or positional relationships are based on the directions or positional relationships shown in the accompanying drawings. This is merely for the convenience of description and does not indicate or imply that the device or element must have a specific orientation, be constructed and operated in a specific orientation. Therefore, it cannot be understood as a limitation on the present invention. In addition, the terms "first" and "second" are used for descriptive purposes only and cannot be understood as indicating or implying relative importance. In the description of the present invention, the meaning of "multiple" is two or more, unless otherwise clearly and specifically defined.

[0034] In the present invention, unless otherwise expressly specified or limited, the terms "mounted," "connected," "connect," "fixed," etc. should be understood broadly. For example, they may refer to fixed connection, detachable connection, or integration; mechanical connection or electrical connection; direct connection or indirect connection through an intermediate medium; internal communication between two components or interaction between two components. Those skilled in the art will understand the specific meanings of the above terms in the present invention based on specific circumstances.

[0035] Example 1, see Figure 1 As shown, this embodiment proposes a polynya inversion method based on multi-channel satellite data, including:

[0036] Step 1, data preprocessing step, includes: obtaining the reflectivity anomaly value RA of the visible light radiation brightness value and the brightness temperature anomaly value BTA of the thermal infrared radiation brightness value, and obtaining the initial value of the reflectivity anomaly threshold and the initial value of the brightness temperature anomaly threshold of the polynya in the satellite data.

[0037] In this embodiment, multiple thermal infrared channels and visible light channels of the MERSI-II sensor of the FY-3D satellite are used to perform polynya inversion.

[0038] Step 2: Obtain pixel credibility step, calculate the credibility of each pixel in the satellite data, and form a pixel credibility matrix .

[0039] Step 3: Use an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification to adjust the brightness temperature anomaly threshold and reflectivity anomaly threshold. The thresholds include upper and lower thresholds. Calculate the brightness temperature anomaly membership and reflectivity anomaly membership. The brightness temperature anomaly membership is the membership of the pixel corresponding to the brightness temperature anomaly to the polynya. The reflectivity anomaly membership is the membership of the pixel corresponding to the reflectivity anomaly to the polynya. The comprehensive probability that the pixel belongs to the polynya is determined based on each membership and the pixel credibility matrix. The iterative dual-threshold adaptive optimization algorithm based on fuzzy classification includes:

[0040] Calculate the brightness temperature anomaly center of the current target class and reflectivity anomaly center , the target class is polynya, the brightness temperature anomaly center is the center value of the brightness temperature anomaly threshold, and the reflectivity anomaly center is the center value of the reflectivity anomaly threshold. According to BTA and Calculate the membership of brightness temperature anomalies , according to RA and Calculate reflectivity outlier membership , k represents the number of iterations, j represents the thermal infrared channel number, and i represents the visible light channel number.

[0041] The pixels are classified according to the upper and lower thresholds of brightness temperature anomaly, and the inter-class mean of each category is calculated as the brightness temperature inter-class mean. The upper and lower thresholds of brightness temperature anomaly are updated according to the brightness temperature inter-class mean. The pixels are classified according to the upper and lower thresholds of reflectivity anomaly, and the inter-class mean of each category is calculated as the reflectivity inter-class mean. The upper and lower thresholds of reflectivity anomaly are updated according to the reflectivity inter-class mean.

[0042] The iteration is terminated when the brightness temperature threshold change and reflectivity threshold change of two adjacent iterations simultaneously meet the following conditions:

[0043] ;

[0044] .

[0045] Based on the membership degree obtained in the last iteration and the pixel credibility matrix, the comprehensive probability that the pixel belongs to the polynya is calculated. Comprehensive probability That is, it reflects the probability that the pixel belongs to the polynya, and realizes the inversion of the polynya.

[0046] The polynya channel inversion method based on multi-channel satellite data in this embodiment is based on the fact that the thermal infrared channel of satellite data can better reflect the thermal radiation characteristics of ground objects, such as Figure 2a 、 Figure 2b As shown in Figure 2, the visible light channel of satellite data can better reflect the details of polynyas, such as Figure 2c As shown in the figure, the characteristics of high brightness temperature and low reflectivity of ice channels in polar regions are utilized. By looking for brightness temperature anomaly data of the thermal infrared channel and reflectivity anomaly values ​​of the visible light channel as processing objects for inverting ice channels, the two are combined to utilize the advantages of visible light data in identifying the boundaries and textures of land objects, and the strengths of thermal infrared data in distinguishing land objects of different temperatures, thereby achieving complementary advantages.

[0047] When judging brightness temperature anomaly data and reflectivity anomaly data, the pixels are classified according to the current brightness temperature anomaly threshold and reflectivity anomaly threshold, and the inter-class mean of each class is calculated. The inter-class mean can reflect the overall distribution of each class of anomalies. The reflectivity anomaly threshold and brightness temperature anomaly threshold are recalculated based on the inter-class mean, so that this scheme can better adapt to the dynamic complexity of the polar environment, reduce the misjudgment rate, and improve the inversion accuracy of the interglacial waterway.

[0048] In some embodiments, in step one, the method for obtaining RA and BTA includes: obtaining visible light channel data and thermal infrared channel data of satellite data, performing radiometric calibration on the visible light channel data to obtain a visible light radiation brightness value, and obtaining a reflectivity anomaly value RA of the visible light radiation brightness value; performing radiometric calibration on the thermal infrared channel data to obtain a thermal infrared radiation brightness value, converting the thermal infrared radiation brightness value into the brightness temperature of an equivalent black body, and calculating the brightness temperature anomaly value BTA; obtaining an initial value of the reflectivity anomaly threshold based on RA, and obtaining an initial value of the brightness temperature anomaly threshold based on BTA.

[0049] The data of the visible light channel is processed. The radiometric calibration technology is used to convert the digital signal DN measured by the sensor into a physically meaningful radiometric brightness value Ref according to the radiation response characteristics of the sensor and the calibration parameters. Then, a mathematical model of the solar zenith angle and the radiometric brightness value Ref is established. The change in the radiometric brightness value Ref caused by the change in the solar zenith angle is corrected by this model, that is, the error caused by the change in the solar zenith angle is eliminated to ensure data consistency. Then, the reflectivity anomaly value RA is calculated based on the FY-3D / MERSI-II radiometric brightness value Ref, and then the reflectivity data anomaly membership is calculated.

[0050] In this embodiment, a quantitative relationship is established between the digital signal measured by the FY-3D / MERSI-II sensor and the actual radiation characteristics of the ground object, so as to accurately obtain the radiation energy information of the ground object.

[0051] right The reflectance data of the three visible light channels are calibrated using a specific formula:

[0052] ;

[0053] .

[0054] Where, Normalize the i-th visible light channel The value is the brightness value of each pixel in the remote sensing image, which records the intensity of electromagnetic waves reflected or radiated by the ground objects and is expressed as a unitless integer value.

[0055] are the calibration coefficients of the corresponding channels (corresponding to the 1st, 2nd, and 3rd columns respectively), and For the corresponding channel datasets EV_1KM_RefSB and EV_250_Aggr.1KM_RefSB attributes, is the reflectivity of the i-th channel.

[0056] Because the solar altitude determines the angle between sunlight and the ground, the smaller the solar altitude, the longer the path of light through the atmosphere, the greater the atmospheric scattering and absorption of light, the less solar radiation reaches the ground, and the less radiation reflected back to the sensor by objects. Conversely, the larger the solar altitude, the more radiation reflected back to the sensor by objects.

[0057] Process the data of the thermal infrared channel and calculate the radiance of the two thermal infrared channel data respectively , and then according to the Planck inverse transform formula and brightness temperature correction coefficient Convert to equivalent blackbody brightness temperature , calculate brightness temperature anomaly , and then calculate the brightness temperature anomaly membership.

[0058] Specifically, the radiance value of the jth thermal infrared channel measured by the FY-3D / MERSI-II sensor is And the calibration parameters given by the sensor and Perform calibration calculations to obtain the radiometry of the jth thermal infrared channel (mW / (m2 cm sr)):

[0059] .

[0060] in, , .

[0061] Then, using the inverse Planck transform formula, the radiance Convert to equivalent blackbody brightness temperature , its specific form is:

[0062] .

[0063] Then use the channel brightness temperature correction coefficient and Will Convert to channel blackbody brightness temperature and , the formula is as follows:

[0064] .

[0065] Among them, when The brightness temperature correction coefficient of the channel is , ,when The brightness temperature correction coefficient of the channel is , .

[0066] Then the brightness temperature of the two thermal infrared channels Convert to brightness temperature anomaly image , its specific form is:

[0067] .

[0068] in, is a sliding window, is the brightness temperature anomaly image, For the corrected thermal infrared channel brightness temperature data, for The median image of is the pixel location.

[0069] In some embodiments, in step 1, the method of obtaining the initial value of the reflectivity anomaly threshold and the initial value of the brightness temperature anomaly threshold includes:

[0070] Sort the data in RA or BTA from small to large;

[0071] Remove the first p1% of the data and the last p2% of the data.

[0072] The minimum value in the retained data is the initial value of the reflectivity anomaly lower threshold or the initial value of the brightness temperature anomaly lower threshold, and the maximum value in the retained data is the initial value of the reflectivity anomaly upper threshold or the initial value of the brightness temperature anomaly upper threshold.

[0073] In some embodiments, p1% and p2% can be set to 5% or other numbers, so that the retained data is the data ranked 5%-95%.

[0074] In some embodiments, step one further includes the step of performing zenith angle correction on the visible light radiation brightness value and the thermal infrared radiation brightness value.

[0075] In the same remote sensing image, if the solar altitude angles in different areas are different, it will cause differences in the grayscale values ​​of different areas on the image, which requires correction.

[0076] First, the solar zenith angle Convert to radians:

[0077] .

[0078] The relationship between the solar altitude angle and the reflectivity of the ground is used to establish the correction formula:

[0079] .

[0080] in is the solar zenith angle, is the zenith angle converted to radians, This is the reflectivity of the i-th visible light channel after solar zenith angle correction. This correction effectively eliminates image reflectivity errors caused by the solar zenith angle, ensuring that the reflectivity data accurately reflects the optical properties of the polynya and its surrounding environment.

[0081] Similarly, the above method can be used to correct the solar zenith angle of the thermal infrared channel.

[0082] The vast snow cover of the Arctic provides excellent conditions for NDSI to demonstrate its superior ability to distinguish between snow and clouds. Due to the relatively stable acquisition of FY3D / MERSI-II data in the Arctic, NDSI data is of high quality. It provides accurate ground feature background information for cloud detection, especially when distinguishing clouds above snow. It clearly delineates cloud boundaries and accurately identifies cloud presence, providing a reliable basis for subsequent cloud detection and analysis.

[0083] To minimize cloud interference with polynya channel detection, this method combines the Normalized Difference Snow Cover Index (NDSI), the Cloud Detection Product (CLW), the High Cloud Cover Product (HCCP), and the Total Cloud Cover Product (TCP) to assess the reliability of pixel retrieval. NDSI data can be used to distinguish snow from other surface features and assist in determining cloud presence; the Cloud Detection Product (CLW) directly identifies cloud-covered areas; and the High Cloud Cover Product (HCCP) and Total Cloud Cover Product (TCP) provide quantitative information on cloud cover. By comprehensively analyzing these data, pixel reliability is assessed to minimize the impact of cloud cover on polynya channel detection. Snow Cover Index data can be used to distinguish snow from other surface features and assist in determining cloud presence; the Cloud Detection Product (CLW) directly identifies cloud-covered areas; and the Total Cloud Cover Product (TCP) provides quantitative information on cloud cover.

[0084] Based on the polynya pixel membership and pixel inversion reliability, the probability of each pixel belonging to the polynya category is calculated, thereby achieving a refined inversion of polynyas and improving the accuracy and reliability of detection. In some embodiments, the reliability of each pixel in the satellite data is calculated based on the Normalized Snow Cover Index (NDSI), the cloud detection product (CLW), and the total cloud cover product (TCP) in step 2, including:

[0085] ;

[0086] Pixel credibility matrix Make corrections:

[0087] when When , the credibility of the corresponding pixel is 0;

[0088] when When , the credibility of the corresponding pixel is 0;

[0089] in, is the preset threshold.

[0090] The calculation formula of the normalized snow cover index data NDSI is as follows:

[0091] .

[0092] in 0.65 for the 3rd channel of the MERSI-II sensor data, 1.64 for channel 6 in the MERSI-II sensor data.

[0093] Due to the omission of total cloud cover products, NDSI and CLM are added for supplementary correction.

[0094] In some embodiments, step three further includes performing data scaling processing on BTA and RA so that their value range is between [-1, 1], including:

[0095] Using sliding window Traverse the brightness temperature anomaly image or reflectance outlier image .

[0096] Sliding window Scaling of brightness temperature anomalies or reflectivity anomalies within the image:

[0097] .

[0098] .

[0099] in, For sliding window Internal brightness temperature anomaly, For sliding window The minimum brightness temperature anomaly value, For sliding window Maximum brightness temperature anomaly within the For sliding window Internal reflectivity anomaly, For sliding window The inner minimum reflectivity outlier, For sliding window The maximum reflectivity outlier within the .

[0100] After traversing the entire brightness temperature anomaly image, all Update brightness temperature anomalies , all Update reflectivity outliers .

[0101] In some embodiments, the brightness temperature anomaly center of the current target class and reflectivity anomaly center The calculation method is:

[0102] .

[0103] .

[0104] In some embodiments, the brightness temperature anomaly membership The calculation method is:

[0105] ;

[0106] Reflectivity outlier membership The calculation method is:

[0107] ;

[0108] in, is the brightness temperature width parameter of the jth thermal infrared channel, is the reflectivity width parameter of the i-th visible light channel, and n is a constant.

[0109] In some embodiments, in step 1, the reflectivity outlier image The methods for obtaining include:

[0110] Using sliding window Traverse the reflectance data of the i-th visible light channel;

[0111] ;

[0112] For sliding window Reflectivity data within for The median image of .

[0113] In some embodiments, in step 3, the method for calculating the brightness temperature inter-class mean is:

[0114] Pixels are divided into the following three categories based on the upper and lower thresholds of brightness temperature anomaly:

[0115] ;

[0116] ;

[0117] .

[0118] represents the brightness temperature value of pixel x of the jth thermal infrared channel at the kth iteration, represents the lower limit threshold of brightness temperature anomaly of the j-th thermal infrared channel at the k-th iteration, represents the upper threshold of brightness temperature anomaly of the j-th thermal infrared channel at the k-th iteration.

[0119] 、 as well as The inter-class means are:

[0120] ;

[0121] ;

[0122] .

[0123] Update the upper and lower thresholds for brightness temperature anomalies:

[0124] ;

[0125] ;

[0126] in, is a dynamic adjustment coefficient used to compensate for the influence of ice-water mixed pixels. In this embodiment .

[0127] renew :

[0128] .

[0129] in, , is the experience value, is a preset constant.

[0130] Similarly, the calculation method of the inter-class mean of reflectivity is:

[0131] Pixels are divided into the following three categories based on the upper and lower reflectivity anomaly thresholds:

[0132] ;

[0133] ;

[0134] .

[0135] Represents the reflectance value of pixel x of the i-th visible light channel at the k-th iteration, represents the lower limit threshold of the reflectivity anomaly of the i-th visible light channel at the k-th iteration, Represents the upper threshold of the reflectivity anomaly of the i-th visible light channel at the k-th iteration.

[0136] 、 as well as The inter-class means are:

[0137] ;

[0138] ;

[0139] .

[0140] Update the upper and lower thresholds for brightness temperature anomalies:

[0141] ;

[0142] .

[0143] renew :

[0144] .

[0145] in, , is the experience value, is a preset constant.

[0146] The iteration is terminated when the brightness temperature threshold change and reflectivity threshold change of two adjacent iterations simultaneously meet the following conditions:

[0147] ;

[0148] .

[0149] In some embodiments, .

[0150] In some embodiments, in step 3, the multi-source membership and the credibility weight are combined to calculate the comprehensive probability that the pixel belongs to the polynya:

[0151] The comprehensive probability that a pixel belongs to a polynya The calculation method is:

[0152] .

[0153] in, is the number of thermal infrared channels, is the number of thermal infrared channels.

[0154] In some embodiments, after step 3, the scale of the sliding window is changed, and steps 1 to 3 are performed again to calculate the comprehensive probability under the current scale sliding window, and obtain ,like Figure 3a-3c As shown, where q represents the scale number.

[0155] Will The probability value of is assigned to g, and the other probability values ​​are assigned to 0.

[0156] All pixels with a comprehensive probability value of g are marked as connected components to find multiple connected areas.

[0157] Calculate the width of each connected region .

[0158] Calculate the mean width of the polynya detected under the sliding window of each scale respectively and standard deviation .

[0159] according to as well as Calculate the weight of sliding windows of each scale :

[0160] .

[0161] Calculate the final probability P that a pixel belongs to a polynya:

[0162] .

[0163] Then The area is considered to be a polynya. Areas considered to be potential polynyas, such as Figure 4 shown.

[0164] like Figure 5 As shown, this is the sea ice density image on April 2, 2021. Figure 6 As shown, the method of this embodiment is used to Figure 5 Inverted image of polynyas in the Arctic.

[0165] In some embodiments, the width of the polynya is calculated , as follows:

[0166] First, the probability value of the polynya is calculated using the 3×3 neighborhood structure S. Mark the connected components and assign a label to each polynya, and we get tag objects, where .

[0167] For each connected object , calculate the number of pixels, the formula is as follows:

[0168] .

[0169] Right now Indicates all tags that meet equal The total number of points.

[0170] in, express The median coordinate is The connected area identifier to which the pixel belongs, for The total number of connected regions in , Indicates the The number of pixels in a connected region.

[0171] For each connected region , define its coordinate set:

[0172] .

[0173] Then split the continuous segments by row, and for each row y, extract the column coordinate sequence The difference between adjacent coordinates can be calculated:

[0174] .

[0175] when When , it indicates that the current segment is disconnected and the sum of all continuous segments of the current row needs to be calculated. .

[0176] Then we can get the average row width:

[0177] .

[0178] Where, Represents a set of consecutive segments The number of elements in Indicates a continuous segment length.

[0179] Similarly, we can get the average column height:

[0180] .

[0181] Taking into account the information in the row and column directions, the minimum value of the average row width and average column height is taken as the connected area. The final width:

[0182] .

[0183] From a physical perspective, different windows correspond to different spatial frequency responses. Small windows preserve high-frequency details and can reflect the fine structure of polynyas, while large windows capture low-frequency trends and reveal the overall characteristics and distribution trends of polynyas.

[0184] Multi-scale weighted feature fusion is based on the Gaussian distribution principle, using a Gaussian weighting function to weight features of different scales. In polynya channel inversion, windows of different sizes have different advantages and disadvantages for detecting polynyas. The Gaussian weighting function can rationally assign weights to features of different scales in the fusion results based on the difference between the theoretical detection width of the window and the median width of the target channel.

[0185] For a sliding window of scale q , can be calculated and the width of the polynya , then when the sliding window The number of scales is When z sliding windows are used, the width of the polynya can be calculated. The mean and standard deviation .

[0186] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by ordinary technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. A method for inverting polynyas based on multi-channel satellite data, characterized in that: include: Step 1: Data preprocessing step, including: obtaining the reflectance anomaly value of the visible light radiation brightness value , and obtain the brightness temperature anomaly value of the thermal infrared radiation brightness value , obtain the initial values ​​of the reflectivity anomaly threshold and the brightness temperature anomaly threshold of the polynya in the satellite data; Step 2: Obtain pixel credibility step, calculate the credibility of each pixel in the satellite data, and form a pixel credibility matrix ; Step three, using an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification, adjust the brightness temperature anomaly threshold and the reflectivity anomaly threshold, the thresholds include an upper threshold and a lower threshold, calculate the brightness temperature anomaly membership and the reflectivity anomaly membership, the brightness temperature anomaly membership is the membership of the pixel corresponding to the brightness temperature anomaly to the polynya, the reflectivity anomaly membership is the membership of the pixel corresponding to the reflectivity anomaly to the polynya, and judge the comprehensive probability that the pixel belongs to the polynya based on each membership and the pixel credibility matrix. The iterative dual-threshold adaptive optimization algorithm based on fuzzy classification includes: Calculate the brightness temperature anomaly center of the current target class and reflectivity anomaly center The target class is polynyas, the brightness temperature anomaly center is the center value of the brightness temperature anomaly threshold, and the reflectivity anomaly center is the center value of the reflectivity anomaly threshold. According to BTA and Calculate the membership of brightness temperature anomalies , according to RA and Calculate reflectivity outlier membership , k represents the number of iterations, j represents the thermal infrared channel number, and i represents the visible light channel number; Classify the pixels according to the brightness temperature anomaly upper threshold and the brightness temperature anomaly lower threshold, calculate the inter-class mean of each class, which is the brightness temperature inter-class mean, and update the brightness temperature anomaly upper threshold and the brightness temperature anomaly lower threshold according to the brightness temperature inter-class mean; classify the pixels according to the reflectivity anomaly upper threshold and the reflectivity anomaly lower threshold, calculate the inter-class mean of each class, which is the reflectivity inter-class mean, and update the reflectivity anomaly upper threshold and the reflectivity anomaly lower threshold according to the reflectivity inter-class mean; The iteration is terminated when the brightness temperature threshold change and reflectivity threshold change of two adjacent iterations simultaneously meet the following conditions: ； ； Based on the membership degree obtained in the last iteration and the pixel credibility matrix, the comprehensive probability that the pixel belongs to the polynya is calculated. .

2. The polynya channel inversion method based on multi-channel satellite data according to claim 1, characterized in that: In step one, as well as The acquisition method includes: obtaining visible light channel data and thermal infrared channel data of satellite data, performing radiation calibration on the visible light channel data to obtain a visible light radiation brightness value, and obtaining a reflectivity anomaly value of the visible light radiation brightness value. ; Perform radiation calibration on the thermal infrared channel data to obtain the thermal infrared radiation brightness value, convert the thermal infrared radiation brightness value into the brightness temperature of the equivalent black body, and calculate the brightness temperature anomaly value ;according to Get the initial value of the reflectivity anomaly threshold, and Get the initial value of the brightness temperature anomaly threshold.

3. The polynya channel inversion method based on multi-channel satellite data according to claim 2, characterized in that: In step 1, the method for obtaining the initial value of the reflectivity anomaly threshold and the initial value of the brightness temperature anomaly threshold includes: Will or The data in is sorted from small to large; Remove the first p1% of the data and the last p2% of the data; The minimum value in the retained data is the initial value of the reflectivity anomaly lower threshold or the initial value of the brightness temperature anomaly lower threshold, and the maximum value in the retained data is the initial value of the reflectivity anomaly upper threshold or the initial value of the brightness temperature anomaly upper threshold.

4. The polynya channel inversion method based on multi-channel satellite data according to claim 2, characterized in that: Step one also includes the step of performing zenith angle correction on the visible light radiation brightness value and the thermal infrared radiation brightness value.

5. The polynya channel inversion method based on multi-channel satellite data according to claim 1, characterized in that: In step 2, the reliability of each pixel in the satellite data is calculated based on the Normalized Snow Cover Index (NDSI), the Cloud Cover Product (CLW), and the Total Cloud Cover Product (TCP), including: ； Pixel credibility matrix Make corrections: when When , the credibility of the corresponding pixel is 0; when When , the corresponding pixel’s credibility is 0; is the preset threshold.

6. The polynya channel inversion method based on multi-channel satellite data according to claim 1, characterized in that: Step 3 also includes data scaling of BTA and RA so that their values ​​are within the range of [-1, 1], including: Using sliding window Traverse the brightness temperature anomaly image or reflectance outlier image ; Sliding window Scaling of brightness temperature anomalies or reflectivity anomalies within the image: ； ； in, For sliding window Internal brightness temperature anomaly, For sliding window The minimum brightness temperature anomaly value, For sliding window Maximum brightness temperature anomaly within the For sliding window Internal reflectivity anomaly, For sliding window The inner minimum reflectivity outlier, For sliding window the maximum reflectivity outlier within; After traversing the entire brightness temperature anomaly image, all Update brightness temperature anomalies , all Update reflectivity outliers ; Brightness temperature anomaly membership The calculation method is: ； Reflectivity outlier membership The calculation method is: ； in, is the brightness temperature width parameter of the jth thermal infrared channel, is the reflectivity width parameter of the i-th visible light channel, and n is a constant.

7. The polynya channel inversion method based on multi-channel satellite data according to claim 6, characterized in that: In step 1, reflectivity outlier image The methods for obtaining include: Using sliding window Traverse the reflectance data of the i-th visible light channel; ; For sliding window Reflectivity data within for The median image of .

8. The polynya channel inversion method based on multi-channel satellite data according to claim 7, characterized in that: In step 3, the calculation method for the brightness temperature inter-class mean is: Pixels are divided into the following three categories based on the upper and lower thresholds of brightness temperature anomaly: ； ； ； represents the brightness temperature value of pixel x of the jth thermal infrared channel at the kth iteration, represents the lower limit threshold of brightness temperature anomaly of the j-th thermal infrared channel at the k-th iteration, represents the upper threshold of brightness temperature anomaly of the j-th thermal infrared channel at the k-th iteration; 、 as well as The inter-class means are: ； ； ； Update the upper and lower thresholds for brightness temperature anomalies: ； ； in, is the dynamic adjustment coefficient; renew : ; in, , is a preset constant.

9. The polynya channel inversion method based on multi-channel satellite data according to claim 8, characterized in that: In step three, The comprehensive probability that a pixel belongs to a polynya The calculation method is: ； in, is the number of thermal infrared channels, is the number of visible light channels.

10. The polynya channel inversion method based on multi-channel satellite data according to any one of claims 1 to 9, characterized in that: After step 3, the scale of the sliding window is changed, and steps 1 to 3 are executed again to calculate the comprehensive probability under the current scale sliding window. , where q represents the scale number; Will The probability value of is assigned to g, and the other probability values ​​are assigned to 0; Mark the connected components of all pixels with a comprehensive probability value of g and find multiple connected areas; Calculate the width of each connected region ; Calculate the mean width of the polynya detected under the sliding window of each scale respectively and standard deviation ; according to as well as Calculate the weight of sliding windows of each scale : ； Calculate the final probability P that a pixel belongs to a polynya: 。

Images