A bias correction scheme for FY-3E/HIRAS-II observation data assimilation

Meteorological satellite data have been extensively utilized in global numerical weather prediction systems and have a positive impact to improve forecast accuracy. In order to correctly assimilate satellite radiance observations in data assimilation systems, the systematic observation biases must be corrected to conform to a Gaussian normal distribution with a mean of 0. By selecting appropriate air-mass predictors through correlation assessment, a two-step bias correction scheme (namely the scan-angle bias correction and the air-mass bias correction) is established in this paper based on radiance observations of FY-3E/ HIRAS-II from 1 to 31 January 2023. The results indicate that FY-3E/HIRAS-II O−B (observation-simulation) bias exhibits scanning angle bias dependence from nadir to limb field of view. Statistics have found that this scanning angle bias does not depend on latitude band. After scan-angle bias correction using statistical scan-angle correction coefficients, the dependence of the O−B biases on the scan angle can be eliminated. The second step is to perform air-mass correction. Our correction scheme is compared with the air-mass bias correction scheme in NCEP-GSI. Although the scan angle influence is also considered in NECP-GSI scheme, it does not account for the water vapor effect in the atmosphere. Consequently, the correction effect is not good for channels with lower peak height of weighting function, resulting in a slightly residual positive bias after correction. The combination of air-mass predictors (model surface skin temperature, model total column water vapor, thickness of 1000–300 hPa, and thickness of 200–50 hPa) selected through importance assessment in this study effectively eliminates the air-mass biases. The systematic biases between observed brightness temperature and background simulated brightness temperature from background atmospheric field for all HIRAS-II channels significantly decrease after bias correction, and the bias distribution essentially follows a Gaussian normal distribution with a mean of 0. The FY-3E/HIRAS-II data assimilation experiments show that the selected air-mass predictors (EXP-2 scheme) is the most effective among the four experiments. The mean O−B and O−A in all channels are the smallest after bias correction. Compared with the independent ERA5 objective analysis fields, the EXP-2 scheme has a significant improvement for the temperature analysis at upper air and near surface. The water vapor profiles of the EXP-2 scheme are the closest to ERA5 at almost all height levels.

1 Introduction

The quality of the numerical weather prediction (NWP) largely depends on the accuracy of the initial atmospheric conditions, provided by the data assimilation system (Auligné and McNally, 2007). Satellite data is an important input observation source to data assimilation systems, it is of great significance to correctly assimilate satellite radiance data to improve the accuracy of the numerical weather prediction (Zhang et al., 2023).

Satellite observations have the advantages of wide coverage and high temporal and spatial resolution, which greatly complement the data gaps in areas lacking conventional observations worldwide. Moreover, satellite radiance in infrared or microwave bands exhibit strong sensitivity to meteorological elements (e.g., temperature, humidity) within the atmospheric structure (including the Earth's surface) (Li et al., 2019). Notably, satellite-borne infrared hyperspectral atmospheric sounders such as AIRS (Atmospheric Infrared Sounder), CrIS (Cross-track Infrared Sounder) and IASI (Infrared Atmospheric Sounding Interferometer) can obtain global meteorological observations with multiple channels and high vertical and spectral resolution. These observations, particularly three-dimensional atmospheric temperature and humidity profiles, have been extensively utilized in global numerical weather prediction with a significant positive impact. In contrast to other infrared hyperspectral atmospheric sounders, HIRAS-II (Hyperspectral Infrared Atmospheric Sounder-II) is carried on FY-3E, the world's first early morning polar-orbiting meteorological satellite. The satellite launched in 2021 effectively fills the data gap of satellite observations within a 6 h assimilation window, ensuring almost 100 % coverage of satellite observations within the assimilation time window. As the world's only infrared atmospheric sounder operating in an early morning polar orbit, it is imperative to assimilate its observations into data assimilation systems (Zhang et al., 2024).

One of the most important steps in data assimilation is bias correction, especially for the bias correction of satellite observation (Yin et al., 2020). Satellite observation currently account for the vast majority of all assimilated atmospheric observations, and satellite observations have a strong influence on the quality of the main output meteorological parameters (e.g., temperature, winds and humidity) from the analysis (including the reanalysis), as well as the additional information generated by the assimilating forecast models (e.g., precipitation, cloudiness and radiative fluxes). During assimilating the meteorological satellite-observed radiance data, the target functional atmospheric data assimilation method based on statistical optimal estimation requires the errors of satellite observations (O) and background simulations (B) all to be an unbiased Gaussian distribution. Therefore, if uncorrected satellite observations are directly absorbed by the data assimilation system, the accuracy of the analyzed fields will be affected, thus affect the forecast accuracy of NWP (Dee, 2004). However, satellite observations (O) and background (B) in fact are always systematically biased. The systematic bias between O and B may originate from observation errors of the instrument itself, simulation errors of the fast radiative transfer model, forecast errors of the NWP system (as the input atmospheric state parameters to radiative transfer models), errors introduced during data preprocessing steps, among others. The sources of systematic bias are complex, but the bias can be quantitatively estimated by statistical O−B. If the unbiased assumption is not valid, the deviation of the observation error ${\mathit{\mu }}_{\mathrm{o}}=\stackrel{\mathrm{‾}}{O-T}$ and the deviation of the simulation error ${\mathit{\mu }}_{\mathrm{b}}=\stackrel{\mathrm{‾}}{B-T}$ must be subtracted from the observation and simulation brightness temperature, respectively. T represents the true value of brightness temperature. In the presence of errors in observation and background, ${\mathit{\mu }}_{O-B}=\stackrel{\mathrm{‾}}{O-B}$, μ_O−B is the deviation between observed and simulated brightness temperatures. This expression provides a basis for bias estimation and bias correction, so the value of μ_O−B can be estimated even in the absence of the true value based on the statistical samples of O−B. In order to properly use satellite observations in the assimilation system, the O−B must first be corrected to produce an unbiased analysis (Zou, 2020).

For a long time, satellite radiance data assimilation has mainly used a bias correction scheme that relies on “air-mass”, which use quantities related to “air state” as predictors to implement bias correction. Its theoretical basis is that there is a linear correlation between the spatial and temporal variations of the satellite observation biases and the predictors, and the higher the correlation, the more obvious the correction effect is (Zou, 2020). McMillin et al. (1989) first proposed a bias correction scheme for the TIROS Operational Vertical Sounder (TOVS) using the observed radiance from MSU (Microwave Sounder Unit) channels 2, 3 and 4 as air-mass predictors. Eyre (1992) adjusted the scheme by incorporating cloud radiance, but the air-mass biases still remained. It is worth mentioning that Harris and Kelly (2001) proposed a revolutionary bias correction scheme that divides the bias into two parts: the unavoidable scan angle bias of cross-track scanning instruments and the air-mass bias caused by the NWP mode and the fast radiative transfer mode. On one hand, the scheme incorporated the latitude dependency into scan angle bias correction, and on the other hand, it replaced observation-based predictors variables with model background-based predictors. Although the scheme has achieved remarkable progress, it is a static off-line scheme that relies on bias correction coefficients calculated in advance based on historical samples and does not account for the evolution of biases with time and with weather systems. Subsequently, many researchers developed variational bias correction schemes. Dee (2004) proposed a variational scheme for adaptive radiance bias estimation and correction based on the ECMWF (European Centre for Medium-Range Weather Forecasts) assimilation system. Zhu et al. (2014) objectively evaluated the effect of the variational correction scheme in NCEP (the National Centers for Environmental Prediction) – GSI (Gridpoint Statistical Interpolation).

Due to the dependence of systematic biases on instruments (or sensors), many scholars have developed specific bias correction schemes suitable for respective satellite instruments. Liu et al. (2007) proposed a bias correction scheme based on the radiance data from the ATOVS (Advanced TIROS Operational Vertical Sounder) instruments onboard NOAA-15/16/17 polar-orbiting meteorological satellites. Li et al. (2016) proposed a bias correction scheme for IASI (Infrared Atmospheric Sounding Interferometer) suitable for the GRAPES assimilation system, based on the approaches of Harris and Kelly. Li et al. (2019) assessed the capability of air-mass predictors for bias correction in the NCEP GSI assimilation system using CrIS radiance data and proposed an improved bias correction scheme based on the periodic characteristics of observation minus background biases. Yin et al. (2020) evaluated the observation quality of the FY-4A/GIIRS (Geostationary Interferometric Infrared Sounder) longwave infrared channels using GRAPES 4D-Var assimilation system and applied an off-line bias correction scheme to correct the O−B biases of these channels. Liu et al. (2024) used the CMA-GFS (the China Meteorological Administration Global Forecast System) model to assimilate FY-3E/HIRAS-II data using a combination of 1000–300 hPa thickness, 200–50 hPa thickness, and 50–10 hPa thickness as predictors.

Most of the above radiance bias correction studies mainly focus on spaceborne microwave radiometers, and there is limited research on observation bias in infrared hyperspectral atmospheric sounders, particularly for the HIRAS-II onboard the early morning polar-orbiting satellite FY-3E launched recently. Liu et al. (2024) use a combination of 1000–300 hPa thickness, 200–50 hPa thickness and 50–10 hPa thickness as predictors to correct FY-3E/HIRAS-II biases, but do not provide the reasons for predictors selection. Therefore, a bias correction scheme suitable for FY-3E/HIRAS-II is established in this paper based on the selection of the optimal air-mass correction predictor combination using its radiance observation from 1 to 31 January 2023. In addition, the scheme is compared with the air-mass correction scheme in NECP-GSI to quantitatively evaluate its bias correction effect. Finally, the potential positive impact of the scheme to enhance the accuracy of data assimilation is check based on one-month assimilation experiments.

2 Data and Model

This article investigates a bias correction study on the radiance observations from the hyperspectral infrared atmospheric sounder HIRAS-II aboard the Fengyun-3E satellite. HIRAS-II is an interferometer Fourier transform spectrometer with continuous coverage of infrared wavelength from 3.92 to 15.38 µm, consisting of 3041 channels (after apodization) with a spectral resolution of 0.625 cm⁻¹. HIRAS-II measures Earth and atmosphere in the conventional mode through a cross-track rotary scan mirror that provides the scan angles vary from −50.4 to +50.4°. Each scan line observes 32 fields of regard (FOR), including 28 contiguous ground targets, 2 cold spaces and 2 onboard blackbody targets. Each FOR consists of a 3 × 3 array of fields of view (FOVs) and the approximate resolution for the nadir FOV is 14 km. A total of 31 d HIRAS-II Level 1 radiance data from 1 to 31 January 2023 is used in our research as the satellite observations (O). The HIRAS-II Level 1 radiance data obtained from the Fengyun Satellite Data Center: http://data.nsmc.org.cn (last access: 18 May 2025).

The National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) 6 h forecast filed valid at 00:00, 06:00, 12:00, and 18:00 UTC, are used as input to the fast radiative transfer model. (Data are available at https://rda.ucar.edu/datasets/d084001/, last access: 18 May 2025). GFS data is a regular latitude-longitude grid data with a spatial resolution of 0.25° × 0.25° and the atmosphere is divided into 41 vertical layers from 1000 to 0.01 hPa. The atmospheric state parameters include the profile of temperature, humidity and ozone, et al.

The GFS data are spatialy-temporaly matched to HIRAS-II FOVs as follows: for each HIRAS-II field of view (FOV), perform bilinear interpolation on GFS forecast value by selecting the 4 closest grid points. Since the samples collected in this study are clear-sky observations over sea, the clear-sky atmosphere does not change much within 0 to 3 h, so based on the observation time of each HIRAS-II FOV, select the two spatially matched GFS data that are closest to HIRAS-II's observation time for linear interpolation.

The fast radiative transfer model employed in this study is RTTOV v12.3. RTTOV is a widely used fast radiative transfer model, suitable for satellite radiometers, spectrometers, and interferometers in visible, infrared and microwave bands. It can simulate satellite-observed brightness temperature based on the atmospheric and surface state vectors input by users (Saunders et al., 2018). The spatially and temporally matched GFS data are used as input to calculate simulated brightness temperature (B) for FY-3E/HIRAS-II assuming clear-sky condition.

3 Quality control

To ensure the rationality and effectiveness of the HIRAS-II data used for statistics, a quality control process is performed before calculating the bias correction coefficients. The following steps are included:

1. Cloud detection. Cloud detection is a crucial step before satellite data assimilation. Since the infrared hyperspectral atmospheric sounder cannot penetrate clouds during the measurement, its observations are highly susceptible to the influence of clouds and precipitation, and simulations of the fast radiative transfer modes under cloudy conditions have considerable uncertainty, so selecting the clear-sky field of view with high confidence can greatly reduce simulation errors in the fast radiative transfer models. In this study, the temporally and spatially matched FY-4A/AGRI (Advanced Geostationary Radiation Imager) cloud mask products are employed for FY-3E/HIRAS-II FOV cloud detection. When the viewing angle of the satellite is large, the field of view becomes distorted and spatial resolution decreases. To minimize the impact of field-of-view deviation, this study selects only the observations from the FY-4A/AGRI with a viewing zenith angle less than 60° (The region is approximately from 55° S to 55° N and from 50 to 155° E) for matching. The HIRAS-II has a coarse spatial resolution of 14 km at nadir, whereas the AGRI cloud mask product has a higher spatial resolution of 4 km. Therefore, approximately 4 × 4 AGRI pixels are co-located within each HIRAS-II FOV. These clear HIRAS-II FOVs are retained during statistics when all the spatially matched AGRI pixels within the HIRAS-II FOV are identified as clear. The FY-4A/AGRI Level 2 4 km cloud mask product with 4 km resolution can be downloaded from the Fengyun Satellite Data Center: http://data.nsmc.org.cn (last access: 18 May 2025).

2. Surface type detection. The surface emissivity calculated by the fast radiative transfer model is relatively accurate when the underlying surface within a satellite field of view is more uniform. The underlying surface over land is relatively complex, leading to some errors during simulation, while ocean surfaces are relatively uniform. The surface type of each FOV is determined based on the FY-3E/HIRAS-II Level 1 LSM (Land Sea Mask) product. Only these ocean satellite-observed scenes are kept during statistics.

3. Data Thinning. Each field of regard (FOR) of FY-3E/HIRAS-II consists of 3 × 3 arranged FOVs with 14 km spatial resolution. Only the observations from the central fifth FOV within each FOR are retained, so the spatial resolution is approximately 45 km after thinning.

4. Clear-sky detection using window channels. To further eliminate the influence of clouds, a threshold test based on the O−B biases at these window channels 925, 970, and 1111 cm⁻¹ is performed for each thinned FOV. The FOV with O−B bias exceeding −4 to 4 K at any channel is discarded.

5. Outlier detection. The FOV with observed brightness temperature in any channel exceeding the value range (150–350 K) will be discarded. If the O−B biases in any channels exceed three times its standard deviation, the FOV is also discarded.

After the above quality control steps, a total of 67 534 samples were counted. 34 844 samples from 1 to 14 January 2023 are used as the training data for fitting bias correction coefficients, and 32 690 samples from 15 to 31 January 2023 are used as the testing data to examine the correction effect.

4 Bias Correction Scheme

HIRAS-II radiance bias correction is divided into two steps referring to the offline static bias correction method of Harris and Kelly (2001) in this study.

4.1 Scan-angle bias correction

The field of view is susceptible to deformation as the scan angle increases when the satellite sensor performs a cross-track scanning to both sides of scan lines, which leads to unavoidable observation biases relative to the nadir FOV. The optical pathlength is longer as the scan angle increases, also leading to scan angle biases. In addition, although the fast radiative transfer models are become increasingly sophisticated and accurate, they can not only simulate the brightness temperature at the nadir, but also simulate the brightness temperature at other scan angles. However, atmospheric inhomogeneity increases with the increase of scan angle. If the model cannot fully simulate atmospheric inhomogeneity, it may also lead to bias in the simulated brightness temperature that depend on the scan angle (Zou et al., 2011). For each channel, the global or regional mean bias of every scan position (angle) relative to the central nadir position is calculated as:

$$\begin{array}{}\text{(1)}& S\left(\mathit{\theta }\right)=\stackrel{\mathrm{‾}}{R}\left(\mathit{\theta }\right)-\stackrel{\mathrm{‾}}{R}(\mathit{\theta }=\mathrm{0})\end{array}$$

Where $\stackrel{\mathrm{‾}}{R}$ refers to the mean radiance at different scan angles, θ represents the scan position (scan angle) and S denotes the scan bias. A significant improvement made by Harris and Kelly (2001) to this scheme is incorporating the dependency of biases on latitude band. The Earth is divided into 18 latitude bands with 10° interval and the correction coefficient is computed respectively. After a long period of lager sample size statistics, it has been shown that the scan biases of HIRAS-II do not exhibit a pronounced dependence on latitude band as found in microwave instruments. Therefore, the influence of latitude is not considered for the scan bias correction in this study.

4.2 Air-mass bias correction

The O−B biases generally exhibit variations associated with the properties of air-mass (and the surface) due to the inaccuracies in the radiance calibration of satellite instruments, the error from the fast radiative transfer models and NWP systems. The air-mass bias correction scheme primarily involves establishing a multivariate linear regression equation that relates the air-mass predictors $x_i(i=\mathrm{1},\mathrm{2},\mathrm{\dots },n)$ to the air-mass biases, the air-mass biases r_j for each channel j can be calculated:

$$\begin{array}{}\text{(2)}& r_j={\sum }_{i=\mathrm{1}}^n{a}_{ji}{x}_i+c_j\end{array}$$

Where a_ji and c_j are calculated by least square fitting using a large number of samples.

$$\begin{array}{}\text{(3)}& a_{ji}={\sum }_{k=\mathrm{1}}^n\langle D_j,x_k\rangle {\left[\mathit{X},\mathit{X}\right]}_{ki}^{-\mathrm{1}}\end{array}$$

Here, 〈…,…〉 represents covariance, k is the sample number, X is the vector of x_k, and D_j denotes the O−B bias in channel j.

The success of air-mass bias correction depends on the selection of air-mass predictors. The commonly used air-mass predictors in the ECMWF and NCEP-GSI assimilation systems are shown in Table 1. Here, p represents the air-mass predictors in the ECMWF assimilation system and p^′ represents the air-mass predictors in GSI. The predictor p₁, p₂ , p₅ and p₆ reflect the mode background-errors at different layers and various dependencies within the forward model. The predictor p₃ represents the systematic error of near-surface channels. In addition, it can compensate somewhat for the different emissivity characteristics of different surface types for these channels (Harris and Kelly, 2001). The water vapor is one of the main input components in the fast radiative transfer mode, the predictor p₄ can reflect the model simulated error to some extent. The predictors p₈–p₁₀ are employed to correct residual biases after scan bias correction (Dee and Uppala, 2009). The Parameter thickness is calculated as

$$\begin{array}{}\text{(4)}& {\mathrm{Pred}}_{\mathrm{thickness}}=kth\times {\sum }_i^{N-\mathrm{1}}tv\left(i\right)\times \ln {\displaystyle \frac{P\left(i\right)}{P(i+\mathrm{1})}}\end{array}$$

Here, $kth={\mathrm{gas}}_{\mathrm{constant}}/\mathrm{gravity}({\mathrm{gas}}_{\mathrm{constant}}=$ 287.0  J K⁻¹ kg⁻¹, gravity = 9.81 N kg⁻¹), N is atmosphere levels, P is atmospheric pressure, tv is the parameter characterizing atmospheric temperature and humidity and calculated as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}tv\left(i\right)& ={\displaystyle \frac{T\left(i\right)+T(i+\mathrm{1})}{\mathrm{2}}}\\ & \times \left[\mathrm{1.0}+\mathrm{0.608}\times {\displaystyle \frac{q\left(i\right)+q(i+\mathrm{1})}{\mathrm{2}}}\right]\end{array}\end{array}$$

where, T and q represent RTM level temperatures and moistures, respectively.

Table 1 Bias predictors implemented in ECMWF and GSI.

In the NECP-GSI assimilation system, the predictor $p_{\mathrm{0}}^{\prime }$ represents a global constant offset. The predictor $p_{\mathrm{1}}^{\prime }$ is a function of the satellite scan angle θ and is primarily used to correct residual scan biases. The predictor $p_{\mathrm{2}}^{\prime }$ is only used for the correction of clear-sky microwave instrument radiance over the ocean to eliminate residual cloud interference. For non-microwave instruments, the value of the predictor is set to 0. $p_{\mathrm{3}}^{\prime }$ and $p_{\mathrm{4}}^{\prime }$ is the predictor of “temperature lapse rate”, τ and T represents the vertical variation rate of transmittance and temperature, respectively. The predictor $p_{\mathrm{4}}^{\prime }$ represents the convolution of τ and T, and $p_{\mathrm{3}}^{\prime }$ is the square of the former. The predictor $p_{\mathrm{3}}^{\prime }$ and $p_{\mathrm{4}}^{\prime }$ reflect instrument and RTM errors. When the frequency of channel is shifted or the spectral settings in the RTM are inaccurate, the calculated transmittance profile will move up/down in the atmosphere (if the atmosphere is not isothermal). If the transmittance is moved up slightly, the weight function of the channel is shifted upwards and the brightness temperature should be decreased if the temperature decreases with height. Conversely, the brightness temperature will increase in the case of a temperature inversion. (Zhu et al., 2014).

Harris and Kelly (2001) as well as Liu et al. (2007) used four predictors (model surface skin temperature, model total column water vapor, thickness of 1000–300 hPa and thickness of 200–50 hPa) as the optimal combination to the TOVS's air-mass bias correction. Auligné et al. (2007) employed a predictor combination of thickness of 1000–300 hPa, thickness of 200–50 hPa, thickness of 50–5 hPa and thickness of 10–1 hPa to correct the air-mass biases in ATOVS. The predictor combinations used in the above-mentioned studies are primarily designed for microwave instruments. The predictors used in the ECMWF for all infrared hyperspectral instruments operating on polar-orbiting platforms are summarized in Table 2 (Auligné et al., 2007; Collard and McNally, 2009; Eresmaa et al., 2017). In order to evaluate the optimal combination for FY-3E/HIRAS-II, several typical channels (737.5, 900, 1040, 1279.375, 1476.25 and 1809.375 cm⁻¹) are taken as examples in our study. The importance of predictors p₁–p₆ for FY-3E/HIRAS-II is evaluated based on the diagnostic scheme proposed by Auligné et al. (2007). Since we use a two-step bias correction scheme, angle-related predictors are not involved in the assessment. This diagnostic scheme evaluates the bias correction effect based on the ability of each predictor to reduce the root mean square error of radiance biases. The results are shown in Fig. 1. Figure 1a and b show the spectral positions and the weighting functions (WF) of each selected channel. The wavenumber of HIRAS-II channel No. 141, 401, 626, 1008, 1323 and 1855 are 737.5, 900, 1040, 1279.375, 1476.25 and 1809.375 cm⁻¹, respectively. The peak heights of weighting function are 535.2, 1070.9, 29.1, 852.8, 300 and 500.2 hPa in sequence. The six channels correspond to the long-wave CO₂ absorption band, the long-wave window channel, the O₃ absorption band, and the water vapor absorption bands at three different heights in turn. Figure 1c indicates the importance assessment of different predictors in different channels, where the horizontal axis indicates the channel No. and the vertical axis is the diagnostic coefficients. A higher coefficient indicates a stronger correlation. The different colored bars represent different predictors. It can be seen from Fig. 1c that the O₃ channel 626 and the water vapor channel 1323 with higher height of weight functions have a higher correlation in predictors thickness of 10–1 and 50–5 hPa. However, other four channels commonly used for data assimilation systems are strongly correlated with the model background surface skin temperature, model total column water vapor, thickness of 1000–300 hPa, and thickness of 200–50 hPa. Data assimilation systems generally do not assimilate strong O₃ absorption channels at present and the water vapor content in upper atmosphere is scarce, so a predictor combination including model surface skin temperature, model total column water vapor, thickness of 1000–300 hPa and thickness of 200–50 hPa is selected to correct the air-mass biases for HIRAS-II in this research. Furthermore, to assess the effectiveness of our bias correction scheme a comparison is made in this study between the correction results of the NCEP-GSI predictor combination ($p_{\mathrm{0}}^{\prime }$–$p_{\mathrm{4}}^{\prime }$) and the proposed method.

Table 2 Predictors used in bias correction of different infrared hyperspectral instruments for ECMWF. The p_x in the table correspond to Table 1.

Figure 1 (a) Simulated HIRAS-II brightness temperature spectrum (black curve) by the RTTOV using the US76 standard atmospheric profile, the selected 6 typical channels (red dots) and the assimilated 485 channels (black dots). (b) weighting function of the selected channels. (c) relevance diagnostics of bias-correction.

As we discussed earlier, model surface skin temperature, model total column water vapor, thickness of 1000–300 hPa and thickness of 200–50 hPa are selected as predictors for bias correction. The multiple linear regression model requires a linear relationship between continuous independent variables and dependent variables. Therefore, we examine the correlation between HIRAS-II O−B bias and predictors p₁–p₆ based on samples from 1 to 31 January 2023. Figure 2a–f show the distribution of a typical HIRAS-II channel 1855 (1809.375 cm⁻¹, the peak height of weight function is 500 hPa) O−B bias with predictors p₁–p₆, respectively. In order to obtain a correct analysis of the relationship between the O−B bias and predictors p₁–p₆, the O−B biases eliminate the influence of scan bias. The horizontal coordinate of each subgraph is the value of each predictor, the vertical coordinate is the O−B bias and the color represents the number of samples. The black solid line in figure shows the first-order polynomial fit of the O−B bias and predictors, the correlation coefficients are given in the subtitles of each figure. From Fig. 2a–f, there is a certain linear correlation between the O−B bias and predictors p₁–p₄, with the highest correlation coefficient of 0.49, and there is no obvious linear relationship between the predictors p₅–p₆.

Figure 2 Scatterplots of (a) skin temperature, (b) total column water vapor, (c) thickness of 1000–300 hPa, (d) thickness of 200–50 hPa, (e) thickness of 50–20 hPa and (f) thickness of 10–1 hPa with respect to HIRAS-II O−B for channel 1855 (1809.375 cm⁻¹, 500 hPa). Black curves show the first-order polynomial fitting. Color shade represents data number.

5 Result

The bias correction coefficients are calculated using above method based on FY-3E/HIRAS-II observations from 1 to 14 January 2023. Subsequently, the FY-3E/HIRAS-II O−B biases in all channels from 15 to 31 January 2023 are corrected based on the statistical coefficients. The correction results of the selected six typical channels mentioned above are analyzed in detail.

5.1 The result of scan bias correction

The variation of mean biases with the scan angles before and after scan-angle correction for the six typical channels is illustrated in Fig. 3. The dot-dash line represents the bias distribution before correction and the dashed and solid lines represent the bias distribution after correction for scan angle and air-mass, respectively. The vertical axis is the mean O−B biases and the horizontal axis is the positions of the field of regards (FORs) for each scan line (i.e., the scan angle) with the nadir between FOR14 and FOR15. It is can be seen from Fig. 3 that the biases change with the scan angles for all channels before correction. The bias increases as the scan angle increases, especially for the lower tropospheric water vapor channel No. 1008 with value up to 1.5 K caused by the larger scan angle with respect to nadir. Additionally, there is an asymmetrical bias distribution on both sides of the nadir in all channels, particularly for the lower height channels 141 and 1008. This is primarily caused by the non-90° inclination angle of the polar orbit satellite, leading to inconsistent latitude of the field of view on both sides of the scan line. The phenomenon resulted in the HIRAS-II observation being higher on one side and a lower on the other side, which in turn causes an asymmetric distribution of O−B bias. After the scan-angle bias correction, the mentioned biases from the limb FOR relative to nadir and the asymmetrical biases on both sides have been eliminated and the mean bias of each scan position is basically consistent. However, there are still some residual biases (dotted line) in the properties of air mass caused by inaccurate simulation (some channels reaching 0.5–1 K). It can be clearly seen from the black solid line in Fig. 3 that the mean biases of all scan positions in all channels approach 0 K after the air-mass correction.

Figure 3 The mean O−B bias of HIRAS- II varies with the scan position before and after the bias correction. The subplots are channel (a) No. 141, (b) 401, (c) 626, (d) 1008, (e) 1323 and (f) 1855 in turn. (The dot-dashed line represents the O−B bias before bias correction, while the dotted line and the solid line are the results after the scan-angle and air-mass correction, respectively.)

5.2 Comparison with GSI's air-mass bias correction scheme

Based on the O−B biases of each HIRAS-II channel after scan-angle bias correction from 1 to 14 January 2023 and corresponding air-mass predictors (model surface skin temperature, model total column water vapor, thickness of 1000–300 hPa and thickness of 200–50 hPa) data, the air-mass correction coefficients are fitted using Eq. (3). Then the coefficients are applied for air-mass bias correction (referred to as EXP-2). Figure 4a to f and g to l show the scatterplots of observed brightness temperature (O) and background simulated brightness temperature (B) for the six representative channels before and after the air-mass bias correction. The dashed line represents the y=x contour and the color shade represent the density of the radiance data number. The values in each subplot are the mean value and standard deviation of O−B. From Fig. 4a, b, and d, it is evident that significant negative biases are exhibited in channels 141, 401, and 1008 (with a lower height of weight function) when the scene temperature is high before the EXP-2 correction, while water vapor channels with a higher height of weight function (Fig. 4e–f) show a relatively warm bias compared to the background simulation (with mean biases ranging from 0.2 to 0.4 K). The scatter of these channels are all concentrated and evenly distributed near the y=x contour after air-mass correction, with mean bias close to 0 K and a significantly reduced standard deviation.

Figure 4 The scatterplots of observed versus RTTOV simulated brightness temperature (a–f) before and (g–l) after air-mass bias correction for HIRAS-II channel No. 141, 401, 626, 1008, 1323 and 1855. Color shade is the data number.

The bias correction scheme proposed in this paper has been compared with the air-mass bias correction scheme of NCEP-GSI to further evaluate its correction effectiveness. The air-mass predictors used in GSI for infrared instruments are $p_{\mathrm{0}}^{\prime }$, $p_{\mathrm{1}}^{\prime }$, $p_{\mathrm{3}}^{\prime }$, and $p_{\mathrm{4}}^{\prime }$ listed in Table 1 (hereafter referred to as EXP-GSI). Figure 5 illustrates the histogram of O−B biases before and after correction using two different schemes. The x axis corresponds to the O−B bias and y axis is the probability density function (PDF) of O−B. The black curves in the figure represents the observation residuals O−B before correction and the red and blue curves are the O−B residuals after correction using the EXP-2 and EXP-GSI schemes, respectively. The EXP-GSI scheme still displays notable positive biases in most channels after correction, especially in the channels with a height of weight function below 500 hPa. Considering that ninety percent of atmospheric moisture is confined below 500 hPa, the omission of water vapor as a predictor in the EXP-GSI scheme could potentially explain the suboptimal correction results. In contrast, the biases in all channels of the EXP-2 scheme are distributed close to a Gaussian distribution centered at zero, with the most significant correction effect in channel 1008. However, the corrected bias distribution of window channel 401 still exhibits a long tail on the left side, which may be attributed to the incomplete removal of cloud contamination scene during quality control.

Figure 5 The probability density function of O−B bias before (the black curves) and after bias correction for HIRAS-2 channel (a) No. 141, (b) 401, (c) 626, (d) 1008, (e) 1323 and (f) 1855. (the red curves for EXP-2 BC and blue curves for EXP-GSI).

In order to examine whether the distribution of O−B biases after correction satisfies a normal distribution with a mean of 0, Table 4 presents the values before and after correction for these two schemes. Mean represents the mean value, STD represents the standard deviation, Kurtosis value indicates the steepness of the sample distribution and Skewness is the asymmetry of the sample distribution. Kurtosis of 3 and skewness of 0 indicate a normal distribution. It can be seen that the kurtosis and skewness values in all channels after EXP-GSI correction have not changed significantly, while both the observation residuals and standard deviations in all channels show a significant reduction after correction by EXP-2. The O−B bias of all other channels expect for channel 141 with higher kurtosis 5.4 are close to a normal distribution with kurtosis value 3 and skewness value 0. This indicates a significant improvement in the correction effect for EXP-2.

Table 3 The statistics before and after bias correction.

5.3 Analysis results from 1-month experiments

It is important to verify the potential impact of different BC methods on HIRAS-II radiance data assimilation. The analyzed fields will deviate from the true state if the model first guess or observation contain systematic errors. This study uses WRFDA V4.4 to validate the effects of different bias correction schemes on NWP. The WRFDA model has been developed by the National Center for Atmospheric Research (NCAR). It provides different methods of data assimilation and can assimilate a wide range of observations. The WRF three-dimensional variational data assimilation system (WRF-3DVar) minimizes the so-called variational cost function J(x) as Eq. (6) (Barker et al., 2004).

$$\begin{array}{}\text{(6)}& \begin{array}{rl}J\left(\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-x_{\mathrm{b}}\right)}^T{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-x_{\mathrm{b}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{y}-H\left(x\right){)}^T{\mathbf{R}}^{-\mathrm{1}}(\mathit{y}-H\left(x\right))\end{array}\end{array}$$

Where x is the atmospheric state vectors, x_b is the first guess (background), B is the background error covariance, y is the observation vector, R is the observation error covariance and H stands for the observation operator by using RTTOV v12.3.

A new HIRAS-II data assimilation module is created in WRFDA by adding the reading and QC (The QC steps as described in Sect. 3, but the MR method (Eyre and Menzel, 1989) was used for cloud detection) interfaces. A total of 485 channels are selected for data assimilation, with specific positions in the spectrum shown as black dots in Fig. 1a. Four parallel experiments are designed to assess the effects of different bias correction on analysis. These experiments differed in predictors and their detailed configurations are shown in Table 3. The initial and boundary conditions for all experiments are provided by the NECP GFS 6 h forecast data. The simulation domain of all experiments is approximately from 0 to 60° N and from 70 to 150° E, with a 9 km grid spacing on 889 × 828 horizontal grids and 60 vertical levels up to 1 hPa. The bias correction coefficients used for experiments are fitted by the predictors and HIRAS-II O−B biases obtained from 17 to 31 July 2023, and the fitted coefficients are used in 1-month assimilation experiments from 1 to 31 August 2023.

Table 4 The setting of three experiments.

This study assesses the effect of different bias correction schemes on NWP based on 1-month assimilation experiments from 1 to 31 August 2023. The mean O−B bias and mean O−A (observation minus analysis) bias at 485 assimilated channels during the 1-month experiments period after QC are plotted in Fig. 6a and b, respectively. The horizontal coordinate is the ordinal numbers of the assimilation channels, arranged from longwave to shortwave. The vertical coordinate is the mean bias. The black curves, blue curves and red curves in figures represent the results of experiment NO BC, experiment EXP-GSI and experiment EXP-2, respectively. As shown in Fig. 6a, the majority of the HIRAS channels without bias correction exhibit significant negative biases up to −2.5 K. The O₃ absorption bands from channel NO. 161 to NO. 182 with an average bias of 5 K due to the fixed default O₃ profiles used in RTTOV. Although the EXP-GSI scheme can reduce the absolute deviation of most channels, it still shows a certain positive deviation. It is consistent with the results shown in Fig. 5. After EXP-2 bias correction, the average biases of almost all channels are around 0 and only channel NO. 183 to NO. 223 and channel NO. 452 to NO. 472 still have a small negative bias (maximum not exceeding −0.32 K). This indicates that EXP-2 scheme is effective in the bias correction for HIRAS-II assimilated channels. From Fig. 6b, there is obvious difference in O−A biases with and without bias correction. The mean O−A bias of the experiment without bias correction (NO BC) shows large deviations in the O₃ absorption band (channel NO. 161 to NO. 182), the near-surface water vapor absorption band (channel NO. 189 to NO. 297) and the shortwave CO₂ absorption band (channel NO. 452 to NO. 485). There is a significant improvement in O−A for all channels relative to O−B after bias correction (both for EXP- GSI and EXP-2). The mean O−A of EXP-2 (red curve) is essentially 0 for all channels except the O₃ absorption band and the shortwave CO₂ absorption band. It shows that the EXP-2 scheme is effective in improving the analysis.

Figure 6 (a) Mean O−B biases for 485 assimilation channels without BC (black curve), with the EXP-GSI BC (blue curve) and with the EXP-2 BC (red curve), (b) Mean O−A biases for 485 assimilation channels without BC (black curve), with the EXP-GSI BC (blue curve) and with the EXP-2 BC (red curve). The results are sampled from the 1-month assimilation experiments from 1 to 31 August 2023. Only data that passed the QC process are used.

In order to specifically analyze the improvement of the analysis by different schemes, three data assimilation experiments at 00:00 UTC 7 August 2023 were selected for analysis. Figure 7 shows the spatial distribution of HIRAS-II O−B biases. The colored dots in figures represent the values of the O−B bias, with colder colors being negative bias and warmer colors being positive bias. The shading indicates the observed FY-4A AGRI brightness temperature of the window channel 13 at a central wavelength of 12 µm. Figure 7a–c show the distribution of O−B biases passed QC in 900 cm⁻¹ channel (1000 hPa) for the NOBC, EXP-GSI and EXP-2, respectively. Figure 7d–f are as above but for 1476 cm⁻¹ channel (300 hPa). It can be seen that the O−B biases without BC have either cold or warm bias (Fig. 7a and d). The most O−B biases after EXP-GSI BC are slightly warmer (Fig. 7b and e), especially in the region from 0 to 30° N. The overall O−B biases after EXP-2 BC is near 0 (Fig. 7c and f) and the correction effect is significant.

Figure 7 Spatial distributions (colored dots) of the HIRAS-II O−B biases without BC (left), with the EXP-GSI BC (middle) and with the EXP-2 BC (right) for (a)–(c) the 900 cm⁻¹ channel (1000 hPa) and (d)–(f) 1476 cm⁻¹ channel (300 hPa) at 00:00 UTC 7 August 2023.

To further validate DA results, the analysis is verified against the ERA5 0.25° × 0.25° data for each experiment. Figure 8a and c show the vertical RMSE (Root Mean Square Error) profiles of temperature and water vapor, respectively. The black solid line represents experiment NOBC, the blue solid line represents experiment EXP-GSI, the red solid line represents experiment EXP-2 and the black dotted line represents experiment EXP-CONV. Compared with experiment EXP-COVN, experiments EXP-2, EXP-GSI and EXP-NOBC assimilated a large amount of HIRAS-II data covering from the ground to the upper atmosphere. Therefore, it can be seen from Fig. 8a that the temperature analysis field accuracy is effectively improved (closest to ERA5), including experiment EXP-NOBC. The RMSE of temperature analysis fields from NOBC, EXP-2 and EXP-GSI significantly decrease relative to EXP-CONV in all atmosphere (Fig. 8a), the RMSE of water vapor analysis fields slightly decrease relative to EXP-CONV below 800 hPa (Fig. 8c). The RMSE profiles from EXP-2 and EXP-GSI normalized by the RMSE from the experiment NOBC are displayed in Fig. 8b and d to validate the influence of different bias correction schemes. The horizontal axis represents the RMSE percentage reduction relative to NOBC. A value less than 100 % indicates a decrease in RMSE after bias correction and values greater than 100 % means the RMSE increase reversely. Figure 8b shows that EXP-2 correction scheme effectively improves the temperature analysis fields accuracy in the upper and near-surface levels. EXP-2 scheme shows a better improvement than EXP-GSI schemes, especially in the upper levels. For the water vapor analysis fields (Fig. 8d), all BC schemes showed less significant improvements relative to NOBC compared to the temperature analysis fields, with changes in RMSE within 0.5 %. It is still the EXP-2 scheme reducing the largest RMSE at 400 to 800 hPa atmosphere. This may be due to the warmer O−B biases after the EXP-GSI correction.

Figure 8 Averaged vertical RMSE profiles of (a) temperature and (c) water vapor compared with ERA5. The vertical RMSE profiles of (b) temperature and (d) water vapor are normalized by the RMSE from the NOBC experiment. Error bars give statistical significance intervals for differences from the NOBC at the 95 % level.

6 Conclusions

This paper establishes a two-step bias correction scheme for the innovation vector (O−B) based on the FY-3E/HIRAS-II radiance data from 1 to 31 January 2023. Furthermore, a cross-comparison is conducted with the NCEP-GSI's air-mass bias correction scheme to objectively evaluate the effectiveness of the scheme. In addition, we briefly investigate the effectiveness of this scheme in the data assimilation system. The main conclusions are as follows:

1. Due to the high sensitivity of the FY-3E/HIRAS-II innovation vector (O−B) to instrument scan angles, it is imperative to perform a scan bias correction. The distribution of scan biases is independent of latitude, so the division of the latitude band is not necessary during scan correction. The biases of limb FOV with respect to nadir and the asymmetry in biases on the two sides of the scan lines have been eliminated after scan-angle correction.

2. The air-mass biases can be effectively eliminated by selecting the optimal combination of air-mass predictors in this study based on correlation evaluation. The systematic biases between the observed brightness temperature and the simulated brightness temperature in all channels are reduced and the standard deviation is also significantly decreased after correction. Additionally, the O−B biases basically follow a Gaussian distribution with a mean of 0.

3. The correction effect of the NCEP-GSI's air-mass bias correction scheme to these channels mid-lower tropospheric layer channels (with a height of weight function below 500 hPa) is unsatisfactory. This could be attributed to the omission of total column water vapor as a predictor in this scheme.

4. Compared the data assimilation analysis fields from different schemes with the independent ERA5 objective analysis fields, the EXP-2 scheme has a significant improvement for the temperature analysis at upper air and near surface. The water vapor profiles of the EXP-2 scheme are the closest to ERA5 at 400 to 800 hPa height levels.

This bias correction scheme is just a preliminary experiment for the FY-3E/HIRAS-II data assimilation. At present, the off-line static correction is adopted and the variation of O−B bias with time and weather system is not considered. The upcoming step will involve implementing a variational bias correction scheme, the bias correction coefficient will change with time and weather system. In addition, this study only validates the effectiveness of the scheme for HIRAS-II and will continue to validate its applicability for similar infrared instruments in the future. Finally, this study only briefly evaluates the impact of the bias correction scheme on data assimilation system. The impact on NWP will be further evaluated in the future in actual extreme weather individual cases (e.g., convective precipitation and typhoons).