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.

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.

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.

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 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.

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.

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):

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

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.

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.

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.

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.