## 3 Methodology
### LEO/GEO-Integrated AOD Fusion Process  
#### 3.2.1 BME Method  
The BME method is based on a knowledge processing framework with cognitive principles that can deal with the general knowledge base G (e.g., stochastic physical laws, empirical relationships, and scientific theories) and a specific knowledge base S (e.g., hard and soft data), which is beyond the reach of spatial statistics and linear geostatistics [57; 58]. Hard data refers to data with accurate measurement or less error, while soft data generally has a certain degree of incompleteness or uncertainty. In this study, AOD observations provided by AERONET sites can be regarded as hard data, while AOD retrieved from satellite sensors could be regarded as soft data. The BME method contains three essential epistemological analysis stages: (1) the prior stage, (2) meta-prior stage, and (3) integration stage [59], which are briefly described as follows.  
In the prior stage, the task is to obtain the joint probability density function (PDF) from G. For a vector of random variables \(x_{map}=[x_{soft1},\ x_{soft2}\ \ldots\ x_{softm},\ x_{k}]\) at points \(x_{soft1}\), \(x_{soft2}\ \ldots\ x_{softm}\) with valid pixels and estimation point \(x_{k}\) with an invalid pixel, the expected information entropy function ([PERSON]'s information criterion) contained in the joint PDF of \(f_{x}(x_{map})\) is defined by Equation (5):  
\[H=-\int dx_{map}f_{x}(x_{map})logf_{x}(x_{map}) \tag{5}\]In general, when entropy is larger, that is, when there is more information, the probability distribution of the targeted random variables is more in tune with the actual situation [60; 61]. The Lagrange multiplier \(\lambda_{\alpha}\) is introduced to maximize entropy under the constraint of a set of functions of statistical spatiotemporal covariance moments (i.e., G), which is expressed as \(g_{\alpha}(x_{map})\) in this study [62]. The constructed function is expressed as follows:  
\[F[x_{map},f_{x}(x_{map}),\lambda_{\alpha}]=-\int dx_{map}f_{x}\big{(}x_{map} \big{)}logf_{x}\big{(}x_{map}\big{)}-\sum_{\alpha=1}^{N}\lambda_{\alpha}\bigg{[} \int g_{\alpha}(x_{map})f_{x}\big{(}x_{map}\big{)}dx_{map}-g_{\alpha}\bigg{]} \tag{6}\]  
where \(\mathrm{g}_{\alpha}\) denotes the mathematical expectation of \(g_{\alpha}(x_{map})\). By setting the first partial derivative with respect to \(x_{map},f_{x}(x_{map})\), and \(\lambda_{\alpha}\) to zero, the joint PDF can then be solved using Equation (7):  
\[f_{x}\big{(}x_{map}\big{)}=\frac{\exp\left[\sum_{\alpha=1}^{N}\lambda_{\alpha }g_{\alpha}(x_{map})\right]}{\int\exp\left[\sum_{\alpha=1}^{N}\lambda_{\alpha }g_{\alpha}(x_{map})\right]dx_{map}} \tag{7}\]  
In the meta-prior stage, the task is to express S appropriately, which refers to the multi-source satellite-based AODs serving as soft data and are expressed as \(x_{soft1}\), \(x_{soft2}\)... \(x_{softn}\) in this study. The construction process is described in detail in the following sections.  
In the final integration stage, the posterior PDF can be obtained after updating the initial prior PDF by integrating G and S through operational Bayesian conditionalization [63], which is expressed as follows:  
\[f\Big{(}x_{k}\Big{|}x_{soft1},x_{soft2}\ldots x_{softm}\Big{)}=\frac{f_{x}\Big{(} x_{soft1},x_{soft2}\ldots x_{softm},x_{k}\Big{)}}{f_{x}\Big{(}x_{soft1},x_{ soft2}\ldots x_{softm}\Big{)}}=\frac{f_{x}\big{(}x_{map}\big{)}}{f_{x}\Big{(}x_{ soft1},x_{soft2}\ldots x_{softm}\Big{)}} \tag{8}\]  
Finally, the nonlinear expected value at estimation point \(x_{k}\) is calculated using Equation (9):  
\[\overline{x_{k}}=\int x_{k}f\Big{(}x_{k}\Big{|}x_{soft1},x_{soft2}\ldots x_{ softm}\Big{)}dx_{k} \tag{9}\]