Denoising Sentinel-1 Extra-Wide mode cross-polarization images over sea ice

—Sentinel-1 (S1) extra-wide (EW) swath data in cross-polarization (horizontal-vertical, HV or vertical-horizontal, VH) are strongly affected by the scalloping effect and thermal noise, particularly over areas with weak backscattered signals, such as sea surfaces. Although noise vectors in both the azimuth and range directions are provided in the standard S1 EW data for subtraction, the residual thermal noise still significantly affects sea ice detection by the EW data. In this paper, we improve the denoising method developed in previous studies to remove the additive noise for the S1 EW data in cross-polarization. Furthermore, we propose a new method for eliminating the residual noise (i.e. multiplicative noise) at the sub-swath boundaries of the EW data, which cannot be well processed by simply subtracting the reconstructed 2-D noise field. The proposed method of removing both the additive and multiplicative noise was applied to EW HV-polarized images processed using different Instrument Processing Facility (IPF) versions. The results suggest that the proposed algorithm significantly improves the quality of EW HV-polarized images under various sea ice conditions and sea states in marginal ice zone (MIZ) of the Arctic. This is of great support for the utilization of cross-polarization SAR images in wide swaths for intensive sea ice monitoring in polar regions.


I. INTRODUCTION
PACEBORNE Synthetic Aperture Radar (SAR) has played important roles in Earth observations since the emergence of satellite remote sensing. As an active microwave sensor, SAR can provide high-resolution images independent of clouds and daylight. The C-band ERS-1/2, Envisat, RADARSAT-2, Sentinel-1(S1) A/B, Gaofen-3, X-band TerraSAR-X/TanDEM-X, Cosmo-SkyMed, L-band JERS-1 and ALOS PALSAR are representatives of operational spaceborne SAR, which have been continuously acquiring images over the Earth for nearly three decades.
Along with the development of SAR techniques, dual and quad-polarizations are increasingly applied in more domains than single-polarization [1], e.g., detection of targets in the sea [2], [3], [4] and retrieval of sea surface dynamics parameters [5], [6], [7]. Compared with the co-polarized (horizontal-horizontal, This  HH or vertical-vertical, VV) radar signal, the cross-polarized (horizontal-vertical, HV or vertical-horizontal, VH) signal is less sensitive to the incidence angle and wind speed [8], [9], which is particularly advantageous in sea ice monitoring. Several studies have shown that the combination of HH and HV polarizations can significantly improve the accuracy of sea ice classification and ice-water discrimination [10], [11]. However, open water and newly formed sea ice often have lower radar backscatter values in HV polarization than in HH polarization channels [11]. These cross-polarization values are closer to the noise floor and therefore often have a lower signal-to-noise ratio (SNR) introducing artifacts produced by thermal noise [12]. Thermal noise is an additive background noise that causes a noise floor, and it can be processed with the same processing gains applied to the true signal [13]. Thermal noise is especially obvious throughout an image in wide swath observation modes [13], [14].
To achieve a large observational area but not significantly reduce spatial resolution, Scanning SAR (ScanSAR) and Terrain Observation by Progressive Scans SAR (TOPSAR) technologies have been developed. The burst mode is used in both ScanSAR and TOPSAR techniques, which is necessary to provide wide swath coverage [15]. In the ScanSAR imaging mode, with the satellite moving ahead, the sensor periodically switches the antenna beam to several "sub-swathes" in the range direction to obtain wide swath coverage. The beam acquires a finite periodical sequence of echoes, i.e., bursts, for each sub-swath [15]. In each burst, the echo intensity at the edges of the antenna beam is much larger than those close to its center, which causes the "scalloping effect" (e.g., referring to Fig. 1). To reduce the scalloping effect and further improve image quality, the TOPSAR acquisition mode was proposed. With the TOPSAR technique, in addition to steering the beam in range as that in ScanSAR mode, the beam is also electronically steered from backward to forward in the azimuth for each burst (e.g., [15]). Although one of the purpose of applying the TOPSAR acquisition by the Sentinel-1 (S1) Interferometry wide (IW) swath and extra-wide (EW) swath modes is to reduce the scalloping, the fact is that the it seems to still exist, which is particularly distinct in the cross-polarization data over areas with low signal-to-noise ratio (SNR). , While the scalloping noise changes along the azimuth direction with a variable burst-wise intensity, another thermal noise, normally called the noise floor, varies as a function of the range position [17]. In the HV polarization image of ScanSAR and TOPSAR wide swath mode, the boundaries of the noise floor are highlighted as discontinuous sharp changes in intensity, which is called the banding effect [14], [18], [19]. Fig. 1(a) shows an S1 HV polarization EW swath image with a swath width of approximately 400 km containing five subswaths, where the scalloping effect along the azimuth direction and the banding effect (bright strips) among sub-swaths are clearly visible. The upper panel in Fig. 1(b) shows the variation in the range-averaged profile of the backscatter values along the azimuth direction in the first sub-swath (EW1, counting from the near-range direction). The peak-trough-peak variations in the radar backscatter signal represent the bursts. The bottom panel in Fig. 1(b) shows the variation in the azimuthal averaged radar backscatter along the range direction (of the whole image), showing upward and downward ("V-shape") features in each sub-swath as well as high-intensity jumps at the inter-swath boundaries, which is consistent with the noise floor variations. Notably, the intensity of the center of EW1 is also particularly large, and this is often observed in the EW image in crosspolarization, particularly when open water or thin ice presents in this sub-swath.
The existence of scalloping in the azimuth direction and the noise floor in the range direction significantly disturb the signal (often weak) in the SAR cross-polarized images acquired over the sea surface. To remove scalloping, the azimuth antenna pattern (AAP) for ScanSAR [14] and the azimuth antenna element pattern (AAEP) for TOPSAR [13] are the appropriate beam pattern corrections. In addition, a band-stop finite impulse response (FIR) filter [17], the deconvolution technique [20], the wavelet multiresolution analysis (MRA) [21] and the Kalman filter of image intensity [14] or image texture [22] are also used to reduce scalloping. In the S1 EW data, the scalloping noise vectors are provided in the product annotation files along with the range noise vectors since the Instrument Processing Facility (IPF) version 2.90.
For removal of the noise floor, the noise equivalent sigma zero (NESZ), i.e., the range noise vector, is provided in the SAR product annotation files as a look-up table to do the subtraction. [8], [9], [11], [12], [16], [17], [23]. However, the image after subtracting the noise field that was interpolated from NESZ vectors still has some problems in its application for denoising. On the one hand, the low-level signal may appear too noisy to be useful after the noise-subtraction [18]. Thus, some studies give a certain threshold for the low signals [11], [24] or just discard the negative values from analysis [25]. On the other hand, the residual noise present in the inter-swath boundaries is still significant in the denoised images [16]. For sea ice monitoring by SAR cross-polarized data, the residual noise in open water can be equivalent to or close to the radar backscatter of thin ice, which leads to misinterpretations between sea ice and open water [18]. In the texture features, such as the graylevel co-occurrence matrix (GLCM) features, the residual noise is even amplified [18], [26]. The reasons can be attributed to two aspects. One is that the NESZ cannot accurately describe the noise floor causing the residual additive noise. The other is that the intensity values in the inter-swath boundaries are of higher scattering entropy [25], caused by the multiplicative noise, which cannot be removed solely by the NESZ subtraction [29].
Using the given noise information (NESZ vectors) to modify and subtract the noise field is a preferable approach to remove the additive noise of the intensity image [13]. There are limited studies about the further correction of NESZ vectors. A technique for SAR denoising based on the NESZ is to compute the concentration of "bad" pixels in a sliding window and then multiply this matrix to the mean value of the noise floor to build the corrected noise field [27]. Another noise scaling approach based on the local SNR estimation is also computed in a sliding window [28]. Recently, Park et al. [13] proposed an effective denoising method for the S1 EW data in cross-polarization, which can find optimal parameters to modify the provided noise vector over homogeneous areas, e.g., open water and dense pack ice, to subtract the additive noise. However, the multiplicative noise in the inter-swath boundaries of EW data still exists after the denoising. Thus, the researchers proposed another method called the "texture noise correction" method in [29], which is built based on the noise equivalent standard deviation (NESD) and signal-plus-noise-to-noise ratio (SNNR) to rescale the local data using a sliding window, thus introducing the coefficients for the NESD model. In the complete noise removal method of [13] and [29], the coefficients are determined by averaging a large number of HVpolarized images over a homogeneous area with a low SNNR.
Our motivation is to improve the quality of S1 EW images in cross-polarization (HV) acquired in polar regions for better sea ice monitoring. . However, as described above, the existing scalloping in the azimuth and varied noise floor in the range direction significantly reduce the quality of EW images in HV polarization and lead to misinterpretation of sea ice or open water, as well as different ice types. Therefore, denoising the EW image in HV polarization is a key step in the entire processing chain of sea ice detection in the Arctic by S1 EW data.
We improve the denoising algorithm proposed in [13] by changing the methods of calculating the scaling and balancing parameters, which are used to modify the ESA-provided NESZ noise vector. Furthermore, we developed a method to remove the multiplicative noise in additive noise-subtracted EW images.
Following the Introduction, the general calibration of noise removal of the S1 EW data in HV polarization is described in Section II. In Section III, the proposed denoising algorithm, including removing additive and multiplicative noise, is described in details. Applications of the algorithm to denoise EW data acquired in the marginal ice zone (MIZ) of the Arctic and processed by different IPF versions are presented in Section IV. We also compared the denoising results using the proposed method with the results using the methods introduced in [13] and [29]. The summary and conclusions are given in the last section.

A. Sentinel-1 EW Data
The S1 SAR constellation consists of two satellites, S1A and 1B. The two C-band (5.4 GHz) SARs have four operational imaging modes of acquiring data: Interferometric-wide (IW), EW, strip map (SM) and wave mode (WV). The EW mode data in HH and HV polarizations are frequently acquired in the Arctic for sea ice monitoring.
The EW mode image covers incidence angles from approximately 18.9° to 47.0° with a swath width of 400 km. The data have a pixel spacing of 40 m × 40 m [8]. The repeat cycle of the S1A and S1B satellites is 6 days, and this cycle can be shorter up to 2 days in the Arctic region. The high spatial and temporal resolutions are particularly conducive to polar observations [23]. In addition, the NESZ of the EW mode is expected to be lower than that of the IW mode data [1], [7], which also highly favors sea ice monitoring. In this paper, the S1 EW data used are the ground range detection (GRD) products, which are acquired via the Copernicus Open Access Hub (scihub.copernicus.eu).
ESA has been improving the estimation of the NESZ noise vectors for removal of the additive noise based on continuously updated IPF versions. However, even with the newest IPF version 3.1 (released in June 2019), slight under-or overcompensation and unbalancing still exist in noise-subtracted images using the provided noise vectors, as shown in Section IV. Moreover, the multiplicative noise that can be largely amplified in texture images has not been well removed in all IPF versions.

B. Standard Radiometric Calibration and Additive Noise Subtraction of S1 EW Data
The standard radiometric calibration and noise removal of the S1 EW data are completed using the following equation [30]: is the image digital number at pixel . is the calibration coefficient, and is the ESA-provided noise vector in the range direction. The right term of the equation indicates that the calibrated sigma naught of the signal is equal to the raw sigma naught value , of pixel minus the NESZ values , . Prior to S1 data processing version IPF 2.90 (available from March 2018), only the thermal noise vectors in the range direction are provided by the ESA, while the descalloping gain needs to be computed using the AAEP look-up where , and , are the scalloping noise in the azimuth direction and thermal noise in the range direction, respectively. Here, , is the ESA-provided range noise vector. The detailed computation of the descalloping gain was given in [13], which is also provided in the Appendix of this paper for reference. We specified several intermediate variables and processes that were not clear to us when we used the method in [13] to calculate the descalloping gain, which might be helpful to others to process the historical S1 EW data. In the new IPFs after version 2.90, the descalloping vectors have also been provided for the radiometric calibration of the S1 EW data. [13] The main concept of the denoising algorithm in [13] proposed for the S1 EW data in cross-polarization is to find an optimal noise scaling factor and an inter-swath power balancing factor that can yield the most flattened backscatter profile along the range direction over homogeneous areas to modify the ESA-provided noise vectors and then reconstruct a complete 2D noise field. The noise removal function is given as follows:

C. The Denoising Method Proposed in
The terms in parentheses are the reconstructed 2-D noise field.
, is the optimal noise scaling factor, and , is the inter-swath balancing power factor.
, and , represent the raw sigma naught and the noise-subtracted sigma naught, respectively.
1, … ,5 is the sub-swath order from near to far range, indicating that denoising is performed in each sub-swath. Referenced values of the two factors of different IPF versions are given in [13]. Fig. 2(a) shows an EW HV-polarized image acquired over the Arctic MIZ, which is radiometrically calibrated using (2), and Fig. 2(b) is the noise-subtracted image using the method proposed in [13]. In Fig. 3(a), the profiles of the original signal, NESZ and the scaled NESZ noises are shown. The of second sub-swath (EW2) is less than 1.0, while the of other sub-swaths are all larger than 1.0. Fig. 3(b) shows the scaled noise-subtracted and balanced range profiles, in which the alignment of red lines (linear fit lines of scaled noise-subtracted signal) is the balancing process described in [13]. The profiles in Fig. 3(b) seem to have good scaling and alignment. However, in the noise-subtracted image shown in Fig. 2(b), excessive scaling appears in the last three sub-swaths, whereas inadequate scaling appears in the second sub-swath, and balancing is conducted inaccurately in all sub-swaths that lead to new nonhomogeneity among open water or sea ice. This inconsistent performance of denoising over all range profiles of the HVpolarized image is because the radar backscatter of the sea ice is much larger than that of the open water, so that ice dominates the average range profile in each sub-swath and the small values in open water are covered. Thus, we found that in the ice-water mixed image, the factors and calculated based on the open water regions are more accurate. In the algorithm proposed in [13], the entire S1 EW images are divided into five azimuthal sub-images and the average range profile is acquired for each sub-image to calculate the optimized noise scaling factor and balancing factor . The authors discarded profiles having powers higher than the estimated noise power by 3 dB, which selects the calm open water sub-images that of low backscatter, i.e., one type of homogeneous surface to acquire the optimized factors. The averaged values of and for each sub-swath from IPF 2.4x to 2.7x were provided. Then, these averaged values of and are used for the NESZ correction of all other images, which is similar to the NESZ vectors that are measured by averaging the power profile from calm water images and then used for thermal noise subtraction for all other images [1], [16], [31]. When one uses these average parameters for the EW images processed by one IPF version, it is evident that the noise vectors are not accurately scaled and balanced in some cases. As shown in Table I in [13], the fluctuation of scaling factors can be up to 0.07, while it is found a slight variation of 0.02 can lead to over-scaling or under-scaling. For the values of the balancing factors in [13], their fluctuations are much greater. Thus, scaling and power balancing parameters derived for each single image should be more reasonable and accurate than the averaged values.

D. Multiplicative Noise Correction
After additive noise subtraction, the multiplicative noise in S1 EW images still exists, which is particularly evident as large deviations in the center of the first sub-swath and inter-swath boundaries, as shown in Fig. 2(b). This multiplicative noise significantly disturbs the signal, and the effect is even "enlarged" in the texture images in some cases (shown later in Section III). Therefore, further removal of multiplicative noise from the additive noise-subtracted EW images should be included in the full processing chain of the denoising process.
We use the sketch map shown in Fig.4 to illustrate the concept of further removing the multiplicative noise in the additive noise-subtracted images. Fig. 4(a) shows an ideal radar signal of the first sub-swath of the additive noise-subtracted HV-polarized image, where the sigma naught values in the boundaries exhibit large variations, indicating high multiplicative noise, while the signals in the low-amplitude zones exhibit little multiplicative noise. If the large amplitude signal is compressed to be consistent with the low amplitude, the amplitudes along the entire range direction will become uniform. One possible process is to identify a profile (e.g., the solid line shown in Fig. 4(b)), which has a contrasting trend to the envelop (dashed line in Fig. 4(a)) of the sigma naught variation. By multiplying the profile of the solid line with the sigma naught profile, a flat profile along the range direction can be achieved, as demonstrated in Fig. 4(c).
In the following section, the improved additive noise subtraction method and the proposed multiplicative noise subtraction approach based on the above description are presented.

III. METHODOLOGY
The proposed method of removing both additive and multiplicative noise includes three major steps: determination of the noise scaling factor, power balancing factor and removal of multiplicative noise, which are all conducted based on the S1 images in linear unit. The detailed processing and explanations are given below.

A. Determination of the Optimal Noise Scaling Factor
Computation of the scaling factor includes four steps, as shown in the flow chart presented in Fig. 5. First, an EW HVpolarized intensity image is divided into numbers of evenly spaced blocks along the azimuth direction in each sub-swath, and each block has a size of 200 azimuth lines. The red rectangles shown in Fig. 2(a) present two examples of such blocks. The noise scaling factor , for each block is calculated using the functions (A-10) to (A-12) given in the Appendix.
Second, the homogeneous blocks are identified using a variance factor, which is calculated using the following formulas: is the variance in the deviation in block . is the range profile of the un-denoised sigma naught, and _ is obtained by shifting the scaled NESZ vectors to the un-denoised sigma naught profile.
where • is the averaging operator of all pixels in block .
_ is the scaled NESZ vector, and is the difference between the global averages of and _ . Fig. 6 (a) and (b) show the range profiles of calm open water and ice-water mixed regions, respectively, corresponding to the two red rectangles in Fig. 2(a). In the plots, the un-denoised sigma naught value profile is shown by the black line, from which the deviation of is obtained (the light gray line at the bottom) by subtracting the approximate mean values _ (red line). Additionally, the dashed line denotes the scaled NESZ vectors acquired by the optimal scaling factor , , and is the  To find the threshold of a variance factor that can discriminate homogeneous and inhomogeneous blocks, we selected 100 S1 EW images (acquired in August and November of 2018) which represent five types of sea ice and open water conditions, i.e., the calm open water, windy/rough open water, thin ice, thick ice, and mixture of thick ice and open water. The variance factors in the blocks of these images were subsequently calculated. Fig. 7(a) shows the combined histograms of the calculated values of the 100 images. Then we fitted the histograms of the values of the five different types of images representing various sea ice and open water conditions, which are presented by the curves in different styles in Fig. 7(b). Note that the histograms of the values of calm (the first peak) and rough sea surface (the second peak) were combined to a double-peak curve, as presented by the solid black curve. It is clear that the values of the blocks in different open water and sea ice conditions are generally discriminated. In this study, the threshold of (denoted ) is chosen to be 10 . (approximating 8 10 ). If the in blocks are below the threshold, the blocks are regarded homogeneous ones (i.e., their histograms are represented by the solid and dashed black lines in (b), which generally corresponding to calm open water, windy open water, young ice, and level first-year ice), and they are used to calculate the scaling and balancing parameters.
Following the identification of the homogeneous blocks, we compute the mean values of , for these homogeneous blocks to obtain one optimal , for each sub-swath as follows: , , | Here, "{}" in equation (7) denotes an assemble that contains all scaling factors acquired in the homogeneous blocks in each subswath.
If the variance factors in one sub-swath are all larger than the threshold, we directly average all the , values to obtain , in a sub-swath. The , values of the five sub-swathes of the case are 1.62, 0.85, 0.96, 1.04, and 1.04, respectively. Fig.8(a) shows the scaled noise subtracted image using (3) and (9), where we can see that the noise floor is cleanly subtracted. However, the intensity of EW2 is slightly higher than others, which can be amplified more in the texture image. Thus, the next step is to  Fig. 2a. The red rectangle in region A is the 25 th image block in the first sub-swath EW1, which has an optimal scaling factor of 1.58, and the variance is 1.4156  10 -8 . Region B is the 25 th image block in the second sub-swath EW2, in which the optimal scaling factor is 0.8 and the variance is 1.1609  10 -6 . determine the power balancing parameter , .

B. Determination of the Power Balancing Factor
The principle of inter-swath balancing is to retain the sigma naught values of the noise-subtracted EW image at one level by aligning the intensity values at the boundaries. In [13], the way the power balance factor , is determined based on the linear fit lines of the range profiles. As shown in Fig. 3(b), the difference between the two endpoints of the two linear fit lines (red lines) in the first inter-swath boundary (i.e. the boundary between EW1 and EW2) is the initial balancing power , that is chosen in [13]. However, the linear fit is not applicable for balancing sub-swaths covering thick ice or both sea ice and open water, and the highly inaccurate balancing can cause more additive disturbances in the open water regions, as shown in Fig.  2(b). To solve this problem, we optimize the method of determining the power balancing factor , by using the fundamental principle of boundary pixels' mean alignment. Similar to the method of determining the noise scaling factor , presented above, computation of the factor , is also based the variance factor computed in each block. Detailed steps of determining the inter-swath power balancing factor , are given below. First, the original balancing factor , of each sub-swath is assigned to 0 to construct the modified 2D noise field _ , i.e., using only the noise scaling factor determined in the previous step to reconstruct the noise field.
According to the S1 Product Specification file [32] and product annotation file, the swath merge type of S1 EW data is 'Varying Optimal Range Cut Line'. By reading the firstRangeSample and lastRangeSample values in the annotation file, we first find the minimum lastRangeSample and the maximum firstRangeSample to confirm the merge boundary (illustrated in Fig. 8(b)). In addition, we found that the pixels on the left side close to the minimum lastRangeSample (dotted line in Fig. 8(b)) have jump deviations that can highly disturb the accuracy of average values balancing. Thus, we further discard the 100 pixels (marked in the Fig. 8(b) by the doublehead arrow) in the sub-swath on the left. The pixels between two solid black lines in Fig. 8(b) represent the actual boundary between two neighboring sub-swaths. To obtain the average values of both sides of pixels, 20 pixels from an image block in the left sub-swath (e.g., EW1) and another 20 pixels from the right sub-swath (e.g., EW2) of the boundary are extracted, which are denoted as 1 and 2 , respectively (two red lines shown in Fig. 9). For 1 and 2 , the mean sigma naught of _ , 1 and _ , 2 are derived. Thus, we have the power balancing factor of the jth image block by using (10). Note here that the sub-swath index number is from 2 to 5. indicates the image pixels along the azimuth direction in block .
We use a threshold with an empirical value of 0.002 to limit the balancing power Similarly, when , calculated by (13) is null or no block has a larger variance factor than the threshold, we calculate the average of , in all blocks to obtain , . After deriving , and , , we can use (9) directly to construct the modified 2D noise field , .
However, because the balance factor , is an average value of the image blocks in one sub-swath, which may include various heterogeneous surfaces, it is not accurate enough to balance the whole EW image by only one iteration. Therefore, the new balanced modified 2D noise (9) is used to repeat the process from (10) to (13) until the determined inter-swath balancing power , is less than 10 . Finally, to maintain the radiometric consistency, the adjustment of total power in the modified noise field is also carried out using the following formula: This formula is similar with (A-15). However, we replaced the single in (A-15) using • including noise vectors in both the azimuth and range directions, which we think is the more complete noise field close to the actual noise level. Moreover, because the noise level of the first sub-swath is especially high, as presented in Fig. 1, the second item on the right side of (14) is based only on the value differences over the last four sub-swaths. Then, the complete noise-subtracted image can be acquired using (3), as shown in Fig. 10(a).

C. Further Removal of Multiplicative Noise
After the noise subtraction and power balancing processes described above are completed, the EW HV-polarized image in Fig. 10(a) seems to have a good quality. However, the corresponding homogeneity texture image calculated based on GLCM shown in Fig. 10(b) suggests that the effect of multiplicative noise is "enlarged".
Multiplicative noise removal processing is briefly described in section II that one needs to determine a trend profile to compress the large variations in the additive noise subtracted signals, which are indeed induced by the multiplicative noise. We found that the variable trend of the additive noisesubtracted signal is similar to the variations in the ESAprovided NESZ vectors ( ). Therefore, we decided to scale the NESZ vectors to an appropriate amplitude to represent the envelop of the additive noise-subtracted radar signal variations (as illustrated by the dashed line shown in Fig. 4(a)). In fact, the multiplicative noise primarily exists in the boundaries of subswaths, with a border zone in the back boundaries. We identify the "open water" blocks based on the local mean SNNR values calculated over a 200  200 pixels window nearest to the back boundaries to conduct the signal correction. The average SNNR threshold of 2.0 (demonstrated in [29]) is adopted to estimate the boundaries of one block being open water or not. If the local mean SNNR is less than 2.0, the local region in the back boundary are open water, taking the block A in Fig. 2(a) as an example, the functions and processing can be described in detail as follows: 1) Fig. 11(a): The deviation (solid dark gray curve, the same in other plots in Fig. 11) is the difference between the additive noise-subtracted sigma naught profile (light gray curve) and the fitted mean profile (dashed dark gray line, denoted as ). The is obtained by applying two times of the box-car smoothing operation. 2) Fig. 11(b): The purpose of this this step is to estimate the relative amplitudes of the NESZ and its deviation, depending on the positions of the boundary and minimum in a sub-swath, as marked in the plot. The two regions of both 100 pixels widths are located based on the profile (brown line). Then, the local averaged values of in the boundary and the minimum pixels (highlighted by the cyan lines) are calculated, which are denoted as Here, is the absolute value operator. Using the local average SNNR, the local regions before the back boundary can be identified as open water or thin ice, but this is only valid in the full block in a sub-swath. Therefore, we adopted additional judgement of the difference factor in (15). If the difference factor is less than 0, this block is regarded as a dense ice region, and it is not necessary to conduct the multiplicative noise removal in the following steps. Otherwise, this block is regarded as an open water region, and further multiplicative noise removal is performed. 3) Fig. 11(c): is multiplied by the difference factors to obtain the scaled noise profile (brown line): Compared to in Fig. 11(b), the amplitude of scaled noise in Fig. 11(c) is in the similar scale with the deviation . We obtained the ratio value of 0.0172 for this block. Fig. 11(d): The magnitude of is normalized by:

4)
It is represented by the brown line in the plot. We take 0.1 as the minimum to avoid the exponent arithmetic values being negative in the next step. 5) Fig. 11(e): The is further modified by using the logarithm of the ratio factors as an exponent to obtain the final-modified noise profile (brown line): After modification, the is on the proper size for removing the deviation floor. The constants 0.1 in (17) and the ratio factors ( , ) jointly determine the shape of . The term 1/ 2 represents the total scaling of deviation . We magnified the range of 400 pixels around the minimum value in Fig. 11(d) and (e) for comparison, suggesting that the shape was further changed (less curved) by the exponent. We obtained the exponential term abs( ⁄ ) in (18) as 0.4839 for this block. 6) Fig. 11(f): The deviation is multiplied by the reciprocal of to obtain a uniform deviation profile (light blue line):

1/
( 1 9 ) Until this step, the scaled deviation becomes homogeneous in the range direction. However, the amplitude of the scaled deviation becomes larger or smaller than the original amplitude. In Fig. 11(f), the amplitude of is approximately 1.5 times the amplitude of . 7) Fig. 11 (g)  at the minimum region (same as that in Fig.  11b), respectively.
is the ratio of two mean values, and we obtained the value of 0.67 for this block. Moreover, is the mean line of the noise-subtracted sigma naught values, same as in Fig. 11(a). The first term on the right side of (21) is the corrected deviation of the normal original amplitude (light blue line in Fig.11(g)), and _ is the final corrected signal (light gray line in Fig. 11(h)). The effectiveness of multiplicative noise removal can be clearly found in comparison of Fig. 11(a) and (h). After conducting the independent processing in each block, the inconsistency between blocks appears; thus, we calculate the average values of the modification parameters (see Table I), i.e. the difference ratio in (16), the logarithm ratio in (18), and the in (20), to apply to all homogeneous blocks again in each sub-swath. Fig. 12 shows the multiplicative noise removed intensity image and the consequent GLCM homogeneity texture image, which look much better than the one shown in Fig. 10, suggesting that the further removal of multiplicative noise performs well.
Balancing factor , Variance factor , or Threshold of the variance factor (7) The difference ratio ⁄ ( 1 6 ) The logarithm ratio ⁄ The regression ratio ( 2 0 ) Table I lists the main parameters used in the proposed method to denoise both additive and multiplicative noise. As the scaling and balancing factors appear a few times in different questions and here only list the equations where they appear in the first time.

IV. RESULT
In this section, we present a few examples showing various sea ice conditions and various sea states, to further demonstrate the effectiveness of the proposed denoising algorithm.

A. Denoising EW HV-polarized Images Processed by Different IPF Versions
To validate the proposed denoising algorithm, hundreds of images acquired in different seasons with varying sea state and sea ice conditions are tested. Here, as shown in Fig. 13, we selected four examples of EW HV-polarized images processed by IPF versions 2.7 to 3.1, which also present various sea ice and wind conditions, to demonstrate the denoising results. Fig.  13 (a), (c), (e) and (g) are the additive noise-subtracted images using the ESA-provided noise vectors. Obviously, the main problem of the images is the large residual noise that affects the image GLCM homogeneity, especially in the first sub-swath. By applying the proposed method, both additive noise and multiplicative noise are all well removed, as presented in Fig.  13 (b), (d), (f) and (h). Fig. 13a shows a scene imaging a large area of floating ice, with a maximum wind speed of 8 m/s (derived from the ERA-5 reanalysis model; the same model is used for other examples). The low backscatter values of sea ice regions in the upper-right panel are close to those of open water having a low SNR. Thus, the additive noise modification and multiplicative noise correction are conducted over nearly the entire image.
The second example shows a distinct boundary between sea ice and open water, where the region has a maximum wind speed of approximately 10 m/s that causes some wind streaks in the image. Under this condition, the maximum variance factor of open water blocks exists in the first sub-swath with a value smaller than the threshold . Therefore, the open water blocks are all detected to be homogeneous. Based on that, the noise scaling and balancing parameters of the denoising process are calculated more accurately.
The third example presents a case of highly mixed sea ice and water, where the brash ice spreads widely over open water. After subtracting the additive noise using the ESA-provided noise vectors, one can find that the residual noise is particularly large in the center and boundary of the first sub-swath because large open water regions are present there. By applying the proposed method, the image is well balanced, and the multiplicative noise is effectively removed.
With respect to the fourth case that is processed by the newest IPF version 3.10, the noise is reduced substantially using the ESA-provided noise vectors, as shown in Fig. 13(g), suggesting that the current IPF version has significantly improved the accuracy of additive noise estimation. However, slight unbalance and distinct multiplicative noise still have some impacts on the sub-swath boundaries. The effect is more evident in the texture image (not shown). Fig. 13(h) shows the final denoised image using the proposed method, in which multiplicative noise as well as the additive noise are removed, though the scaled and balanced noise vectors are slightly different from the original noise vector that provided by ESA in annotation files. The scaling parameter is 1.06, 0.99, 0.93, 1.18, and 1.11, and the power balancing is -1.99, -0.00, -0.65, -2.47, and 2.94 (all with scale of 10 -4 ). Moreover, the wind speed in this case exhibits large spatial variations, ranging from 8 m/s to 24 m/s, causing a rough sea surface of high radar backscatter, but the algorithm performs well at reducing the noise.
In summary, the presented four cases suggest that the proposed algorithm performs well and is highly adaptability for different circumstances of sea ice and open water.

B. Comparison of the Proposed Denoising Method with the
Methods in [13] and [29] In this sub-section, we present one case to compare the denoising using our proposed method and that using the methods in [13] and [29]. Because the scaling and balancing parameters provided in [13] are for the EW images processed by the IPF prior to version 2.7x, and the NESD model coefficients provided in [29] are derived based on the EW data processed by IPF versions 2.5~2.9, we chose one image processed by IPF version 2.72 (Fig. 14) that can be denoised by directly using the parameters given in [13] and [29]. Fig. 14(a) and (b) are denoised results (both additive and multiplicative noise removal) by using the proposed method in this paper and combining the methods proposed in [13] and [29], respectively. The corresponding GLCM homogeneity texture images are shown in Fig. 14(c) and 14(d), respectively. A visual inspection suggests that the denoised result using the proposed method here better solves the problems of unbalance and large deviations of boundaries between sub-swathes, which is further indicated in the texture images. Nevertheless, both methods indeed significantly improve the quality of the HV-polarized images in the EW mode. The proposed method has high computing efficiency because its calculation is based on the divided blocks, while the method proposed in [29] relies on much smaller sliding windows with a size of tens of pixels. Furthermore, high adaptivity is implemented by the processing conducted in any single image with independently calculated parameters. Meanwhile, it should be noted that our method does not handle the residual scalloping noise along the azimuth direction, although this has smaller impacts on the ice classification methods or other ice detection applications.
From the previous images, the macroscopic noise removal of  [13] and [29]. (c) and (d) Corresponding GLCM homogeneity texture images. The black regions are land masks, and the red rectangle is used in Fig.  16.
(Image ID: S1A_EW_GRDM_1SDH_20170125T184255_20170125T184359_014992_0 187AF_9340.SAFE (IPF 2.72) the entire image has been experimentally demonstrated. However, at local scales, the changes in the sigma naught and standard deviation values also need to be analyzed. According to the illustration in [29], under a low SNNR, the image contrast is strongly decreased in the denoised image, causing a blurry appearance. To better illustrate the changes during the multiplicative denoising process, we analyzed a block marked by the red rectangle in Fig. 14(a). Fig. 15(a), (b) and (c) show the intensity values after additive noise subtraction, both additive noise and multiplicative noise subtraction using the proposed method and that using the method proposed in [29] in this block. To highlight the differences among different noise subtraction methods applied to the block, the corresponding standard deviation values are shown in Fig. 15(e), (d) and (f), respectively. From the intensity images of the single block, it is hard to observe the differences between the multiplicative noise-removed images in Fig. 15(b) and (c) and the only additive noise-subtracted image in (a), which, however, is better highlighted in the standard deviation images. Comparing Fig.15(e) (result based on the proposed denoising method) to (d) and (f), it is found that the multiplicative noise is well removed not only at the front boundary (right part the rectangle) but also at the back boundary (left part of the rectangle). This single case comparison may suggest that the denoising algorithm proposed in this study has better performance than others for some conditions.

C. The Variation in Sigma Naught between Neighboring Denoised S1 EW Images
Demonstration of the denoising effect of the proposed method is based on a single image. In this subsection, we used an example to verify the variation in sigma naught between the neighboring images covering a large area. Fig. 16 shows a mosaic images consisting of 4 EW scenes in HV polarization, acquired in the Barents Sea on December 13, 2018. All the images are denoised using the proposed method. The case shows a scenario of mixture of open water and sea ice in the MIZ. The visual inspection suggests that the noise is well subtracted for various sea ice and open water conditions in a large coverage, and both fine features of the sea ice and open water are well reserved. The gray line (in fact, a rectangle with a width of 100 pixels) in the mosaic image is across two neighboring scenes. The corresponding mean profile of the denoised sigma naught along the line is shown in Fig.16(b). Note that the two scenes are overlapped between approximately 6500 to 10329 pixels counting from west to east. The red line in (b) marks the far range boundary of the EW scene in the east. The profile suggests that there is no signal jump of sigma naught in the frame boundaries between the two neighboring scenes. The difference of the sigma naught in a small area (200 100 pixels) at the right side of the red line between the two scenes is 0.0024 in linear unit. Note that the difference of the local incidence angles in the small area between the two scenes is 15.88 degree (30.43 degree versus 46. 31degree). This large difference of local incidence angle naturally contributes to the discrepancy of sigma naught in the same area, though to some extent in cross-polarization data.
The case further indicates that the denoising processing, though conducted independently from scene from scene, does not induce additional artificial effects on the original sigma naught values.

V. SUMMARY AND CONCLUSION
The S1A and S1B EW mode data in combination of the HH and HV polarizations have been widely and frequently acquired in the Arctic for intensive sea ice monitoring since these satellites were launched. Although the data in HV polarization have great advantages in sea ice detection, they are strongly affected by various noise. While the scalloping effect in the azimuth direction can be removed using the AAEP method or the ESA-provided noise vector (since the processing version IPF 2.90), removal of the thermal noise floor presenting in the range direction of the HV-polarized EW image still faces many challenges. In this paper, we present an effective denoising algorithm based on previous studies. Some of the conclusions are discussed below.
The denoised algorithm proposed in [13] represents important progress on improving the image quality of EW data in HV polarization. The major concept involves using a scaling factor and power balancing factor to modify the ESAprovided noise vectors. However, we found that this algorithm does not perform effectively in many cases when using the average parameters or using the algorithm globally for a single image. This is because circumstances of sea ice and open water are complicated, and EW HV-polarized images often present a high mixture of sea ice and open water, while our aim is to detect sea ice cover at a high spatial resolution in MIZ. Different from their algorithm used to calculate the average values of the two factors from segmented azimuthal blocks, which is onefifth of all azimuth lines, we segment the image into more azimuthal blocks (approximately 50 blocks in each sub-swath) and introduce a variance factor to discriminate between homogeneous and inhomogeneous blocks to derive a more accurate scaling factor and balancing factor based on the local homogeneous regions in each sub-swath. This method enables the denoising process to be carried out accurately in each image alone.
Although the EW HV-polarized image quality has been significantly improved based on the improved noise subtraction process, the residual multiplicative noise still has great impacts on the noise-subtracted images, which is particularly evident in the image texture features; these features are informative information for sea ice detection in many scenarios. This should be the reason that the authors of [13] published another paper in [29]. In our study, we found that the deviation of the multiplicative noise along the range direction has a similar trend to the ESA-provided NESZ profile. Therefore, the deviation floor of residual noise is corrected by using the modified NESZ profile. The results prove that the method of multiplicative noise removal is effective, as is particularly evident in texture images.
We applied the denoising algorithm to hundreds of EW images to test its performance. Some examples presented in this paper demonstrate that the proposed algorithm can handle data not only processed by different IPF versions but also representing various open water and sea ice states over the Arctic MIZ. We further compared the denoising results using the proposed method with those from the method in [13] and [29]. In the overall intensity and texture images, both methods significantly improve HV-polarized image quality. Our method better removes the additive noise based on more accurately calculating the scaling and power balancing factors. With respect to the multiplicative noise removal, the proposed method has better performance in the inter-swath boundaries of EW images. Furthermore, we conducted the denoising process in segmented blocks rather than sliding windows with a size of tens of pixels, as was performed in [29], which can improve the computing efficiency to some extent.
In the full processing of denoising, few parameters, e.g., the threshold value of for screening homogeneous regions, is determined based on the overall evaluation and statistical analysis of hundreds of EW images. Some exceptional cases may exist that do not satisfy these empirical parameters. Due to the assumption that multiplicative noise removal is conducted in separated blocks along the entire range direction in each subswath, there may be some multiplicative noise remaining in icewater mixed blocks. This requires further studies.
In summary, the proposed denoising method generally has good performance on removing both additive and multiplicative noise presenting in the EW HV-polarized images. This can benefit many applications on monitoring sea ice in polar regions based on S1 cross-polarization data in wide swath. One of such application on extracting sea ice cover by EW HV-polarized data is introduced in [33].

ACKNOWLEDGEMENT
We would like to thank ESA for providing the S1 data available for worldwide users. We particularly would like to thank the anonymous reviewers for their constructive comments and suggestions which are of great help for improvement of this manuscript.

APPENDIX
The denoise method proposed in [13] for the EW HVpolarization image first reconstructs a complete 2D noise field using two parameters: a scaling parameter and a balancing parameter, and then subtracts the noise field from the raw sigma naught values. This method can be divided into three steps: 1) azimuth descalloping; 2) noise power scaling and inter-swath power balancing; and 3) local residual noise power compensation.
Although the parameters needed in this process can mostly be found in product annotation files and user manuals, there are still some ambiguous parameters that are difficult to find. To reduce the workload on reproducing this method, we specified several intermediate variables and processes to better understand the method described in [13].,

A. Azimuth Descalloping
The azimuth descalloping function is defined by a two-way AAEP, which changes with the antenna steering angle [30]: where is the antenna element length, is the radar's wave length, and is the antenna steering angle. However, instead using the AAEP function, the authors in [13] used the AAEP vectors in auxiliary data provided by ESA through the Sentinel quality control webpage (https://qc.sentinel1.eo.esa.int). The AAEP lookup table values are provided corresponding to antenna steering angle increasing from negative to positive values at a certain interval, of which at the middle the steering angle is zero.
First, the angle for each focused burst azimuth time must be calculated as follows [30]: where is the focused burst azimuth time, is the satellite velocity, is the azimuth frequency modulation rate, and is the Doppler Centroid rate introduced by the antenna scanning.
can be calculated as follows [15], [30]: where is the antenna steering rate, with a unit of rad/s [15]. is computed by the function in the annotation file. Note that there is usually confusion between these two parameters [30], but here, they have been clarified after our experiments.
The focused burst azimuth time relies on the burst center and boundary positions, which can be computed from the times in the "Antenna Elevation Pattern Data Set Record", i.e., "antennaPatternList" in the Level 1 Product Annotation file of the GRD product. Detailed functions can be found in [13]. Note that the parameter is the sampling rate in the azimuth direction, i.e., "azimuthFrequency" in the annotation file. is the number of azimuth lines for each sub-swath of the input image, i.e., "numberOfInputLines".
is the number of bursts for each sub-swath of the input image, i.e., length of "azimuthFmRate List".
After acquiring the antenna steering angle , the descalloping gain can be calculated one burst at a time as follows [13], [30]: 1) Identifying the two boundary angles in the AAEP LUT, and are determined as follows ( A -5 ) 2) The AAEP values (in dB) are computed for the focused burst azimuth times of each azimuth line , , … , , from the AAEP LUT values corresponding to angle samples between and at a certain interval by linear interpolation: : , , is converted into the linear scale: The descalloping gain , which varies only along the azimuth lines in one sub-swath, should multiply the noise field obtained from the ESA-provided range noise vectors to remove thermal noise in both the azimuth and range directions.
The complete noise removal function is given as follows: The terms in the parentheses are the reconstructed complete 2D noise field, and the parameters have been interpreted in (3) in Section II.

B. Noise Power Scaling and Inter-swath Power Balancing
For reconstructing the complete 2D noise field, two factors are introduced in [13], namely, the optimal noise scaling factor for the optimal noise scaling factor, and for the inter-swath balancing power. Here, is the weight factor acquired by the absolute values of the gradient of , i is the range cell number, ̂ is the linear fit model for the denoised signal , and and are the slope and intercept of the linear fit model for the sub-swath index number = {2, …,5}. In the selection of the optimized scaling factor , the range of varies from 0.5 to 2.0. Another detailed explanation of each parameter can be found in [13].
The authors in [13] use the flatness of range profiles to acquire the optimal noise scaling factor and the inter-swath balancing power in each sub-swath. In functions (A-10) to (A-15), the linear fit model limits the profiles as homogeneous (open water or dense pack ice). Therefore, the two parameters cannot be accurately obtained for each image, including inhomogeneous (ice-water) images.

C. local residual noise power compensation
To remove the negative power values generated by the noise subtraction process and retain the total power balance locally, the local SNR-dependent radiometric correction developed by Balss et al. [28] was adopted in [13].