Ice Concentration Estimation Method with Satellite Based Microwave Radiometer by Means of Inversion Theory

Ice concentration estimation method with satellitebased microwave radiometer by means of inversion theory is proposed. Through experiments, it is found that the proposed methods are superior to the existing methods, the NASA Team algorithm and the Comiso's Bootstrap algorithm with up to 45% of improvement on ice concentration estimation accuracy based on the simulation study. Also 1.5 to 2.1% of improvement was achieved for the proposed method compared to the NASA Team and Comiso's Bootstrap algorithms for the actual The Special Sensor Microwave Imager (SSM/I) data of Okhotsk using Japanese Earth Resources Satellite: JERS-1/Synthetic Aperture Radar: SAR data as a truth data for estimating ice concentration. Keywords—Ice concentration; Microwave radiometer; Inversion Theory; Comiso’s Bootstrap algorithm; The Special Sensor Microwave Imager (SSM/I); Japanese Earth Resources Satellite: JERS-1; Synthetic Aperture Radar: SAR


I. INTRODUCTION
Observation of sea ice is important in considering the global environment. This is because the sea ice areas on both poles of the earth cover 10% of the sea surface and not only have a great influence on the heat balance of the earth, the movement of the ocean and the atmosphere, but are also the places where the effects of global warming phenomena are likely to appear [1]. In addition to scientific observation purposes, it also has the purpose of preventing marine accidents, such as ensuring navigation and sea work safety in polar regions and ice floes.
Global maps of sea ice concentration, age and surface temperature derived from NIMBUS-7 satellite onboard Special Sensor Microwave Radiometer: SSMR: A case study is conducted and well reported together with importance of ice concentration estimation for global warming [2]. A microwave technique for mapping thin ice is investigated [3].
Special Sensor Microwave Imager: SSM/I concentrations using the bootstrap algorithm as the NASA standard product is well documented [4]. Also, temperature corrected bootstrap algorithm is proposed [5].
Estimating occupancy of in-pixel covering class by solving inverse problem, is conducted and well reported [6]. Area ratio estimation (mixing ratio estimation) by pixel category decomposition is proposed and validated [7]. Category decomposition based on maximum likelihood estimation is proposed [8].
Advanced Microwave Scanning Radiometer: AMSR is well reported in terms of requirements and preliminary design study [9]. Method for proportion estimation of mixed pixels by means of inversion problem solving is proposed for mixing ratio estimation [10].
On the other hand, inversion techniques for proportion estimation of Mixed Pixel: Mixels in high spatial resolution of satellite image is proposed [11]. Inversion for emissivitytemperature separation with Advanced Spaceborne based Sensor for Thermal Emission and Radiation: ASTER data is proposed [12]. Meanwhile, method for ice concentration estimation with microwave scanning radiometer data by means of inversion is proposed [13].
Estimation accuracy of ice concentration is not good enough for the global warming problem-solving. Therefore, strong demands of improvement of ice concentration estimation accuracy are raised among the global change research community. In this paper, a microwave radiometer installed on an artificial satellite is used to observe sea ice, and the method of category decomposition is used based on the brightness temperature data of multiple frequencies observed by the microwave radiometer. The author proposes a method to estimate ice concentration.
In the following section, related research works are described. Then, the proposed method is described followed by experimental set-up together with experimental results. After that, concluding remarks and some discussions are described.

II. RELATED RESEARCH WORKS
A method for ice concentration estimation with microwave radiometer data by means of inversion techniques is proposed [14]. An inversion for emissivity-temperature separation with ASTER data is proposed [15]. Spatial resolution enhancement by means of inversion is proposed [16]. Inversion techniques for proportion estimation of Mixels in high spatial resolution of satellite image analysis is also proposed [17].
Ice concentration estimation based on local inversion is proposed [18]. Application of inversion theory for image 88 | P a g e www.ijacsa.thesai.org analysis and classification is investigated [19]. Sea Surface Temperature: SST estimation method with linearized inversion of Radiative Transfer Equation: RTE code for Advanced Earth Observing Satellite: ADEOS / Ocean Color and Temperature Scanner: OCTS is proposed [20].
Estimation of SST, wind speed and water vapor with microwave radiometer data based on simulated annealing is proposed as one of the microwave radiometer data applications [21]. Nonlinear optimization-based SST estimation methods with remote sensing satellite-based Microwave Scanning Radiometer: MSR data is proposed. [22].
Simultaneous estimation of geophysical parameters with microwave radiometer data based on accelerated Simulated Annealing: SA is proposed [23]. Meanwhile, sensitivity analysis for water vapor profile estimation with infrared sounder data based on inversion is proposed [24].
Data fusion between microwave and thermal infrared radiometer data and its application to skin sea surface temperature, wind speed and salinity retrievals is proposed [25]. Comparative study of optimization methods for estimation of SST and ocean wind with microwave radiometer data is proposed [26].

A. Traditional Method for Estimation of Sea Ice Concentration
The estimation of sea ice concentration using a microwave radiometer has been actively performed since 1972, when the satellite NIMBUS-5 equipped with the microwave radiometer Electrically Scanning Microwave Radiometer: ESMR was launched. Gloersen where, C is the sea ice concentration, Tb is the brightness temperature observed by the sensor, Ts is the physical temperature of the 1-year ice, ε is the emissivity of the ice, and is 0.92 in the case of the 1-year ice in the nadir. The constant 135 is the sum of the brightness temperature of the open water surface (120K) and the atmospheric radiation (15K).
U.S. Navy uses the same NIMBUS-5 / ESMR, and as shown in Fig. 1, the brightness temperature when sea ice is 100% is 240K, and the brightness temperature when seawater is 100% is 135K, and the observed brightness temperature during that period is 135K.
The sea ice concentration was expressed as a linear relationship with the brightness temperature, and the sea ice concentration map was created constantly. In Japan, MOS-1 / MSR is often used to estimate the sea ice concentration of one-year ice. For example, a formula using a band of 31 GHz has been proposed as follows. Here, is the output (digital number) D of the 31 GHz channel. C = 4.17D-220.83 (2)

B. Comiso's Bootstrap Algorithm
Furthermore, recently, NASA Team algorithm and Comiso's Bootstrap algorithm have been proposed to improve the estimation accuracy [3], [4]. In the NASA Team algorithm, a linear combination function of the sum / difference ratio (PR) of the vertical and horizontal polarization channels at 19 GHz and the sum / difference ratio (GR) of the vertically polarized waves at 37 GHz and 19 GHz is used by regression analysis beforehand. It is estimated using the obtained coefficients. At this time, it is considered that it is possible to estimate the composition ratio separately for one-year ice and perennial ice, and it is also called a weather filter to remove cloud and water vapor. It is also being done.
f=(C 1 +C 2 PR+C 3 GR+C 4 PR*GR)/D C m =(C 9 +C 10 PR+C 11 GR+C 12 PR*GR)/D (6) D=C 5 +C 6 PR+C 7 GR+C 8 PR*GR (7) where, the Weather Filter means that if the GR calculated from 19 and 37 GHz is 0.05 and the GR calculated from 19 and 22 GHz is 0.045 or more, the sea ice concentration is 0.
Here, C f and C m are the composition ratios of one-year and perennial ice. In Comiso's Bootstrap algorithm, training samples are pre-instructed by a human in advance, and they are referred to obtain the sea ice concentration by the following formula.
where, IC is the sea ice density, T bx is the brightness temperature of different frequency channels (19 and 37 GHz are often used), and T bx w is the open surface and ice brightness temperature, respectively. Represents temperature. Furthermore, in the Comiso's Bootstrap algorithm, a sea surface mask algorithm is also proposed, and correction of sea ice temperature is also considered [5].
Thus, several formulas for estimating sea ice concentration have been proposed, but it is said that the antenna brightness temperature of a microwave radiometer is expressed as a 89 | P a g e www.ijacsa.thesai.org linear combination by the mixing ratio of the brightness temperature of sea ice and sea water. It is based on assumptions. However, as can be seen from the twodimensional scatter diagram of vertically polarized waves of 19 and 37 GHz of SSM / I in the Arctic Ocean region (image shown in Fig. 3) on January 1, 1989, as shown in Fig. 2.
The distributions of 1-year ice and perennial ice have a large overlap and show a large dispersion, and further, they also largely overlap with those of cloud or water vapor-rich regions. Therefore, it is difficult to improve the estimation accuracy of sea ice concentration. In addition, the antenna brightness temperature of the microwave radiometer changes depending on other factors such as water vapor and cloud water content. Therefore, an algorithm that considers those influences is necessary. In this research, we tried to use the method of category decomposition for sea ice concentration estimation. In doing so, the effects of cloud and water vapor are added by adding them to the category.

A. Linear Model of Pixel Spectral Vector
The author has already proposed a method of category decomposition based on the solution of the inverse problem6). According to this, when k types of categories are included in the pixels of the multiple spectral image, if the area ratio of each category is a j , (j = 1, ..., k), the spectral vector P is It is represented by the following equation by a typical spectral vector M j .
The area ratio can be estimated by the following equation, where n is the number of observation channels and k is the number of decomposed categories.
P is a vector that aligns the brightness temperatures observed at different frequencies, and M j is a vector that consists of representative brightness temperatures of categories such as sea ice and open water surface in each observation channel. A is a vector showing the area ratio of the category. Hereinafter, P, A, M j , and M will be referred to as an observation vector, an area ratio vector, an average vector, and an average matrix, respectively.
It is already known that M can be estimated by extracting training samples from the image, and P is an observed value. Based on this basic model, the area ratio of each category (one of them is sea ice concentration) is obtained by using the inverse problem-solving method. This method includes the Moore-Penrose type generalized inverse matrix, the observed least squares method, the area ratio least squares method, the maximum likelihood grid search method, etc., as shown in [6].

B. Moore-Penrose Type Generalized Inverse Matrix
In the above linear model, if the mean matrix M is regular, there is an inverse matrix: Therefore, the area ratio can be easily calculated. However, since M is rarely regular, the Moore-Penrose generalized inverse matrix shown below is needed.
However, since this method does not consider the measurement error of the data, the estimation error may be large, or the total of the mixing ratios may not be 1. Therefore, the method using the least squares method with the following constraint conditions is used.

C. Estimation using the Least Squares Method based on Observed Value
If the observed value contains an error, the accuracy is improved by minimizing the error between the observed vector P and the estimated value P '= MA of the true value. The restraint conditions at that time are as follows: If the non-negative condition a j =0 is removed, it can be analytically solved by using Lagrange's undetermined multiplier method as the constrained least squares method. Therefore, let λ be an undetermined multiplier and consider the following function F.
A can be obtained by solving the following simultaneous equations. Eq. (17) Substituting this into Eq. (18) gives and λ is determined.
Therefore, from Eq. (19), the estimated value of the area ratio vector A is obtained as.
By the way, in the actual measurement data, the matrix M is not always equal to the true representative value vector M o . Let this difference be E. If the true area ratio vector A 0 and the observation vector does not include an error, the following equation holds.
Therefore, when M and P are given by Eq. (22), the area ratio vector A 1 obtained by Eq. (21) is which is an estimated error of the area ratio vector can be estimated by the following Eq. (24).
where α is the angle between u and (M t M) -1 u, and β is the angle between u and M + EA 0 . On the other hand, the area ratio vector A' calculated algebraically for the same data and its estimation error are as follows: and from Eq. (24), in some cases a much larger error can occur than in the algebraic solution.
That is, the estimation accuracy depends on the accuracy of the matrix M composed of the representative values of each category.

D. Estimation by Area Ratio Least Squares Method [7]
In the least-squares method of observed values, the area ratio was obtained because the problem of the representativeness of M can be expressed as the error of the observation vector. However, when using the actual measurement data, the representativeness of M is more important than P, and it has a great influence on the estimation accuracy of the area ratio. Therefore, the estimated value N of the generalized inverse matrix M 0 + of the true representative value matrix and the area ratio vector A are set as unvalued, and A that minimizes the residual difference between M + and N is obtained from the following formula.
The area ratio vector obtained by this method agrees with the one obtained by the following least squares method with constraints.

|A-M
where, if M and P are given in the same way as in the case of the observed least squares method, the estimated area ratio vector A 2 and the estimation error | A 2 -A 0 | are given by (37) and (38).
The estimation accuracy of the area ratio least squares method is as good as the second term on the right side of (38), and the problem that the estimation accuracy is extremely poor depending on the property of the matrix M as in the conventional method is solved. It seems to be that. u t A = 1 is an equation of a plane in k-dimensional vector space, where A is a variable, and its normal vector is u. M + P is also a point in the k-dimensional vector space. This is to find the point on the plane u t A 1 with the minimum distance from the point M + P. In other words, it is enough to find the foot of the perpendicular line from the point M + P to the plane u t A = 1. Since the perpendicular is parallel to the vector u, the equation of the perpendicular can be written as v with λ as a parameter.

A=M + P+
(39) Eq. (36) is obtained by finding the intersection of this and the plane.
then eliminating A and finding λ.
E. Maximum Likelihood Grid Search Method [8] As in the case of the least square method, the observation value vector is P, the area ratio vector is A, and the representative value matrix is M. At that time, it is assumed that the observation value is given in a form in which the observation error ε is added to the linear combination of M and A.
Then, assume that the element m ij protection of M follows the normal distribution: ( * , 2 ) with mean m ij and variance σ ij 2 , and ε i follows ( , 2 ) , and consider the observation value vector P as a random variable. At that time, the observed value of the i-th band pi follows the normal distribution ( * , 2 ) of mean mi and variance σ ij 2 expressed by the following equations (42) and (43). However, the representative values of each category are assumed to be independent of each other.
where diag represents a diagonal matrix. The probability Q (p i ) that the observed value p i of the i-th band is observed is expressed by the following equation.
where the probability that the observed value vector P is observed is Q (P), it is expressed by the following equation.
The calculated area ratio A is the value when Q (P) is at its maximum. Therefore, the area ratio A that minimizes the following equation R (P) must be obtained.

R(P)=-ln{Q(P)}
(46) In practice, the ratio of each category is changed every 1%, R (P) is calculated for all the combinations, and A is determined. Therefore, it is called the maximum likelihood grid search method.

V. EXPERIMENTS
The most important thing in valid intercomparing between methods is correct data. Here, we tried the evaluation by the simulation data that can give the correct answer and the evaluation by the correct answer data created from the synthetic aperture radar image. 92 | P a g e www.ijacsa.thesai.org

A. Simulation Data Used
The SSM / 1 data of 19GHz vertical, horizontal polarization, 22GHz vertical polarization, 37GHz vertical, horizontal polarization of the above-mentioned January 1, 1989 Arctic were used. Simulation data is created based on the training sample data extracted from SSM / I data, and the proposed method is applied to compare the accuracy. The procedure for creating the simulation data was as follows.
1) Create mixture ratio data for multiple categories.
2) Extract training samples for the category to be classified from the SSM / I image.
3) Considering the variance of the training sample in (2), output the data given the error according to the normal distribution according to the mixing ratio in (1) Table I is a training sample of each category extracted from the actual SSM / I image. Table II shows the Root Mean Square: RMS error and CPU Time (Elapsed time) from the correct data.

B. Simulation Result
In this table, LSQ is the least squares and MLH is the maximum likelihood grid search method, respectively. As for the inverse problem-solving method, the method with the constraint condition is more accurate. It is also slightly better than the NASA Team and Comiso's Bootstrap algorithms. Furthermore, the maximum likelihood grid search method was confirmed to improve the accuracy by 45% with respect to the Comiso's Bootstrap algorithm. However, as the number of categories increases, the amount of calculation increases exponentially, so it takes considerably longer than the linear solution method.

C. Experimental Method with Real SAR Data
SAR has a much higher spatial resolution than SSM / I. In addition, in the case of L-band SAR, sea ice generally has a higher backscattering cross section than the sea discharge surface, so it is relatively easy to identify it if effects such as sea surface wind speed are taken into consideration. Furthermore, because of its long wavelength, it is hardly affected by clouds. Utilizing this, sea ice concentration data corresponding to SSM / I images was created from SAR images, and the accuracy of the inverse problem-solving method was evaluated using the data as correct answer data.
First, for each pixel of the SAR image, the categories of sea ice and open water surface are classified using the difference in backscattering cross section. Due to the difference in resolution, there are about 7,000 SAR pixels corresponding to SSM / I pixels. Since the latitude and longitude are known in both images, the sea ice concentration was estimated from the SAR image by matching both images based on them.
JERS-1 / SAR images of the Sea of Okhotsk facing the coastline of northern Hokkaido were used. The date and time are from February 4 to 9, 1994, and the exact location is 143~146 degrees east longitude and 43.5~46 degrees north latitude. One pixel of JERS-1 is 12.5mx 12.5m, but this time we used the data processed to 300mx 300m. Each pixel is data corresponding to −20 dB to 5.5 dB. Here, the important point is the threshold that distinguishes sea ice from open water. Fig.  4 is a histogram of pixel values of the SAR image used this time.
The higher backscattering cross section is sea ice, and the lower backscatter cross section is the open water surface distribution. From this figure, it is appropriate to set it to 60 to 70 (actual value is -14dB to -13dB). In the experiment, -13.5 dB was used as the threshold.

D. Experimental Result with Real SAR Data
Looking at the RMS error of each method, Moore-Penrose generalized inverse matrix (13.63%), observed least squares (12.96%), area ratio least squares (12.80%), NASA algorithm (12.66%), Comiso's Bootstrap algorithm (12.58%) and the maximum likelihood grid search method (12.40%) were confirmed, and it was confirmed that the accuracy was high in almost the same order as the comparative evaluation results by simulation. However, a low estimation accuracy was obtained overall, which is considered to be due to the following reasons.

1) Error in alignment
2) Error of estimated data from SAR image 3) Moving sea ice Due to the large difference in the resolution of the images used this time, 25 km square and 0.3 km square, and the physical characteristics of sea ice concentration, it is possible that a slight error in the registration caused a large error in the estimation results.
The author artificially generated images that were shifted by 0.5 pixels in SMM / I and in the north, south, east, and west. From the result, there is a possibility that the maximum is shifted by 0.5 pixels. If SSM / I pixel shifts by 0.5 pixel, it is 12.5km, which is quite large. In addition, since the sea ice present at this site is drift ice, it is considered that the difference between the observation times of both data is also relevant. Furthermore, as can be seen from the histogram, the distribution of backscatter values of sea ice and the distribution of objects on the sea surface largely overlap. An estimation error occurs from here. The threshold that minimizes the estimation error can be set but cannot be 0. The reason for this distribution is that the surface condition of the drift ice is complicated and has various forms.

VI. CONCLUSION
From the result of the simulation experiment and the experiment using the sea ice concentration estimated from the synthetic aperture radar, it was confirmed that the inverse problem solution is effective for the sea ice concentration. In particular, the least-squares method with the constraint of the square of the estimation error of the observed value and the estimation error of the mixture ratio was found to be effective. However, the accuracy of estimation is slightly inferior to that of the maximum likelihood grid search method because the calculated area ratio does not include a non-negative condition and the variance of the brightness temperature of the training samples in the category is not taken into consideration.
Compared to these, the NASA Team and Comiso's Bootstrap algorithms have almost the same estimation accuracy as the least-squares method that constrains the squares of the estimation error of the observed value and the estimation error of the mixture ratio, and the maximum likelihood grid search method. It was also confirmed that the estimation accuracy was slightly inferior.
The maximum likelihood grid search method was 2.1% and 1.5% higher than that of NASA Team and Comiso's Bootstrap algorithm, respectively, based on the results of experiments using the correct solution of sea ice concentration estimated from the synthetic aperture radar in the Okhotsk sea area. It was confirmed that the estimation accuracy of was improved. The reason is that only the maximum likelihood grid search method considers not only the average value of the brightness temperature of each category but also the variance.

VII. FUTURE RESEARCH WORKS
The proposed method must be tested with the other microwave radiometer data as well as synthetic aperture data. Also, influences of sea ice concentration on global warming based on the estimated concentrations. 94 | P a g e www.ijacsa.thesai.org