Inter-ice water channel inversion method based on multi-channel satellite data

Abstract

The invention discloses an inter-ice water channel inversion method based on multi-channel satellite data, and belongs to the technical field of satellite image processing, and the method comprises the steps: obtaining a reflectivity abnormal value of a visible light radiation brightness value and an initial value of a reflectivity abnormal threshold value, and obtaining a brightness temperature abnormal value BTA of a thermal infrared radiation brightness value and an initial value of a brightness temperature abnormal threshold value; step 2, acquiring a pixel credibility matrix; and step 3, adjusting a brightness temperature abnormal threshold value and a reflectivity abnormal threshold value by adopting an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification, calculating a brightness temperature abnormal value membership degree and a reflectivity abnormal value membership degree, and judging the comprehensive probability that the pixel belongs to the inter-ice water channel according to each membership degree and the pixel credibility matrix. According to the method, the dynamic complexity of the polar region environment can be better adapted by utilizing the advantages of the visible light data in the aspect of identifying the boundary and texture of the ground features and the advantages of the thermal infrared data in the aspect of distinguishing the ground features at different temperatures, and the inversion precision of the inter-ice water channel is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of image processing, and more specifically, relates to a method for retrieving leads based on multi-channel satellite data. Background Art

[0002] Polar (Antarctic and Arctic) leads are key elements in the polar sea ice cover and play extremely important roles in aspects such as the polar climate system, ecological environment, and human activities. With the continuous change of polar sea ice, the importance of its research has become increasingly prominent. A lead is a linear or narrow crack area in the sea ice cover, and its formation is mainly due to external forces such as wind and ocean currents causing uneven stress on the sea ice, resulting in its rupture. These external forces cause the sea ice to displace and deform, and when the stress exceeds the strength of the sea ice itself, leads are formed. In winter, due to the large temperature difference between the ocean and the atmosphere in the newly formed leads, the heat of the seawater is quickly dissipated, and thin ice will quickly form, but there are still obvious differences from the surrounding thick ice. This difference is not only reflected in the physical structure, such as the thickness and density of the thin ice being different from those of the thick ice, but also in the thermodynamic properties. The temperature conduction and heat capacity of the thin ice are also different from those of the thick ice, thus having an important impact on the energy exchange in the polar region.

[0003] Since the late 1970s, remote sensing images obtained by satellite sensors have been used to detect polar leads. Willmes (2015) et al. detected leads based on MODIS thermal infrared images using the method of subtracting the background threshold. Hoffman et al. (2019) used the thermal infrared band of MODIS to design a multi-step method to effectively detect and characterize leads. Qu (2021) et al. proposed an improved algorithm for retrieving Beaufort leads using Terra / MODIS thermal infrared temperature anomaly images. Hoffman et al. (2021) combined the thermal infrared data of MODIS and VIIRS satellites and used artificial intelligence (AI) technology to explore the potential and challenges of AI technology in sea ice research and proposed solutions. Qiu (2023) et al. detected Arctic leads using high-resolution infrared images provided by the Thermal Infrared Spectrometer (TIS) on the Sustainable Development Science Satellite 1 (SDGSAT-1). The detection method based on a single near-infrared channel has certain limitations. The emissivity of newly formed thin ice leads (thickness < 5 cm) in the near-infrared band has a low difference from that of the leads not covered by thin pancakes, and the feature separability is poor. In addition, the anisotropic reflection caused by snow cover will mask the edge features of the leads, resulting in an increased missed detection rate for leads with a width < 100 m. These problems make it difficult for a single near-infrared channel to meet the detection requirements of the full scene.

[0004] Key et al. (1993) used the thermal infrared and visible images of the Advanced Very High Resolution Radiometer (AVHRR) to analyze the occurrence of leads, attempting to obtain lead width information and provide basic data and method references for subsequent research. They performed single-parameter segmentation by setting a fixed threshold method, which was limited by the dynamic complexity of the polar environment. Especially in the ice-water transition zone, there was a problem of high misjudgment rate and it could not adapt to the dynamically changing polar environment.

[0005] In addition, the existing technology only simulated the atmospheric effect through the radiative transfer model (LOWTRAN), there was a problem of cloud interference, and it was even more impossible to suppress cloud interference in real time. Summary of the Invention

[0006] In order to solve the technical problem that when inverting leads in the existing polar region, single-parameter segmentation is performed by setting a fixed threshold method, which is limited by the dynamic complexity of the polar environment and has a high misjudgment rate, the present invention proposes a method for inverting leads based on multi-channel satellite data, which can solve the above problems.

[0007] In order to solve the above technical problems, the present invention is implemented by adopting the following technical solutions: A method for inverting leads based on multi-channel satellite data, characterized by comprising: Step 1, a data preprocessing step, including: obtaining the reflectance anomaly value of the visible light radiance value , and obtaining the brightness temperature anomaly value of the thermal infrared radiance value , obtaining the initial value of the reflectance anomaly threshold and the initial value of the brightness temperature anomaly threshold of the leads in the satellite data; Step 2, a step of obtaining the pixel credibility, calculating the credibility of each pixel in the satellite data, and forming a pixel credibility matrix ; Step 3, adopting an iterative double-threshold adaptive optimization algorithm based on fuzzy classification to adjust the brightness temperature anomaly threshold and the reflectance anomaly threshold. The threshold includes an upper threshold and a lower threshold, calculating the membership degree of the brightness temperature anomaly value and the membership degree of the reflectance anomaly value. The membership degree of the brightness temperature anomaly value is the membership degree of the pixel corresponding to the brightness temperature anomaly value to the leads, and the membership degree of the reflectance anomaly value is the membership degree of the pixel corresponding to the reflectance anomaly value to the leads. According to each membership degree and the pixel credibility matrix, judge the comprehensive probability that the pixel belongs to the leads. The iterative double-threshold adaptive optimization algorithm based on fuzzy classification includes: Calculating the brightness temperature anomaly center of the current target class and the reflectance anomaly center , the target class is lead, the center of the brightness temperature anomaly is the central value of the brightness temperature anomaly threshold, and the center of the reflectance anomaly is the central value of the reflectance anomaly threshold. According to BTA and Calculate the membership degree of the brightness temperature anomaly value , according to RA and Calculate the membership degree of the reflectance anomaly value , 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 upper threshold and lower threshold of the brightness temperature anomaly, calculate the between-class mean of each class, which is the between-class mean of the brightness temperature. Update the upper threshold and lower threshold of the brightness temperature anomaly according to the between-class mean of the brightness temperature; Classify the pixels according to the upper threshold and lower threshold of the reflectance anomaly, calculate the between-class mean of each class, which is the between-class mean of the reflectance. Update the upper threshold and lower threshold of the reflectance anomaly according to the between-class mean of the reflectance; When the change amounts of the brightness temperature threshold and the reflectance threshold in two adjacent iterations simultaneously satisfy the following conditions, terminate the iteration: ; ; Calculate the comprehensive probability that the pixel belongs to the lead according to the membership degree obtained in the last iteration and the pixel credibility matrix .

[0008] Compared with the prior art, the advantages and positive effects of the present invention are as follows: The lead inversion method based on multi-channel satellite data of the present invention, since the thermal infrared channel of satellite data can better reflect the thermal radiation characteristics of ground objects, and the visible light channel of satellite data can better reflect the detailed information of leads. It utilizes the characteristics of high brightness temperature and low reflectance existing in leads in the polar region, searches for the brightness temperature anomaly data of the thermal infrared channel and the reflectance anomaly value of the visible light channel as the processing objects for lead inversion, combines the two, and utilizes the advantages of visible light data in identifying ground object boundaries and textures, as well as the advantages of thermal infrared data in distinguishing ground objects with different temperatures, to achieve complementary advantages.

[0009] When determining the brightness temperature anomaly data and the reflectance anomaly data, classify the pixels according to the current brightness temperature anomaly threshold and reflectance anomaly threshold, calculate the between-class mean of each class. The between-class mean can reflect the overall distribution of each class of anomaly values, and recalculate the reflectance anomaly threshold and brightness temperature anomaly threshold according to the between-class mean, so that the present solution can better adapt to the dynamic complexity of the polar environment, reduce the misjudgment rate, and improve the lead inversion accuracy.

[0010] After reading the detailed description of the embodiments of the present invention in conjunction with the drawings, other features and advantages of the present invention will become clearer. Brief Description of the Drawings

[0011] Figure 1 is a flowchart of an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 2a is an image of 10.8μm visible light channel data in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 2b is 12.0 in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention visible light channel data image; Figure 2c is an image of visible light channel data in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 3a is an image of retrieved leads using a small-scale sliding window in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 3b is an image of retrieved leads using a medium-scale sliding window in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 3c is an image of retrieved leads using a large-scale sliding window in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 4 is an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention, where Figures 3a - 3c the fused image of retrieved leads; Figure 5 is an image of sea ice concentration in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention; Figure 6 is an image of retrieved leads in the Arctic region in an embodiment of the method for retrieving leads based on multi-channel satellite data proposed by the present invention. Detailed Description of the Invention

[0012] The following further describes in detail the specific embodiments of the present invention with reference to the accompanying drawings.

[0013] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Apparently, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0014] It should be noted that in the description of the present invention, the terms indicating directions or positional relationships such as "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc. are based on the directions or positional relationships shown in the accompanying drawings. This is only for 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 should not be construed as a limitation to the present invention. In addition, the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance. In the description of the present invention, "a plurality of" means two or more unless otherwise specifically defined.

[0015] In the present invention, unless otherwise clearly specified and limited, the terms "installed", "connected", "connected", "fixed", etc. should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or integrated; it can be a mechanical connection or an electrical connection; it can be directly connected or indirectly connected through an intermediate medium, and it can be the communication inside two elements or the interaction relationship between two elements. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific situations.

[0016] Example 1, see Figure 1 As shown, this embodiment proposes an ice floe lead inversion method based on multi-channel satellite data, including: Step 1, a data preprocessing step, including: obtaining the reflectance anomaly value RA of the visible light radiance value and obtaining the brightness temperature anomaly value BTA of the thermal infrared radiance value, and obtaining the initial value of the reflectance anomaly threshold and the initial value of the brightness temperature anomaly threshold of the ice floe lead in the satellite data.

[0017] In this embodiment, multiple thermal infrared channels and visible light channels in the MERSI-II sensor of the FY-3D satellite are used for ice floe lead inversion.

[0018] Step 2, a step of obtaining the pixel credibility, calculating the credibility of each pixel in the satellite data, and forming a pixel credibility matrix .

[0019] Step 3: Use the iterative dual-threshold adaptive optimization algorithm based on fuzzy classification to adjust the brightness temperature anomaly threshold and the reflectivity anomaly threshold. The thresholds include the upper threshold and the lower threshold. Calculate the membership degree of the brightness temperature anomaly value and the membership degree of the reflectivity anomaly value. The membership degree of the brightness temperature anomaly value is the membership degree of the pixel corresponding to the brightness temperature anomaly value to the polynya. The membership degree of the reflectivity anomaly value is the membership degree of the pixel corresponding to the reflectivity anomaly value to the polynya. Judge the comprehensive probability that the pixel belongs to the polynya according to each membership degree 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 the reflectivity anomaly center , the target class is the polynya, the brightness temperature anomaly center is the central value of the brightness temperature anomaly threshold, and the reflectivity anomaly center is the central value of the reflectivity anomaly threshold. According to BTA and calculate the membership degree of the brightness temperature anomaly value , according to RA and calculate the membership degree of the reflectivity anomaly value , k represents the number of iterations, j represents the thermal infrared channel number, and i represents the visible light channel number.

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

[0021] When the change amounts of the brightness temperature threshold and the reflectivity threshold in two adjacent iterations simultaneously meet the following conditions, terminate the iteration: ; .

[0022] According to the membership degree and the pixel credibility matrix obtained in the last iteration, calculate the comprehensive probability that the pixel belongs to the polynya . The comprehensive probability reflects the probability that the pixel belongs to the polynya, and realizes the inversion of the polynya.

[0023] In the polynya inversion method based on multi-channel satellite data in this embodiment, since the thermal infrared channel of the satellite data can better reflect the thermal radiation characteristics of the ground object, as Figure 2a , Figure 2b shown, the visible light channel of the satellite data can better reflect the detailed information of the polynya, as Figure 2cAs shown in the figure, the characteristics of high brightness temperature and low reflectivity in the polynyas in the polar region are utilized. By searching for the brightness temperature anomaly data in the thermal infrared channel and the reflectivity anomaly values in the visible light channel as the processing objects for inverting polynyas, and combining the two, the advantages of visible light data in identifying object boundaries and textures and the advantages of thermal infrared data in distinguishing objects with different temperatures are utilized to achieve complementary advantages.

[0024] When determining the brightness temperature anomaly data and the reflectivity anomaly data, the pixels are classified according to the current brightness temperature anomaly threshold and the reflectivity anomaly threshold, and the between-class means of each class are calculated. The between-class means can reflect the overall distribution of each class of anomaly values, and the reflectivity anomaly threshold and the brightness temperature anomaly threshold are recalculated based on the between-class means, so that this solution can better adapt to the dynamic complexity of the polar environment, reduce the misjudgment rate, and improve the inversion accuracy of polynyas.

[0025] In some embodiments, in step one, the methods for obtaining RA and BTA include: obtaining the visible light channel data and the thermal infrared channel data of satellite data, performing radiometric calibration on the visible light channel data to obtain the visible light radiance value, and obtaining the reflectivity anomaly value RA of the visible light radiance value; performing radiometric calibration on the thermal infrared channel data to obtain the thermal infrared radiance value, converting the thermal infrared radiance value into the brightness temperature of an equivalent blackbody, and calculating the brightness temperature anomaly value BTA; obtaining the initial value of the reflectivity anomaly threshold according to RA, and obtaining the initial value of the brightness temperature anomaly threshold according to BTA.

[0026] The data of the visible light channel is processed. Using the radiometric calibration technology, based on the radiation response characteristics and calibration parameters of the sensor, the digital signal DN measured by the sensor is converted into a radiance value Ref with physical significance. Then, a mathematical model of the solar zenith angle and the radiance value Ref is established, and the change value of the radiance value Ref caused by the change of the solar zenith angle is corrected through this model, that is, the error caused by the change of the solar zenith angle is eliminated to ensure the consistency of the data. Then, the reflectivity anomaly value RA is calculated based on the FY-3D / MERSI-II radiance value Ref, and then the membership degree of the reflectivity data anomaly value is calculated.

[0027] 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 object, so as to accurately obtain the radiation energy information of the object.

[0028] For the reflectivity data of the three visible light channels, a specific formula is used for calibration: ; .

[0029] In the formula, Normalize for the i-th visible light channel The value is the brightness value of each pixel in the remote sensing image, recording the intensity of the electromagnetic wave reflected or radiated by the ground object, and is expressed as an integer value without units.

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

[0031] Since the solar altitude angle determines the angle between the sun's rays and the ground, the smaller the solar altitude angle, the longer the path of the light through the atmosphere, the stronger the scattering and absorption of the light by the atmosphere, the less solar radiation reaches the ground, and the less radiation reflected by the ground object back to the sensor. Conversely, the larger the solar altitude angle, the more radiation reflected by the ground object back to the sensor.

[0032] Process the data of the thermal infrared channels, and calculate the radiance of the data of the two thermal infrared channels respectively , and then convert to the equivalent blackbody brightness temperature according to the Planck inverse transformation formula and the brightness temperature correction coefficient, calculate the brightness temperature anomaly value , and then calculate the membership degree of the brightness temperature anomaly value.

[0033] Specifically, the radiance value measured by the FY-3D / MERSI-II sensor for the j-th thermal infrared channel, as well as the correction parameters and given by the sensor are used for calibration calculation to obtain the radiance (mW / (m2 cm sr)) of the j-th thermal infrared channel: .

[0034] Among them, , .

[0035] Then, using the Planck inverse transformation formula, convert the radiance to the equivalent blackbody brightness temperature , and its specific form is: .

[0036] Then, with the help of the channel brightness temperature correction coefficients and convert Convert to channel blackbody brightness temperature and , the formula is as follows: .

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

[0038] Then convert the brightness temperatures of the two thermal infrared channels into a brightness temperature anomaly value image , and its specific form is: .

[0039] Among them, is the sliding window, is the brightness temperature anomaly value image, is the brightness temperature data of the corrected thermal infrared channel, is the median image of is the pixel position.

[0040] In some embodiments, in step one, the method for obtaining the initial value of the reflectivity anomaly threshold and the initial value of the brightness temperature anomaly threshold includes: Sort the data in RA or BTA from smallest to largest; Remove the data in the first p1% and the data in the last p2%.

[0041] The minimum value of the remaining 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 of the remaining data is the initial value of the reflectivity anomaly upper threshold or the initial value of the brightness temperature anomaly upper threshold.

[0042] In some embodiments, p1% and p2% can be set to 5% or other numbers. In this way, the remaining data is the data sorted from 5% to 95%.

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

[0044] In the same remote sensing image, if the solar elevation angles in different regions are different, it will cause differences in the gray values of different regions in the image, and correction is required.

[0045] First, convert the solar zenith angle to radian system: .

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

[0047] in is the solar zenith angle, is the zenith angle converted to radians, is the reflectivity of the ith visible light channel after solar zenith angle correction. After correction, the image reflectivity error caused by the solar zenith angle can be effectively eliminated, ensuring that the reflectivity data can accurately reflect the optical characteristics of the polynya and its surrounding environment.

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

[0049] The vast snow cover in the Arctic provides an excellent condition for NDSI to play its advantage in distinguishing snow from clouds. Since the acquisition of FY3D / MERSI-II data in the Arctic is relatively stable, the NDSI data quality is high. It can provide accurate ground object background information for cloud detection, especially when distinguishing clouds above the snow surface, it can clearly outline the boundaries of clouds and accurately identify the existence of clouds, providing a reliable basis for subsequent cloud detection analysis.

[0050] In order to eliminate the interference of clouds on the detection of polynyas as much as possible, this method combines the normalized snow cover index NDSI, cloud detection product CLW, high cloud cover product and total cloud cover product TCP to obtain the credibility of pixel inversion. NDSI data can be used to distinguish snow from other objects and assist in judging the existence of clouds; cloud detection products can directly identify cloud-covered areas; high cloud cover products and total cloud cover products provide quantitative information on cloud cover. By comprehensively analyzing these data, the credibility of pixels is evaluated to reduce the impact of cloud-covered areas on the detection of polynyas. Snow cover index data can be used to distinguish snow from other objects and assist in judging the existence of clouds; cloud detection products can directly identify cloud-covered areas, while total cloud cover products provide quantitative information on cloud cover.

[0051] Based on the membership of polynya pixel and the credibility of pixel inversion, the probability of each pixel belonging to the polynya category is calculated to achieve a refined inversion of polynyas and improve the accuracy and reliability of detection. In some embodiments, the credibility of each pixel in the satellite data is calculated according to the normalized snow cover index data NDSI, the cloud detection product CLW and the total cloud cover product TCP in step 2, including: ; Pixel Credibility Matrix To make corrections: when When it is [a certain condition], the confidence of the corresponding pixel is 0; When it is [a certain condition], the confidence of the corresponding pixel is 0; Among them, is a preset threshold.

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

[0053] Among them is the data of the 3rd channel 0.65 of the MERSI-II sensor data, is the data of the 6th channel 1.64 of the MERSI-II sensor data.

[0054] Due to the missed detection phenomenon in the total cloud cover product, NDSI and CLM are added for supplementary correction.

[0055] In some embodiments, step three further includes performing data scaling processing on BTA and RA so that their value ranges are between [-1, 1], including: Using a sliding window to traverse the brightness temperature anomaly value image or the reflectivity anomaly value image .

[0056] Performing scaling processing on the brightness temperature anomaly value or the reflectivity anomaly value within the sliding window : .

[0057] .

[0058] Among them, is the brightness temperature anomaly value within the sliding window , is the minimum brightness temperature anomaly value within the sliding window , is the maximum brightness temperature anomaly value within the sliding window ; is the reflectivity anomaly value within the sliding window , is the minimum reflectivity anomaly value within the sliding window , is the maximum reflectivity anomaly value within the sliding window .

[0059] After traversing the entire brightness temperature anomaly value image, all update the brightness temperature anomaly value , all Update reflectivity anomaly value 。

[0060] In some embodiments, the center of the brightness temperature anomaly of the current target class and the center of the reflectivity anomaly are calculated as follows: 。

[0061] 。

[0062] In some embodiments, the membership degree of the brightness temperature anomaly value is calculated as follows: ; The membership degree of the reflectivity anomaly value is calculated as follows: ; wherein, is the brightness temperature width parameter of the j-th thermal infrared channel, is the reflectivity width parameter of the i-th visible light channel, and n is a constant.

[0063] In some embodiments, in step one, the method for obtaining the reflectivity anomaly value image includes: Using a sliding window to traverse the reflectivity data of the i-th visible light channel; ; is the reflectivity data within the sliding window , is the median image of.

[0064] In some embodiments, in step three, the method for calculating the between-class mean of the brightness temperature is: Dividing the pixels into the following three categories according to the upper threshold of the brightness temperature anomaly and the lower threshold of the brightness temperature anomaly: ; ; 。

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

[0066] , and The between-class means of are respectively: ; ; .

[0067] Update the upper threshold and lower threshold of the brightness temperature anomaly: ; ; wherein, is a dynamic adjustment coefficient used to compensate for the influence of ice-water mixed pixels. In this embodiment, .

[0068] Update : .

[0069] wherein, , is an empirical value, is a preset constant.

[0070] Similarly, the calculation method of the between-class mean of reflectivity is: Divide the pixels into the following three categories according to the upper threshold and lower threshold of the reflectivity anomaly: ; ; .

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

[0072] , and The between-class means of are respectively: ; ; .

[0073] Update the upper threshold and lower threshold of the brightness temperature anomaly: ;

[0074] 。

[0075] Update : 。

[0076] Among them, is the experience value, is the preset constant.

[0077] When the change amount of the brightness temperature threshold and the change amount of the reflectivity threshold between two adjacent iterations simultaneously satisfy the following conditions, the iteration is terminated: ; 。

[0078] In some embodiments, 。

[0079] In some embodiments, in step three, the comprehensive probability that a pixel belongs to an open water lane is calculated by combining the multi-source membership degree and the credibility weight:

[0080] The comprehensive probability that a pixel belongs to an open water lane is calculated as follows: 。

[0081] Among them, is the number of thermal infrared channels, is the number of thermal infrared channels.

[0082] In some embodiments, after step three, it further includes changing the scale of the sliding window, and then executing steps one to three again to calculate the comprehensive probability under the sliding window of the current scale, and obtaining , as Figures 3a - 3c shown, where q represents the scale number.

[0083] Assign the probability value of to g, and assign other probability values to 0.

[0084] Perform connected component labeling on all pixels with a comprehensive probability value of g, and find multiple connected regions.

[0085] Calculate the width of each connected region 。

[0086] Calculate the mean value and the standard deviation of the widths of the open water lanes detected under the sliding windows of each scale respectively.

[0087] According to and Calculate the weights of the sliding windows at each scale :

[0088] 。

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

[0090] Then consider the area of as a polynya, and consider the area of as a potential polynya, as shown in Figure 4 。

[0091] As Figure 5 shown, it is the sea ice concentration image on April 2, 2021. As shown in Figure 6 ,it is the polynya image in the Arctic region inverted by using the method of this embodiment for Figure 5 。

[0092] In some embodiments, calculate the width of the polynya , specifically as follows: First, use the 3×3 neighborhood structure S to perform connected component labeling on the polynya probability value , and assign a label to each polynya, then obtain labeled objects, where 。

[0093] For each connected object , calculate its number of pixels, and the formula is as follows: 。

[0094] That is, represents the total number of all points that satisfy the label equal to 。

[0095] Among them, represents the identification of the connected region to which the pixel with coordinates in belongs, is the total number of connected regions in , and represents the number of pixels in the

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

[0097] Then split the continuous segments by row. For each row y, extract the column coordinate sequence The calculation of the difference between adjacent coordinates can then be obtained: .

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

[0099] Then the average row width can be obtained: .

[0100] In the formula, represents the number of elements in the set of continuous segments , represents the length of the continuous segment .

[0101] Similarly, the average column height can be obtained: .

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

[0103] From a physical mechanism perspective, different windows correspond to different spatial frequency responses. Small windows retain high-frequency details and can reflect the fine structure of leads; large windows capture low-frequency trends and show the overall characteristics and distribution trends of leads.

[0104] The multi-scale feature weighted fusion is based on the Gaussian distribution principle, and the features at different scales are weighted by the Gaussian weight function. In the inversion of leads, windows at different scales have their own advantages and disadvantages in detecting leads. The Gaussian weight function can reasonably allocate the weights of features at different scales in the fusion result according to the difference between the theoretical detection width corresponding to the window and the median value of the target lead width.

[0105] For the sliding window at the q scale , the following can be calculated and the width of the lead . Then when the number of scales of the sliding window is , the mean value and the standard deviation of the lead widths under z sliding windows can be calculated .

[0106] Certainly, 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 those of ordinary skill in the art within the scope of the essence of the present invention shall also fall within the protection scope of the present invention.

Claims

1. An ice lead inversion method based on multi-channel satellite data, characterized in that, Including: Step 1, data preprocessing step, including: obtaining the reflectance anomaly value of the visible light radiance value , and obtaining the brightness temperature anomaly value of the thermal infrared radiance value , obtaining the initial values of the reflectance anomaly threshold and the brightness temperature anomaly threshold of the ice channel in the satellite data; Step 2, the step of obtaining the pixel credibility, calculates the credibility of each pixel in the satellite data and forms a pixel credibility matrix ; Step 3: Adopt an iterative dual-threshold adaptive optimization algorithm based on fuzzy classification to adjust the brightness temperature anomaly threshold and the reflectivity anomaly threshold. The threshold includes an upper threshold and a lower threshold, calculate the membership degree of the brightness temperature anomaly value and the membership degree of the reflectivity anomaly value. The membership degree of the brightness temperature anomaly value is the membership degree of the pixel corresponding to the brightness temperature anomaly value to the polynya. The membership degree of the reflectivity anomaly value is the membership degree of the pixel corresponding to the reflectivity anomaly value to the polynya. Judge the comprehensive probability that the pixel belongs to the polynya according to each membership degree and the pixel credibility matrix. The iterative dual-threshold adaptive optimization algorithm based on fuzzy classification includes: Calculate the center of the brightness temperature anomaly of the current target class and the center of the reflectivity anomaly , where the target class is an ice channel, the center of the brightness temperature anomaly is the central value of the brightness temperature anomaly threshold, and the center of the reflectivity anomaly is the central value of the reflectivity anomaly threshold. According to BTA and calculate the membership degree of the brightness temperature anomaly value , according to RA and calculate the membership degree of the reflectivity anomaly value , 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 upper limit threshold of the brightness temperature anomaly and the lower limit threshold of the brightness temperature anomaly, calculate the between-class mean of each class, which is the between-class mean of the brightness temperature, and update the upper limit threshold of the brightness temperature anomaly and the lower limit threshold of the brightness temperature anomaly according to the between-class mean of the brightness temperature; Classify the pixels according to the upper limit threshold of the reflectivity anomaly and the lower limit threshold of the reflectivity anomaly, calculate the between-class mean of each class, which is the between-class mean of the reflectivity, and update the upper limit threshold of the reflectivity anomaly and the lower limit threshold of the reflectivity anomaly according to the between-class mean of the reflectivity; When the change amounts of the brightness temperature threshold and the reflectivity threshold in two adjacent iterations simultaneously meet the following conditions, terminate the iteration: ； ； Calculate the comprehensive probability of a pixel belonging to an open lead based on the membership degree obtained from the last iteration and the pixel credibility matrix. .

2. The method for retrieving polynya based on multi-channel satellite data according to claim 1, characterized in that In Step 1, and the acquisition method includes: acquiring the visible light channel data and the thermal infrared channel data of satellite data, performing radiometric calibration on the visible light channel data to obtain the visible light radiance value, and acquiring the reflectance anomaly value of the visible light radiance value ; performing radiometric calibration on the thermal infrared channel data to obtain the thermal infrared radiance value, converting the thermal infrared radiance value into the brightness temperature of an equivalent blackbody, and calculating the brightness temperature anomaly value ; according to acquiring the initial value of the reflectance anomaly threshold, and according to acquiring the initial value of the brightness temperature anomaly threshold.

3. The method for inverting leads 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: Sort the data in or in ascending order; Remove the first p1% of the data and the last p2% of the data; The minimum value in the remaining data is the initial value of the lower limit threshold of the reflectivity anomaly or the initial value of the lower limit threshold of the brightness temperature anomaly, and the maximum value in the remaining data is the initial value of the upper limit threshold of the reflectivity anomaly or the initial value of the upper limit threshold of the brightness temperature anomaly.

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

5. The ice-strewn lane inversion method based on multi-channel satellite data according to claim 1, wherein In step 2, calculating the credibility of each pixel in the satellite data according to the normalized difference snow index data NDSI, the cloud detection product CLW, and the total cloud cover product TCP includes: ； Correct the pixel credibility matrix as follows: When the confidence level of the corresponding pixel is 0; When the confidence level of the corresponding pixel is 0; where is a preset threshold value.

6. The method for retrieving leads based on multi-channel satellite data according to claim 1, wherein Step 3 also includes performing data scaling processing on BTA and RA so that their value ranges are between [-1, 1], including: Adopt a sliding window Traverse the brightness temperature anomaly image Or the reflectivity anomaly image ; For the sliding window Perform scaling processing on the brightness temperature outliers or reflectivity outliers within it: ； ； Among them, is the brightness temperature anomaly value within the sliding window ; is the minimum brightness temperature anomaly value within the sliding window ; is the maximum brightness temperature anomaly value within the sliding window ; is the reflectivity anomaly value within the sliding window ; is the minimum reflectivity anomaly value within the sliding window ; is the maximum reflectivity anomaly value within the sliding window ; After traversing the entire brightness temperature anomaly image, all Update the brightness temperature anomaly , all Update the reflectivity anomaly ; Membership degree of brightness temperature anomaly value The calculation method is as follows: ； Membership degree of reflectivity anomaly value The calculation method is as follows: ； wherein, is the brightness temperature width of the j-th thermal infrared channel parameter, is the reflectivity width parameter of the i-th visible light channel, and n is a constant.

7. The method for retrieving leads based on multi-channel satellite data according to claim 6, wherein In Step 1, the method for obtaining the reflectivity anomaly value image includes: Adopt a sliding window Traverse the reflectance data of the i-th visible light channel; ; is the reflectance data within the sliding window, and is the median image of 8. The method for retrieving leads based on multi-channel satellite data according to claim 7, wherein In step 3, the calculation method of the between-class mean of the brightness temperature is: Classify the pixels into the following three categories according to the upper limit threshold of the brightness temperature anomaly and the lower limit threshold of the brightness temperature anomaly: ； ； ； represents the brightness temperature value of the pixel x in the j-th thermal infrared channel at the k-th iteration, represents the lower limit threshold of the brightness temperature anomaly in the j-th thermal infrared channel at the k-th iteration, represents the upper limit threshold of the brightness temperature anomaly in the j-th thermal infrared channel at the k-th iteration; , and The between-class means of are respectively: ； ； ； Update the upper limit threshold of the brightness temperature anomaly and the lower limit threshold of the brightness temperature anomaly: ； ； Among them, is a dynamic adjustment coefficient; Update : ; Among them, , is a preset constant.

9. The method for retrieving leads based on multi-channel satellite data according to claim 8, wherein In step 3, Comprehensive probability that a pixel belongs to a polynya is calculated as follows: ； Among them, is the number of thermal infrared channels, is the number of thermal infrared channels.

10. The method for retrieving leads based on multi-channel satellite data according to any one of claims 1-9, characterized in that After Step 3, it further includes changing the scale of the sliding window, and executing Steps 1 to 3 again to calculate the comprehensive probability under the sliding window of the current scale, obtaining , where q represents the scale number; Assign the probability value to g and assign other probability values to 0; Perform connected component labeling on all pixels with a comprehensive probability value of g, and find multiple connected regions; Calculate the width of each connected region ; Calculate the mean of the widths of the detected leads under the sliding windows at each scale and the standard deviation ; According to and calculate the weights of the sliding windows at each scale : ； Calculate the final probability P that the pixel belongs to the polynya: 。

Images