A multi-scale second-order autoregressive recursive filter approach for the sea ice concentration analysis

1. Introduction

In the field of ocean data assimilation, the three-dimensional variational method is one of the commonly used methods. In variational methods, the background error covariance matrix plays an important role in the quality of the analyzed field. It affects the propagation of observed information among different variables and within the internal space of variables (Cao et al., 2008). However, this matrix has a huge scale and is often ill-conditioned. Therefore, explicit creation based on selected relevant scales is usually not adopted (Li et al., 2011). More studies focus on using recursive filter methods to approximately simulate the background error covariance matrix (Lorenc, 1992; Hayden and Purser, 1995; Huang, 2000; Purser et al., 2003a, b).

Recursive filter is a continuous approximation method. It modifies the information of the background field through changes in spatial characteristic scales in a series of iterations, thus reducing the regional influence (He et al., 2011). In the variational method, a Gaussian filter can be approximately simulated by iterating a first-order filter several times (Purser et al., 2003a). The characteristic scale of Gaussian filter is a low-pass filter parameter. When the given characteristic scale is small, there are short waves in the analyzed field. If the characteristic scale is large, it may occur excessive smoothing. The determination of filter parameter needs to ensure that the effect of the filter is similar to that of convolution on a field using a covariance function (Vandenberghe and Kuo, 1999; Zhang et al., 2004). Gaussian or second-order auto-regressive (SOAR) functions are generally used to describe the background error horizontal correlation function. Research has shown that the Gaussian correlation function do not have sufficient power at small scales (Vandenberghe and Kuo, 1999; Zeng, 2006). To address this issue, Purser et al. (2003b) replaced the single Gaussian function by superimposing Gaussian functions of different scales. He et al. (2011) fitted three Gaussian correlation functions with different characteristic scales by linear superposition to obtain a fourth-order recursive filter. It not only preserves large-scale information but also displays many meso- and small-scale information. The regional three-dimensional variational (3D-VAR) assimilation based on multi-characteristic scale recursive filters has also been verified to obtain more mesoscale information in meteorological element forecasting applications (Wu et al., 2018; Zhu and Li, 2021). In addition, in recent years, other studies on scale separation or multi-scale information extraction have also achieved successful applications in the marine field, such as multigrid data assimilation method (Li et al., 2008), 3D-VAR based on sequential filter method (Wu et al., 2011), sequential 3D-VAR method (He et al., 2008), scale-selective data assimilation method (Peng et al., 2010), and multi-scale 3D-VAR based on diffusion filter method (Li et al., 2011).

It should be noted that according to Zhang et al. (2004), in one-dimensional recursive filter experiments, the above two background error horizontal correlation functions can approximate the analytic solution of the original function with the increase of filtering times. However, the Gaussian function requires 10 filters, while the SOAR function requires only three filters. Due to the limited selection of fixed filter parameters, the filter results based on the SOAR function still cannot obtain information at all scales. To take the advantage of SOAR function and solve the problem of multi-scale extraction, a reconstruction method for gridding observation based on variational optimization technology, called multi-scale second-order auto-regressive recursive filter (MSRF) scheme, is designed in this study. The MSRF scheme is based on the spatial multi-scale recursive filter (SMRF) method firstly proposed by Chen (2011), then applied in sea ice concentration (SIC) analysis by Zhang et al. (2020). Here MSRF uses a SOAR filter to replace the cascade of multiple first-order recursive filters in the SMRF scheme. The MSRF scheme only needs to perform two times of recursive filter to be equivalent to a SOAR correlation function (Lorenc, 1992), instead of several times (8 times; Zhang et al., 2020) of recursive filter as in the SMRF scheme. At the same time, the spectral response function of both the SOAR correlation function and its Laplace operator has a larger power spectrum at meso- and small-scale (Zhuang et al., 2021), which is beneficial for obtaining meso- and small-scale information in analysis. In addition, the MSRF scheme also gains additional benefits in multi-scale processes: (1) it has convergence; (2) in the minimization process, the linear search algorithm automatically determines the step size and analyzed results are adjusted without manual intervention; (3) observations at all spatial scales can be processed in one iteration.

The paper is organized as follows. In Section 2, the theory of SOAR correlation function is firstly introduced and the SOAR filter concept is proposed. Secondly, based on the SMRF method, the SOAR filter is applied to multi-scale 3D-VAR analysis, and the MSRF scheme is designed. In Section 3, two-dimensional ideal experiments are conducted on the MSRF scheme based on single point observations. Section 4 further applies the MSRF scheme to an example of Arctic SIC, and compares the MSRF scheme with the SMRF scheme in terms of analyzed results and computational efficiency. The conclusions are summarized and discussion is presented in Section 5.

2. Methods

2.1 The SOAR function

Assuming the input vector is ${\boldsymbol{a}}_{M\times 1} $ and the output vector is $ {\boldsymbol{c}}_{M \times 1} $, a one-dimensional recursive filter can usually be represented as follows (Lorenc, 1992; Hayden and Purser, 1995):

where $ \alpha $ is the filter parameter $ (0<\alpha< 1) $, setting $ b_{0}=0 $ and $ c_{M+1}=0 $. By combining the above two equations, it can be concluded that:

Let $ a_{i}={\mathrm{exp}}({{\mathrm{j}} ki \Delta x}) $ and $ c_{i}=\phi(k) {\mathrm{exp}}({{\mathrm{j}} k i \Delta x}) $ (j is the imaginary notation, k represents the wave number and $ \Delta x$ is the grid interval) in the above equation, and the frequency response function is obtained as

When $ k \Delta x $ is infinitesimal, the above equation after N iterations can be approximated as

The determination of filter parameter needs to ensure that the effect of the filter is similar to that of convolution on a field using a covariance function. A SOAR function is considered as the correlation function, which in the one-dimensional case can be expressed as

where $ c> 0 $, $ \sigma_{b}^{2} $ is the variance. Its spectral response is as follows (Zeng, 2006):

Comparing Eq. (5) and Eq. (7), it can be found that, in addition to satisfying the following two equations,

it needs to multiply the final result of the filter by the following equation:

Let (Zhang et al., 2004; Zeng, 2006)

according to Eqs (8), (9), and (11), it can be calculated that:

By using the scale factor of Eq. (10), the filter parameter N can be calculated from Eqs (11) and (12), then the recursive filter process of Eqs (1) and (2) can be performed. Moreover, the result of only performing two recursive filters (satisfying Eq. (8)) is equivalent to the SOAR correlation function of Eq. (6) (Lorenc, 1992), which is named the SOAR filter in this study.

2.2 The MSRF scheme

Observation data provide information on various wavelength scales of ocean (or sea ice) elements. Traditional 3D-VAR methods only perform a single 3D-VAR analysis, which is difficult to extract multiple wavelengths simultaneously. The analyzed results may be contaminated by noise such as observation errors or irregular data distribution. Compared to the multi-grid 3D-VAR method (Li et al., 2008; Xie et al., 2011), which separately solves 3D-VAR analysis on grids of different resolutions, the SMRF method can process all spatial scale observations in once analysis process (Zhang et al., 2020).

The main idea of the SMRF method is to minimize the difference between the estimated field and the observed field through variational optimization techniques. The unconstrained minimization problem below is considered:

where x is the analyzed field, x^o is the observed field. H is the observation operator from the analysis space to the observation space. R is the covariance matrix of observation errors. To suppress observed noise, a recursive filter operator D with a filter parameter $ \eta$ (using a smaller value) is considered to be used. Let $ {\boldsymbol{x}}=\mathrm{D} {\boldsymbol{w}} $ (it represents the filter operator is applied to the raw data) and Eq. (13) is replaced with:

The gradient of the cost function $ {\mathrm{J}}(w) $in the above equation is as

For the unconstrained minimization problem of Eq. (14), minimization algorithms can be used to solve it (such as steepest descent method, L-BFGS method and conjugate gradient method). The specific calculation steps are as follows.

(1) Given the initial guess value $ {\boldsymbol{w}}={\boldsymbol{w}}_{0} $ and the constant value of the filter parameter $ \eta $ in the [0, 1] interval;

(2) Perform recursive filter on $ {\boldsymbol{w}} $ using filter operator $ \text{D}$ with filter parameter $ \eta $, and calculate $\mathrm{D} {\boldsymbol{w}}$;

(3) Calculate $\mathrm{J}({\boldsymbol{w}})$ and $ \nabla \mathrm{J}({\boldsymbol{w}}) $;

(4) Use minimization algorithms to optimize $ {\boldsymbol{w}} $;

(5) Go to Step (3) until the convergence condition is met;

(6) The final analyzed value is $ {\boldsymbol{x}}_{{\rm{a}}}=\mathrm{D} {\boldsymbol{w}} $.

As mentioned by Zhang et al. (2020), when observations are scarce or irregularly distributed, the minimization problem solved according to the above steps is usually ill posed. In this case, the gradient values of the cost function are spatially discontinuous. If the “flawed” gradient is introduced into a general gradient-based minimization algorithm, the updated estimates will contain erroneous small-scale signals. In the SMRF scheme, in addition to applying recursive filters to the variable $ {w} $, the recursive filter process is also applied to the gradient of the cost function. Moreover, in each iteration of the minimization process, the filter parameter also changes to extract wavelengths of different scales. Considering that the SOAR filter has certain advantages in capturing meso- and small-scale signals and computational efficiency, the SOAR filter is applied to both the variable $ {w} $ and the gradient of the cost function in the MSRF scheme. The basic steps of the MSRF scheme are briefly described as follows.

(1) Define the cost function and give the initial guess value $ {\boldsymbol{w}}={\boldsymbol{w}}_{{{0}}} $.

(2) Filter $ {w} $ using a recursive filter operator D with a fixed small scale factor $ \beta' $, and calculate $ {\mathrm{D}} {\boldsymbol{w}} $.

(3) Calculate observational residuals $ {\boldsymbol{x}}^{\mathrm {o}}-\mathrm{HD} {\boldsymbol{w}} $ and the gradient $\nabla \mathrm{J}({\boldsymbol{w}})$.

(4) Given a larger scale factor $ \beta $, calculate another filter parameter based on Eqs (10)–(12). Apply it to the negative direction of the gradient, yielding $\mathrm{E}(-\nabla \mathrm{J}({\boldsymbol{w}}))$ ($ \mathrm{E} $ is another recursive filter operator). The result is used as the descent direction.

(5) Using the line search algorithm (Moré and Thuente, 1994) to determine the descent step size and adjust ${\boldsymbol{w}}$.

(6) Decrease the scale factor $ \beta $ according to the given reduction formula.

(7) Repeat Steps (2)–(6) until the convergence condition is met.

(8) The final analyzed value $ {\boldsymbol{x}}_{{\bf{a}}}=\mathrm{D} {\boldsymbol{w}} $.

The basic principle of the MSRF scheme is analyzed below, and the corresponding algorithm flow chart is shown in Fig. 1. Similar to the SMRF scheme, the initial guess value $ {\boldsymbol{w}}_{0} $ is firstly given to the independent variable $ {w} $ of the cost function and a smaller scale factor $ \beta' $ is also given ($ \beta' $ has the same meaning as $ \beta $ of Eq. (10), but $ \beta' $ is a fixed value). The recursive filter operator D with parameter $ \beta' $ is applied to the initial estimation. Then the gradient $ \nabla \mathrm{J}\left({\boldsymbol{w}}_{0}\right) $ can be calculated. All scales of information can pass through the filter because of the smaller scale factor, meaning that the gradient $ \nabla \mathrm{J}\left({\boldsymbol{w}}_{0}\right) $ retains almost all information (only filtering out noise). Then, a large enough scale factor $ \beta $ is given and the SOAR filter (Eqs (10)–(12)) is applied to the negative gradient direction $ -\nabla {\rm{J}}\left({\boldsymbol{w}}_{0}\right) $, yielding ${{\rm{E}}}\left(-{\rm{J}}\left({\boldsymbol{w}}_{0}\right)\right) $. Since the SOAR filter is used here, instead of filtering multiple times based on basic first-order recursive filter used by the SMRF scheme, it only needs to perform the filtering process twice. Because the scale factor $ \beta $ is directly proportional to the filter parameters in Eq. (12), a sufficiently large filter parameter has been selected at the initial iteration. Therefore, $ \mathrm{E}\left(-\nabla \mathrm{J}\left({\boldsymbol{w}}_{0}\right)\right) $ represent the “longest” wavelength information of the observational residuals at the initial guess. To obtain the information of amplitude loss in the filtering process, $ \mathrm{E}\left(-\mathrm{J}\left({\boldsymbol{w}}_{0}\right)\right) $ is used as the search direction. The positive definite operator $ \mathrm{E} $ ensures that this is a descending direction. The line search algorithm is used to find the appropriate step size l along this direction, and $ {w} $ is adjusted. The new estimated value $ {\boldsymbol{w}}_{1}={\boldsymbol{w}}_{0}+l \times\left(-\mathrm{J}\left({\boldsymbol{w}}_{0}\right)\right) $. In fact, at this time, the “maximum” scale of the observational residual at the initial guess has been “completely” extracted. The scale factor $ \beta $ is decreased according to the given reduction formula so that the “maximum” scale of the current observational residual can be extracted in the second iteration process. It is merged with $ {\boldsymbol{w}}_{1} $ to generate a new estimated value $ {\boldsymbol{w}}_{2} $. As $ \beta $ gradually decreases, the information of each scale from long wavelength to short wavelength is extracted successively. Finally, the iteration is stopped when the filter parameter $ \alpha $ calculated based on $ \beta $ and Eqs (10)–(12) is less than 0. The final analyzed field is obtained.

Figure  1.  Flow chart of the MSRF scheme.

DownLoad: Full-Size Img PowerPoint

3. Single-observation ideal experiment

To investigate the transmission effect of SOAR filter on observation information in multi-scale 3D-VAR analysis, two-dimensional single point experiments are carried out. A recursive filter operator with varying scale factor $ \beta $ is selected to filter $ -\nabla \mathrm{J}({\boldsymbol{w}}) $, and the Eq. (14) is solved through a step-by-step iterative. In addition, for the convenience of comparison, the filter operator with different scale factor $ \beta $ are applied to ${\boldsymbol{w}}$ separately, and the Eq. (14) is solved directly using the minimization algorithm Limited Memory-Broyden Fletcher Goldfarb Shanno (L-BFGS) algorithm. Furthermore, the propagation capability of the MSRF and SMRF schemes is compared visually.

3.1 Data and parameters

The analysis area is given as a square area, with a latitude range of 0°–10°N and a longitude range of 0°–10°E. The resolution of the analysis grid is 0.25° × 0.25°, that is, $ \Delta x $ in Eq. (11) is taken as 0.25. A single point observation with a value of 1.0 is placed in the center of the area at (5°N, 5°E). No value is assigned to other positions. The number of filtering passes $ N $ is set to 2. In the MSRF experiment, the varying scale factor $ \beta $ is set as follows:

where $ M $ represents the total number of iterations. To ensure that a sufficient long wavelength is extracted at the initial iteration step, $ \beta_{1} $ is set to 7.3. When i = 1, the filter parameter $ \alpha $ is approximately 0.87 according to Eqs (10) –(12) (this is close to the initial $ \alpha $ setting of 0.9 by Zhang et al. (2020)); when $ i=M $, $ \alpha \rightarrow 0 $. When $ \alpha< 0 $, the iteration stops. The observation operator $ \mathrm{H} $ is obtained by bilinear interpolation, the variance $ \sigma_{b}^{2} $ is set to 1.0 and the initial guess $ \psi_{0} $ is set to 0.

3.2 Results

Figure 2 shows the results of solving Eq. (14) based on the L-BFGS algorithm when the scale factor $ \beta $ of the SOAR filter is set to 0.6 (Fig. 2a), 1.2 (Fig. 2b), 1.8 (Fig. 2c), and 2.6 (Fig. 2d), respectively. It can be seen that different scale factors have different performance on the transmission ability of observation information. In Fig. 2a, when $ \beta $ is set to a smaller value, small-scale information contained in the observation is maximally extracted, and the analyzed result is close to the observation. However, its propagation range is very limited, and in practical applications, the signals may not able to propagate in sparse observations or vacant areas. When $ \beta $ is relatively large (Fig. 2d), the observed signal is obviously propagated over a wider area, but at the cost of sacrificing some necessary shortwave observation information. From the above results, it can be found that the selection of a single scale factor cannot meet all requirements. An expected mechanism for transmitting observation information should both ensure the accuracy of information transmission and maximize the transmission of observation signals. This requires simultaneous selection of different scale factors in one mechanism to extract information at different scales. The MSRF proposed in this study can achieve this to a certain extent.

Figure  2.  Spread of observational information using the L-BFGS scheme when $ \beta $ = 0.6 (a), 1.2 (b), 1.8 (c) and 2.6 (d), respectively.

DownLoad: Full-Size Img PowerPoint

Figure 3 shows the propagation results of observed signals under different iteration steps when the scale factor $ \beta'$ of the MSRF scheme is set to a smaller value of 0.35. The corresponding $ \beta $ values for Figs 3a, b, c, and d are approximately 4.5, 2.5, 1.8, and 1.2, respectively. As $ \beta $ gradually decreases, the signals from large-scale (long wave) signals to small-scale (short wave) are corrected. Compared to the single scale observation signals extracted by different scale factors in Fig. 2, the MSRF scheme extracts signals of all scales in one analysis. This not only ensures that the accuracy of single point observation signals in analyzed results does not reduce, but also spreads to a wider area.

Figure  3.  Spread of observational information in the MSRF scheme when $ \beta' $ = 0.35. a−d are the results at iteration 70, 120, 140, and 160, respectively.

DownLoad: Full-Size Img PowerPoint

To further analyze the signal propagation ability of the MSRF scheme in two-dimensional ideal experiments, the observed signal propagation results of the SMRF scheme in the same analysis domain are visually compared, as shown in Fig. 4. Figure 4a shows the surface diagram corresponding to the propagation results of the MSRF scheme when total number of iterations $ M $ of Eq. (16) is 210. To ensure maximum detail information, the initial fixed value of the filter parameter for the SMRF scheme in Fig. 4b is set to 0.1. The total number of iterations is also selected as 210. It is not difficult to see from Fig. 4 that both of them have successfully extracted information at various scales. The difference between the two is that the SMRF scheme clearly exhibits isotropic characteristics, while the MSRF scheme tends to be anisotropic at the outer edge of the propagation range. At the same time, under the same number of iteration steps, the MSRF scheme has a wider signal propagation range, with 645 grid points in analysis domain where the signal propagation value is greater than or equal to 0.05. This number is only 374 in the SMRF scheme. Furthermore, the iteration completion time of the MSRF scheme is only about 1/3 of that of the SMRF scheme. It can be concluded that the MSRF scheme has a certain degree of improvement over the SMRF scheme in signal propagation range and computational efficiency.

Figure  4.  Surface plot of the MSRF scheme (a) and the SMRF scheme (b). The results at iteration 210.

DownLoad: Full-Size Img PowerPoint

4. Application to SIC experiment

To explore the ability of the MSRF scheme in extracting spatial multi-scale information in the actual ocean (or sea ice) environment, two-dimensional data assimilation experiments are conducted using Arctic SIC observation data. By comparing with the L-BFGS scheme, the principle of the MSRF scheme to extract information in the data-void region is further analyzed. Moreover, the MSRF and SMRF schemes are qualitatively and quantitatively analyzed in terms of assimilation error and computational efficiency.

4.1 Data and parameters

The SIC used in the experiment is derived from the dataset named Neal-Real-Time Defense Meteorological Satellite Program (DMSP) Special Sensor Microwave Imager/Sounder (SSMIS) Daily Polar Gridded Sea Ice Concentrations, Version 2 provided by the National Snow and Ice Data Center (Meier et al., 2021). The dataset is obtained by using the NASA team algorithm (Cavalieri, 2012) to process data received from SSMIS on DMSP satellites. The data range is from November 1, 2021 to present, with a spatial resolution of 25 km × 25 km and a temporal resolution of days. The data covers the entire Arctic.

The “true” field from SSMIS Arctic SIC observation data on August 14, 2023 is shown in Fig. 5a. Since satellite data typically have higher spatial resolution than the model (Yang et al., 2022), one grid point is retained for every four grid points. Moreover, some observations are removed from areas with dramatic changes in SIC to verify the ability of the MSRF scheme in extracting multi-scale information, as shown in the data-void region in Fig. 5b. A total of 1382 observations are remained to reconstruct the “true” field (Fig. 5b). In this experiment, the grid interval $ \Delta x $ in Eq. (11) is taken as 25. Because it is expanded by 100 times compared to the single point observation experiment in Section 3, the formula for the varying scale factor needs to be expanded by 100 times compared to Eq. (16):

Figure  5.  SSMIS Arctic sea ice concentration observation on August 14, 2023 (a) and the location of selected observation points for data assimilation (b).

DownLoad: Full-Size Img PowerPoint

where $ \beta_{1} $ is set to 730. The other parameter settings are the same as in Section 3.

4.2 Results of the L-BFGS scheme

Figure 6 shows the SIC analyzed results obtained by the L-BFGS scheme based on the SOAR filter when the scale factors are set to 110 (Fig. 6a), 200 (Fig. 6b), 450 (Fig. 6c), and 700 (Fig. 6d), respectively. From the comparison with the “true” field (Fig. 5a), it is not difficult to see that the selection of $ \beta $ is not an easy task. When $ \beta $ is selected as a smaller value (Fig. 6a), it shows good ability to capture shortwave information and can reflect more detailed information in observations. However, due to the small propagation range of the observed signal (as shown in Fig. 2a), it is difficult to extract the observed longwave information in the data-void region. In contrast, when $ \beta $ is set to a large value (Fig. 6d), there are no significant data void regions. Nonetheless, compared with the “true” field (Fig. 5a), the analyzed results are too smooth and the detailed information is lost because the long wavelength information is mostly extracted. In addition, the analyzed results of SIC (Figs 7a, c, and e) and the descent direction (Figs 7b, d, and f) corresponding to $ \beta $ = 110 at iteration steps of 3, 5, and 7 are also presented. From the right column of Fig. 7, it can be seen that the gradient varies dramatically near the data-void region due to the irregular distribution of observations, and the descent direction is spatially incoherent. This contributes to the analyzed field (left column of Fig. 7) updated in this direction not being able to achieve data filling in this region.

Figure  6.  Sea ice concentration analyzed field on August 14, 2023 solved by using the L-BFGS scheme when $ \beta $ = 110 (a), $\beta $ = 200 (b), $\beta $ = 450 (c), and $\beta $ = 700 (d).

DownLoad: Full-Size Img PowerPoint

Figure  7.  Sea ice concentration analyzed field (left column) and the descent direction ($ -\nabla \mathrm{J} $) (right column) from the L-BFGS scheme ($ \beta $ = 110) at iteration 3, 5, and 7, respectively.

DownLoad: Full-Size Img PowerPoint

4.3 Results of the MSRF scheme

When the initial fixed value of the scale factor $ \beta' $ is set to 35, the SIC analyzed field on August 14, 2023 constructed by the MSRF scheme is shown in Fig. 8a. Compared with the “true” field (Fig. 5a), it can be found that the construction results basically reflect the variation trend of SIC in different regions. Especially near the Beaufort Sea, three contour lines of different SIC are well reproduced. The fundamental reason for the improvement of the construction results compared to the L-BFGS scheme is that in the MSRF scheme, the gradient is filtered and long wavelength of observational residuals is extracted to establish the descent direction. It smooths the gradient with sharp changes in the data-void region in the L-BFGS scheme (Figs 7b, d and f). Figures 9b, d, and f demonstrate that with the increase of iterations increases, the scale factor $ \beta $ gradually decreases, and the wavelengths from long wave to short wave in the descent direction are sequentially extracted. Therefore, the updated analyzed field in this direction also achieves the sequential acquisition of observed information from long wave to short wave (Figs 9a, c, and e). The final analyzed field contains all wavelength information, effectively reproducing the SIC field.

Figure  8.  Sea ice concentration analyzed field from the MSRF scheme with $ \beta' $ = 35 and $ M $ = 215 (a) and SIC analyzed field from the SMRF scheme (b) on August 14, 2023.

DownLoad: Full-Size Img PowerPoint

Figure  9.  Sea ice concentration analyzed field (left column) and the descent direction (right column) of the MSRF scheme ($ \beta' $ = 35) at iteration 7, 80, and 140.

DownLoad: Full-Size Img PowerPoint

4.4 Comparison between the MSRF scheme and the SMRF scheme

The SIC on August 14, 2023 is also taken as an example to conduct two-dimensional assimilation experiments on the SMRF scheme. The fixed and small filter parameter $ {\alpha} $ in the SMRF scheme is set to 0.1, while the initial value of another variable filter parameter is set to 0.9 and the number of iterations is set to 500. The SIC analyzed field constructed by the SMRF scheme (Fig. 8b) is also very similar to the “true” field (Fig. 5a), demonstrating good multi-scale information extraction ability.

Further qualitative and quantitative analysis are conducted on the SIC analyzed results constructed by the MSRF and SMRF schemes. Figure 10 shows the distribution of differences between the two schemes and observations. As can be seen from the figure, the analyzed results of the two schemes show significant deviations in different degrees in the data-void region, except for the sea ice margin region. The MSRF scheme tends to underestimate the SIC, as shown in the data-void region around (80.5°N, 150°W) marked in the pink pentagrams in Fig. 10a. On the contrary, the SMRF scheme exhibits smaller bias in this region (Fig. 10b), but it has a larger positive bias in the data-void region near (85°N, 88°E) marked in the blue pentagrams than the MSRF scheme. The similar results are also reflected in the histogram of SIC deviation in Fig. 11. Positive deviation indicates that the analyzed value is greater than the observed value, otherwise the reverse. The deviation value of the MSRF scheme has a higher frequency in the interval [–0.3, –0.1) (0.0461) than the SMRF scheme (0.0413), while the frequency in the interval [0.1, 0.3) (0.0288) is lower than the SMRF scheme (0.0340). The main reason for the different results of the two schemes is that different types of filters have different filling structures for the regions of missing observations. From single point observation experiments, it was found that the MSRF scheme has a wider observation propagation range and exhibits anisotropic characteristics in the outer edge region of the range. This makes it more susceptible to the influence of a wider range of observation data in the iterative process of analyzed field. Nevertheless, the SMRF scheme tends to smoothly transition the data-void region based on available observations around it.

Figure  10.  Differences between the analyzed field and the true sea ice concentration field in the MSRF scheme (a) and the SMRF scheme (b). The blue and pink pentagrams are located at (87°N, 88°E) and (80.5°N, 150°W), respectively.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Histogram of the deviation between the analyzed field and the true sea ice concentration field in the MSRF scheme (a) and SMRF scheme (b).

DownLoad: Full-Size Img PowerPoint

The overall performance of the two schemes is further analyzed from the perspectives of regional overall error and computational efficiency. The root mean square error (RMSE) and mean absolute deviation (MAD) are calculated as follows:

where i is the label of the grid point, $ x_{i}^{{\mathrm{a}}} $ represents the analyzed value of the i-th grid point, $ x^{\mathrm{o}}_i$ represents the observed value of the i-th grid point and n represents the total number of non-NaN grid points.

Although both schemes have certain deviations in filling the data-void region, from the perspective of overall regional error (Table 1), the MSRF scheme is significantly smaller than the SMRF scheme. The RMSE and MAD of the MSRF scheme are reduced by approximately 0.6 % and 3.2 % compared to the SMRF scheme, respectively. Moreover, the frequency of the deviation between the analyzed results of the MSRF scheme and the observations within the interval [–0.1, 0.1) (0.9182) is also greater than that of the SMRF scheme (0.9175) (Fig. 11). Further, the model performance of the two schemes over a period of time (August 11–31, 2023) is evaluated. Figure 12 shows the overall regional RMSE and MAD for the two schemes, respectively. It can be seen that in the whole period, the deviations of analyzed results of the MSRF scheme are smaller than those of SMRF scheme. In addition, under the same experiment configuration, the MSRF scheme also obtained smaller errors in the SIC analyzed results of the two schemes on August 8, 2022 (bold values in Table 1). It is worth mentioning that the MSRF scheme has higher computational efficiency while achieving analyzed results closer to the “true” value. The iteration steps of the MSRF scheme are less than half of the SMRF scheme, and the calculation time is only about 1/7 of the SMRF scheme.

Table  1.  Comparison of RMSE, MAD and CPU computation time between the MSRF scheme and SMRF scheme on August 14, 2023 and August 8, 2022 (bold)

Scheme	RMSE	MAD	Iteration step	CPU time/s
MSRF scheme	0.0652	0.0297	215	12.137
	0.0655	0.0289	215	12.558
SMRF scheme	0.0656	0.0307	500	88.936
	0.0660	0.0301	500	85.161

 | Show Table

DownLoad: CSV

Figure  12.  Comparison of RMSE (a) and MAD (b) in the MSRF scheme and the SMRF scheme from August 11 to August 31, 2023.

DownLoad: Full-Size Img PowerPoint

5. Conclusions and discussion

In this study, the SOAR correlation function is introduced into the multi-scale 3D-VAR analysis. The SOAR filter, which only requires two filters, is used to replace the cascade of multiple first-order recursive filters in the SMRF scheme, and the MSRF scheme is proposed. In addition to filtering the independent variables of the cost function, the MSRF scheme also applies SOAR filter to the gradient of the cost function. The descent direction is established by smoothing the sharp changes of the gradient to extract the long wave of the observational residual. As the scale factor gradually decreases, wavelengths of various scales from long wave to short wave are extracted successively. In the two-dimensional SIC experiment, the MSRF scheme has a smaller deviation between the analyzed results and observations compared to the SMRF scheme, and the calculation time is only 1/7 of that of the SMRF scheme.

The SOAR filter combines SOAR function with recursive filters, and only two filtering processes are performed in each direction, greatly improving computational efficiency. Meanwhile, the filter demonstrates good propagation ability for observed signals in multi-scale 3D-VAR applications. In the supplementary experiment, the MSRF scheme also shows good performance when applied to the observation of sea surface temperature with uneven distribution. These show the promising potential of the MSRF scheme for further application in ocean or sea ice forecast systems in the future. However, the MSRF scheme has some parameters that require manual adjustment, such as the selection of fixed-value scale factor and the determination of scale factor decreasing formula. In addition, the selection of initial scale factor $ \beta $ in the scale factor decreasing formula affects the calculation result and efficiency. When the initial $ \beta $ is twice or larger than the pre-set value in the study, more than 60% of the filter parameter $ \alpha $ calculated in all iteration steps is greater than or equal to 0.85. This makes too many iterative steps used to extract long wavelengths, while the extraction of short wavelengths is insufficient. Meanwhile, larger initial $ \beta $ also means higher computational costs. The initial value selected in this study has been verified by a large number of experiments, which can not only extract almost all scale wavelengths, but also have high computational efficiency. In the future, we will attempt to optimize the algorithm to automatically find the optimal initial scale factor and reduce labor costs.

Furthermore, according to the study of Yang et al. (2022), the low-order recursive filter in the SMRF scheme was replaced by a Van Vliet fourth-order recursive filter, and the proposed multi-scale high-order recursive filter (MHRF) scheme has good multi-scale information extraction ability. The MSRF scheme proposed in this study tends to be more isotropic in single point signal transmission performance compared to the MHRF scheme. Moreover, the MAD of the MSRF scheme between the SIC analyzed results and observations on August 14, 2023 reduced by approximately 1.66% compared to the MHRF scheme. It is worth noting that the MHRF scheme has relatively independent filtering processes in various directions and has certain advantages in computational efficiency. Compared with the MSRF scheme, its calculation time is shortened by about 2.62 s. In the future, we will attempt to further improve the MSRF scheme in this study based on high-order recursive filters, in order to achieve better results in terms of computational efficiency and accuracy.