High-precision ASI sea ice concentration inversion method based on CGAN for data correction

Abstract

The present invention provides a high-precision ASI sea ice concentration inversion method based on CGAN for data correction, belonging to the field of satellite remote sensing technology, and comprising the following steps: S1. Data screening: find a relatively stable relationship between the 89 GHz data and the 36 GHz data that are not disturbed or slightly disturbed by external environments such as clouds and water vapor, and screen out the 89 GHz data with large interference; S2. Data correction: the generation network and the discriminant network of CGAN are trained adversarially, and the stable relationship between the 36 GHz data with high reliability and the 89 GHz data that are not disturbed or slightly disturbed by external environments such as clouds and water vapor is used to correct the 89 GHz data with large interference; S3. Based on the corrected 89 GHz data, the ASI sea ice concentration algorithm is used to invert the sea ice concentration. The present invention can effectively correct the mixed pixel sea ice concentration and greatly reduce the error caused by the atmosphere.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite remote sensing, and particularly relates to a high-precision ASI sea ice concentration inversion method based on CGAN for data correction. Background Art

[0002] The changes of polar sea ice are closely related to global atmospheric changes, ocean circulation changes, etc., and play an important role in the study of global climate change. Sea ice concentration (SIC, sea ice concentration) is one of the important parameters for studying the spatio-temporal variation characteristics of sea ice.

[0003] In recent decades, many inversion algorithms have been proposed to estimate more accurate SIC values. Comiso (1986, 1995) proposed the bootstrap algorithm, which mainly uses the variation of sea ice emissivity with frequency band and different physical properties of seawater to invert sea ice concentration. Cavelieri et al. (1984, 1991) proposed the NT algorithm, which establishes an inversion equation for sea ice concentration based on the difference in radiometric brightness temperature and combined with the polarization gradient rate and spectral gradient rate. In 2015, Liu et al. (2015) proposed an algorithm for calculating Antarctic sea ice concentration using the total constraint least squares algorithm based on the NT algorithm. Cavalieri et al. (1991) used the brightness temperature data of SSM / I to propose the NASA Team algorithm, which can be used for the inversion of the concentration of first-year ice and multi-year ice. Markus et al. (2000) added the 89GHz vertical and horizontal polarization brightness temperature data on the basis of NASA Team and proposed the NASA Team2 algorithm. Lomax et al. (1995) used the SSM / I 85.5GHz data to propose the Lomax algorithm to calculate sea ice concentration. Hao (2015) improved the existing NASA Team algorithm by introducing the AMSR-E 6.9GHz data to improve the calculation accuracy of multi-year ice results. Kern and Heygster (2001) proposed the SEA LION algorithm, which uses the polarization value in the 85.5GHz frequency band to invert sea ice concentration, and uses a radiative transfer model and atmospheric data obtained from numerical weather prediction models and / or SSM / I measurements to correct the polarization difference. The ASI algorithm (ARTIST sea ice algorithm) was developed in the 1998 project "Research on Arctic Radiation and Turbulent Exchange". Svendsen et al. based on the concept of "polarization-corrected temperature" and used data near the 90GHz band to invert sea ice concentration (Spencer et al., 1989; Svendsen et al., 1987). Kaleschke et al. improved the algorithm proposed by Svendsen et al. and used higher-resolution SSM / I 85GHz data to conduct a mesoscale numerical simulation of the atmospheric boundary layer at the Arctic sea ice edge (Kaleschke et al., 2001). One advantage of this algorithm is that compared with other algorithms using 85GHz band data, it does not require additional data input, can directly invert the results according to the measured brightness temperature data, and has similar results to the sea ice concentration algorithms using other channels (Kern, 2004). Spreen et al. (2008) applied the ASI algorithm to AMSR-E data and obtained the corresponding sea ice concentration calculation formula. Wang (2009) proposed a method for calculating the concentration of multi-year ice based on the different characteristics between first-year ice, multi-year ice and seawater in the 89GHz band.Su et al. (2013) conducted a series of experiments, including the interpolation algorithm "join points" based on the ASI algorithm and weather filters. Zhang et al. (2012) proposed a method for calculating sea ice concentration using multi-band and dual polarization based on the radiation characteristics of sea ice and sea water.

[0004] The resolution of satellite data and the inversion algorithm are crucial for accurately providing sea ice concentration. Currently, the microwave data product with the highest resolution is the 6.25 km resolution sea ice concentration gridded product retrieved by the University of Bremen, Germany, using the ASI algorithm based on the 89 GHz high-frequency data of the Advanced Microwave Scanning Radiometer - Earth Observing System (AMSR-E) (Spreen et al., 2008). Although the ASI algorithm has advantages, compared with low-frequency data, the 89 GHz band data is more affected by atmospheric cloud liquid water, water droplets, water vapor, and the density of snow particles on the ice surface. Especially in thin ice areas, when the liquid water content in the cloud is high or a cyclone passes by, the inversion results have large errors; at the same time, it is also very sensitive to the density of snow particles on the sea ice surface area; therefore, the ASI algorithm requires weather filter processing (Spreen et al., 2008). Although using weather filters can eliminate some errors, the concentration of pixels cannot be modified during mixing. To obtain more accurate sea ice concentration results, it is necessary to screen and correct the 89 GHz data that is greatly disturbed by external environments such as clouds and water vapor. Summary of the Invention

[0005] In view of this, the present invention provides a high-precision ASI sea ice concentration inversion method based on CGAN for data correction, which can greatly reduce the errors caused by the atmosphere by effectively correcting the sea ice concentration of mixed pixels.

[0006] To solve the above technical problems, the present invention provides a high-precision ASI sea ice concentration inversion method based on CGAN for data correction, including the following steps:

[0007] S1. Data screening: Find a relatively stable relationship between the 89 GHz data that is not affected or slightly affected by the external environment and the 36 GHz data, and screen out the 89 GHz data that is greatly disturbed.

[0008] S2. Data correction: The generator network and discriminator network of CGAN are trained adversarially, and the stable relationship between the highly reliable 36 GHz data and the 89 GHz data that is not affected or slightly affected by the external environment is used to correct the 89 GHz data that is greatly disturbed.

[0009] The model function of CGAN is shown in formula (3):

[0010] (3)

[0011] Among them, is the data affected by the external environment, is the additional information, is the input random noise, is the data not affected by the external environment output after inputting the random noise, additional information, and the data affected by the external environment into the generation network of the CGAN, is the probability that the discriminant network determines that the input data is fake; since the goal of the generation network of the CGAN is to make the generated data as close as possible to the data not affected by the external environment, the loss function is set to to ensure that the probability of the discriminant network outputting a fake image is as small as possible, and the goal of the discriminant network is to improve the ability to judge the difference of the input data. Therefore, the larger the better, and at the same time, it is hoped that the noise influence is as small as possible, so the loss function is set to , using to represent this game process;

[0012] The training process of the CGAN model is as follows:

[0013] (1) Before training, rotate and translate the dataset to increase the number of datasets, and then normalize the training set and test set;

[0014] (2) Input the training set into the generation network, and then perform continuous batch normalization + convolution + ReLU + pooling operations to complete the downsampling operation;

[0015] (3) Perform continuous operations of deconvolution, BN, and activation function on the feature map obtained by downsampling to complete the upsampling;

[0016] (4) In the same network layer, connect the output feature maps of downsampling and upsampling; in different network layers, from the top neural network layer to the bottom neural network layer, fuse the output feature maps of downsampling and the upsampling of the next neural network layer. The fused output feature map continues to be connected by the upsampling feature map of the next neural network layer, and this iteration continues until there is no corresponding upsampling in the next layer, and then obtain the mapping relationship between 89 GHz data and highly reliable 36 GHz data;

[0017] (5) Then input the test set data, that is, the relationship between the undisturbed 89 GHz data and 36 GHz data, and the mapping relationship between 89 GHz data and highly reliable 36 GHz data obtained in step (4) into the discriminant network;

[0018] (6) Then, perform the operations of continuous batch normalization + convolution + ReLU + pooling to complete the downsampling operation;

[0019] (7) Finally, use cross-entropy to discriminate the results obtained by the discriminant network. If the loss function reaches the minimum value, output the corrected 89 GHz data; otherwise, return to step (2) and repeat the above steps until the loss function reaches the minimum.

[0020] S3. Based on the corrected 89 GHz data, use the ASI sea ice concentration algorithm to retrieve the sea ice concentration.

[0021] Further, the method for screening 89 GHz data with large interference is specifically as follows.

[0022] (1) Under clear weather conditions, plot the polarization ratio scatter plot with the polarization ratio of 36 GHz data as the abscissa and the polarization ratio of 89 GHz data as the ordinate.

[0023] The polarization ratio formula is as shown in (1).

[0024] (1)

[0025] Among them, , are the vertical polarization brightness temperature and the horizontal polarization brightness temperature respectively.

[0026] (2) In the polarization ratio scatter plot, equally divide the abscissa into several intervals, calculate the average value and standard deviation of the ordinate in each interval. Then, use twice the difference between the average value and the standard deviation in each interval to plot the least squares best fit curve, and this curve is approximately a quadratic equation, as shown in (2).

[0027] (2)

[0028] Among them, is the polarization ratio of the 89 GHz data, is the polarization ratio of the 36 GHz data, and a, b, c are constants.

[0029] Further, since the Unet network can integrate the features of different network layers through skip connections, improving the denoising performance, and the Unet network has strong self-adaptability and can effectively retain the structural information of the image, so the Unet is used as the generation network G.

[0030] Among them, the Unet includes an input layer, a convolutional layer, a pooling layer, an activation layer, and an output layer.

[0031] Specifically, set the kernel of the pooling layer to 2×2, set the size of the convolutional filter to 3×3, and select ReLU as the activation function.

[0032] The function of the discriminant network D is to discriminate two sets of mapping relationships. D uses a CNN network. After the image pair is input into the discriminant network, feature extraction is first performed through consecutive downsampling layers. The convolutional kernel size in the downsampling layers is 4×4, and the stride is 2. To prevent the occurrence of gradient disappearance when updating parameters using the gradient descent method, normalization processing is performed on the convolutional operations of the downsampling layers, and the ReLU activation function is used for non-linear mapping. Finally, the Sigmoid function is used in the discriminant layer to output the value of each pixel, and the cross-entropy is used to calculate the final loss.

[0033] The core idea of the CGAN model: Achieve Nash equilibrium through the game between the generator network and the discriminant network, that is, both the discriminant network and the generator network achieve good results. The purpose of the generator network is to generate data that is close to being unaffected by the external environment, improve the generation ability and reduce the discrimination ability of the discriminant network.

[0034] The beneficial effects of the above technical solutions of the present invention are as follows:

[0035] The present invention proposes a method for correcting remote sensing image data based on an improved Conditional Generative Adversarial Network (CGAN). Through the adversarial training of the generator network and the discriminant network of CGAN, and using the stable relationship between the highly reliable 36GHz data and the 89GHz data that is less affected or not affected by external environments such as clouds and water vapor, the 89GHz data with large interference is corrected. This method can greatly reduce the error caused by the atmosphere by effectively correcting the mixed pixel sea ice data.

[0036] On this basis, the ASI sea ice concentration algorithm is used to invert the sea ice concentration. This method makes full use of the highly reliable 36GHz brightness temperature information, which can make up for the confusion of the brightness temperature information of sea ice and sea water caused during the observation process, improve the inversion accuracy of SIC, and at the same time achieve high spatial resolution.

[0037] The results show that compared with the results of the ASI sea ice concentration inversion algorithm, the sea ice concentration inversion method based on CGAN is feasible. Compared with the sea ice concentration obtained from Landsat data, the sea ice concentration inversion method based on CGAN significantly improves the inversion accuracy of sea ice concentration. Brief Description of the Drawings

[0038] Figure 1 It is a flowchart for correcting the 89GHz affected data based on the improved CGAN model of the present invention;

[0039] Figure 2 It is the polarization ratio scatter plot in the embodiment of the present invention;

[0040] Figure 3 The sea ice concentration results obtained in the embodiments of the present invention;

[0041] Figure 4 The sea ice concentration results obtained in Comparative Example 1 of the present invention;

[0042] Figure 5 In the embodiments of the present invention Figure 3 The sea ice concentration results of the selected area;

[0043] Figure 6 In Comparative Example 1 of the present invention Figure 4 The sea ice concentration results of the selected area;

[0044] Figure 7 The sea ice distribution results obtained from Landsat 8 OLI data using the reflectivity threshold method in Comparative Example 3 of the present invention. Detailed implementation manners

[0045] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, 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 described embodiments of the present invention fall within the scope of protection of the present invention.

[0046] As Figures 1-3 shown in FIGS. 4 and 5: This embodiment provides a high-precision ASI sea ice concentration inversion method based on CGAN for data correction, including the following steps:

[0047] S1. Data screening: Find a relatively stable relationship between the 89 GHz data and the 36 GHz data that are not affected by the external environment or are slightly affected, and screen out the 89 GHz data with large interference;

[0048] Among them, the method for screening the 89 GHz data with large interference is specifically

[0049] (1) Under clear weather conditions, draw a polarization ratio scatter plot with the polarization ratio of the 36 GHz data as the abscissa and the polarization ratio of the 89 GHz data as the ordinate;

[0050] The polarization ratio formula is as shown in (1),

[0051] (1)

[0052] Among them, , are the vertical polarization brightness temperature and the horizontal polarization brightness temperature respectively;

[0053] (2) In the polarization ratio scatter plot, divide the abscissa into several intervals, calculate the mean and standard deviation of the ordinate in each interval, and then use twice the difference between the mean and the standard deviation in each interval to draw the least squares best fit curve, which is approximated by a quadratic equation as shown in (2).

[0054] (2)

[0055] Among them, is the polarization ratio of the 89 GHz data, is the polarization ratio of the 36 GHz data, and a, b, and c are constants.

[0056] The polarization ratio scatter plot is as Figure 2 shown, specifically obtained, and a = 3.5504, b = 0.1876, c = 0.0061.

[0057] S2. Data correction: The generator network and discriminator network of CGAN are trained adversarially. Using the stable relationship between the highly reliable 36 GHz data and the 89 GHz data that is less affected or not affected by the external environment, the 89 GHz data with large interference is corrected.

[0058] The model function of CGAN is as shown in formula (3):

[0059] (3)

[0060] Among them, is the data affected by the external environment, is the additional information, is the input random noise, is the data not affected by the external environment output after inputting the random noise, additional information, and the data affected by the external environment into the generator network of CGAN, is the probability that the discriminator network determines that the input data is fake; since the goal of the generator network of CGAN is to make the generated data as close as possible to the data not affected by the external environment, the loss function is set to to ensure that the probability of the discriminator network outputting a fake image is as small as possible. And the goal of the discriminator network is to improve the ability to distinguish the differences in the input data. Therefore, the larger the better, and at the same time, it is hoped that the influence of the noise is as small as possible. Then the loss function is set to , and is used to represent this game process;

[0061] The training process of the CGAN model is as follows:

[0062] (1) Before training, rotate and translate the dataset to increase the number of datasets, and then normalize the training set and test set.

[0063] (2) Input the training set into the generation network, and then perform consecutive operations of batch normalization + convolution + ReLU + pooling to complete the downsampling operation;

[0064] (3) Perform consecutive operations of deconvolution, BN, and activation function on the feature map obtained by downsampling to complete the upsampling;

[0065] (4) In the same network layer, connect the output feature maps of downsampling and upsampling; in different network layers, from the top neural network layer to the bottom neural network layer, fuse the output feature maps of downsampling and the upsampling of the next neural network layer. The fused output feature map continues to be connected by the upsampling feature map of the next neural network layer, and this iteration continues until there is no corresponding upsampling in the next layer, and then obtain the mapping relationship between 89 GHz data and highly reliable 36 GHz data;

[0066] (5) Then input the test set data, that is, the relationship between 89 GHz data and 36 GHz data without interference, and the mapping relationship between 89 GHz data and highly reliable 36 GHz data obtained in step (4) into the discriminant network;

[0067] (6) Then perform consecutive operations of batch normalization + convolution + ReLU + pooling to complete the downsampling operation;

[0068] (7) Finally, use cross-entropy to discriminate the result obtained by the discriminant network. If the loss function reaches the minimum value, output the corrected 89 GHz data, otherwise return to step (2) and repeat the above steps until the loss function reaches the minimum;

[0069] S3. Based on the corrected 89 GHz data, adopt the ASI sea ice concentration algorithm to retrieve the sea ice concentration.

[0070] Among them, since the Unet network can integrate the features of different network layers through skip connections, improving the denoising performance, and the Unet network has strong self-adaptability and can effectively retain the structural information of the image, so the Unet is adopted as the generation network G;

[0071] Among them, the Unet includes an input layer, a convolutional layer, a pooling layer, an activation layer, and an output layer;

[0072] Specifically, set the kernel of the pooling layer to 2×2, set the size of the convolutional filter to 3×3, and select ReLU as the activation function;

[0073] The function of the discriminant network D is to discriminate two sets of mapping relationships. D uses a CNN network. After the image pair is input into the discriminant network, feature extraction is first performed through consecutive downsampling layers. The convolution kernels in the downsampling layers are all 4×4 in size, and the strides are all 2. To prevent the occurrence of gradient vanishing when updating parameters using the gradient descent method, the convolution operations in the downsampling layers are normalized, and the ReLU activation function is used for non-linear mapping. Finally, the Sigmoid function is used in the discriminant layer to output the value of each pixel, and the cross-entropy is used to calculate the final loss.

[0074] Specifically, the Antarctic sea ice concentration was retrieved as 98 by the ASI algorithm using the AMSR-2 brightness temperature data corrected on February 1, 2021.

[0075] As Figure 3 shown are the sea ice concentration results obtained by the method of the present invention.

[0076] Figure 5 For Figure 3 the sea ice concentration results obtained by the method of the present invention in the selected area in

[0077] Comparative Example 1

[0078] The Antarctic sea ice concentration was retrieved using the uncorrected AMSR-2 brightness temperature data on February 1, 2021 by the ASI algorithm, i.e., ASI algorithm + weather filter.

[0079] As Figure 4 shown are the sea ice concentration results obtained by the method of this comparative example.

[0080] Figure 6 For Figure 4 the sea ice concentration results obtained by the method of this comparative example in the selected area in

[0081] Comparative Example 2

[0082] The Antarctic sea ice concentration was retrieved using the uncorrected AMSR-2 brightness temperature data on February 1, 2021 by the ASI algorithm. The Antarctic sea ice concentration retrieved by the ASI algorithm was 88.

[0083] Comparative Example 3:

[0084] The results were verified using the high-resolution optical remote sensing data Landsat8 OLI.

[0085] Based on Landsat 8 data (resolution: 30 m), the area between 160° - 175° E and 74° - 79° S near the Ross Sea was selected for further verification. Based on Landsat 8 data, we used the NSDI (Normalized Difference Snow Index) calculated by the green band and shortwave infrared band threshold method to detect sea ice distribution. This method can identify ice and water according to the reflectance difference between sea ice and sea water in the red and near-infrared regions (Perovich, 1996; Riggs et al, 1999; Riggs and Hall, 2015; Hallet al, 2001; Liu et al, 2016). Then, the proportion of the number of sea ice pixels in the corresponding AMSR-2 pixel grid was counted, and this proportion was used as the sea ice concentration result of Landsat 8 OLI, specifically 99.

[0086] As Figure 7 shown is the sea ice distribution result obtained from Landsat 8 OLI data using the reflectance threshold method. Among them, Figure 7 in the figure, the white area is sea ice, the black area is sea water, and the gray area is land.

[0087] Verification analysis:

[0088] Combining the examples and the results of Comparative Examples 1-2 shows that the sea ice concentration results inverted by the ASI algorithm based on the corrected data are relatively close to the sea ice concentration results obtained by the ASI algorithm plus the weather filter. To further verify the sea ice concentration results inverted by the ASI algorithm based on the corrected data, we used the sea ice concentration results obtained from high-resolution Landsat 8 data (resolution 30 m) for verification. The sea ice concentration calculated using the ASI algorithm based on the corrected data is 98, which is closer to the sea ice concentration of 99 obtained from the high-resolution data than the sea ice concentration of 88 calculated using the ASI algorithm.

[0089] The ASI algorithm based on the corrected data significantly changes the sea ice concentration of mixed pixels, thereby reducing the impact of weather on high-frequency data. Compared with the sea ice concentration obtained by the ASI algorithm plus the weather filter, the sea ice concentration results inverted by the ASI algorithm based on the GAN-corrected data have higher accuracy.

[0090] The method for correcting sea ice data can also be applied to other data sources, providing new method support for the correction of sea ice data based on microwave radiometer data.

[0091] The above is the preferred implementation manner of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as the protection scope of the present invention.

Claims

1. A high-precision ASI sea ice concentration inversion method for data correction based on CGAN, characterized in that: It includes the following steps: S1. Data screening: Find the relatively stable relationship between the 89 GHz data and the 36 GHz data that is not affected by or has little interference from the external environment, and screen out the 89 GHz data with large interference; S2. Data correction: The generation network and discriminant network of CGAN are trained adversarially. Using the stable relationship between the highly reliable 36 GHz data and the 89 GHz data that is not affected by or has little interference from the external environment, correct the 89 GHz data with large interference; The model function of CGAN is shown in formula (3): (3) Among them, is the data affected by the external environment, is the additional information, is the input random noise, is the data not affected by the external environment output after inputting the random noise, additional information, and the data affected by the external environment into the generation network of the CGAN, is the probability that the discriminator network determines that the input data is fake; since the goal of the generation network of the CGAN is to make the generated data as close as possible to the data not affected by the external environment, the loss function is set to to ensure that the probability of the discriminator network outputting a fake image is as small as possible, and the goal of the discriminator network is to improve the ability to judge the difference of the input data. Therefore, the larger the better, and at the same time, it is hoped that the influence of the noise is as small as possible, so the loss function is set to using to represent the process of this game; The training process of the CGAN model is as follows: (1) Before training, rotate and translate the dataset to increase the number of datasets, and then normalize the training set and test set; (2) Input the training set into the generation network, and then perform continuous batch normalization + convolution + ReLU + pooling operations to complete the downsampling operation; (3) Perform continuous operations of deconvolution, BN, and activation function on the feature map obtained by downsampling to complete the upsampling; (4) In the same network layer, connect the output feature maps of downsampling and upsampling; in different network layers, from the top neural network layer to the bottom neural network layer, fuse the output feature maps of downsampling and the upsampling of the next neural network layer. The fused output feature map continues to be connected by the upsampling feature map of the next neural network layer, and this iteration continues until there is no corresponding upsampling in the next layer, and then obtain the mapping relationship between the 89 GHz data and the highly reliable 36 GHz data; (5) Then input the test set data, that is, the relationship between the 89 GHz data and 36 GHz data without interference, and the mapping relationship between the 89 GHz data and the highly reliable 36 GHz data obtained in step (4) into the discriminant network; (6) Then perform continuous batch normalization + convolution + ReLU + pooling operations to complete the downsampling operation; (7) Finally, use cross-entropy to discriminate the results obtained by the discriminant network. If the loss function reaches the minimum value, output the corrected 89 GHz data, otherwise return to step (2) and repeat the above steps until the loss function reaches the minimum; S3. Based on the corrected 89 GHz data, use the ASI sea ice concentration algorithm to invert the sea ice concentration.

2. The high-precision ASI sea ice concentration inversion method for data correction based on CGAN according to claim 1, wherein: Since the Unet network can integrate the features of different network layers through skip connections, improving the denoising performance, and the Unet network has strong self-adaptability and can effectively retain the structural information of the image, so Unet is used as the generation network G; Among them, Unet includes an input layer, a convolutional layer, a pooling layer, an activation layer, and an output layer; Specifically, set the kernel of the pooling layer to 2×2, set the size of the convolutional filter to 3×3, and select ReLU as the activation function; The function of the discriminative network D is to discriminate two sets of mapping relationships. D uses a CNN network. After the image pair is input into the discriminative network, feature extraction is first performed through consecutive downsampling layers. The convolutional kernels in the downsampling layers are all 4×4 in size, and the strides are all 2. To prevent the vanishing gradient problem when updating parameters using the gradient descent method, the convolutional operations in the downsampling layers are normalized, and the ReLU activation function is used for non-linear mapping. Finally, the Sigmoid function is used in the discriminative layer to output the value of each pixel, and the cross-entropy is used to calculate the final loss.

Images