Irradiance Restoration Based Shadow Compensation Approach for High Resolution Multispectral Satellite Remote Sensing Images

Numerous applications are hindered by shadows in high resolution satellite remote sensing images, like image classification, target recognition and change detection. In order to improve remote sensing image utilization, significant importance appears for restoring surface feature information under shadow regions. Problems inevitably occur for current shadow compensation methods in processing high resolution multispectral satellite remote sensing images, such as color distortion of compensated shadow and interference of non-shadow. In this study, to further settle these problems, we analyzed the surface irradiance of both shadow and non-shadow areas based on a satellite sensor imaging mechanism and radiative transfer theory, and finally develop an irradiance restoration based (IRB) shadow compensation approach under the assumption that the shadow area owns the same irradiance to the nearby non-shadow area containing the same type features. To validate the performance of the proposed IRB approach for shadow compensation, we tested numerous images of WorldView-2 and WorldView-3 acquired at different sites and times. We particularly evaluated the shadow compensation performance of the proposed IRB approach by qualitative visual sense comparison and quantitative assessment with two WorldView-3 test images of Tripoli, Libya. The resulting images automatically produced by our IRB method deliver a good visual sense and relatively low relative root mean square error (rRMSE) values. Experimental results show that the proposed IRB shadow compensation approach can not only compensate information of surface features in shadow areas both effectively and automatically, but can also well preserve information of objects in non-shadow regions for high resolution multispectral satellite remote sensing images.


Introduction
Optical satellite remote sensing images play an increasingly important role in a number of applications like building reconstruction, height evaluation, resource assessment, change detection and precision agriculture, along with the increasingly spatial resolution of optical satellite remote sensing images provided by high spatial resolution (HSR) satellites, such as QuickBird, WorldView-2, WorldView-3, and Jilin-1 [1][2][3][4][5][6][7][8]. However, shadows inevitably formed by clouds and land surface features interfere more seriously in these HSR optical satellite remote sensing images [9]. Notably, the shadow often weakens feature information due to the absence of direct light compared with the corresponding non-shadow [9][10][11]. Besides, significant information not only exists in most non-shadow ρ is the surface reflectance of the non-shadow area. Specifically, the surface irradiance (E nshw ) of the non-shadow aera includes three parts: direct irradiance (E d ), scattered irradiance (E s ) (also named diffuse irradiance), and ambient (or reflected) irradiance (E a ). However, the surface irradiance (E shw ) of the shadow area mainly consists of E s and E a , because E d is often inevitably occluded by clouds and surface features for the shadow area. The surface irradiance compositions of the non-shadow area and the shadow area are expressed in the following equations [25,35,38]: Generally, the radiant existence (M nshw ) of the Earth's surface in the non-shadow area is expressed with the corresponding E nshw : where ρ nshw is the surface reflectance of the non-shadow area.
In addition, the radiance from the Earth's surface of the non-shadow area (L earth nshw ) can be expressed with the corresponding radiant existence (M nshw ) of the Earth surface in the non-shadow area: where π is the solid angle. Therefore, L earth nshw can be directly described with the corresponding E nshw by substituting Equation (3) for M nshw in Equation (4), as follows [39]: Our purpose is to restore information of shadow regions in the acquired images. Therefore, we mainly concentrate on the light propagation path from the Earth surface to the satellite sensor. In the light propagation process, atmospheric influence should be taken into consideration, because the radiance of land-surface features propagates through the atmosphere before arriving at the satellite sensor. The atmospheric influence mainly results in an attenuated radiance from the Earth surface, and a path radiance (L p ) that is directly reflected to the satellite sensor by the atmosphere, but does not arrive at the Earth surface [40]. Hence, the total radiance of the non-shadow area arriving at the satellite sensor (L sensor nshw ) can be expressed with the attenuated L earth nshw and L p : L sensor nshw = τL earth nshw + L p (6) where τ is the atmospheric transmittance from the Earth surface to the satellite sensor.
Here, the image radiance of the non-shadow area is calculated in accordance with Equations (1) and (5) as the following equation: Similarly, the image radiance of the shadow area is also expressed in a similar way: L sensor shw = τρ shw π E shw + L p = τρ shw π (E s + E a ) + L p (8) where ρ shw is the surface reflectance of the shadow area.
Thus, image radiance can be described in a general way combining Equations (7) and (8), as follows [11]: where ρ is the surface reflectance, and k is the shadow approximation parameter, in which k = 0 refers to shadow and k = 1 denotes non-shadow. Additionally, the absence of E d in the shadow area contributes to the difference between non-shadow and shadow areas. For retrieving information of the shadow area containing the same type features to the nearby non-shadow area, an assumption can be arranged in which the same surface irradiance of the non-shadow area is owned by the nearby shadow area. Namely, the shadow area receives E s , E a , and a hypothesized direct irradiance (E hyp d ) under the assumption above. In this case, the corresponding hypothesized radiance of the shadow area (L hyp shw ) can be considered as the compensated radiance, which is received by the satellite sensor hypothetically but not actually, as shown in Equation (10).   In the shadow compensation process, because E hyp d is developed under the assumption that the shadow area contains the same type features as the nearby non-shadow area, E hyp d can be regarded approximately equal to E d in the corresponding nearby non-shadow area (i.e., E hyp d = E d ). Therefore, Equation (10) can be rewritten as in the following: where r is the irradiance coefficient defined by Equation (12).
In order to calculate r with image radiance, L p is first moved from the right-hand side to the left-hand side in both Equations (7) and (8), then items of Equation (8) are divided by the corresponding items of Equation (7), as shown in Equation (13): In addition, the surface reflectance of the non-shadow area is reasonably equal to that of the nearby shadow area containing the same type features (i.e., ρ nshw = ρ shw ) [9]. Thus, Equation (13) can be further rewritten as Equation (14): Hence, the irradiance coefficient r can be obtained in the image radiance from with the following equation: For clarity, the compensated image radiance of the shadow area (L comp shw ) can be represented by the following equation: Note that Equation (16) only needs to estimate L p avoiding the uncertainty of the estimation of parameters τ and ρ. Thus, the solution of Equation (16) is promising.
For pixels of both shadow and non-shadow areas in the target image, Equation (16) can be further expressed as follows:

Workflow of the IRB Approach
In our study, the target image is first divided into shadow and non-shadow regions with the logarithmic shadow index (LSI) shadow detection method and further refined with a certain manual operation. Consequently, shadow regions are compensated with the proposed IRB method pixel-by-pixel and non-shadow regions are left unchanged. The developed IRB approach is mainly accomplished in five steps described as in the following, as shown in Figure 2.

Workflow of the IRB Approach
In our study, the target image is first divided into shadow and non-shadow regions with the logarithmic shadow index (LSI) shadow detection method and further refined with a certain manual operation. Consequently, shadow regions are compensated with the proposed IRB method pixel-by-pixel and non-shadow regions are left unchanged. The developed IRB approach is mainly accomplished in five steps described as in the following, as shown in Figure 2. Workflow chart for the IRB approach.

•
Step 1: Radiance calibration Optical satellite remote sensing images are usually saved in a digital number (DN) value form. For instance, the Worldview-3 data is saved in the 11-bit DN value form [41]. Radiance calibration is often applied prior to shadow correction, in which the image DN value is converted into radiance pixel by pixel with the image metadata in the acquired image files through a certain calibration algorithm or directly over a professional software. In this study, the radiance calibration process is carried out with the radiometric calibration module in the ENVI 5.2.

•
Step 2: Shadow detection Workflow chart for the IRB approach.

•
Step 1: Radiance calibration Optical satellite remote sensing images are usually saved in a digital number (DN) value form. For instance, the Worldview-3 data is saved in the 11-bit DN value form [41]. Radiance calibration is often applied prior to shadow correction, in which the image DN value is converted into radiance pixel by pixel with the image metadata in the acquired image files through a certain calibration algorithm or directly over a professional software. In this study, the radiance calibration process is carried out with the radiometric calibration module in the ENVI 5.2.

•
Step 2: Shadow detection In the shadow correction process, shadow detection is a pre-process prior to shadow compensation, which extracts shadow regions in the target image and classifies the target image into non-shadow and shadow parts [42][43][44][45]. In this step, the LSI shadow detection method [46] is first utilized to produce the shadow mask image. Problems of non-shadow misclassification and shadow omission are still inevitable for the current shadow detection process. Consequently, we manually refine the shadow detection result of the LSI shadow detection method to further alleviate shadow detection problems.

•
Step 3: Path radiance estimation The path radiance L p is an important element for the whole shadow compensation process, which can be estimated with various algorithms, like the dark object subtraction (DOS), the histogram method, and the improved dark object subtraction (IDOS) [47,48]. In this study, we estimate the L p value per band of the target image with the IDOS algorithm.

•
Step 4: Irradiance coefficient computation As shown in Equation (15), the irradiance coefficient r is derived with L sensor nshw , L sensor shw and L p . As mentioned in Step 3 above, L p values are estimated with the IDOS algorithm. In this step, Sensors 2020, 20, 6053 7 of 23 we compute L sensor nshw and L sensor shw with the Minkowski norm on the basis of the color constancy theory [33,34], respectively, as shown in Equation (18): where a is the gain factor, p is an adjustable parameter, f (x, y) is the radiance value of the target image with the coordination (x, y), and M × N refers to the amount of pixels in either the non-shadow area or the shadow area of the target image.

•
Step 5: Shadow compensation and optimization The shadow compensation task is initially accomplished with the shadow detection result (Step 2), L p (Step 3) and r (Step 4) in accordance with the solution of Equation (17).
For pixels belonging to the shadow area (i.e., k = 0), refining parameters α and β are additionally applied to optimize the compensation result by the solution of Equation (17). The optimized compensated radiance L opt shw is shown with the following equation: where α and β are refining parameters. As shown in Equation (19), refining parameters α and β are additionally applied to optimize the shadow compensation result by the initial solution of Equation (17), which contributes to the refined resulting image after shadow compensation. Particularly, parameter α is mainly employed to optimize the impact of the original shadow on the compensation result, and parameter β is utilized to refine the influence of the path radiance and the irradiance coefficient on the compensation result. However, it is unrealistic to form an explicit mathematical model for the refining parameters used in the image interpretation of high resolution remote sensing images, because complex and diverse ground features are usually caught in these images. Therefore, we refer to the experimental solution with similar optimization for defining the refining parameters for shadow detection of high resolution remote sensing images by Zhang et al. [49]. We finally decide to choose the experimental solution to defining the refining parameters α and β. As for the restriction of the refining parameters α and β, the additional experimental results reveals that the restriction of α + β is both simple and valid for optimizing the compensation result, although there are many ways to restrict the refining parameters α and β. With these considerations, we further explore the restriction of α + β by ranging α + β from 1 to 6 in experiments with both the quantitative measurement, relative root mean square error (rRMSE), and the visual sense of compensation results. We finally determine to restrict the refining parameters α and β with α + β = 3. Subsequently, under the restriction of α + β = 3, we traverse the refining parameters α and β in the interval of 0.1 from 0.1 to 3, and define values of α and β. A related discussion and analysis on the parameter sensitivity of compensation results are provided in additional experiments in Section 4.

Test Images
The proposed IRB shadow compensation approach was accomplished on a DELL personal computer with 64-bit Windows 7 operation system equipped with a 3.2 GHz CPU and a 4 GB RAM. To verify the shadow compensation performance of the proposed IRB approach, we run comparative experiments with many images from WorldView-3 of Tripoli, Libya, and Rio de Janeiro, Brazil, and WorldView-2 of Washington DC, USA, which are respectively called WV3-Tripoli, WV3-Rio, and WV2-WDC. The corresponding discussion is provided in the next section (Section 4: Discussion). Particularly, qualitative and quantitative assessments are performed in this section to evaluate the shadow compensation performance of the proposed IRB approach against several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM) with two test images from WorldView-3 of Tripoli, Libya [41], as shown in Figure 3a,b (respectively called Tripoli-1 and Tripoli-2). Specifically, the test image Tripoli-1 in Figure 3a is a 400 × 300 pixel image, which covers typical ground objects, such as shadow, various scale urban buildings, asphalt roads, bare land, and grass. The test image Tripoli-2 in Figure 3b is a 260 × 195 pixel image mainly consisting of shadow, buildings, asphalt roads, grass, playground, and parks.

Test Images
The proposed IRB shadow compensation approach was accomplished on a DELL personal computer with 64-bit Windows 7 operation system equipped with a 3.2 GHz CPU and a 4 GB RAM. To verify the shadow compensation performance of the proposed IRB approach, we run comparative experiments with many images from WorldView-3 of Tripoli, Libya, and Rio de Janeiro, Brazil, and WorldView-2 of Washington DC, USA, which are respectively called WV3-Tripoli, WV3-Rio, and WV2-WDC. The corresponding discussion is provided in the next section (Section 4: Discussion). Particularly, qualitative and quantitative assessments are performed in this section to evaluate the shadow compensation performance of the proposed IRB approach against several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM) with two test images from WorldView-3 of Tripoli, Libya [41], as shown in Figure 3a,b (respectively called Tripoli-1 and Tripoli-2). Specifically, the test image Tripoli-1 in Figure 3a is a 400 × 300 pixel image, which covers typical ground objects, such as shadow, various scale urban buildings, asphalt roads, bare land, and grass. The test image Tripoli-2 in Figure 3b is a 260 × 195 pixel image mainly consisting of shadow, buildings, asphalt roads, grass, playground, and parks.

Qualititave Evaluation
For evaluating shadow compensation results, the subjective visual sense evaluation method intuitively reflects the compensation effect of specific regions in target images and is conducive to analyzing the performance among various shadow compensation methods for a specific shadow [14,16]. In this part, the shadow compensation effect is first analyzed for test images Tripoli-1 and Tripoli-2 from the perspective of qualitative visual sense comparison. Shadow compensation results of test images by the proposed IRB shadow compensation approach are respectively discussed in terms of R, G and B components. Then, a comprehensive visual sense comparison is also made between shadow compensation results by various shadow compensation algorithms (i.e., IRB, LCC, LSCR, MSR, HF, and DELM). Figure 4 lists the resulting images of the test image Tripoli-1 compensated before and after by the presented IRB approach in terms of R, G, and B components. As shown in Figure 4b sense improvement, which is consistent with the compensation result of R component in Figure 4b, when compared with those in Figure 4c,e. Consequently, not only is obvious improvement in visual sense acquired by the IRB shadow compensation approach for the test image Tripoli-1, but good consistency is also shown among the compensation results of R, G, and B components.
Sensors 2020, 20, x FOR PEER REVIEW 10 of 24   Figure 5a,b, the information of typical ground targets is clearly restored by the IRB method in Figure 5b, like grass under building shadow in region E, asphalt roads in region F, tops of buildings in region G, and   Figure 5a,b, the information of typical ground targets is clearly restored by the IRB method in Figure 5b, like grass under building shadow in region E, asphalt roads in region F, tops of buildings in region G, and grass under small shadow in region H, which are mostly covered by shadow in Figure 5a. Therefore, it can be observed that the R component of Tripoli-2 is well compensated by the proposed IRB method, and the visual sense of the compensated R component has been significantly improved compared with the original one. Analogous to the compensation result of the R component of Tripoli-2, G and B components also acquire a dramatical improvement in terms of the gray value of shadow regions, as shown in Figure 5d,f. Moreover, the compensation results of R, G, and B components are relatively consistent. Hence, not only is meaningful improvement of the visual sense achieved by the IRB shadow compensation approach, but strong consistency is also shown in compensation results in terms of R, G, and B components of Tripoli-2.
Sensors 2020, 20, x FOR PEER REVIEW 11 of 24 grass under small shadow in region H, which are mostly covered by shadow in Figure 5a. Therefore, it can be observed that the R component of Tripoli-2 is well compensated by the proposed IRB method, and the visual sense of the compensated R component has been significantly improved compared with the original one. Analogous to the compensation result of the R component of Tripoli-2, G and B components also acquire a dramatical improvement in terms of the gray value of shadow regions, as shown in Figure 5d,f. Moreover, the compensation results of R, G, and B components are relatively consistent. Hence, not only is meaningful improvement of the visual sense achieved by the IRB shadow compensation approach, but strong consistency is also shown in compensation results in terms of R, G, and B components of Tripoli-2. In addition to visual sense comparison regarding compensation results by the IRB shadow compensation approach in terms of R, G, and B components of test images Tripoli-1 and Tripoli-2, In addition to visual sense comparison regarding compensation results by the IRB shadow compensation approach in terms of R, G, and B components of test images Tripoli-1 and Tripoli-2, related analysis is also carried out to compare compensation results by the IRB approach and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). Figure 6 illustrates the shadow compensation results of the test image Tripoli-1 by the IRB approach and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). As presented in Figure 6a, obviously, the resulting image by the IRB method reveals more details about ground targets under shadow and significant improvement is also shown in a visual sense by comparing the shadow compensation result in Figure 6a with the original test image Tripoli-1 in Figure 3a. A similar color effect of shadow regions is also achieved in the compensation result to that of the adjacent similar objects, such as small buildings and grass in region A, blue houses in region B, and nearby shadow of buildings in regions C and D shown in Figure 6a. As illustrated in Figure 6b, the shadow compensation result by LCC also improves the visibility of objects in the shadow regions. Similarly, the shadow compensation result by LSCR and DELM in Figure 6c,f also shows great improvement of shadow areas in visual sense, even though the improvement effect is not as obvious as the results by IRB and LCC listed in Figure 6a,b. However, the shadow effect is still obvious in the resulting image by the MSR shadow compensation method shown in Figure 6, like almost similar shadow regions A-D appearing in Figure 3a. Additionally, color distortion and non-shadow interference occur obviously in the resulting image by HF, as shown in Figure 6e. Generally speaking, the proposed IRB algorithm well restores the information of objects in shadow regions of the test image Tripoli-1, and achieves a relatively good visual sense compared with the comparative shadow compensation methods.
Sensors 2020, 20, x FOR PEER REVIEW 12 of 24 related analysis is also carried out to compare compensation results by the IRB approach and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). Figure 6 illustrates the shadow compensation results of the test image Tripoli-1 by the IRB approach and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). As presented in Figure 6a, obviously, the resulting image by the IRB method reveals more details about ground targets under shadow and significant improvement is also shown in a visual sense by comparing the shadow compensation result in Figure 6a with the original test image Tripoli-1 in Figure 3a. A similar color effect of shadow regions is also achieved in the compensation result to that of the adjacent similar objects, such as small buildings and grass in region A, blue houses in region B, and nearby shadow of buildings in regions C and D shown in Figure 6a. As illustrated in Figure 6b, the shadow compensation result by LCC also improves the visibility of objects in the shadow regions. Similarly, the shadow compensation result by LSCR and DELM in Figure 6c,f also shows great improvement of shadow areas in visual sense, even though the improvement effect is not as obvious as the results by IRB and LCC listed in Figure 6a,b. However, the shadow effect is still obvious in the resulting image by the MSR shadow compensation method shown in Figure 6, like almost similar shadow regions A-D appearing in Figure 3a. Additionally, color distortion and non-shadow interference occur obviously in the resulting image by HF, as shown in Figure 6e. Generally speaking, the proposed IRB algorithm well restores the information of objects in shadow regions of the test image Tripoli-1, and achieves a relatively good visual sense compared with the comparative shadow compensation methods. The linear correlation correction (LCC) method [16,21]. (c) The light source color rate (LSCR) method [33]. (d) The multi-scale Retinex (MSR) method [24,36]. (e) The homomorphic filtering (HF) method [25]. (f) The direct and environment light based method (DELM) [11,37].  [16,21]. (c) The light source color rate (LSCR) method [33]. (d) The multi-scale Retinex (MSR) method [24,36]. (e) The homomorphic filtering (HF) method [25]. (f) The direct and environment light based method (DELM) [11,37].
Additionally, Figure 7 also depicts shadow compensation results of the test image Tripoli-2 by the IRB approach and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). Similar to the discussion of Tripoli-1, IRB, LCC, LSCR, and DELM also perform well in shadow compensation of Tripoli-2 and achieve a dramatical improvement of shadow regions in visual sense. For instance, as shown in Figure 7a-c and f, most ground objects under shadow are well restored, like grass under building shadow in region E, asphalt roads in region F, tops of buildings in region G, and grass under small shadow in region H. However, as presented in Figure 7d, shadow still exists obviously in the resulting image by MSR, like large building shadow in regions E, F and G, and small house shadow in region H. In addition, color distortion and non-shadow interference still appear in the resulting image by HF shown in Figure 7e. In general, obvious improvement of shadow in the test image Tripoli-2 in terms of visual sense is achieved in resulting images by IRB, LCC, LSCR, and DELM.
All in all, by comprehensively comparing shadow compensation results of test images Tripoli-1 and Tripoli-2 by IRB and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM), it can be seen that the proposed IRB shadow compensation approach can well restore features of shadow regions and achieve a relatively good overall visual sense. the test image Tripoli-1. (a) The irradiance restoration based (IRB) method. (b) The linear correlation correction (LCC) method [16,21]. (c) The light source color rate (LSCR) method [33]. (d) The multi-scale Retinex (MSR) method [24,36]. (e) The homomorphic filtering (HF) method [25]. (f) The direct and environment light based method (DELM) [11,37]. All in all, by comprehensively comparing shadow compensation results of test images Tripoli-1 and Tripoli-2 by IRB and several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM), it can be seen that the proposed IRB shadow compensation approach can well restore features of shadow regions and achieve a relatively good overall visual sense.

Quantitative Assessment
In addition to the qualitative evaluation above, the quantitative assessment method is also an effective assessment way. Meanwhile, the quantitative assessment method is also used to evaluate the performance of shadow compensation methods for further quantitatively assessing the shadow compensation results. The root mean square error (RMSE) is a widely used indicator for quantitative assessment [9]. However, considering the difference of the absolute RMSE values for shadow compensation results in terms of R, G, and B components, the relative RMSE (rRMSE) is more reasonable instead of the absolute RMSE to assess shadow compensation results in terms of R, G, and B components. Accordingly, the rRMSE is defined by the following equation [9]:

Quantitative Assessment
In addition to the qualitative evaluation above, the quantitative assessment method is also an effective assessment way. Meanwhile, the quantitative assessment method is also used to evaluate the performance of shadow compensation methods for further quantitatively assessing the shadow compensation results. The root mean square error (RMSE) is a widely used indicator for quantitative assessment [9]. However, considering the difference of the absolute RMSE values for shadow compensation results in terms of R, G, and B components, the relative RMSE (rRMSE) is more reasonable instead of the absolute RMSE to assess shadow compensation results in terms of R, G, and B components. Accordingly, the rRMSE is defined by the following equation [9]: where f result (x, y) and f re f (x, y) are radiance values of the shadow compensation resulting image and the corresponding reference image with the coordination (x, y), and M × N refers to the amount of pixels. Specifically, the performance of the IRB approach is evaluated by comparing rRMSE of the resulting image against the rRMSE of the original image and the rRMSE of resulting images by several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM), where a smaller rRMSE value delivers a better shadow compensation performance. In the following section, rRMSE is mainly utilized to quantify the performance of the proposed IRB shadow compensation approach for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components. Figure 8a,b depict rRMSE of shadow compensation results for test images compensated before and after by various shadow compensation algorithms (i.e., IRB, LCC, LSCR, MSR, HF, and DELM) in terms of R, G, and B components, respectively. As shown in Figure 8a,b, obviously, rRMSE values of resulting images by MSR and HF are relatively high, which shows the poor shadow compensation ability of MSR and HF. However, relatively small rRMSE values of resulting images are acquired by IRB and LCC compared with the rRMSE of the original images for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, which reveals the good shadow compensation performance of IRB and LCC. Similarly, small rRMSE values of resulting images are also achieved by LCSR and DELM compared with the rRMSE of the original images. Since the rRMSE values of resulting images by IRB and LCC are smaller than those by LCSR and DELM, IRB and LCC can deliver a better shadow compensation performance. Additionally, Figure 9a,b illustrate ∆ rRMSE between non-shadow regions in the original image and those in the shadow compensation resulting images by various shadow compensation algorithms for test images in terms of R, G, and B components. As depicted in Figure 9a,b, ∆RMSE values are more than zero for resulting images by MSR and HF, which reveals that feature information in non-shadow regions are corrupted to some extent. However, ∆RMSE values are zero or approximately zero for resulting images by IRB, LCC, LSCR, and DELM, which delivers that IRB, LCC, LSCR, and DELM can well preserve feature information in non-shadow regions, while compensating for shadow. Additionally, Figure 9a,b illustrate ∆rRMSE between non-shadow regions in the original image and those in the shadow compensation resulting images by various shadow compensation algorithms for test images in terms of R, G, and B components. As depicted in Figure 9a,b, ∆RMSE values are more than zero for resulting images by MSR and HF, which reveals that feature information in non-shadow regions are corrupted to some extent. However, ∆RMSE values are zero or approximately zero for resulting images by IRB, LCC, LSCR, and DELM, which delivers that IRB, LCC, LSCR, and DELM can well preserve feature information in non-shadow regions, while compensating for shadow.
Generally speaking, the proposed IRB method not only further alleviates the color distortion problem compared with typical image enhancement based shadow compensation methods (like MSR and HF), but also acquires a relatively good visual sense similar to that by LCC. Besides, the relatively small rRMSE and zero or approximately zero ∆RMSE also reveal the relatively good performance of the proposed IRB method for shadow compensation and non-shadow preservation. Better agreement is also shown for R,G, and B components of results in terms of both visual sense and quantitative measurement (rRMSE) by the proposed IRB method against the compared methods.
Additionally, Figure 9a,b illustrate ∆ rRMSE between non-shadow regions in the original image and those in the shadow compensation resulting images by various shadow compensation algorithms for test images in terms of R, G, and B components. As depicted in Figure 9a,b, ∆RMSE values are more than zero for resulting images by MSR and HF, which reveals that feature information in non-shadow regions are corrupted to some extent. However, ∆RMSE values are zero or approximately zero for resulting images by IRB, LCC, LSCR, and DELM, which delivers that IRB, LCC, LSCR, and DELM can well preserve feature information in non-shadow regions, while compensating for shadow. Generally speaking, the proposed IRB method not only further alleviates the color distortion problem compared with typical image enhancement based shadow compensation methods (like MSR and HF), but also acquires a relatively good visual sense similar to that by LCC. Besides, the relatively small rRMSE and zero or approximately zero ∆ RMSE also reveal the relatively good performance of the proposed IRB method for shadow compensation and non-shadow preservation. Better agreement is also shown for R,G, and B components of results in terms of both visual sense and quantitative measurement (rRMSE) by the proposed IRB method against the compared methods.

Discussion
The presented IRB algorithm achieves a good shadow compensation performance in the comparative experiments above with test images Tripoli-1 and Tripoli-2. However, as described in the workflow of the proposed IRB algorithm in Section 2, uncertainties may appear in Step 3 (path radiance estimation) and Step 4 (irradiance coefficient computation). Meanwhile, in Step 5 (shadow compensation), the shadow compensation result is also sensitive to refining parameters α and β. In this section, a corresponding discussion is given to evaluate the effect of these elements. Thus, additional experiments are performed to discuss potential influence elements of the IRB approach on the shadow compensation results with test images Tripoli-1 and Tripoli-2.

Influence Analysis of Path Radiance Estimation
The IDOS algorithm is utilized to estimate the path radiance value for each band of the target image. The relative scattering model and the starting haze value (SHV) band, which are used in the estimation of the path radiance with the IDOS algorithm, usually maintain the estimation result. Available bands of WorldView-3 [41] and regular five relative scattering models [47] are listed in Tables 1 and 2, respectively. Because the test image is captured under sunny and cloudless atmospheric conditions, the corresponding atmospheric condition can be regarded as very clear. Namely, the λ −4 relative scattering model [47] is applied during the estimation of path radiance. Hence, properly selecting the SHV band becomes the main influential element in the path radiance estimation. Near-IR2 857-1039 Table 2. Five relative scattering models [47].

Atmospheric Conditions Relative Scattering Model
Very clear Consequently, in order to study the influence of the path radiance estimation on the shadow compensation results by the IRB algorithm, the available bands of test images (i.e., B1, B2, B3, B4, B5, B6, B7, and B8) are respectively set as the SHV band in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 10a,b illustrate the rRMSE of shadow compensation results with various SHV bands for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, respectively. As illustrated in Figure 10a,b, it can be observed that the SHV band slightly impacts the shadow compensation performance of the IRB algorithm. Specifically, when B1, B2, or B3 is set as the SHV band for the IRB algorithm, relatively lower rRMSE values of R, G, and B components are almost acquired at the same time for test images Tripoli-1 and Tripoli-2, which shows that a better shadow compensation effect is achieved for test images Tripoli-1 and Tripoli-2 with B1, B2, or B3 as the SHV band. Accordingly, in this study, B2 is finally utilized as the optimized SHV band for test images.

Influence Analysis of Irradiance Coefficient Computation
The Minkowski norm is utilized in the irradiance coefficient computation together with the path radiance. As for the calculation of the Minkowski norm, related studies [33,34] noted that uncertainties of shadow correction effect in natural and aerial images are mainly sensitive to the parameter p. In order to further explore the influence of the parameter p on the shadow compensation effect of HSR optical satellite remote sensing images, the parameter p is respectively set as from 1 to 30 with an integral interval of 1 in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 11a,b respectively present the rRMSE of shadow compensation results with various p values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components. As illustrated in Figure 10a,b, it can be observed that the SHV band slightly impacts the shadow compensation performance of the IRB algorithm. Specifically, when B1, B2, or B3 is set as the SHV band for the IRB algorithm, relatively lower rRMSE values of R, G, and B components are almost acquired at the same time for test images Tripoli-1 and Tripoli-2, which shows that a better shadow compensation effect is achieved for test images Tripoli-1 and Tripoli-2 with B1, B2, or B3 as the SHV band. Accordingly, in this study, B2 is finally utilized as the optimized SHV band for test images.

Influence Analysis of Irradiance Coefficient Computation
The Minkowski norm is utilized in the irradiance coefficient computation together with the path radiance. As for the calculation of the Minkowski norm, related studies [33,34] noted that uncertainties of shadow correction effect in natural and aerial images are mainly sensitive to the parameter p. In order to further explore the influence of the parameter p on the shadow compensation effect of HSR optical satellite remote sensing images, the parameter p is respectively set as from 1 to 30 with an integral interval of 1 in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 11a,b respectively present the rRMSE of shadow compensation results with various p values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components.
As can be observed in Figure 11a,b, obviously, relatively lower rRMSE values are acquired for R, G, and B components with p value of 4, 5, 6, or 7 for test images Tripoli-1 and Tripoli-2. Additionally, the rRMSE curves in Figure 11a,b show an increasing trend for test images Tripoli-1 and Triploli-2 with the p value from 1 to 30. Consequently, in this study, we handle test images with parameter p value of 5.
parameter p. In order to further explore the influence of the parameter p on the shadow compensation effect of HSR optical satellite remote sensing images, the parameter p is respectively set as from 1 to 30 with an integral interval of 1 in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 11a,b respectively present the rRMSE of shadow compensation results with various p values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components. As can be observed in Figure 11a,b, obviously, relatively lower rRMSE values are acquired for R, G, and B components with p value of 4, 5, 6, or 7 for test images Tripoli-1 and Tripoli-2. Additionally, the rRMSE curves in Figure 11a,b show an increasing trend for test images Tripoli-1 and Triploli-2 with the p value from 1 to 30. Consequently, in this study, we handle test images with parameter p value of 5.

Sensitivity Analysis of Refining Parameters
In this study, the initially derived shadow compensation solution is further optimized with parameters α and β. According to the optimized solution given in Equation (19), the resulting images are sensitive to refining parameters α and β. In particular, the refining parameter α optimizes the partial impact of the original shadow on the compensation results, and β mainly affects the influence of both the path radiance and irradiance coefficient on the compensation results. In this study, the relationship of α and β is roughly restricted with α + β. Subsequently, the influence of different α + β is analyzed by compensating test images with various α + β ranging from 1 to 6. Figure 12a,b illustrate rRMSE of shadow compensation results with various α + β values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components. As presented in Figure 12a,b, relatively low rRMSE values are acquired for test images Tripoli-1 and Tripoli-2 when α + β = 2 or α + β = 3. In this study, refining parameters α and β are restricted with α + β = 3.

Sensitivity Analysis of Refining Parameters
In this study, the initially derived shadow compensation solution is further optimized with parameters α and β . According to the optimized solution given in Equation (19), the resulting images are sensitive to refining parameters α and β . In particular, the refining parameter α optimizes the partial impact of the original shadow on the compensation results, and β mainly affects the influence of both the path radiance and irradiance coefficient on the compensation results. In this study, the relationship of α and β is roughly restricted with α β + .
Subsequently, the influence of different α β + is analyzed by compensating test images with various α β + ranging from 1 to 6. Figure 12a Additionally, the sensitivity of refining parameters α and β are further analyzed by valuing α and β from 0.1 to 3 with an interval of 0.1 in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 13a,b depict rRMSE of shadow compensation results with various α and β values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, respectively. As presented in Figure 13a  Additionally, the sensitivity of refining parameters α and β are further analyzed by valuing α and β from 0.1 to 3 with an interval of 0.1 in the additional experiments with test images Tripoli-1 and Tripoli-2. Figure 13a,b depict rRMSE of shadow compensation results with various α and β values for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, respectively. As presented in Figure 13a,b, lower rRMSE values are acquired for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, when parameters (α, β) are set as (2.4, 0.6), (2.5, 0.5), (2.6, 0.4), (2.7, 0.3), or (2.8, 0.2), which shows that better shadow compensation results are achieved with these parameter values of (α, β). Hence, in this study, we develop the IRB shadow compensation approach over test images with α = 2.6 and β = 0.4. components, respectively. As presented in Figure 13a,b, lower rRMSE values are acquired for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components, when parameters ( α , β ) are set as (2.4, 0.6), (2.5, 0.5), (2.6, 0.4), (2.7, 0.3), or (2.8, 0.2), which shows that better shadow compensation results are achieved with these parameter values of ( α , β ). Hence, in this study, we develop the IRB shadow compensation approach over test images with 2.6 α = and 0.4 (a) (b)

Shadow Difference Analysis
In order to discuss shadow difference between different shadow areas in the same image, we particularly labelled several shadow regions with different shadows and analyzed the compensation results with rRMSE for test images Tripoli-1 and Tripoli-2. Specifically, as shown in Figure 3a

Shadow Difference Analysis
In order to discuss shadow difference between different shadow areas in the same image, we particularly labelled several shadow regions with different shadows and analyzed the compensation results with rRMSE for test images Tripoli-1 and Tripoli-2. Specifically, as shown in Figure 3a As shown in Figures 14a-c and 15a-c, smaller rRMSE values of various shadow areas of the resulting images are acquired for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components compared with the corresponding original ones, which reveals the fact that the IRB

Shadow Difference Analysis
In order to discuss shadow difference between different shadow areas in the same image, we particularly labelled several shadow regions with different shadows and analyzed the compensation results with rRMSE for test images Tripoli-1 and Tripoli-2. Specifically, as shown in Figure 3a, regions A-D are labelled in the test image Tripoli-1 (i.e., small buildings and grass in region A, blue houses in region B, and nearby shadow of buildings in regions C and D). Similarly, regions E-H are also labelled in the test image Tripoli-2 shown in Figure 3b (i.e., grass under building shadow in region E, asphalt roads in region F, tops of buildings in region G, and grass under small shadow in region H). Figures 14a-c and 15a-c respectively depict the rRMSE of various shadow areas of images compensated before and after by the IRB approach for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components. As shown in Figures 14a-c and 15a-c, smaller rRMSE values of various shadow areas of the resulting images are acquired for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components compared with the corresponding original ones, which reveals the fact that the IRB approach can well compensate feature information of shadow regions. As for the distinction of the  Figures 14a-c and 15a-c, smaller rRMSE values of various shadow areas of the resulting images are acquired for test images Tripoli-1 and Tripoli-2 in terms of R, G, and B components compared with the corresponding original ones, which reveals the fact that the IRB approach can well compensate feature information of shadow regions. As for the distinction of the rRMSE value for different shadow regions (like regions A-D in the test image Tripoli-1 and regions E-H in the test image Tripoli-2), the reflectivity of different ground features is the main reason. Therefore, external work will be carried out in the future for exploring the reflectivity of typical ground features.

IRB Method Generalization Analysis
As described in Section 2, the IRB approach is evaluated with many images (i.e., WV3-Tripoli, WV3-Rio, and WV2-WDC) to validate the IRB method generalization. Figure 16a-c respectively illustrate rRMSE values of both the original WV3-Tripoli images and the shadow compensated WV3-Tripoli images in terms of R, G, and B components. Similarly, Figures 17a-c and 18a-c also depict the corresponding rRMSE values for WV3-Rio and WV2-WDC images compensated before and after by the IRB approach in terms of R, G, and B components.
As can be observed in Figure 16a-c, smaller rRMSE values are achieved by the resulting WV3-Tripoli images in terms of R, G, and B components compared with the original WV3-Tripoli images. Figures 17a-c and 18a-c also show similar phenomena for WV3-Rio and WV2-WDC images, which together deliver the IRB approach and is able to well compensate shadow information for HSR multispectral satellite remote sensing images. Provided this situation, two test images of WV3-Tripoli are particularly selected in this paper to assessing the IRB approach both qualitatively and quantitively against several shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM), as previously described in Section 2.
Sensors 2020, 20, x FOR PEER REVIEW 20 of 24 As described in Section 2, the IRB approach is evaluated with many images (i.e., WV3-Tripoli, WV3-Rio, and WV2-WDC) to validate the IRB method generalization.   Sensors 2020, 20, x FOR PEER REVIEW 20 of 24 As described in Section 2, the IRB approach is evaluated with many images (i.e., WV3-Tripoli, WV3-Rio, and WV2-WDC) to validate the IRB method generalization.    As can be observed in Figure 16a-c, smaller rRMSE values are achieved by the resulting WV3-Tripoli images in terms of R, G, and B components compared with the original WV3-Tripoli images. Figures 17a-c and 18a-c also show similar phenomena for WV3-Rio and WV2-WDC images, which together deliver the IRB approach and is able to well compensate shadow information for HSR multispectral satellite remote sensing images. Provided this situation, two test images of WV3-Tripoli are particularly selected in this paper to assessing the IRB approach both

Conclusions
In this paper, we developed and validated a simple but effective shadow compensation approach for high spatial resolution optical multispectral satellite remote sensing images on the basis of the satellite sensor imaging mechanism and radiative transfer theory. According to the distinct surface irradiance difference between non-shadow and nearby shadow areas, an assumption is first arranged for shadow areas under which shadow areas acquire the absent direct irradiance compared to the nearby non-shadow areas containing the same type features. Then, the solution for shadow compensation is further developed by estimating the path radiance and computing irradiance coefficient. Finally, refining operations are applied to optimize the initial shadow compensation solution. Additionally, comparative experiments are carried out to validate the proposed irradiance restoration based (IRB) shadow compensation. Shadow compensation performance of the proposed IRB method is assessed by both qualitative visual sense comparison and quantitative analysis against shadow compensation results by other comparative shadow compensation methods (i.e., LCC, LSCR, MSR, HF, and DELM). The shadow compensation results and the corresponding small rRMSE of the resulting images together reveal that the proposed IRB approach not only preserves information of features in non-shadow areas, but also performs well in restoring information of objects in shadow regions. Additionally, the proposed IRB method is more reasonable because it is developed on the basis of the satellite sensor imaging mechanism and radiative transfer theory, rather than only image processing (like LCC, MSR, and HF). The result images are more suitable for the subsequent image interpretation of high-resolution remote sensing images (such as change detection and surface feature reflectivity reversion). In the future, we will attempt to improve shadow compensation performance by further settling shadow detection problems and considering more influential details for the path radiance estimation combined with our current study.