A Survey of Multimodal Data Fusion in Earth Observation-Remote Sensing

1. Introduction

Definitions of data fusion have long varied across application domains, data type, and even among individual scientists [1,2], challenging the possibility of a single universally accepted definition. The first attempt at a formal definition of data fusion was in 1987 by the Data Fusion Sub-panel of JDL, which defined it with a military undertone, yet from a comprehensive viewpoint as

a process dealing with the association, correlation, and combination of data and information from single and multiple sources to achieve refined position and identity estimates, and complete and timely assessments of situations and threats, and their significance. The process is characterized by continuous refinements of its estimates and assessments, and the evaluation of the need for additional sources, or modification of the process itself, to achieve improved result [3].

This definition has over the years evolved into concision as "the process of combining data to estimate entity states" [4] and while it has been influential in certain domains, various fields have developed their own definitions and frameworks, which may or may not directly build upon the JDL definition. As noted in [5] from the field of classical Metrology, data fusion is defined as the combination of different information sources for new knowledge about physical quantities, events, and situations to achieve higher precision. [6] and [7] restricted the definition to the fusion of sensory data, while [8] and [9] presented it from a generic and holistic perspective as multidisciplinary, and from multiple sensors, information sources, and decision-makers. Although these definitions may be attuned to their respective fields of application, the common underlying tenet is the unification of complementary information from multiple sources into an improved single-view perspective. Similarly, from the Special Interest Group of EARSeL in 1996 emerged a formal definition of data fusion in the domain of EO-RS [10]. Data fusion in this regard was adopted as "a formal framework in which are expressed means and tools for the alliance of data originating from different sources" [11]. This aimed to establish a structured conduct of data fusion in Earth observation remote sensing and has become a widely accepted definition among researchers and practitioners. This review focuses on data fusion in EO-RS, set against a backdrop of increasing preference for remote sensing data in fields even outside the earth sciences for various applications.

The modality of a sensor is informed by internal (technical) specifications such as the resolutions (spatial, spectral, temporal, and radiometric) and imaging mechanism, as well as external factors like observation angle and acquisition time [12]. Typical earth observation remote sensing modalities include Panchromatic, Multi-spectral, Hyper-spectral, Radar, LiDAR, thermal, and microwave imaging. Each of these modalities is designed and tailored for a specific range of the EM spectrum, view angle, and application(s); hence, none can singularly provide a true and accurate approximation of an observed scene [13]. Multimodal data fusion mitigates this constraint by combining the near-precision of various modalities in a complementary fashion for improved data representation.

The early 2010s marked the beginning of a dramatic increase in earth observation satellites, mainly due to improved technology, reduced launch cost, and the increasing demand for high-resolution and real-time data for various fields of application; which includes agriculture, urban planning, environmental monitoring, and disaster management [14,15]. This, in turn, translates to a vast repository of disparate datasets that require data fusion to fully exploit and leverage the complementary aspects of their respective modalities [16,17]. Consequently, several multimodal data fusion techniques have emerged over the past years and have in recent times, gained significant traction, particularly at the turn of the third decade of the 21st century. This is largely driven by the increasing influx of EO data and improved computational power.

To put these fusion techniques into perspective, several taxonomies have been developed since the 1980s. [6] classified fusion techniques into four groups based on the level at which the fusion is performed; they are signal-level, pixel-level, feature-level, and decision-level fusions. [18] incorporated the transforms of disparate spatiotemporal resolutions onto a common grid and developed a three-level taxonomy, which includes sub-feature level, feature level, and decision level fusion techniques. A most recent taxonomy was proposed by [19], which classifies existing techniques based on the underlying mathematical principle. They inferred that data fusion algorithms are fundamentally designed with either of the following: a governing equation that describes the calculation process and geometric features of a surface (Gaussian Process algorithms), a design to reduce noise by assigning weights to outputs (Weighted Least Square Algorithms), or a design to learn patterns and relationships in datasets (Machine learning Algorithms). Some notable fusion techniques include Principal Component Analysis, sparce encoding, Hue Intensity Saturation, gradient descent, Bayesian inference, Wavelet transform, Probabilistic modeling, Simple Block Average, Laplacian Pyramid Decomposition, Discrete Cosine Transform, and Machine Learning (which has significantly contributed to the surge in the development of multimodal fusion methods today) [20,21]. The variations and growing number of classifications can be attributed to the fact that even within individual domains such as EO-RS, there isn’t a universal theoretical basis/framework upon which a comprehensive classification scheme can be developed for multimodal data fusion techniques [2,22].

While existing taxonomies offer valuable perspectives, none consider the underlying analytical paradigm driving fusion methods. All methods, regardless of what defines them or what taxonomy they constitute, are fundamentally model- or data-driven. In addition to an overview of the state of fusion methods, we present a classification based on Machine learning and Traditional methods of fusion in EO-RS, which should provide an understanding and distinction between established techniques and emerging approaches. Traditional methods of data fusion are based on predefined mathematical or statistical models underlined by fixed assumptions [20] as opposed to machine learning models that learn and adapt to their respective input data. We then examine several case studies that illustrate the practical applications of these approaches. But first, we highlight the main groups of pre-processing techniques employed as precursors to data fusion. Although not within the scope of this review, it is important to highlight the major prior operations that are equally fundamental to the fidelity of multimodal data fusion. Regardless of how sophisticated a fusion algorithm is, an inaccurate and inadequate pre-processing (including failure to preserve modal integrity) can significantly compromise results and render it ineffective. This paper is organized as follows: Section 1 introduces the research and provides a concise background on pre-processing techniques essential for fusion accuracy. Section 2 presents a classification based on the analytical paradigm, categorizing techniques into Machine Learning and Traditional approaches. Section 3 explores case studies demonstrating practical applications. Section 4 presents the challenges and limitations of data fusion in EO-RS. Section 5 presents the Discussion and outlook.

1.1. Pre-Processing Techniques

Although recent advancements in satellite and other remote sensing technologies have improved the spatial, temporal, and spectral resolutions of scene imaging, EO data still lack perfect continuity in space and time. This makes the alignment of the resultant disparate observations into a cohesive unit imperative and one that requires efficient strategies to create parity and extrapolate pseudo-continuous scenes, making the data appear more seamless and useful for analysis. Hence, the need and importance of pre-processing. Pre-processing refers to all activities of image rectification and restoration; thus, it aims to create a nearly true representation of an observed scene. It is so crucial that its accuracy is a major factor in determining the quality of the fused image [23]. To the best of our knowledge, very few or none of the publications on EO remote sensing data fusion have highlighted and acknowledged the pivotal role of prior or pre-processing techniques. We present a highlight of the most common groups of these techniques. Before these state-of-the-art data fusion algorithms can be effectively implemented, the tasks of projecting data onto a common grid or geographic space, assimilating observations of the same timeline, matching feature points between multiple images, and atmospheric corrections, among many others, are crucial. These techniques are generally grouped under geometric and radiometric corrections/pre-processing [24].

• Geometric Pre-processing: This refers to all operations aimed at correcting spatial distortions in images and inconsistencies between images due to sensor specifications/mechanics, platform instability or movement, terrain relief, curvature, and rotation of the earth, etc. [25]. Examples include Radial and Tangential distortions, Skew, along-track scan error, scale error, etc. These distortions are unique to each imaging system’s self-geometric characteristics and, hence, require tailored models or approaches for correction. Even among images from the same sensor, a one-model correction solution is not always guaranteed, as sensor properties may change over time due to age, temperature, mechanical, and orbital stress. The sources of geometric distortions in remote sensing imagery can be broadly categorized into two: Observer-related distortions, which stem from the imaging system or sensor, and Observed-related distortions, which arise from the environment or objects being observed (e.g., atmosphere, Earth’s surface) [26]. Observer-related distortions, being predominantly systematic, are relatively easy to correct due to their consistency and predictability. In contrast, errors related to the observed scene can be more challenging to rectify, as they depend on factors such as time, position, and atmospheric conditions at the moment of image acquisition that may not always be readily available [27]. Geometric corrections have become increasingly critical in modern remote sensing due to the shift toward computerized data processing and the growing complexities of multimodal data fusion. Commonly used geometric correction methods include image registration, orthorectification, resampling, and georeferencing.

• Radiometric Pre-processing: For a given pixel, the value recorded is a representation of energy emitted or reflected from both an observed target on the surface and its atmospheric interactions. The ideal imaging system should be able to measure the reflected energy from a scene accurately and uniformly. However, atmospheric perturbations affect the propagation of incident and reflected light, while factors such as platform instability and sensor properties (e.g., viewing geometry, noise, and response drift) further compromise the recorded digital numbers, introducing radiometric inconsistencies that degrade signal accuracy. Radiometric pre-processing aims to correct these errors associated with scene illumination through the process of sensor calibration, atmospheric correction, noise reduction[28], etc, for an improved estimation of an observation. Transforming the TOA spectral radiance to earth surface reflectance is the most crucial step in radiometric correction[29] as without this, signals would primarily reflect unintended observation of atmospheric conditions rather than actual surface characteristics.

For clarity and brevity, this paper uses abbreviations throughout the main text. A complete list of abbreviations and their corresponding full terms is provided in the Abbreviations section.

2. Taxonomy of Multimodal Data Fusion Techniques

2.1. Fusion Level

We highlight an underlying concept of all fusion approaches. Across all disciplines, fusion problems are a priori defined by the questions "when", "where", and "how" to fuse [18,30]. This concept is fundamental, as it forms the basis and building blocks for a wide range of techniques and taxonomies of multimodal data fusion in EO-RS. On this basis, data fusion techniques can generally be grouped under the following themes:

• Signal-level fusion: This aims to create an image with a better signal-to-noise ratio than the source images. Fusion at this level combines the raw electrical signals measured from different sensors to create an enhanced representation of the observed scene.

• Pixel-level fusion: Pixel-level fusion is fusion at the lowest level (pixels) where bits of measured quantities are merged[31]. This approach combines images by applying mathematical operations, such as averaging, to create an enhanced image that would otherwise not be obtained from the individual source images.

• sub-feature level fusion: sub-feature involves fusing different spatial and temporal scales into a common multidimensional grid[18]. Fusion is at a lower level and involves integrating features within the extracted features, such as specific spectral bands or texture patterns from hyperspectral imagery.

• feature-level fusion: This is relatively a high-level fusion technique as opposed to sub-feature level fusion. This involves combining features such as edges, textures, or shapes directly extracted from multiple images into a new and improved representation of all the source images [32].

• decision-level fusion: The decision or results of independent processing algorithms on the disparate modalities are fused using specific techniques at a higher level to produce a final decision[33,34]. This approach aims to improve the accuracy of the final output via a thorough and robust decision-making process.

• Hybrid-level fusion: Combines multiple fusion levels to produce a final output. This technique leverages the complementary advantage of each fusion level to improve the final output.

2.2. Machine Learning Methods

2.2.1. Deep Learning

Deep Learning, a subset of Machine Learning, is computation through multilayer neural networks and processing [35]. As concise as it is, deep learning is in fact a complex network of multilayered artificial neurons that simulate the complex decision-making process of the human brain. This complex architecture allows the capture of intricate and subtle patterns and relationships present in data. Hence, deep learning can provide finer-grained and detailed characterisation of EO imaging, and has garnered increasing attention in the field of EO-RS [36]. Accurate estimation of land use and land cover in an image scene is crucial for precision environmental monitoring and biodiversity conservation [37], agricultural management and urban planning [38], climate change studies [39], and disaster management [40]. [41] in a study in Indonesia developed ForestNet; a semantic segmentation-based Convolutional Neural Network that fuses the visible bands of Landsat 8 imagery and several variables of disparate modalities including NCEP-CFSv2 and USGS-SRTM data, to identify direct drivers of primary forest loss. A technique that can detect multiple clusters, land use, and land cover within the same region to ascertain drivers of the change. ForestNet, unlike canonical classification, takes a semantic segmentation approach, which conducts classification at the pixel level to identify intra-classes of land use or cover classes within a single image. [42] proposed a hybrid CNN which fuses a low resolution 3-meter PlanetScope imagery with a 3cm DJI Phantom 4 Pro UAV imagery to measure the areal extent of mangrove sites in Baja California. Using a modified efficientnet-b0 feature extractor (added a Perceptron), the network was able to combine the coarse optical image with the finer RGB UAV image, which was then fed into a fully connected layer for the task of classification (which could not be otherwise achieved with an efficientnet alone). Active cloud and aerosol profilers lack the wide swath their passive optical counterparts possess; hence, both are often combined for a detailed and yet extensive coverage of atmospheric profiles. A commonly fused data-pair is CALIOP and VIIRS, for which [43] developed an unsupervised DMSDA model to learn feature representations between the two. This approach works by mapping VIIRS into the feature space of the CALIOP cloud profile for a finer cloud characterization within the VIIRS image. This complements VIIRS with improved cloud profiling by providing detailed vertical cloud structure information. [44] developed a fully convolutional Siamese network for multimodal image registration, a generalized approach at the time, where numerous studies focused on specialized networks for registration of unimodal images. Each branch of the network has seven convolutional layers with filters ranging from 32 to 128. Specific to this design, batch normalization is applied only to the seventh layer. This is to preserve information encoded in negative values.

2.2.2. Classical Machine Learning

The backbone of today’s revolutionary sub-field of deep learning is Classical machine learning. It has laid the foundations through algorithms such as decision trees, SVM, and KNN, among other data-driven methods for deep learning today. Yet, Classical Machine Learning continues to excel in EO-RS data modeling. For example, the Random Forest machine learning algorithm, in spite of recent advances in neural networks, continues to achieve great results in land use and land cover classifications [45,46]. For structured and non-complex data patterns and relationships, CML yields excellent results with high efficiency and reduced strain on computing resources, and it remains a preferred algorithm in this field today [47]. [48] developed a gradient-boosting machine learning model, FI-GBM, to estimate the full coverage spatio-temporal distribution of hourly ground-level nitrogen oxide ($N{O}_2$) over Tangshan in China. By combining in-situ data from ground stations with the spatial advantage of TROPOMI and LCSs, the proposed FI-GBM algorithm augments AQM $N{O}_2$ by assimilating TROPOMI and LCS data $N{O}_2$ to infer hourly ground $N{O}_2$ across time and space. A novel approach to integrate machine learning algorithms for the fusion of several multimodal datasets was conducted by [49] for the retrieval of VWC. A feature-level fusion of satellite and reanalysis data by employing five ensemble decision tree algorithms, including Random Forest, Extreme and Light Gradient Boosting Machines, Gradient Boosted Decision Trees, and Bagging Tree algorithms. Although the integrated approach in itself is nothing new, its application for VWC quantitative retrieval modeling is pioneering. [50,51] employed a random forest fusion machine learning strategy to downscale 25km course resolution AMSR-E soil moisture data using 1km MODIS data products (LST, NDVI, EVI, Albedo, LAI, and ET). The latter further extended the resulting AMSR-E product to create a high-resolution soil moisture drought index over Tangshan in China.

The inherent strengths and limitations of individual machine learning algorithms are some of the key drivers for data fusion in EO-RS. Integrating multiple modeling algorithms, particularly when dealing with highly complex datasets, creates a more robust workflow and significantly improves the accuracy of the results [52]. [53] proposed a decision-level machine learning fusion approach to extract urban buildings from Landsat 5, THEOS, and WorldView-3 satellite images. After the respective classifications of buildings using a SVM, KNN, and an ANN on computed spectral indices, the Majority Voting ensemble method was used to find the intersections of propositions from the various classifiers for the final output. Although decision-level fusion offers promising accuracy and robustness, it can be computationally intensive and time-consuming, especially for complex tasks. Therefore, it may not be the most efficient choice for simpler or less data-intensive applications.

2.3. Traditional Methods

2.3.1. Geostatistical Methods

Originally developed for mining problems, Geostatistics provides a set of techniques for estimating a given quantity by exploiting the spatial dependence between two or more statistically heterogeneous datasets. It leverages the existing spatial dependence between a scene (sensor data) and its true representation to estimate the relationship between. Some of the most commonly used methods include kriging-based methods, Bayesian methods, and geostatistical simulation methods. These methods provide the advantage of better spatial autocorrelation, computational efficiency, and better uncertainty quantification, among many others. [54] extended the Fixed-Rank Krieging geostatistical modeling technique for the fusion of AIRS and OCO-2 to estimate carbon dioxide in the PBL. In spite of using synthetic OCO data (due to OCO-2 launch failure coupled with insufficient ground truth data), the estimates returned were accurate representations of the global distribution of PBL $C{0}_2$ except for a few places where synthetic data was sparse. Similarly, [55] combined AERONET, MISR and MODIS datasets for improved estimates of AOT using universal Kriging (which can take multisensor and spatial features as predictors). This was a novel study at the time because the majority of prior studies in this regard only fused passive sensor data, which are characteristically coarse. The estimated AOTs were significant improvements from their original versions, especially in areas with high correlation between the input data. A major challenge with traditional kriging in remote sensing is the intensive computation of inverting large datasets. With a time complexity of $O\left(N^3\right)$, applying this technique to massive matrices is mostly infeasible. However, reducing the dimensions of such matrices promises a way around this conundrum. Given this, [56] proposed an approach of using SSDF that employs spatial random effects modeling (a variant of kriging for dimensionality reduction) to fuse data from AIRS and CrIMSS in the creation of enhanced near-surface air temperature data over the conterminous United States. This generated a refined dataset that exhibited consistent accuracy and errors with NOAA’s Integrated Surface Database.

Moving away from traditional kriging with variograms, [57] employed MPS to integrate ground points with coarse raster elevation data (SRTM and GMTED2010), producing higher a spatial resolution output. Their approach built upon the FILTERSIM MPS algorithm by [58], with modifications to better account for spatial structures in digital elevation model fusion. The modified MPS method demonstrated superior performance compared to traditional and MPS-based kriging across all tested metrics (ME, SDE, RMSE, MAE, and SSIM).

2.3.2. Other Traditional Methods

Traditional methods, also known as Standard methods as indicated earlier, are techniques developed over the years that rely on statistical and mathematical models. Geostatistics, although a traditional method in this sense, is an umbrella term encompassing various methods that focus on modeling the spatial relationships between variables to estimate a quantity. In this context, we refer to other traditional methods as simply non-geostatistical methods of EO-RS data fusion, yet underlined with not-so-flexible mathematical assumptions of a phenomenon or the relationship between quantities. In a passive-active sensor fusion, [59] employed a Wavelet-PCA fusion approach in a comparative study of the performance of pixel and object-based image classification over the alluvial Sanjing Plain in China. Prior performances of the wavelet-PCA fusion method [60] inspired its adoption for the integration of multispectral GF-1 and ALOS-PALSAR imagery for land cover classification. In a similar study, [61] evaluated the accuracy of the Random Forest and Support Vector Machine algorithms for land cover classification by using DWT to fuse HJ1B and ALOS-PALSAR coregistered images over arid and semi-arid landscapes. Although the focus of this work was on the performance of the classification algorithms, the wavelet transform fusion considerably mitigates information loss caused by modality-specific variations in the source images. Wavelet transform continues to demonstrate superior performance in remote sensing image analysis, particularly for classification tasks, by reducing both spectral degradations and spatial distortions in fused images across multiple scales [62,63]. Primarily, estimating AOD via the fusion of disparate remote sensing sources has been with passive-passive sensor fusion. This is due, among other things, to their widespread availability and large coverage, allowing for broader regional aerosol studies. However, this advantage is mired with the presence of cloud, haze, and illumination [64,65], which significantly reduces the spatial completeness of the derived or estimated AODs. In light of this shortcoming, [66] in retrieving AODs over mainland China combined retrievals from MODIS and CALIOP using BME. By combining the principles of maximum entropy and Bayesian inference, BME leverages the use of prior information such as physical laws, empirical models, and expert opinions in fusing the complementary attributes of both inputs for an improved result. The fine-grained CALIOP retrievals combined with MODIS’s wide scanning swath provided detailed AOD characterization and enhanced spatial continuity across large areas, demonstrating clear advantages over passive-sensor-only approaches that lack CALIOP’s active sensing capabilities. Intensity Hue Saturation (IHS) is a traditional fusion technique that is easy to implement but prone to spectral distortion. GTF is an approach proposed by [67] for fusing infrared and gray visible images. GTF avoids designing rules for pixel-level fusion while retaining their radiation and spatial information. To address these shortfalls and exploit the full advantage of both techniques, [68] developed the IHS-GTF hybrid approach for fusing TerraSAR and Sentinel-2A optical data in urban impervious surface extraction over Wuhan, China. The IHS-GTF method significantly outperformed individual techniques (IHS, GTF, and DWT) with R-squared values averaging 0.45 points higher.

A Comprehensive list of these fusion techniques detailing their fusion paradigms, fusion levels and modalites involved can be found in

3. Case Studies

In this Section, we explore the application of multimodal data fusion in EO-RS analysis for effective model performance and improved results. Areas to explore include Aerosol and Cloud detection, land cover change, and Shadow Detection and Image Restoration.

3.1. Aerosol and Cloud Detection

Optical sensors are characteristically designed with a wide scanning swath, making them suited for regional and near-global studies. However, they lack the quality of detailed vertical and horizontal differentiation of clouds and aerosols within an atmospheric column. Radar and LiDAR sensors, on the other hand, provide fine and detailed atmospheric columnar characterizations but lack a wide horizontal scanning range. This trade-off presents an opportunity to harmonize the complementary advantages of both sensors for detailed cloud and aerosol profiles on a wide scale. We present two use cases where multimodal data fusion significantly improves dust-aerosol and cloud machine learning detection algorithms.

Mineral dust aerosols play crucial roles in the bio-geochemical cycle of the Earth. They significantly influence cloud formation, radiative energy transfer, and nutrient supply. In addition, they are critical indicators of air quality and have important implications for human health. Earth observation remote sensing provides a means for monitoring the distribution and transport trajectory of dust aerosols on both local and global scales. Therefore, an efficient and effective dust detection algorithm becomes paramount for reliable retrievals. [69] developed a global dust detection algorithm; an FFNN with three hidden layers to detect dust pixels in cloud-free conditons from daytime VIIRS image. Given the limited scope of many prior studies which often focused on specific regions, restricted dust events, specific viewing angles, and use physical-based methods (eg. threshold and spectral signature-based methods) for model training and evaluation, the main objective of this study is to develop a global detection ML-based algorithm that is independent of study area or dust event. Data for the study is NASA’s SIP global collocated daytime VIIRS and CALIOP product. Although MODIS-CALIOP collocation is very common in this domain, the study opted for VIIRS for two reasons. First, is the extended operational lifetime of VIIRS on JPSS, and second, is the number of diverging VIIRS pixels off CALIOP nadir track. The collocated product from SIPS initially included only VIIRS-on-CALIOP (VIIRS-CALIOP intersections) track pixels. To enhance the robustness of the detection algorithm, four VIIRS pixels adjacent to the CALIOP track were added to the collocated product for improved generalization. To evaluate the algorithm’s performance, two dust aerosol products were employed. The first was the ADP generated using NOAA’s physics-based VIIRS aerosol detection model [70], which identifies dust based on spectral properties in the ultraviolet, blue, and thermal VIIRS bands. The second was the CALIOP aerosol product, which utilizes the V3 algorithm to distinguish between pure dust and polluted dust. Over the ocean, the low performance of FFNN can be attributed to the fact that the model finds it difficult to separate dust from other aerosols. However, it outperforms both ADP and CALIOP by a large margin over land where dust events originate and are predominant. The combined approach leverages the vertical discernment property of CALIOP and the extensive scanning swath of VIIRS for wide area retrievals. Table 1, Table 2 and Table 3 show the detection rate of each approach under clear skies (clouds masked using ADP cloud Mask) over land and the ocean.

The radiative forcings of clouds create feedback uncertainty, which makes it challenging to model cloud properties and actions for accurate climate modeling and projection. In a time where climate change threatens irreversible consequences [71], it is imperative to develop a more comprehensive and robust understanding of cloud dynamics and distribution, as both significantly affect the distribution of heat in the atmosphere [72,73,74]. Cloud profiles are critical for the accurate characterization of cloud albedo and radiative forcings [75], cloud-aerosol interactions [76], climate model validation and evaluation [77] and precipitation formation [78] emphasizing the integral role they play in Earth’s energy balance and climate system. [79] developed a deep CNN that estimates and maps cloud profiles from an active sensor unto a passive sensor product on the same grid. The model learns from labeled narrow nadir cloud pixels and extrapolates to off-nadir unlabeled cloud pixels on several kilometers of swath (approximately 2200km) for the estimation of full vertical cloud profile. Data for the study is the ATCS Dataset; which is a fusion (via alignment) of Polarization and Directionality of the Earth’s Reflectances data (POLDER level 1B) and 2B-CLDCLASS data from cloudSat’s CPR. Each POLDER pixel contains 16 co-registered multi-angle viewpoints from which the observation is made. The model returns a binary cloud/non-cloud binary output, regardless of the cloud profiles preserved in the ATCS, and was assessed using cloud estimations from 2B-CLDCLASS. The study allows for a series of ablation experiments focusing on sensor properties such as view angle, spectral band, and sensor polarimetry to evaluate their respective contributions in identifying cloud pixels. Cloud structure estimations improved with the ablation of two view angles while retaining all spectral bands. The ablation of at least a spectral band reduced model performance, indicating that all spectral bands informed the estimation of cloud structure. The approach achieved a dice score of 73.0% which by current segmentation standards may seem satisfactory at best. However, the absence of a combined human-algorithm annotation approach made accurate and precise pixel prediction difficult since the input image is multimodal, with each constituent having different sensitivities to cloud physics. For instance, the model struggled to identify optically thick clouds and lower layers of multi-layer clouds. Polarimetry had a negligible impact on model performance, responsible for a 0.4% reduction in dice score. This somewhat suggests the Stokes vector parameterization employed may not be ideal for describing the polarization state of reflected energy in this context. Nonetheless, the distribution of identified cloud pixels follows estimations in existing scientific literature and observed atmospheric patterns, with exact markedness of the vertical structure indicating cloud extinction from a height of approximately 14km upward. [80,81,82].

3.2. Land Cover Change

Land cover is the biophysical cover of the earth [83], and it is therefore a direct reflection of shifts in land use, environmental conditions, and human activities on spatio-temporal scales. These shifts, whether driven by natural processes or anthropogenic influences, have far-reaching consequences. Changes in land cover affect climate, the water cycle, biodiversity, urban heat, and food security, among other things. [84,85,86,87,88]. These highlight the importance of monitoring land cover as a key indicator of environmental change, and it reiterates the need for accurate land cover mapping for timely management and conservation efforts [37,89,90]. [36] sort to answer the question, can multimodal data fusion improve classification accuracy? More specifically, can the information gain from the fusion of disparate modalities be linked to improved classification accuracy? To ascertain this, a land cover type classification was conducted over Onondoga County in the central New York region using the Random Forest (RF) ensemble algorithm. Several Random Forest models were developed to assess the impact of different feature combinations on classification accuracy. Models were grouped under single-sensor, two-sensor, and three-sensor feature combinations. Data for the study is coregistered multi-temporal Landsat 5 Thematic Mapper (TM), ALOS-1/PALSAR, and LVIS LiDAR images. Except for the latter, optical and SAR images were sampled during leaf-on, leaf-off, and transition periods in between to evaluate the effect of spectral and scattering variations of some vegetation classes on the task. For each Landsat scene, surface reflectance (SR) for bands 1-5 and 7 was derived, from which the variance for each SR band and the NDVI using the near-infrared (band 5) and red bands (band 3) were computed. ALOS-1/PALSAR features extracted were categorized under intensity, interferometry, polarimetry, and texture, whereas relative height (rh) metrics for all flight lines and measured footprints were derived from LVIS Level 2 Geolocated Surface Elevation and Canopy Height Product (LGE). Table 4 shows the total features extracted from each sensor product.

All model groups showed varying levels of variability in classification accuracies. However, the most significant was observed among single-sensor feature models, with a standard deviation of 9.17, a mean of 61.1% and a statistical range of [74%, 48%], suggesting that this approach rarely makes consistent and reliable approximations of pixels to their true land cover classes. Classifications improved and became more consistent with multi-sensor feature models. Two and three-sensor feature combination models showed minimal variabilities of 4.9 and 3.6 standard deviations away from mean accuracies of 71% and 78%, respectively, indicating pixel assignments and discrimination of different land cover classes improved and remained consistent with the addition of information from multiple sensors. The best performance came from the fusion of LVIS and multitemporal Landsat and PALSAR features, which had an overall accuracy of 83%. It is also worth noting that Landsat-derived features consistently improved model accuracy, regardless of combining them with features from different modalities. This study pioneers the optical-SAR-LiDAR fusion approach in multimodal data fusion, and given the underpinning tenet of this practice, the established positive association between multimodal data fusion and model accuracy in this study does indeed suggest a relationship between the two. Table 5 shows the overall accuracy statistics for each model.

3.3. Shadow Detection and Restoration

Shadow is one of the most common errors in remote sensing data [91]. It occurs when there is an occlusion of direct incident light on a target object [92,93], thereby degrading the quality of imaging, which makes the identification and extraction of useful information very difficult, if not impossible. The financial cost of developing imaging systems, the economic losses incurred by degraded imagery, and the high potential for misinformation from remote sensing image analysis make the detection and restoration of faulty pixels a crucial and urgent task [94,95]. This effort is crucial for enhancing the precision and interpretability of images. The imaging of geologic outcrop formations is often affected by illumination constraints, such as shadows, due to their rugged edges and surfaces. Even HSI, which enables the detailed identification and mapping of minerals and rock types, is insufficient to fully compensate for the loss of information due to the presence of shadow pixels citep9874905. [96] combined SWIR-HSI (970-2500nm) with TLS data for the detection and restoration of shadowed pixels in vertical rock outcrop HSI imagery. HSI was fused with TLS point cloud data for shadow detection via geometric ray tracing. For restoration, TLC intensity information was directly integrated with HSI spectra and then assessed by matching regions in both sunlit and shadowed areas, which is unlike existing methods, which indirectly fuse intensity with the spectra of shadowed pixels. TLS point cloud data was acquired using a 1550nm laser wavelength Reiegl VZ-400 at three scan positions, each about 16m away from the rock face. The points were then coregistered, georeferenced to the WGS84 datum, and transformed to a local geodetic coordinate system in that order to enable shadow computations that utilize local solar zenith and azimuth angles. For occlusion analysis, a triangular mesh was created from the point data to identify locations of shadow pixels in the HSI and ascertain the visibility of each point to the HSI camera. TLS point intensity data were also converted to reflectance for the restoration effort. SWIR-HSI rock outcrop images were captured during morning (SWIR-shade), midday (SWIR-partial), and afternoon (SWIR-sun) hours to reflect different sunlit conditions of the day and as well as for the validation of the restoration technique on the SWIR-partial (which is the image for detection and restoration). All images were coregistered into the pixel space of the designated master image (SWIR-sun). The following corrections/pre-processing were conducted on the HSI image prior to fusion and subsequent analysis.

• Dark current removal.
• Removal of linear image artifact pixels.
• Removal of strong water and carbon dioxide absorption wavelengths and extreme wavelengths in the 896-922nm and 2447-2504nm range of the camera spectrum.
• Conversion of digital numbers (DN) to relative reflectance.

A camera model was developed to facilitate the fusion (registration) of HSI with TLS point-cloud data and intensity information via projection onto a common 2D image space or coordinate system. By first solving the interior (IO) and exterior orientation (EO) parameters of the camera, 3D point-cloud data points were projected from their respective object spaces (x,y,z) onto a 2D image space (x,y). All HSI images were re-projected onto image space following an initial projection into object space by tracing a ray from image pixel and determining its intersection with the outcrop mesh derived from the TLS point cloud. TLS intensity information on the other hand, was rasterized into 8-bit grayscale (active reflectance) images by projecting each point cloud point through the camera model using the IO and EO parameters of the master image. Shadow detection accuracy can be challenging due to the effect of different lightning conditions and their complex interactions with target materials, which sometimes create shadow composites and make edge and boundary detection more difficult [97]. In this study, instead of numerical accuracy metrics, a visual comparison was performed between true shadow pixels (well-defined with distinct boundaries) in the SwIR-shade and the ray tracing detections in the SWIR-partial. By this qualitative approach, an accurate pixel-level differentiation between shadowed and illuminated pixels was established. To achieve direct fusion of TLS information with HSI spectra for pixel restoration, a scale factor is applied to the entire HSI spectrum. This scale factor is calculated from the ratio of active sensor reflectance to HSI reflectance at the active sensor’s operating wavelength. To assess the accuracy of restored pixels, the spectral shape, spectral scale, and band correlation metrics were used to ascertain the preservation of the spectral profile, the effect of the restoration on the magnitude (intensity/brightness) of the signal across the spectra, and the numerical relationship between corresponding bands in the sunlit and shadowed matching material regions, respectively. The approach returned a band correlation of 0.6, about a 0.3 improvement from that of the unrestored SWIR-partial. The Spectral Scale measurements across pixels showed a nearly normal distribution centered at 1 (about 1.2 reduction from the unrestored SWIR-partial), with a variability of 27%. A spectral scale between 0.95 and 1.05 suggests matching materials across the spectrum with very minimal distortions. The statistics from both metrics validate the visual evaluation approach used for shadow detection. Spectral shape metrics for both the restored and original unrestored SWIR-partial images ranged between 0.97 and 1, an indication of quite a preservation of the spectra in both images; nonetheless, the application of the proposed direct method effected no change or improvement. A combination of both the direct and indirect spectra fusion approach did significantly improve the spectral shape of the image, suggesting that indirect shadow restoration, regardless of the computational constraint, does a better job at maintaining spectral profile across pixels. In relation to [98], the unimodal input data approach, coupled with the USRT model employed was able to improve the illumination of shadowed pixels while maintaining the spectral profile of the image. However, via careful visual inspection, it did not achieve the superior perceptual quality of the combined HSI-TLS data approach. This qualitative comparative assessment was chosen because of the differences in evaluation strategies between the two studies.

4. Challenges and Limitations of Data Fusion in Remote Sensing

Heterogeneity is the concept of distinctiveness and peculiarity between entities, which is informed by overt and/or covert qualities fundamental to their respective nature. It is, therefore, only natural, yet an inconvenience, to experience an initial friction or barrier in an attempt to harmonize such disparities. Remote sensing data fusion, although technically feasible and actively practiced, is not seamless by default and, hence, comes with certain challenges which can be attributed to but not limited to platform and sensor design, and externalities such as cost of geodata, limited validation data, computational infrastructure requirements, etc. Needless to say, these constraints can be broadly categorized into internal and External factors.

4.1. Internal Factors

4.1.1. Spectral Resolution

The spectral resolution of a sensor describes its property of detecting radiance from objects in specific wavelengths of the electromagnetic spectrum. It is the sampling rate and bandwidth at which a sensor collects information about a scene [99,100]. Monochromatic sensors detect radiance in a single wavelength band, while multispectral sensors capture data across several discrete spectral bands, typically ranging from 3 to 10 channels [101]. Hyperspectral sensors, as their name suggests, record scene information across hundreds of narrow, contiguous wavelength bands, enabling detailed material characterization and more comprehensive scene analysis through fine spectral resolution [102]. A major challenge in the fusion of mismatching spectral resolutions is spectral distortion [103], which describes alterations or deviations in the true spectral signature or profile of objects in the resulting image. Component Substitution methods (CS) such as IHS, PCA, and Brovey transform introduce spectral distortions when fusing imagery. These techniques operate on the premise that spatial details can be isolated and replaced between images [104]. However, spectral distortions arise from significant differences in spectral profile between the source images when sharp outlines from a single band image with high spatial resolution, e.g. panchromatic image, is injected into a multiband color image with low spatial differentiation and gradual color transitions between feature edges. Speckle noise is another spectral distortion in SAR-optical image fusion. This is caused by the random interference of many elementary reflectors within one resolution cell during the coherent processing of the backscatter signals in SAR imaging [105]. This speckle pattern or artifact can propagate into the fused SAR-optical fusion, bringing distortions into the spectral features of the optical component of the fused product. In spite of recent advances in fusion algorithms in the remote sensing domain, most correction methods like linear filtering, non-linear filtering, and partial differential equation filtering ignore the degradation induced by single bands images such as SAR [106], poorly correct degradations, or trade edge preservation for performance [107].

4.1.2. Spatial Resolution

The ground distance between two adjacent pixels, also known as Ground Sampling Interval, informs the pixel size or the spatial resolution of a terrestrial image [108]. Spatial resolution is therefore the scale or size of the smallest unit of an image or sensor capable of distinguishing objects [109]. One of the most common sources of spatial distortion in the fusion images with discordant spatial resolutions is pan-sharpening, a process that enhances the spatial details of an image using a panchromatic band (with a higher spatial differentiation of features). However, the effect of pan-sharpening is intertwined with the spatial and spectral resolution of the fused product, and as such, the discussion on spectral resolution above gives a glimmer of that. These spatial distortions can manifest in various ways, such as the block effect.

The Block effect is a condition that occurs when pixels of the same land cover class exhibit different reflectance values, making a naturally progressive and continuous phenomenon or class contrastive within by sharp discontinuities between pixels [110]. Block effect is a widely acknowledged problem in spatial unmixing-based methods in remote sensing, and it arises from the use of different local windows in the unmixing of different coarse pixels [111]. Spatial unmixing-based methods are employed when a high spatial resolution image is used to estimate what mix of surface features is likely present in each pixel of a spatially coarse image especially for applications in land cover classification, urban studies, and geological studies. Several studies have attempted to improve the performance of spatial unmixing methods [112,113,114] but have not achieved a significant resolution for blocky artifacts in images. This can be attributed to the following reasons: They are not specifically designed to account for that, they operate on inaccurate assumptions about the temporal dynamics of landcover change, and they ignore the compensation of model residuals [111,115]. The few that significantly minimize the occurrence of blocky artifacts do so at the expense of computational efficiency. Similarly, weight function-based methods, despite receiving considerable attention and extensive development in remote sensing literature [116,117,118,119], also remain susceptible to block artifacts. This vulnerability stems from the fact that they operate on pixel-level processing for the synthesis of fine and coarse-spatial resolution images for spatio-temporal analysis. They also fail to exploit the inherent and valuable shape information of features, which can help identify clusters and maintain the integrity of the natural shapes of features.

4.1.3. Temporal Resolution

Temporal resolution is a measure of the frequency of revisit of a satellite to a scene. More specifically, this is the revisit time of a sensor. The temporal constraint on high spatial resolution imaging, cloudy and shadowy conditions, scanning swath limitations and more specifically, the failure of the scan-line corrector of the enhanced thematic mapper+ on the legacy Landsat 7; suspended on January 19, 2024 [120,121,122], necessitate the fusion of images across different timelines, particularly for spatiotemporal analysis of the environment. Although platforms like Sentinel 2 and 3, and the development of CubeSats have significantly improved mitigating the conundrum of the high spatial and high temporal mutual exclusivity of images [123,124], multitemporal image fusion is still crucial for temporal gap-filling of historical satellite images [125,126,127]. This is, however, not without challenges, as multitemporal images are potentially affected by misregistration errors, varying climate and atmospheric conditions, spectral width and spectral response function, the BRDF effect from changing sensor observation and illumination angles, and phenological changes (which is a major challenge in ecological and urban change detection studies) [128,129]. To minimize these inconsistencies for a somewhat seamless integration, linear transformation models are employed to establish correspondence between these images [130,131,132,133,134], with which deficient and lost pixels can be approximated using derived coefficients [135].

A fundamental misconception underlining many methods of temporal and spatiotemporal fusion is the erroneous assumption of temporal and spatial invariability of surface features over time [114,126,129]. Consequently, these methods are challenged in heterogeneous landscapes and time series imagery with pronounced and abrupt changes in surface cover [128,136]. An example is the Fit-FC method, which blurs fine and sharp edges in highly multi-feature environments [137]. Hence, linear transformations cannot always accurately model the complex spatio-temporal dynamics of real-world scenarios. Non-linear and learning-based modeling approaches that are capable of modeling nuanced and complex relationships are therefore necessary to mitigate these constraints and produce more realistic and physically meaningful fusion results across heterogeneous and rapidly changing landscapes.

4.2. External Factors

4.2.1. Computational Infrastructure Requirements

According to the 2024 Earth Observation Satellite Systems report by Novaspace, EO satellites are expected to triple in number between 2024 and 2033, a 190% increase from the 1,864 in previous decades. This, coupled with the growing usage of UAV, ground instruments, and the continuous improvements in sensor resolutions, means a significant increase in the sheer size (file size and count) of data in an already substantial geospatial data repository [138]. Thus, EO-RS data is considered big data [139] and therefore requires very powerful and highly sophisticated hardware and algorithms for data processing and storage. This data barrage has necessitated the introduction of various processing and storage paradigms as well as the democratization of access to state-of-the-art computational resources [140] such as cloud-based computing (GEE, Google Collab, Microsoft Azure, AWS) and HPC. In spite of these and other existing solutions, the high dimensionality of geodata [141], the high bandwidth demands for data transfer [139,142,143], the intensive energy and data requirements along with specialized hardware demands of AI-driven algorithms, and the cost and geographical marginalization of existing infrastructure challenge the wide-spread adoption and effective utilization of these resources [144,145,146,147,148,149].

4.2.2. Cost of Geodata

The free and open data policy adopted by USGS and ESA for Landsat and Sentinel data products, respectively, has contributed significantly to the use of geospatial data in a wide range of environmental applications [150,151]. However, real-time and very high-precision environmental monitoring and analysis require very high-resolution (spatial, spectral, and temporal) data, which can be very expensive to acquire, especially in remote regions and the developing world for large-scale applications [152]. Although high spatial resolution aerial imagery such as that from the NAIP can be freely available, they are geographically restricted to the United States and severely lack the swath and spectral complements of space-borne imaging systems, which is crucial for large-scale detailed object-property characterization. Commercial data products such as those from PlanetScope and UAVs tend to have higher temporal and spatial resolutions [153], and are therefore often combined with open-access wide-swath products for high spatiotemporal analysis [154]. Furthermore, ground data is also crucial for model and data validation, data gap-filling, and scaling from point measurements to broader areas. However, the high cost of acquisition and maintenance, along with the challenge of achieving adequate spatial distribution of ground stations and equipment, present major barriers to accessing these niche datasets and their eventual integration with disparate modalities [155,156,157,158,159].

4.2.3. Limited Validation Data

Validation data in EO-RS takes several forms, including high-resolution raster, in situ ground measurements, model outputs, and reanalysis data. Validation, as defined by the CEOS WGCV, is the process of independently assessing the quality of derived data products [160], a crucial step in geoanalysis, as proper validation establishes rigor and confidence in remote sensing products and algorithms. Validation assesses the viability of multimodal data fusion and ensures reliable estimations of geophysical quantities through the quantification and propagation of uncertainties. Yet, this indispensable component is challenged by limited available data [161,162] despite mitigation efforts, especially for large-scale analysis. This is due to factors such as cost, terrain fragility and surface heterogeneity which limit data acquisition, untrustworthy reference dataset [163], and spatial scale mismatch between ground-based measurements and remote sensing products [164]. Some suggested and applied mitigation strategies include fusion of multiple validation data sources [165], establishing more in situ ground stations, improving ground stations to expand observation scale, improvements in upscaling transfer functions [163], and resampling [166]. Although these solutions offer potential pathways to address data scarcity, their successful implementation requires careful consideration of their respective inherent challenges. For instance, while resampling interpolates new pixel values, it assumes spatial continuity of features, which can potentially introduce errors that may affect accurate feature representation across heterogeneous landscapes.

Other challenges of multimodal data fusion include selection of appropriate algorithm [167,168], lack of standardized evaluation methods [169,170], misinterpretation of fused product [171], unfamiliarity with remote sensing data and methods especially among non-technical and untrained professionals, and on-demand coverage or real-time data [172].

5. Discussion

In this paper, we have undertaken a comprehensive exploration of data fusion: 1. Its ubiquity and how various fields and their respective applications have shaped its definitions 2). Fusion methods informed by the underlying analytical paradigms 3). Case Studies and 4). Challenges. We have also demonstrated how multimodal data fusion enables comprehensive object characterization by amplifying the information necessary for high-precision environmental analysis. As indicated earlier, improvements in sensor designs and specifications, coupled with the exponential growth in Earth observation data have catalyzed novel analytical techniques and driven significant improvements in existing infrastructure. These have widened the scope of data fusion beyond traditional modalities, allowing for more comprehensive and robust analysis with depth. However, these innovations, as groundbreaking and complementary as they are (especially with advances in fusion techniques, computational and storage infrastructure), either present new sets of challenges or compound existing ones, effectively rendering some solutions potentially inapplicable or highly marginalized to certain domains.

Besides computational inefficiency, many challenges in fusion algorithms stem from very little or no consideration of context and integration of domain knowledge in the learning or modeling process. Context-aware and domain-specific methods, although constrained to particular applications as their names suggest, prioritize precision within their respective sub-domains rather than making granular estimates across multiple fields. Understanding the nature, governing principles and physical laws, and the underlying dynamics within the focus of investigation considerably enhances modeling effectiveness and learning outcomes. Typical examples of contextual and domain-aware approaches are PIML and PINN, which maintain fidelity with physical laws by embedding them directly into data and the learning architecture. This provides a robust, efficient, and improved modeling accuracy in addressing various challenges across settings from mechanics to geophysics and even in the bio-chemical sciences [173,174]. And although working principles such as spatial continuity implemented through advanced fusion techniques like object-based unmixing methods and neural networks (e.g., CNNs) have improved capturing the spatial relationships between features, their accuracies are still constrained by feature heterogeneity and sharp discontinuities [175,176]. The emergence of vision transformers [177] with their self-attention mechanism effectively identify spatial nonstationarity between adjacent pixels in images, which is a limitation of CNNs in remote sensing image segmentation. This capability enables transformers to more effectively mitigate the mixed-pixel problem, leading to more accurate segmentation in heterogeneous landscapes.

While these technical advances in fusion algorithms have improved modeling accuracies, quantifying their comparative effectiveness remains a challenge due to the lack of a standard evaluation framework. Our analysis highlights this problem, particularly in the context of shadow detection and restoration, where metrics remain inconsistent across multiple studies. Shadow detection evaluation typically employs metrics such as precision, recall, F1-score, and balanced error rate, along with domain-specific measures like shadow detection rate and intersection over union. For shadow restoration, image quality metrics such as Peak Signal-to-Noise Ratio, Structural Similarity Index, and mean squared error, as well as spectral fidelity metrics of spectral shape and scale metrics for the restored regions are used. These variations make it more challenging to quantitatively compare the efficacy of fusion techniques and determine their appropriateness under specific conditions. Hence, the qualitative approach of visual inspection has assumed the de facto method when comparing accuracies between studies with different evaluation metrics. While it excels at perceptual evaluation, such as detecting unrealistic textures or contextually inappropriate spectral characteristics, visual inspection is highly subjective and influenced by context and observer perception. This evaluation inconsistency is so pervasive that it persists not only within shadow detection and restoration but also across multiple application domains, such as cloud detection. This fragmented evaluation landscape, as a result, not only hampers cross-study comparisons but also impedes the identification of optimal fusion strategies for specific use cases. This review strongly advocates for prioritizing future research in developing a unified benchmarking framework that establishes standardized evaluation protocols and metrics for assessing both the efficiency and efficacy of fusion algorithms across diverse application domains. Such a framework would enhance interdisciplinary research and collaboration while enabling robust comparative analysis.