A Random Forest Algorithm Combined with Bayesian Optimization for Atmospheric Duct Estimation

Abstract

Inversion of atmospheric ducts is of great importance in the field of performance evaluation for radar and communication systems. Since the model parameters in machine learning play a crucial role in prediction performance, this paper develops a random forest (RF) model integrated with Bayesian optimization (BO) called BO-RF for atmospheric duct prediction, and the BO is adopted to determine appropriate model parameters during the training process. In addition, the K-fold cross-validation (CV) method is also incorporated into the model to obtain the best model partition and overcome the overfitting problem. To test the performance of the proposed model, the results obtained by the BO-RF are compared with other commonly used methods, such as classical RF, extreme gradient boosting (XGBoost) with/without BO, and K-nearest neighbor (KNN) with/without BO. Comparisons demonstrate that BO-RF has the best accuracy and anti-noise ability for the estimation of duct parameters.

1. Introduction

An atmospheric duct is an anomalous electromagnetic propagation phenomenon that commonly appears in the maritime environment and can have a significant impact on radar and communication systems, such as communication blind areas, increased radar sea surface clutter, etc. Therefore, accurate analysis of radar echoes and efficient estimation of atmospheric duct parameters have become active research fields. Atmospheric ducts generally consist of three types, including evaporation ducts, surface-based ducts (SBDs), and elevated ducts. Evaporation ducts and SBDs are two common types in maritime environments due to the rapid changes in vertical atmospheric and humidity profiles with altitude above the sea surface [1]. It should be noted that SBDs exhibit a lower probability of occurrence than evaporation ducts, but their influence on radar and communication systems is much greater than that of evaporation ducts.

Estimation of atmospheric duct refractivity parameters with the refractivity from clutter (RFC) technique has attracted considerable interest in the field of optimization theory. Numerous methods have been proposed to retrieve atmospheric refractivity parameters. Gerstoft et al. [2,3] published simulation results for an eleven-parameter refractivity profile with five parameters in the vertical direction and six parameters in the horizontal direction, solved using a genetic algorithm, and presented the key steps and the latest research progress regarding the RFC technique. Douvenot et al. [4] first employed the least squares support vector, which belongs to a class of machine learning based on the pregenerated radar sea clutter power (RSCP) database produced by the parabolic wave equation (PWE), to estimate SBD. Wang et al. [5] inverted the evaporation duct using particle swarm optimization with RFC. Zhang et al. [6] proposed a four-parameter evaporation duct model and retrieved refractivity parameters from it. Zhao et al. [7] applied classical simulated annealing to the remote sensing of an atmospheric duct. Yang et al. [8,9,10,11] performed the inversion of an atmospheric duct with standard and improved swarm intelligence optimization and presented a comparison of the machine learning algorithm for evaporation duct estimation. Tepecik et al. [12] introduced a novel hybrid model based on an artificial neural network and a genetic algorithm for the inversion problem of refractivity estimation. Lentini et al. [13] calculated the global sensitivity of radar wave propagation to environmental parameters at the maritime atmospheric boundary layer. Penton et al. [14] studied the influence of rough ocean surfaces on evaporation duct atmospheric refractivity inversion with a genetic algorithm. Pozderac et al. [15] adopted a novel transmit–receive array system named the X-band beacon-receiver array to determine the maritime evaporation duct. Yan et al. [16] proposed an estimation method for evaporation ducts using a neural network. Zhu et al. [17] predicted evaporation duct height based on a multilayer perceptron network. Sit et al. [18,19] showed that deep learning and artificial neural networks can be efficiently used to classify and predict atmospheric ducts. It is of great significance to obtain instantaneous environmental information for assessing the performance of radar and communication systems. It can be distinctly seen from the above study that more and more attention has been paid to the field of monitoring the atmospheric environment in the maritime boundary layer with machine learning due to its superiority in real-time compared with swarm intelligence optimization. Random forest (RF) [20,21] is a powerful machine learning algorithm inspired by the ensemble learning method that works by incorporating many decision trees into one algorithm model to enhance its performance. However, the hyper-parameter in RF that directly relates to the performance of machine learning is difficult to determine. In this paper, a combination of RF with Bayesian optimization (BO) is presented to predict atmospheric duct parameters, where BO is adopted to search for the most appropriate hyper-parameter during the training process. To test the performance of the proposed BO-RF model, the RSCP versus propagation distance is numerically simulated by the PWE method with the combinations of refractivity parameters generated with an orthogonal experiment design (OED) method; the results obtained by the BO-RF are compared with the commonly used methods, such as classical RF, extreme gradient boosting (XGBoost) [22] with/without BO, and K-nearest neighbor (KNN) [23] with/without BO. Comparisons demonstrate that BO-RF has the best accuracy and noise immunity for the estimation of atmospheric duct parameters.

2. Forward Model

The forward model commonly includes two parts: the refractivity model and the radar propagation model. In the subsequent section, we will introduce the basic concept in the field of atmospheric ducts.

2.1. Refractivity Model

As the radio wave travels in the electromagnetic duct, the propagation features of the radio wave are controlled by the index of refraction:

$$n=c/v,$$

(1)

where $c$ is the speed of light in a vacuum and $v$ is the speed of light in the medium. The value of $n$ is close to unity, typically 1.00035; however, even a tiny variation in the refractive index may cause a huge change in the propagation features. Hence, a useful dimensionless physical variable $N$, which is the part-per-million change in the index of refraction, named refractivity, is introduced:

$$N=\left(n-1\right)\times {10}^6.$$

(2)

To further account for the curvature of the earth, another variable $M$ called modified refractivity is provided by

$$M\left(z\right)=N+\frac{z}{R_e}\times {10}^6=N+0.157z,$$

(3)

where $R_e$ is the radius of the earth and $z$ is the height above the earth’s surface. As $dM/dz<0$, an electromagnetic duct that can confine radio wave propagation and act like a waveguide is formed [24,25]. In this paper, our emphasis is placed on the three parameters of SBD shown in Figure 1, according to [26].

$$M\left(z\right)=M\left(0\right)+\{\begin{array}{ll}{c}_{1}z\hfill & z\le z_b\hfill \\ c_1{z}_b-M_d\frac{z-z_b}{z_{thick}}\hfill & z_b<z<z_{thick}\hfill \\ c_1{z}_b-M_d+c_2(z-z_b-z_{thick})\hfill & z\ge z_{thick}\hfill \end{array},$$

(4)

where $M\left(0\right)$ is the modified refractivity at the mean sea level, typically taken as $330 \mathrm{M}-\mathrm{units}$, $z$ describes the height from the mean sea level, $z_b$ and $z_{thick}$ represent the thicknesses of the base layer and trap** layer, $c_1$ and $c_2$ denote the slopes of the base layer and top layer, and $M_d$ is the duct strength. Obviously, there are three parameters in the SBD refractivity profile model, and we will employ the parametrization refractivity model of SBD via $\mathit{M}=\left(z_b,z_{thick},M_d\right)$. The tri-linear refractivity profile for SBD is shown in Figure 1.

2.2. PWE Method

PWE can be divided into a narrow-angle and a wide-angle PWE. The propagation elevation of the narrow-angle PWE has good accuracy when it does not exceed $15°$. An atmospheric duct usually occurs when the trap** angle does not exceed $1°$. Therefore, the narrow-angle PWE is selected for studying the radio wave for long-range propagation in atmospheric ducts [27]. In order to analyze the two-dimensional tropospheric propagation problem, the scalar wave equation taking the index of refraction $n$ into account is provided by [28]:

$$\frac{{\partial }^2\phi }{\partial x^2}+\frac{{\partial }^2\phi }{\partial z^2}+k_0^2{n}^2\phi =0,$$

(5)

where $\phi (x,z)$ represents the electric or magnetic field component, $k_0$ is the wavenumber in a vacuum, and $x$ and $z$ are the range and height, respectively. It is a good approximation if the index of refraction varies slowly with the propagation direction. With the help of the reduced function $\phi (x,z)=u(x,z)e^{-i{k}_{0}x}$ and the paraxial approximation, the PWE is expressed by

$$\frac{{\partial }^{2}u}{\partial z^2}+2i{k}_0\frac{\partial u}{\partial x}+k_0^2(n^2-1)u=0.$$

(6)

Then, the split-step Fourier-transform solution to the PWE can be written as [29]

$$u(x_0+\Delta x,z)=e^{i{k}_{0}M\left(z\right)10^{-6}\Delta x}{ℱ}^{-1}\left\{e^{{-i{p}^2\Delta x/2k}_0}ℱ\left\{u\left(x_0,z\right)\right\}\right\},$$

(7)

where $p$ is the transform variable, and $ℱ$ and ${ℱ}^{-1}$ denote the forward Fourier transform and inverse Fourier transform, respectively. If the initial field $u(x_0,z)$ is provided, the field distribution along the vertical direction at each range interval can be easily marched in a horizontal range.

Once the field distribution of the entire computational region is obtained with Equation (7), the RSCP can be easily expressed in dB at range $x$ by radar theory [2],

$$P_c(x,\mathit{M})=-2L(x,\mathit{M})+10 log(x)+\sigma °+C_0,$$

(8)

where $\mathit{M}$ is the parameterized vector of the atmospheric environment, $\sigma °$ is the normalized scattering coefficient from the sea surface, $C_0$ is a constant that contains the parameters of the radar system, and $L(x,\mathit{M})$ is the propagation loss in the atmospheric duct environment that can be conveniently computed from the field distribution.

3. Optimization Algorithm

3.1. Principle of RF Regression

Decision tree is the primary component of RF and belongs to a class of machine learning that generally includes three main steps: feature selection, decision tree generation, and pruning of the decision tree. RF commonly adopts the decision tree that has not completed the pruning optimization operation. Each decision tree makes a prediction for the sample, and the final prediction can be achieved by taking the average of the combined predictions of all decision trees [30]. It has strong robustness and flexibility in dealing with regression problems. However, it has the shortcomings of falling into the local optimum, leading to high complexity and overfitting.

RF is an integrated learning algorithm using decision tree modeling, consisting of multiple independent decision trees. Different decision trees are constructed from random samples, but there is no relationship between them, and the output of the model is determined by each decision tree. In the process of training decision trees, randomly selected features further improve the model’s accuracy and avoid insufficient generalization ability. The detailed steps are given as follows:

• A bootstrap method is adopted to randomly selected $M$ sample points from the original sample dataset $S$ to construct a decision tree and its training subset.
• The training subset is used to train each decision tree. Assume that the number of input features to the dataset is $N$. The split rule for each node is to randomly select $k^{\prime }$ features from the $N$ features as alternative branch features, and then select the best split node from the $k^{\prime }$ features to divide the left and right subspaces. Each split node performs a binary test on each subset, and the test result is sent to the left or right sub-node. The test randomly chooses a subset of features and finds a value with the lowest mean square error $\kappa$ to group and determine its optimal branch, and it can be expressed as

  $$\kappa =\frac{1}{D_1}{\displaystyle {\sum }_{i\in D_1}{\left(E_i-C_1\right)}^2+\frac{1}{D_2}{\displaystyle {\sum }_{i\in D_2}{\left(E_i-C_2\right)}^2}},$$

  (9)

  where $E_i$ is the true value of the ith sample, and $C_1$ and $C_2$ are the sample predictions for the left and right subspaces $D_1$ and $D_2$, respectively.
• Multiple decision trees can be generated via the first two steps above. The final prediction of each decision tree depends on the average of the leaf nodes where the sample points are located. The accuracy of RF is further improved by optimizing parameters such as the number of decision trees, the maximum depth of the decision tree, the minimum number of samples required for node division, the minimum sample number of the leaf nodes, and the number of features in the tree.
• Finally, a RF model is produced by taking the average of the multiple decision trees as follows:

  $$f(x)=\frac{1}{k^{\prime }}{\displaystyle \sum _{i=1}^{k^{\prime }}{E}_i}(x),$$

  (10)

  where $f(x)$ denotes the combined regression model.

3.2. BO Algorithm

BO is a novel and efficient hyper-parameter tuning algorithm [31,32]. In BO, the objective function distribution is first evaluated by Bayes’ theorem, and the new hyper-parameter for the algorithm are then selected according to the optimal sampling point distribution of the objective function. In addition, the update of the hyper-parameter for the next time depends on the objective function value of the optimal sample point. BO is mainly composed of the probabilistic proxy model and the acquisition function. In the following, we will briefly introduce the two parts.

3.2.1. Probabilistic Proxy Model

The updated formula in the probabilistic proxy model is given by

$$p(f|D_{1:t})=\frac{p(D_{1:t}|f)p(f)}{p(D_{1:t})},$$

(11)

where $f$ represents the objective function, $D_{1:t}=\left\{(x_1,f_1),(x_2,f_2),L(x_t,f_t)\right\}$ denotes the collection of sampling point, $p(D_{1:t}|f)$ is the likelihood distribution of observation set, $p(f)$ is the prior probability distribution of $f$, and $p(f|D_{1:t})$ represents the posterior probability distribution of $f$ and is used to modify and test the unknown target function.

The probabilistic proxy model can be divided into parametric and non-parametric models. The Gaussian process is a commonly adopted non-parametric model that has a strong ability to fit functions and is widely used in the field of regression.

A Gaussian process is entirely defined by the mean function $m(x)$ and the covariance function $K(x,x^{\prime })$ as

$$f(x)~GP\left(m\left(x\right),K\left(x,x^{\prime }\right)\right),$$

(12)

where the mean function is a vector and the covariance function is a matrix.

If we assume that there is a collection of sample points $D=\left\{(x_{1:t},y_{1:t})\right\}$, its covariance matrix is given by

$$K=\left[\begin{array}{cccc}k(x_1,x_1)& k(x_1,x_2)& \cdots & k(x_1,x_t)\\ k(x_2,x_1)& k(x_2,x_2)& \cdots & k(x_2,x_t)\\ ⋮& ⋮& \ddots & ⋮\\ k(x_t,x_1)& k(x_t,x_2)& \cdots & k(x_t,x_t)\end{array}\right].$$

(13)

When a set of sample $(x_{t+1},f_{t+1})$ is added to the collection of historical evaluation points, supposing that $k=\left[\begin{array}{cccc}k\left(x_{t+1},x_1\right)& k\left(x_{t+1},x_2\right)& \dots & k\left(x_{t+1},x_t\right)\end{array}\right]$, the corresponding covariance matrix can be updated by

$$K=\left[\begin{array}{cc}K& k^T\\ k& k(x_{t+1},x_{t+1})\end{array}\right].$$

(14)

Using the updated covariance matrix, the posterior probability distribution of $f_{t+1}$ is computed from the first $t$ samples

$$\{\begin{array}{l}p(f_{t+1}|D_{1:t},x_{t+1})~N(\mu ,{\delta }^2)\hfill \\ \mu =k^T{K}^{-1}{f}_{1:t}\hfill \\ {\delta }^2=k(x_{t+1},x_{t+1})-k^T{K}^{-1}k\hfill \end{array}.$$

(15)

3.2.2. Acquisition Function

In the Gaussian process, the acquisition function is adopted to balance the ability of exploration and exploitation and search for the optimal objective function value for the next iteration to obtain the optimal hyper-parameter. Currently, the probability of improvement (PI), expected improvement (EI), and upper confidence bound (UCB) are widely used for acquisition functions. In this paper, PI is selected as the acquisition function, and it can be expressed as follows:

$${\partial }_t(x;D_{1:t})=p(f(x)\le v^{*}-\xi )=\varphi \left(\frac{v^{*}-\xi -{\mu }_t(x)}{{\delta }_t(x)}\right),$$

(16)

where $v^{*}$ represents the current optimal function value, $\varphi (\cdot )$ is the cumulative density function of the standard normal distribution, and $\xi$ is the balance parameter. The problem of local optimum can be overcome to a certain extent and the optimal value can be globally searched by adjusting the parameter $\xi$.

4. A Hybrid BO-RF Model

4.1. Orthogonal Experimental Design Method

In order to train the RF model, an OED method is used to generate the combinations of refractivity parameters. It is well known that OED is a robust way to effectively handle multi-dimension and multi-factor problems because it can select a few representative combinations in the experimental database to obtain the optimal result. The key in OED is the orthogonal array (OA) that can reflect the characteristics of each dimension. Thus, we first need to construct an appropriate OA. Here, we use $L_M(Q^N)$ to represent an OA with N factors and $Q$ levels per factor, where M is the number of combinations of levels. Then, we take a three-factor problem as an example by setting $Q=3$ and $N=3$ to obtain an orthogonal table $L_9(3^3)$ [33,34,35]:

$$L_9(3^3)=\left[\begin{array}{ccc}1& 1& 1\\ 1& 2& 2\\ 1& 3& 3\\ 2& 1& 2\\ 2& 2& 3\\ 2& 3& 1\\ 3& 1& 3\\ 3& 2& 1\\ 3& 3& 2\end{array}\right].$$

(17)

In $L_9(3^3)$, there are three columns, showing that it is suitable for the problem with three parameters at most, and nine rows, demonstrating that one can do nine experiments for a combination with different levels. Each row in (17) represents an experiment with different levels. Normally, we should carry out $3^3=27$ combinations of experiments for a three-factor and three-levels-per-factor problem. However, this task can be reduced to $3^2=9$ combinations of experiments with the help of OA.

An orthogonal expansion is made based on OED. Assuming the lower and upper boundaries of a D-dimension parameter are $u=\left(u_1,u_2,\dots u_D\right)$ and $t=\left(t_1,t_2,\dots t_D\right)$, the range of the corresponding cross-solutions are determined by

$$low=\left(\min \left(u_1,t_1\right),\min \left(u_2,t_2\right),\dots \min \left(u_D,t_D\right)\right),$$

(18)

$$up=\left(\max \left(u_1,t_1\right),\max \left(u_2,t_2\right),\dots \max \left(u_D,t_D\right)\right).$$

(19)

In the cross-solution, the range of the ith dimension variable $x_i$ is $\left[lo{w}_i,u{p}_i\right]=\left[\min (u_i,t_i),\max (u_i,t_i)\right]$. The ith dimension variable $x_i$ can be quantized into $Q$ levels according to the following formula:

$$l_{i,j}=\{\begin{array}{ll}lo{w}_i\hfill & j=1\hfill \\ lo{w}_i+\frac{j-1}{Q-1}(u{p}_i-lo{w}_i)\hfill & 2\le j\le Q-1\hfill \\ u{p}_i\hfill & j=Q\hfill \end{array}.$$

(20)

Then, we use $L_M(Q^N)$ to construct $M$ feasible combinations. More detailed information can be found in [33,34,35].

The inversion of SBD is a three-parameter problem, and the lower and upper boundaries are given by $t=\left(0,20,20\right)$ and $u=\left(100,100,100\right)$, respectively [4]. Let us first construct an appropriate OA $L_M(Q^N)$ with $Q=100$, $N=3$, and $M=\mathrm{10,000}$ for the purpose of generating $\mathrm{10,000}$ training datasets [4].

4.2. Range of the Hyper-Parameter in RF

The hyper-parameter in RF mainly includes the number of decision trees, maximum depth, minimum samples of split, minimum samples of leaf, and maximum features in the tree. The range of hyper-parameter is shown in

Table 1. Range of hyper-parameter in RF.

Description	Hyper-Parameter	Search Range
number of decision trees	n_estimator	(10, 500)
maximum depth of trees	max_depth	(5, 50)
minimum samples of split	min_sample_split	(1, 50)
minimum samples of leaf	min_sample_leaf	(1, 50)
maximum number of features	max_features	(4, 100)

.

4.3. BO-RF Model

The core of the RF model lies in the selection of the hyper-parameter. In this paper, a BO-RF is adopted to predict the refractivity parameter of an atmospheric duct, and the BO is incorporated into the RF to find the best hyper-parameter. In addition, we apply K-fold cross-validation (CV) technology [36] to the BO-RF, which is commonly used to evaluate the training effect and avoid overfitting phenomenon. The CV makes full use of limited data, and the evaluation results may be as close as possible to the performance of the model on the test set, which can be used as an indicator for model optimization.

In BO-RF, the input and output datasets are first used to train the RF model and optimize the RF parameters. Then, the probabilistic proxy model for the estimation of atmospheric duct is constructed and optimized to obtain the best hyper-parameter in RF. Up to this point, the BO-RF model has been successfully established. The flowchart of the BO-RF is shown in Figure 2 for atmospheric duct estimation. The corresponding optimization process is given as follows:

• An OED method is used to generate the combinations of refractivity parameters.
• The split-step Fourier-transform solution to the PWE is adopted to generate the RSCP samples according to the combinations of refractivity parameters.
• A bootstrap method is used to produce decision trees from the training set and then form a random forest.
• A group of hyper-parameter is randomly selected within the scope of the parameter to produce the corresponding initialization point of the sample, then substituted into the RF model to obtain the corresponding objective function $f$ and the initial sample set $D$.
• The probabilistic proxy model in BO is applied to fit $\left(x,f\right)$ and find the most possible evaluation point $x_t$ in the sample, which can be used to obtain the optimal acquisition function. Furthermore, the hyper-parameter is then assigned to RF model to find the objective function value $f_t$ of the optimal sample point and use it as the basis for selecting the hyper-parameter next time.
• A new set of sampling points $\left(x_t,f_t\right)$ is added to the historical sample set $D_{t-1}$, and the Gaussian process is adjusted to optimize the objective function in it.
• The BO stops updating as soon as it reaches the maximum number of iterations. That is to say, the best hyper-parameter, the optimal value of the objective function, and the corresponding BO-RF model are determined.
• The test set is then used to examine the BO-RF model.
• The results are evaluated and analyzed.

5. Results and Discussion

5.1. Simulations of the RSCP in SBD

In this section, the RSCP versus propagation distance are numerically simulated by Equation (8) via the PWE method. In order to estimate the refractivity profile with the RFC technique, the RSCP should be sensitive to the refractivity parameters of SBD. In the simulations, the radar system parameters are given in

Table 2. Radar system parameters.

Parameters	Value
Frequency	2.84 GHz
Transmitting power	91.4 dBm
Beamwidth	0.4°
Antenna height	30.78 m
Transmitting antenna gain	52.8 dB
Polarization	VV

, and the parameterized vector of refractivity profile for SBD is $\mathit{M}=\left(30,42,71\right)$.

Firstly, the curve of RSCP versus the propagation distance at different receive heights are shown in Figure 3. It can be clearly seen that the RSCP in SBD is distinctly sensitive to the receive height.

Furthermore, the coverage diagrams are presented in Figure 4, Figure 5 and Figure 6 to examine the change in RSCP in SBD with every dimension of the parameterized vector. It is obvious that the spatial distributions of RSCP change as each parameter in refractivity profile changes. In other words, the spatial distributions of RSCP vary obviously with the parameters in SBD.

5.2. Estimation of SBD

To test the performance of the BO-RF, the estimation study of three parameters SBD with RFC is carried out, and the results obtained by the BO-RF are compared with the results obtained from commonly used methods, such as classical RF, XGBoost with/without BO, and KNN with/without BO. The radar system parameters are identical to those given in Table 2. In order to fully consider the propagation characteristics of the SBD and reduce the number of training samples, the OED method is adopted to randomly select 10,000 sets of samples of the refractivity parameter $\left(z_b,z_{thick},M_d\right)$, namely a matrix with ten thousand rows and three columns, and 10,000 sets of RSCP from 10.2 km to 50.0 km at an interval of 200.0 m, namely a matrix with ten thousand rows and two hundred columns, generated by the PWE method. It should be noted that we randomly selected 80% of the samples used for training the model, and the other 20% were used to test the regression model. Furthermore, 5% Gaussian noise is also added to 20% of the samples to examine the anti-noise ability. Finally, the chosen RSCP matrix with eight thousand rows and two hundred columns and the refractivity parameter matrix with eight thousand rows and three columns are treated as the input and output matrix to train the regression model, repsectively. In this paper, the default and BO-optimized hyper-parameter are given in

Table 3. Default and BO-optimized hyper-parameter.

Hyper-Parameter	Default Value	Optimized Value
n_estimator	30	382
max_depth	6	42
min_sample_split	4	2
min_sample_leaf	4	1
max_features	auto	50

.

In addition, the following three evaluation indicators, namely the coefficient of determination R², mean absolute error (MAE), and mean square error (MSE), are introduced to test the performance of the BO-RF model for atmospheric duct estimation. Smaller MAE, MSE, and closer R² value to one mean that the algorithm has better performance. The evaluation indicators are defined as follows:

$${\mathrm{R}}^2=1-\frac{{\displaystyle \sum _{i=1}^{m^{\prime }}{\left(X_i-Y_i\right)}^2}}{{\displaystyle \sum _{i=1}^{m^{\prime }}{\left(\overline{Y}-Y_i\right)}^2}},$$

(21)

$$\mathrm{MAE}=\frac{1}{m^{\prime }}{\displaystyle \sum _{i=1}^{m^{\prime }}\left|X_i-Y_i\right|},$$

(22)

$$\mathrm{MSE}=\frac{1}{m^{\prime }}{\displaystyle \sum _{i=1}^{m^{\prime }}{\left(X_i-Y_i\right)}^2},$$

(23)

where $X_i$ is the predicted ith value, $Y$ is the actual ith value, $\overline{Y}$ is the mean of the true values, and $m^{\prime }$ is the number of datasets.

Firstly, a more informative and truthful evaluation indicator, the coefficient of determination R² is adopted to assess the generalization ability of the regression model [37] based on the test sets with different regression models with/without noise, and the values of the coefficient of determination R² are shown in

Table 4. Coefficient of determination R² for the regression model.

Algorithms	R2 without Noise	R2 with Noise
RF	98.92%	98.46%
BO-RF	99.93%	99.82%
KNN	90.67%	87.51%
BO-KNN	91.76%	88.72%
XGBoost	94.68%	93.73%
BO-XGBoost	95.96%	94.83%

, in which the best results are highlighted in bold. From Table 4, we can see that the BO-RF model obtains the best results in the with/without noise case. That is to say, the BO-RF model exhibits good generalization performance and anti-noise ability compared with the other regression models.

Then, the corresponding statistical analysis with respect to MAE and MSE is used to further test the accuracy of the regression model. It can be observed from

Table 5. Comparison of the MAE and MSE for the test results without noise and with 5% Gaussian noise.

Algorithms	Parameters	MAE		MSE
		Without Noise	With 5% Gaussian Noise	Without Noise	With 5% Gaussian Noise
RF	$z_b$/m	0.48	0.49	0.44	0.82
	$z_{thick}$/m	0.35	0.67	0.71	1.18
	$M_d$/M-units	0.48	0.50	0.40	0.95
BO-RF	$z_b$/m	0.14	0.39	0.05	0.62
	$z_{thick}$/m	0.09	0.36	0.06	0.51
	$M_d$/M-units	0.16	0.40	0.09	0.84
KNN	$z_b$/m	1.57	1.91	3.63	7.24
	$z_{thick}$/m	1.96	2.09	8.44	8.37
	$M_d$/M-units	2.61	2.85	9.57	9.91
BO-KNN	$z_b$/m	1.47	1.85	3.45	5.30
	$z_{thick}$/m	1.76	1.84	5.64	6.03
	$M_d$/M-units	2.46	2.70	6.30	8.80
XGBoost	$z_b$/m	0.93	1.06	2.78	4.00
	$z_{thick}$/m	1.42	1.69	4.07	4.15
	$M_d$/M-units	1.12	1.67	3.55	5.06
BO-XGBoost	$z_b$/m	0.64	0.91	1.16	3.14
	$z_{thick}$/m	0.75	0.90	2.11	2.79
	$M_d$/M-units	0.98	0.97	2.08	2.97

that the BO-RF outperforms the other models, such as RF, KNN, BO-KNN, XGBoost, and BO-XGBoost, for the estimation of SBD, and the rank order of the prediction accuracy is followed by BO-XGBoost and BO-KNN from greatest to least. Overall, the regression models optimized with BO are obviously better than those without BO regardless of whether the 5% Gaussian noise is taken into consideration or not. These results suggest that BO is an effective way to find the best hyper-parameter in the regression model.

Furthermore, we randomly select 100 sets of predicted vectors of the refractivity parameter with 5% Gaussian noise from 2000 sets of predicted results to show the distributions of each parameter in SBD. From Figure 7, Figure 8 and Figure 9, it can be roughly observed that the predicted parameters, such as the thickness of the base layer, the thickness of the trap** layer, and the duct strength obtained with BO-RF, BO-KNN, and BO-XGBoost, are in good agreement with the true values. In order to carefully examine the performance of the three methods, the partial enlarged drawing in each figure is also provided. It is clearly shown that the BO-RF obtains the best performance compared with BO-KNN and BO-XGBoost.

Another analysis focused on the accuracy and stability based on 50 independent runs for the inversion of SBD with/without 5% Gaussian noise is provided. Figure 10 and Figure 11 demonstrate the comparison of the histograms of the estimation results achieved by BO-RF, BO-KNN, and BO-XGBoost. It should be pointed out that the red lines denote the position of the true parameter in SBD. It can be seen that the estimation parameters of BO-RF distribute closely to the true parameter compared with the other two models; these mean that the performance of BO-RF performs better than those of BO-KNN and BO-XGBoost in the case of with/without noise. Additionally, the difference between the true and predicted parameters, named inversion error, for the thickness of the base layer, the thickness of the trap** layer, and the duct strength based on the above 50 independent runs is given in Figure 12, Figure 13 and Figure 14. It is clearly shown that the inversion error of BO-RF is also much lower than those of the two models; these also suggest that the accuracy of BO-RF is superior to those of BO-KNN and BO-XGBoost.

In addition, a comparative study focused on the refractivity profile and propagation loss between the true and predicted parameters in SBD is given in the subsequent section. The average inversion parameters of SBD with BO-RF, BO-KNN, and BO-XGBoost are given in

Table 6. Average inversion results of SBD.

Parameters	BO-RF		BO-KNN		BO-XGBoost
	Without Noise	Noise	Without Noise	Noise	Without Noise	Noise
$z_b$/m	32.28	32.43	33.77	34.26	32.92	33.37
$z_{thick}$/m	40.30	40.54	41.95	42.44	41.10	41.62
$M_d$/M-units	66.03	66.33	67.69	68.11	66.76	67.28

.

$z_b$$z_{thick}$$M_d$Figure 15, Figure 16 and Figure 17 give comparisons of the refractivity profile and propagation loss between the true and predicted parameters. Furthermore, a partial enlarged drawing of each figure is presented to exhibit the details. Compared with BO-KNN and BO-XGBoost, we can conclude that the predicted refractivity profiles and their propagation loss obtained with BO-RF exhibit good agreement with the results simulated with true parameters with/without noise. In addition, we can see from Figure 16 and Figure 17 that there are somewhat large deviations between the propagation loss simulated by the true and predicted parameters. This is due to the fact that a small variation in the index of refraction can cause a relatively large change in propagation loss.

6. Conclusions

The accurate inversion of atmospheric ducts is of great significance to the performance evaluation and design of radar and communication systems. In this paper, a BO-RF model is presented to predict the atmospheric duct parameters. Since the hyper-parameter plays a crucial role in the prediction performance of machine learning, BO is adopted to find the appropriate hyper-parameter in RF, and the results suggest that BO is an effective way to find the best hyper-parameter in RF. In addition, OED is used to sample the combinations of refractivity parameters from a range of hyper-parameter since it is a robust way to effectively handle multi-dimension sampling problems. To test the performance of the proposed BO-RF model, the results obtained by the BO-RF are compared with the results obtained by commonly used methods, such as classical RF, XGBoost with/without BO, and KNN with/without BO, with respect to the coefficient of determination R², MAE, and MSE, as well as inversion error. The statistical and comparative analysis indicates that the BO-RF obtains the best results in the case with/without noise and outperforms the other algorithms adopted in this paper for the estimation of atmospheric duct parameters. Future work will mainly focus on the estimation study based on real data with machine learning and the improvement of the existing refractivity model.