Spatio-Temporal Error Sources Analysis and Accuracy Improvement in Landsat 8 Image Ground Displacement Measurements

Because of the advantages of low cost, large coverage and short revisit cycle, Landsat 8 images have been widely applied to monitor earth surface movements. However, there are few systematic studies considering the error source characteristics or the improvement of the deformation field accuracy obtained by Landsat 8 image. In this study, we utilize the 2013 Mw 7.7 Balochistan, Pakistan earthquake to analyze error spatio-temporal characteristics and elaborate how to mitigate error sources in the deformation field extracted from multi-temporal Landsat 8 images. We found that the stripe artifacts and the topographic shadowing artifacts are two major error components in the deformation field, which currently lack overall understanding and an effective mitigation strategy. For the stripe artifacts, we propose a small spatial baseline (<200 m) method to avoid the stripe artifacts effect on the deformation field. We also propose a small radiometric baseline method to reduce the topographic shadowing artifacts and radiometric decorrelation noises. Those performances and accuracy evaluation show that these two methods are effective in improving the precision of deformation field. This study provides the possibility to detect subtle ground movement with higher precision caused by earthquake, melting glaciers, landslides, etc., with Landsat 8 images. It is also a good reference for error source analysis and corrections in deformation field extracted from other optical satellite images.


Introduction
With the continuous improvements of optical satellites in both hardware and software, moderate and high-resolution optical images are increasingly applied in ground surface displacement measurement [1].Bindschadler et al. [2] took the lead in extracting the velocity field of an Antarctic ice stream by matching multi-temporal Landsat 5 TM images.Crippen et al. [3] preliminary investigated sand dunes migration in the southernmost California, USA with SPOT (Systeme Probatoire d'Observation de la Tarre) images.The image matching theory was subsequently expanded as feature tracking in SAR (Synthetic Aperture Radar) imagery [4][5][6] and firstly applied to obtain coseismic deformation of the 1992 Mw 7.3 Landers earthquake [4].Puymbroeck et al. [7] captured surface ruptures and horizontal coseismic deformation field of the 1992 Mw 7.3 Landers earthquake using the SPOT images with cross-correlation technique.After the software Coregistration of Optically Sensed Images and Correlation (COSI-Corr) had been released, the theory and algorithms of the optical images cross-correlation technique were greatly improved and matured, especially in accuracy and computation efficiency [8].The COSI-Corr is originally used to extract crustal deformation induced by earthquake released by Tectonics Observatory at Caltech in 2008 [8,9].The technique is now widely used to monitor the earth surface changes and dynamics [8][9][10][11][12][13][14][15].As richer and constantly updated satellite observation data become available, the application of optical images cross-correlation technique has been expanded to fields such as glacier flow velocity extraction [10][11][12], quantification of sand dune migration [13][14][15], terrain-deformation measurements of slow landslides [16], monitoring process of significant rift events [17], volcanic monitoring [18] and particularly extensive coseismic deformation extraction [19][20][21].
As a continuation of the NASA Landsat program, the Landsat 8 satellite was launched in February 2013 and is operating in a good condition now.Equipped with OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) sensors, Landsat 8 can obtain image data of 11 bands, wherein the panchromatic band (Band 8) has a ground resolution of 15 m, equivalent to ASTER (Advanced Spaceborne Thermal Emission And Reflection Radiometer) VNIR bands but with a larger spatial coverage [22][23][24][25].Since its operation, the Landsat 8 data have been applied to monitor earth surface movements.Many studies demonstrate that the Landsat 8 image cross-correlation technique has advantages in measuring surface deformation, such as mapping ice sheet velocity [26], coseismic deformation extraction [25,[27][28][29] and volcano deformation monitoring [30].Many researchers have done a lot of analyses and error source correction of displacement measurement, but they mainly focused on a few data sources, such as SPOT, ASTER, HiRISE (High Resolution Imaging Science Experiment) and aerial imagery [3,8,11,[13][14][15][16][17][19][20][21]31,32].In general, the deformation field obtained by different types of optical sensors with cross-correlation technique has disparities in error structure and magnitude level.As one of the few free available in-orbit operational data sources, Landsat 8 imagery has a relatively higher resolution, high imaging quality (12 bit) and short revisit cycle (16 days).The Landsat 8 surpasses most satellites in both scientific performances and economic cost.No doubt, earth surface deformation monitoring based on the Landsat 8 image will be more widely used in the future.Currently, there are few detailed studies on the error sources of deformation fields obtained by the Landsat 8 image.Thus, it is necessary and important to systematically study the error source components, characteristics and corresponding correction methods of the Landsat 8 image deformation field.
Based on this, we firstly analyze the error source compositions, distribution patterns and magnitude variation of the Landsat 8 image deformation field.Then, we give comprehensive explanations on the reason of error generation and propose corresponding correction methods.We also compare the results obtained by our method and traditional method, and validate the results' precision.Finally, we discuss the potential application of our proposed methods in time series analysis with the Landsat 8 image.

Experimental Design and Data Processing
Balochistan, Pakistan is selected as the study area because it has dry climate and little cloud cover, which is favorable for optical images acquisition (see Figure 1).In addition, a Mw 7.7 earthquake happened there on 23 September 2013, leaving obvious surface ruptures [25,[27][28][29].Totally, 25 images covering the study area (Path: 154, Row: 42, see Figure 1, Table A1) are selected, pairing freely with a temporal baseline of 16 days (see Figure A1) and calculating mutually image deformation fields.
Considering the computational efficiency, we uses a phase correlation algorithm in frequency domain in the experiments, which can quickly obtain the translation between correlated images and easily separate frequency-dependent noises in images.The Fast Fourier Transformation (FFT) correlation engine embedded in the algorithm of COSI-Corr [8][9][10][11][12][13][14][15] improves the traditional peak correlation algorithm, and its accuracy can reach 1/50 pixel [8].Thus, we select COSI-Corr as the data processing platform.We set both the initial and final search window size as 32 pixel × 32 pixel, step as 4 pixel × 4 pixel and the corresponding ground resolution is 60 m × 60 m.For mitigating the effects of high-frequency noises on the accuracy of correlation results, we set the frequency mask threshold as 0.9 and robust iteration as 2. To maintain a good consistency of results correlated by different temporal images, these configuration parameters are used in cross-correlation calculation proceeding for the whole study.processing platform.We set both the initial and final search window size as 32 pixel × 32 pixel, step as 4 pixel × 4 pixel and the corresponding ground resolution is 60 m × 60 m.For mitigating the effects of high-frequency noises on the accuracy of correlation results, we set the frequency mask threshold as 0.9 and robust iteration as 2. To maintain a good consistency of results correlated by different temporal images, these configuration parameters are used in cross-correlation calculation proceeding for the whole study.In order to obtain sufficient correlation samples and weaken the long-time baseline influence, we select 16 days as the temporal baseline to pair freely the 25 images (Path: 154, Row: 42) in Table A1.To maintain the uniform standard of temporal baseline (16 days) of pair images, we acquire redundant observations containing enough variables, such as different orbit overlap and shadowing differences.This can facilitate the analysis of error source including stripe artifacts and topography-dependent artifacts relevant to these variables.We obtain 16 pairs of image deformation fields totally (East-West and North-South) with a ground resolution of 60 m × 60 m, as shown in Figure 2. It should be noted that these results have not been post-processed, so there are a variety of error sources in deformation fields.In order to obtain sufficient correlation samples and weaken the long-time baseline influence, we select 16 days as the temporal baseline to pair freely the 25 images (Path: 154, Row: 42) in Table A1.To maintain the uniform standard of temporal baseline (16 days) of pair images, we acquire redundant observations containing enough variables, such as different orbit overlap and shadowing differences.This can facilitate the analysis of error source including stripe artifacts and topography-dependent artifacts relevant to these variables.We obtain 16 pairs of image deformation fields totally (East-West and North-South) with a ground resolution of 60 m × 60 m, as shown in Figure 2. It should be noted that these results have not been post-processed, so there are a variety of error sources in deformation fields.

Error Source Analysis and Correction Methods
As shown in Figure 2, by analyzing results in both E/W and N/S components, we can find that there are following types of error sources affecting the measurement accuracy, including temporal decorrelation noises, long wavelength orbital error, stripe artifacts (SA), attitude jitter distortions and topography-dependent artifacts.Figure 3 shows patterns of these errors, but excludes the attitude jitter distortions with a small magnitude.In the following, we will analyze these error sources, and propose some mitigation methods.

Error Source Analysis and Correction Methods
As shown in Figure 2, by analyzing results in both E/W and N/S components, we can find that there are following types of error sources affecting the measurement accuracy, including temporal decorrelation noises, long wavelength orbital error, stripe artifacts (SA), attitude jitter distortions and topography-dependent artifacts.Figure 3 shows patterns of these errors, but excludes the attitude jitter distortions with a small magnitude.In the following, we will analyze these error sources, and propose some mitigation methods.

Error Source Analysis and Correction Methods
As shown in Figure 2, by analyzing results in both E/W and N/S components, we can find that there are following types of error sources affecting the measurement accuracy, including temporal decorrelation noises, long wavelength orbital error, stripe artifacts (SA), attitude jitter distortions and topography-dependent artifacts.Figure 3 shows patterns of these errors, but excludes the attitude jitter distortions with a small magnitude.In the following, we will analyze these error sources, and propose some mitigation methods.

Decorrelation Noises
Decorrelation noises are closely correlated to the radiometric properties of earth surface.Too high, too low or invariant radiation can lead to inconspicuous texture features, and drastic radiometric changes commonly can cause texture defect.Thus, it is difficult to achieve the best correlation based on window matching, showing inaccurate random measurements and low signal-to-noise ratio (SNR) values in corresponding regions.To reduce the effect of decorrelation noises, this study masks out measurements with low SNR (<0.9) and null values (Nan) in image deformation fields, and manually removes the large range of decorrelation noises caused by clouds, water regions, etc.

Long Wavelength Orbital Error
If Landsat 8 images downloaded from the USGS website are not orthorectified enough, they will cause significant signals of linear ramp in image deformation fields.In the image pre-processing, this type of error can be eliminated by rigorous geometric correction to original images if the state vector of the Landsat 8 image is available.In this study, we remove this error using the first-order polynomial fitness [25] because Landsat 8 images do not provide the state vector in these products.The results are detailed in Figure 4.The polynomial curve fitting method is also used to remove the orbital error in InSAR (Interferometric Synthetic Aperture Radar) interferogram [33][34][35][36][37][38][39][40][41][42][43][44].

Decorrelation Noises
Decorrelation noises are closely correlated to the radiometric properties of earth surface.Too high, too low or invariant radiation can lead to inconspicuous texture features, and drastic radiometric changes commonly can cause texture defect.Thus, it is difficult to achieve the best correlation based on window matching, showing inaccurate random measurements and low signalto-noise ratio (SNR) values in corresponding regions.To reduce the effect of decorrelation noises, this study masks out measurements with low SNR (<0.9) and null values (Nan) in image deformation fields, and manually removes the large range of decorrelation noises caused by clouds, water regions, etc.

Long Wavelength Orbital Error
If Landsat 8 images downloaded from the USGS website are not orthorectified enough, they will cause significant signals of linear ramp in image deformation fields.In the image pre-processing, this type of error can be eliminated by rigorous geometric correction to original images if the state vector of the Landsat 8 image is available.In this study, we remove this error using the first-order polynomial fitness [25] because Landsat 8 images do not provide the state vector in these products.The results are detailed in Figure 4.The polynomial curve fitting method is also used to remove the orbital error in InSAR (Interferometric Synthetic Aperture Radar) interferogram [33][34][35][36][37][38][39][40][41][42][43][44].

Stripe Artifacts
Misalignment of Charge Coupled Device (CCD) arrays is a common problem in push-broom satellites with mosaic CCD imaging cameras.The general orthorectification cannot completely eliminate the CCD distortions, which perform as stripe artifacts (SA) in image deformation fields [31].In Figure 5, there are many obvious evenly spaced bands along the flight direction in the crustal deformation field.We find that those SA are correlated with the corresponding spatial positions of CCD line-arrays imaging of pre-and post-event images.A method called "mean subtracting" was proposed to measure the patterns of SA in deformation fields [8,11,31].This method can partly estimate and remove the SA using the mean value along the satellite flight direction in the stable area away from tectonic signals and serious decorrelation regions.

Stripe Artifacts
Misalignment of Charge Coupled Device (CCD) arrays is a common problem in push-broom satellites with mosaic CCD imaging cameras.The general orthorectification cannot completely eliminate the CCD distortions, which perform as stripe artifacts (SA) in image deformation fields [31].In Figure 5, there are many obvious evenly spaced bands along the flight direction in the crustal deformation field.We find that those SA are correlated with the corresponding spatial positions of CCD line-arrays imaging of pre-and post-event images.A method called "mean subtracting" was proposed to measure the patterns of SA in deformation fields [8,11,31].This method can partly estimate and remove the SA using the mean value along the satellite flight direction in the stable area away from tectonic signals and serious decorrelation regions.In this study, we try to understand the relationship between the SA and CCD arrays' relative locations between correlated images.We acquire patterns of the SA for those 16 pairs of images in Figure A2.It can be found that most patterns of the SA jag significantly, and only three pairs present approximate straight lines with the mean value closing to zero.From Table 1, we find that those three pairs are different from the rest pairs in the spatial distances between the correlated image centers.Spatial distances of these three pairs are shorter than 180 m, while that of other pairs are up to several thousands (see Table 1).Here we use "spatial baseline" to define this spatial distance between two correlated images' centers, which can be calculated based on the orthoimages centers' coordinates in master and salve images.To quantify the magnitude level of the SA in different image pairs, we define a value γ as a factor for quantifying the average magnitude of the SA.
where, as shown in Figure A2, i represents the transverse-axis parameter with a total number of n, and T represents the measurement values on the vertical-axis in a single pattern of the SA.Using these data, we perform a regression analysis to examine the relationship between the spatial baseline and γ in Figure 6.We find that variant γ and spatial baseline have an approximately positive linear correlation, meaning the smaller the spatial baseline, the lower the magnitude level of the SA.In this study, we try to understand the relationship between the SA and CCD arrays' relative locations between correlated images.We acquire patterns of the SA for those 16 pairs of images in Figure A2.It can be found that most patterns of the SA jag significantly, and only three pairs ( 1, we find that those three pairs are different from the rest pairs in the spatial distances between the correlated image centers.Spatial distances of these three pairs are shorter than 180 m, while that of other pairs are up to several thousands (see Table 1).Here we use "spatial baseline" to define this spatial distance between two correlated images' centers, which can be calculated based on the orthoimages centers' coordinates in master and salve images.To quantify the magnitude level of the SA in different image pairs, we define a value γ as a factor for quantifying the average magnitude of the SA.
where, as shown in Figure A2, i represents the transverse-axis parameter with a total number of n, and T represents the measurement values on the vertical-axis in a single pattern of the SA.Using these data, we perform a regression analysis to examine the relationship between the spatial baseline and γ in Figure 6.We find that variant γ and spatial baseline have an approximately positive linear correlation, meaning the smaller the spatial baseline, the lower the magnitude level of the SA.In fact, the SA are common in horizontal deformation fields extracted from the current mainstream optical images (e.g., ASTER, SPOT, and Landsat 8).The "mean subtracting" method is a conventional approach to mitigate similar linear artifacts [11,25,27].In this study, we proposed a new method to reduce the SA by selecting correlation image pairs with small spatial baseline.In order to compare the performance of these two methods, we select two image pairs, 21 May 2013-11 May 2015 with a small spatial baseline (<100 m) and 10-26 September 2013 with a large spatial baseline 2324 m, for mapping coseismic deformation of this event.As shown in Figure 7, it has no significant SA in the E/W and N/S components of crustal deformation field produced by image pair 21 May 2013-11 In fact, the SA are common in horizontal deformation fields extracted from the current mainstream optical images (e.g., ASTER, SPOT, and Landsat 8).The "mean subtracting" method is a conventional approach to mitigate similar linear artifacts [11,25,27].In this study, we proposed a new method to reduce the SA by selecting correlation image pairs with small spatial baseline.In order to compare the performance of these two methods, we select two image pairs, 21 May 2013-11 May 2015 with a small spatial baseline (<100 m) and 10-26 September 2013 with a large spatial baseline 2324 m, for mapping coseismic deformation of this event.As shown in Figure 7, it has no significant SA in the E/W and N/S components of crustal deformation field produced by image pair 21 May 2013-11 May 2015.In contrast, the correlation results of image pair 10-26 September 2013 processed by the "mean subtracting" method still have some unmolded residuals, suggesting the "mean subtracting" method cannot completely remove the SA, especially when the stable regions are hard to choose for the calibration.Furthermore, the magnitude of the SA along the flight direction would be variable, not constant, if the topography along the satellite flight direction is changeable.In this study, the statistical results indicate that the image pair with the spatial baseline 200 m can lead to a γ less than 0.095 m, which can be neglected in surface deformation.Thus, it can be a recommended threshold value for spatial baseline selection in other cases.
May 2015.In contrast, the correlation results of image pair 10-26 September 2013 processed by the "mean subtracting" method still have some unmolded residuals, suggesting the "mean subtracting" method cannot completely remove the SA, especially when the stable regions are hard to choose for the calibration.Furthermore, the magnitude of the SA along the flight direction would be variable, not constant, if the topography along the satellite flight direction is changeable.In this study, the statistical results indicate that the image pair with the spatial baseline 200 m can lead to a γ less than 0.095 m, which can be neglected in surface deformation.Thus, it can be a recommended threshold value for spatial baseline selection in other cases.

Satellite Attitude Distortions
As the satellite attitude variations are sparsely sampled during flying, deformation fields (see Figure 8a) have cyclical distortions along the satellite flight direction, which are roll variations on the East-West component and pitch variations on the North-South component.Such error is distributed approximately constant along the satellite across-track direction.Thus, for the deformation fields with a serious attitude jitter distortions, the destriping procedure can be used to reduce cyclical distortions by subtracting the average along azimuth vertical direction based on stable calibration regions [11].

Satellite Attitude Distortions
As the satellite attitude variations are sparsely sampled during flying, deformation fields (see Figure 8a) have cyclical distortions along the satellite flight direction, which are roll variations on the East-West component and pitch variations on the North-South component.Such error is distributed approximately constant along the satellite across-track direction.Thus, for the deformation fields with a serious attitude jitter distortions, the destriping procedure can be used to reduce cyclical distortions by subtracting the average along azimuth vertical direction based on stable calibration regions [11].

Topographic Shadowing Artifacts
Solar position change, described by variation of sun azimuth and sun elevation angle, will cause shadowing difference in orientation and length (see Figure 9).The COSI-Corr algorithm will calculate related "offset" when the shadow texture dominates the correlation window.The experimental results in Figure A3 show that there are obvious topography-dependent biases in the Hinglaj Mountains area of the southeast image deformation fields, especially in the N/S component.Comprehensive analysis of Figure A3 and Table A2 shows that topography-dependent artifacts change with the sun angle disparity.The larger the difference in sun angle, the more obvious the artifacts.Generally, the artifacts are about 3-4 m, and the maximum is up to 10 m in the steep mountain area.Furthermore, the positive and negative variations of topography-dependent artifacts have a good correlation with the shadowing orientation changes (see Figure 9).Thus, we draw a conclusion that shadowing difference can be a source of the topography-dependent artifacts in image deformation fields, which is referred to "topographic shadowing artifacts" (TSA) in this study.
Based on above analysis, we can see the TSA is usually fixed in the spatial location and strongly related to the disparity of sun elevation and azimuth of the image.A factor "radiometric baseline" can be defined to evaluate the influence level of the TSA using the orientation and length of vector offset in Figure 9 and Table 2.The smaller is the radiometric baseline, the lower the level of the TSA and radiometric decorrelation noises becomes.As shown in Figure 7c,d, the radiometric baseline (0.14 h) in N/S component of image pair 10-26 September 2013 is much larger than that (0.05 h) of image pair 21 May 2013-11 May 2015 (h is the height of a reference datum in Figure 9).Thus, the TSA in the former pair is more significant than that of later pair in Figure 7.Moreover, the offset bias caused by the sun angle difference in 10 September 2013 and 26 September 2013 is mainly along N3° orientation (see Equation ( 3)), so the positive values TSA mainly distribute in the N/S component displacement map (see Figure 10a,d).

Topographic Shadowing Artifacts
Solar position change, described by variation of sun azimuth and sun elevation angle, will cause shadowing difference in orientation and length (see Figure 9).The COSI-Corr algorithm will calculate related "offset" when the shadow texture dominates the correlation window.The experimental results in Figure A3 show that there are obvious topography-dependent biases in the Hinglaj Mountains area of the southeast image deformation fields, especially in the N/S component.Comprehensive analysis of Figure A3 and Table A2 shows that topography-dependent artifacts change with the sun angle disparity.The larger the difference in sun angle, the more obvious the artifacts.Generally, the artifacts are about 3-4 m, and the maximum is up to 10 m in the steep mountain area.Furthermore, the positive and negative variations of topography-dependent artifacts have a good correlation with the shadowing orientation changes (see Figure 9).Thus, we draw a conclusion that shadowing difference can be a source of the topography-dependent artifacts in image deformation fields, which is referred to "topographic shadowing artifacts" (TSA) in this study.
Based on above analysis, we can see the TSA is usually fixed in the spatial location and strongly related to the disparity of sun elevation and azimuth of the image.A factor "radiometric baseline" can be defined to evaluate the influence level of the TSA using the orientation and length of vector offset in Figure 9 and Table 2.The smaller is the radiometric baseline, the lower the level of the TSA and radiometric decorrelation noises becomes.As shown in Figure 7c,d, the radiometric baseline (0.14 h) in N/S component of image pair 10-26 September 2013 is much larger than that (0.05 h) of image pair 21 May 2013-11 May 2015 (h is the height of a reference datum in Figure 9).Thus, the TSA in the former pair is more significant than that of later pair in Figure 7.Moreover, the offset bias caused by the sun angle difference in 10 September 2013 and 26 September 2013 is mainly along N3 • orientation (see Equation ( 3)), so the positive values TSA mainly distribute in the N/S component displacement map (see Figure 10a,d).Here we can reduce the TSA in the deformation field by choosing the image pairs with a small radiometric baseline.The good performance of this strategy has been demonstrated by image pair 21 May 2013-11 May 2015 of Figure 7.We can also utilize the combined radiometric baseline to reduce TSA using multi-pair of images, if the single small radiometric baseline image pair is unavailable due to few archived images in the study area.Based on these analyses, we can summarize the strategy of reducing the TSA and radiometric decorrelation noises as the small radiometric baseline method corresponding to two situations of single and combined radiometric baselines.In the situation of single radiometric baseline, there is only one correlated image pair with a very small radiometric baseline.The TSA and radiometric decorrelation noises in the correlation result are very small in this situation, so it can be ignored.Furthermore, in the corresponding situation of combined radiometric baseline, there are two or more correlated image pairs with very similar radiometric baselines.The TSA and radiometric decorrelation noises in those correlation results are very close each other, so it can be cancelled by combination of image pairs.In the 2013 Mw 7.7 Balochistan, Pakistan earthquake study, we find the radiometric baseline of image pair 13-29 September 2014 (0.142 h, N 3.1°) is very similar to that of image pair 10-26 September 2013 (0.139 h, N 2.9°) (see Table 2).Thus, we firstly acquire the image deformation field associated with the TSA and other seasonal errors except the earthquake event, using the correlation image pair 13-29 September 2014.The noise level of the TSA Here we can reduce the TSA in the deformation field by choosing the image pairs with a small radiometric baseline.The good performance of this strategy has been demonstrated by image pair 21 May 2013-11 May 2015 of Figure 7.We can also utilize the combined radiometric baseline to reduce TSA using multi-pair of images, if the single small radiometric baseline image pair is unavailable due to few archived images in the study area.Based on these analyses, we can summarize the strategy of reducing the TSA and radiometric decorrelation noises as the small radiometric baseline method corresponding to two situations of single and combined radiometric baselines.In the situation of single radiometric baseline, there is only one correlated image pair with a very small radiometric baseline.The TSA and radiometric decorrelation noises in the correlation result are very small in this situation, so it can be ignored.Furthermore, in the corresponding situation of combined radiometric baseline, there are two or more correlated image pairs with very similar radiometric baselines.The TSA and radiometric decorrelation noises in those correlation results are very close each other, so it can be cancelled by combination of image pairs.In the 2013 Mw 7.7 Balochistan, Pakistan earthquake study, we find the radiometric baseline of image pair 13-29 September 2014 (0.142 h, N 3.1 • ) is very similar to that of image pair 10-26 September 2013 (0.139 h, N 2.9 • ) (see Table 2).Thus, we firstly acquire the image deformation field associated with the TSA and other seasonal errors except the earthquake event, using the correlation image pair 13-29 September 2014.The noise level of the TSA in the image deformation field extracted from image pair 13-29 September 2014 is very similar to the counterpart in the coseismic deformation field extracted from image pair 10-26 September 2013, due to the similar radiometric baseline in Table 2.After removing the decorrelation noises, orbital error, SA and attitude jitter distortions, we can obtain a mask image from image pair 13-29 September 2014 by setting its deformation field greater than a threshold both in E/W and N/S components (see Figure 10b,e).We choose the non-seismic deformation threshold value as 0.8 m because the surface displacements with 1/10-1/20 pixel magnitude accuracy can be usually detected by COSI-Corr [8].Finally, we mask out corresponding TSA in the coseismic deformation field (10-26 September 2013) using the mask image obtained.This method can also be used to identify other non-targeted deformation regions with a high magnitude level, such as melting glaciers, migratory sand dunes and creeping landslides.Figure 10c, f indicate the E/W and N/S components of coseismic displacement field processed by our small radiometric baseline method in situation of combined radiometric baseline.Comparing between Figure 10d,f, we can find that TSA in the southeast and the moving sand dune area located in the northwest has been masked out effectively.We also use a non-local means filter tool in COSI-Corr to further reduce the noises in the deformation field.

Error Analysis and Accuracy Assessment
Using the methods described in Section 3, we can remove most error sources in Landsat 8 image deformation fields.In order to analyze the error characteristics and validate the performance of our new methods in deformation fields, we have analyzed statistically the error level of a stable area far from the crustal deformation field in image pair 10-26 September 2013 (see the polygon area in Figure 11).The selected area contains a variety of land features, including rolling hills, flat plains and vegetation area.We divide the whole data process into four steps and calculate the means and standard deviations after those procedures in Table 3 and Figure 11.Those four steps' results (R0, R1, R2, and R3) represent the process results of the original COSI-Corr process, after removing orbital error, after removing SA and attitude jitter distortions, and after removing TSA, respectively, in Figure 11.It should be noted that the procedures of reducing SA and attitude jitter distortions have been combined in R2, because they have similar error patterns.in the image deformation field extracted from image pair 13-29 September 2014 is very similar to the counterpart in the coseismic deformation field extracted from image pair 10-26 September 2013, due to the similar radiometric baseline in Table 2.After removing the decorrelation noises, orbital error, SA and attitude jitter distortions, we can obtain a mask image from image pair 13-29 September 2014 by setting its deformation field greater than a threshold both in E/W and N/S components (see Figure 10b,e).We choose the non-seismic deformation threshold value as 0.8 m because the surface displacements with 1/10-1/20 pixel magnitude accuracy can be usually detected by COSI-Corr [8].Finally, we mask out corresponding TSA in the coseismic deformation field (10-26 September 2013) using the mask image obtained.This method can also be used to identify other non-targeted deformation regions with a high magnitude level, such as melting glaciers, migratory sand dunes and creeping landslides.Figure 10c, f indicate the E/W and N/S components of coseismic displacement field processed by our small radiometric baseline method in situation of combined radiometric baseline.Comparing between Figure 10d,f, we can find that TSA in the southeast and the moving sand dune area located in the northwest has been masked out effectively.We also use a non-local means filter tool in COSI-Corr to further reduce the noises in the deformation field.

Error Analysis and Accuracy Assessment
Using the methods described in Section 3, we can remove most error sources in Landsat 8 image deformation fields.In order to analyze the error characteristics and validate the performance of our new methods in deformation fields, we have analyzed statistically the error level of a stable area far from the crustal deformation field in image pair 10-26 September 2013 (see the polygon area in Figure 11).The selected area contains a variety of land features, including rolling hills, flat plains and vegetation area.We divide the whole data process into four steps and calculate the means and standard deviations after those procedures in Table 3 and Figure 11.Those four steps' results (R0, R1, R2, and R3) represent the process results of the original COSI-Corr process, after removing orbital error, after removing SA and attitude jitter distortions, and after removing TSA, respectively, in Figure 11.It should be noted that the procedures of reducing SA and attitude jitter distortions have been combined in R2, because they have similar error patterns.In Figures 11 and 12, we use image pair 10-26 September 2013 instead of the small spatial baseline pair 21 May 2013-11 May 2015 for accuracy analysis because the former is more effective for the 2013 Mw 7.7 Balochistan, Pakistan earthquake study.It has been widely used for geophysical interpretation of this earthquake and also can be used for different error source validation in our study.The comparison in Figure 7 has explicitly proved the advantage of the small spatial baseline method.From Figure 12, both the mean and standard deviation of R1 and R2 become smaller compared to R0.It indicates that these two corresponding methods used in this paper are effective in reducing noises in coseismic deformation fields.After removing the TSA by the small radiometric baseline method, the noise level of R3 is further reduced, in which the standard deviation of E/W component reduces about 3.7% (around 0.02 m) and N/S component reduces about 15.2% (around 0.09 m).It is worth noting that, after subtracting the TSA from the N/S component of the coseismic deformation field (represented by R3), the absolute mean value increased from 0.012 m to 0.06 m (0.012 m vs. 0.06 m).As the TSA in the N/S component of the deformation field mainly have positive magnitudes, masking operation to reduce the TSA will mask amount of positive values.Thus, the mean of process result R3 shows a slight increase in the selected small area.We believe that the process result R3 would not increase the absolute mean value in large areas.As shown in Table 3, the final standard deviations of E/W and N/S components are 50.4cm and 51.8 cm, respectively.From this perspective, the selected threshold (0.8 m) in the process of removing TSA is just about |µ ± 2σ| (µ-mean, σ-standard deviation, [µ − 2σ, µ + 2σ] account for 95.45% measured values), which further validates the statistical reasonability of the threshold.

Comparison with Current Researches
In the process of error analysis, we could know the magnitude level of every type of error sources in Landsat 8 image deformation field, such that the orbital error with around −5-5 m, the SA with around −2-2 m, the attitude distortions with around −0.5-0.5 m, and the TSA with around −10-10 m.

Comparison with Current Researches
In the process of error analysis, we could know the magnitude level of every type of error sources in Landsat 8 image deformation field, such that the orbital error with around −5-5 m, the SA with around −2-2 m, the attitude distortions with around −0.5-0.5 m, and the TSA with around −10-10 m.In our study, we could easily remove the orbital error and attitude distortions.However, we can find little detailed studies with respect to the SA and TSA in Landsat 8 image deformation field.Thus, we endeavor to give comprehensive explanations for the reasons of error sources and propose corresponding weaken methods to reduce the SA and TSA from the perspective of time series analysis in our study.This is a work discriminating from other researches.
Due to differences in CCD sensors, the SA has different patterns and correction methods in different optical satellite images [31].For example, the rigorous CCD calibration model proposed by Leprince et al. [31] has a good applicability for the SPOT4 image, but it is not suitable for the Landsat 8 image due to lacking reference images and auxiliary state vector.In this study, we propose the small spatial baseline method for the Landsat 8 image to reduce the SA.Compared with the traditional "mean subtracting" method, it can obtain much cleaner and more accurate results if suitable spatial baseline image pair is available.Basing on the experimental results, we suggest that the spatial baseline of 200 m can be used as a reference threshold to select correlation image pairs.However, the "mean subtracting" method can be a complementary approach when the archived Landsat 8 image is very limited in the study area.Furthermore, for sudden events study, such as earthquakes, the "mean subtracting" method can also be a complementary approach, because the larger temporal baseline may introduce much larger signals and decorrelation noises into the crustal deformation field.
Topographic shadow effect is a common problem in optical images, its orientation and length variations can lead to the TSA, especially in monitoring glacier and landslide [1,29].Due to no effective reduction method, the TSA still remain in the coseismic deformation field of the 2013 Mw 7.7 Balochistan, Pakistan earthquake [25,27].In this study, we introduce a small radiometric baseline method to reduce the TSA and radiometric decorrelation noises.The TSA and other non-targeted displacements, crustal deformation excluding the coseismic deformation, are effectively masked out by our method, with much smaller mean and standard deviation both in E/W and N/S components.However, it is worth noting that although we could reduce the influence of TSA and radiometric decorrelation noises by our proposed method based on combined radiometric baselines, the crustal signals sometimes are also masked out when the signals and noises simultaneously occur in masked gaps.In addition, the TSA and stereoscopic parallax effects are easily confused, as they both having characteristics associated with topographic features.However, the latter also belongs to the geometric distortions and can be mitigated by a rigorous geometric model [1].We can discriminate them by analyzing base-height ratio or the temporal variation characteristics of these two types of errors.
The Landsat 8 satellite has been continuously monitoring the surface change since the first half year of 2013 and it has acquired sufficient images worldwide.Thus, it enables exploring the time series surface deformation by the multi-temporal Landsat 8 images.The key issue is how to construct the optimal image pairs to accurately acquire the crustal deformation fields.Our small spatial baseline and small radiometric baseline methods based on various error analyses and interpretations can be verified as effective time series analysis methods.The richer the Landsat 8 archived images, the more effective these methods are able to extract the earth surface deformation, with higher precision, which is very helpful for investigating the slow but steady dynamic phenomenon of earth surface, such as the melting glaciers and migratory sand dunes.Furthermore, based on detailed studies on error sources in this paper, we can acquire more precise fault rupture geometry and coseismic displacement field, which may be helpful for us to strengthen the understanding of fault displacement mechanism and improve the fault slip modeling.

Conclusions
In this study, we summarize spatio-temporal distribution characteristics and completely systematic removal processes of error sources in correlation results of Landsat 8 images.Through error analysis, we propose two new conceptual terms of "spatial baseline" and "radiometric baseline", as influential factors with relevance to the SA and TSA, respectively.Then, innovatively, we introduce the small spatial baseline method and small radiometric baseline method to reduce the SA and TSA in deformation field.Confirmatory experiments verify that the spatial baseline of correlated images with below 200 m is effective in reducing the influence of the SA in correlation results.Furthermore, the small radiometric baseline method, using combined radiometric baseline, improves the precision of the deformation field with apparent TSA by around 15% in our study.Meanwhile, the small radiometric baseline method, using single radiometric baseline, can also be helpful for reducing the radiometric decorrelation noises in deformation field.We also expect the content of this study can provide theoretical and experimental guidance for further utilizing the Landsat 8 image to monitor earthquakes, glaciers and so on.Notes: It is difficult to quantify the magnitude of the TSA in the correlation results due to their distribution characteristics, so we only make an orientation statistics of the shadowing offsets.The symbol (+) indicates eastward and northward offsets, correspondingly the symbol (−) indicates westward and southward offsets.These blank represent that the positive or negative characteristics of the measured errors are so tiny to identify.The theoretical calculated values, in which the offset orientation and offset length are equal to the orientation and magnitude of the "radiometric baseline", are calculated using a trigonometry based on the sun angle information as shown in Figure 9.We found that the measured orientations of TSA have a good consistency with the theoretical calculated offset orientations.

Figure 1 .
Figure 1.Coverage of Landsat 8 images over the study area.The red star indicates the epicenter of the 2013 Mw 7.7 Balochistan, Pakistan earthquake.The black boxes are the footprint of the Landsat 8 images in this study.The red rectangle frame in the inset represents the approximate location of the study area in Pakistan.

Figure 1 .
Figure 1.Coverage of Landsat 8 images over the study area.The red star indicates the epicenter of the 2013 Mw 7.7 Balochistan, Pakistan earthquake.The black boxes are the footprint of the Landsat 8 images in this study.The red rectangle frame in the inset represents the approximate location of the study area in Pakistan.

Figure 2 .
Figure 2. E/W and N/S components of image deformation fields.The serial numbers of all sub-graphs correspond to correlation image pairs in Table 1.Displacements are positive toward the East and North in this and following figures.

Figure 3 .
Figure 3. Error sources in E/W and N/S components of image deformation field, selecting the correlation image pair 29 November 2013-15 December 2013 as an example.We can see there are obvious orbital error, decorrelation noises, topography-dependent artifacts, and stripe artifacts in image deformation field.

Figure 2 .
Figure 2. E/W and N/S components of image deformation fields.The serial numbers of all sub-graphs correspond to correlation image pairs in Table 1.Displacements are positive toward the East and North in this and following figures.

Figure 2 .
Figure 2. E/W and N/S components of image deformation fields.The serial numbers of all sub-graphs correspond to correlation image pairs in Table 1.Displacements are positive toward the East and North in this and following figures.

Figure 3 .
Figure 3. Error sources in E/W and N/S components of image deformation field, selecting the correlation image pair 29 November 2013-15 December 2013 as an example.We can see there are obvious orbital error, decorrelation noises, topography-dependent artifacts, and stripe artifacts in image deformation field.

Figure 3 .
Figure 3. Error sources in E/W and N/S components of image deformation field, selecting the correlation image pair 29 November 2013-15 December 2013 as an example.We can see there are obvious orbital error, decorrelation noises, topography-dependent artifacts, and stripe artifacts in image deformation field.

Figure 4 .
Figure 4. E/W and N/S components of image deformation field before and after removing the orbital error, taking the correlation image pair 29 November 2013-15 December 2013 as an example.

Figure 4 .
Figure 4. E/W and N/S components of image deformation field before and after removing the orbital error, taking the correlation image pair 29 November 2013-15 December 2013 as an example.

Figure 5 .
Figure 5. Relevance of the SA (stripe artifacts) and the CCD (Charge Coupled Device) arrays' relative locations superimposed on the E/W component of the crustal deformation field.The yellow rectangles represent CCD arrays' locations acquired on the images from: 10 September 2013 (top); and 26 September 2013 (bottom).The dislocation between the two CCD arrays presents inconsistency in scene center (spatial baseline), which is the main reason causing the SA.The misalignment between CCD sub-arrays in a single CCD chip is identified as the essential reason.Blue profile line represents the average along flight direction of gray block area with a magnitude range of −1.5 m-2 m (−0.1-0.13 pixels).

Figure 5 .
Figure 5. Relevance of the SA (stripe artifacts) and the CCD (Charge Coupled Device) arrays' relative locations superimposed on the E/W component of the crustal deformation field.The yellow rectangles represent CCD arrays' locations acquired on the images from: 10 September 2013 (top); and 26 September 2013 (bottom).The dislocation between the two CCD arrays presents inconsistency in scene center (spatial baseline), which is the main reason causing the SA.The misalignment between CCD sub-arrays in a single CCD chip is identified as the essential reason.Blue profile line represents the average along flight direction of gray block area with a magnitude range of −1.5 m-2 m (−0.1-0.13 pixels).

Figure 6 .
Figure 6.Correlation analysis of γ and spatial baselines of the correlated image pairs.Note that we deleted the first dataset due to high measurement bias, and fit the function using the other 15 datasets.

Figure 6 .
Figure 6.Correlation analysis of γ and spatial baselines of the correlated image pairs.Note that we deleted the first dataset due to high measurement bias, and fit the function using the other 15 datasets.

Figure 7 .
Figure 7.Comparison between the "mean subtracting" method and small spatial baseline method in terms of the SA reduction: (a,b) the E/W components of 10-26 September 2013 using the "mean subtracting" method and 21 May 2013-11 May 2015 using the small spatial baseline method; and (c,d) the N/S components of 10-26 September 2013 using the "mean subtracting" method and 21 May 2013-11 May 2015 using the small spatial baseline method.This comparison further demonstrates that the magnitude level of SA has a strong correlation with spatial baseline of the correlated image pair.In addition, the ellipse lines circle the topography-dependent artifacts area, in which the red arrows in (c) indicate the specific locations.

Figure 7 .
Figure 7.Comparison between the "mean subtracting" method and small spatial baseline method in terms of the SA reduction: (a,b) the E/W components of 10-26 September 2013 using the "mean subtracting" method and 21 May 2013-11 May 2015 using the small spatial baseline method; and (c,d) the N/S components of 10-26 September 2013 using the "mean subtracting" method and 21 May 2013-11 May 2015 using the small spatial baseline method.This comparison further demonstrates that the magnitude level of SA has a strong correlation with spatial baseline of the correlated image pair.In addition, the ellipse lines circle the topography-dependent artifacts area, in which the red arrows in (c) indicate the specific locations.

Figure 8 .
Figure 8. E/W and N/S components before (a) and after (b) subtracting the attitude distortions (red line) along the across-track direction in correlation image pair 16 November 2014-2 December 2014.

Figure 8 .
Figure 8. E/W and N/S components before (a) and after (b) subtracting the attitude distortions (red line) along the across-track direction in correlation image pair 16 November 2014-2 December 2014.

Figure 9 .
Figure 9. Mechanism of topographic shadowing artifacts.Equations below the figures show how to calculate the magnitude and orientation of the offset caused by shadowing variation.(α1, β1) and (α2, β2) represent the sun elevation and azimuth of pre-and post-event images, respectively.L1 and L2 represent the shadow lengths on a flattish surface of a reference datum named O with height h at the sun elevation of α1 and α2, respectively.Equations (2) and (3) show how to calculate the magnitude and orientation θ of vector offset due to shadowing variation, in which ∆α and ∆β denote the sun elevation and sun azimuth disparity, respectively.Equations (4) and (5) represent how to calculate the E/W offset and N/S offset decomposed from the vector offset, and the results are positive toward East and North.

Figure 9 .
Figure 9. Mechanism of topographic shadowing artifacts.Equations below the figures show how to calculate the magnitude and orientation of the offset caused by shadowing variation.(α1, β1) and (α2, β2) represent the sun elevation and azimuth of pre-and post-event images, respectively.L1 and L2 represent the shadow lengths on a flattish surface of a reference datum named O with height h at the sun elevation of α1 and α2, respectively.Equations (2) and (3) show how to calculate the magnitude and orientation θ of vector offset due to shadowing variation, in which ∆α and ∆β denote the sun elevation and sun azimuth disparity, respectively.Equations (4) and (5) represent how to calculate the E/W offset and N/S offset decomposed from the vector offset, and the results are positive toward East and North.

Figure 10 .
Figure 10.Horizontal coseismic deformation of the 2013 Mw 7.7 Balochistan, Pakistan earthquake obtained by our small radiometric baseline method, before and after masking out TSA and other nontargeted offsets: (a) E/W displacement field before masking; (b) E/W mask image, in which the white denotes the region that need to be masked; (c) E/W displacement field after masking; (d) N/S displacement field before masking; (e) N/S mask image; and (f) N/S displacement field after masking.

Figure 10 .
Figure 10.Horizontal coseismic deformation of the 2013 Mw 7.7 Balochistan, Pakistan earthquake obtained by our small radiometric baseline method, before and after masking out TSA and other non-targeted offsets: (a) E/W displacement field before masking; (b) E/W mask image, in which the white denotes the region that need to be masked; (c) E/W displacement field after masking; (d) N/S displacement field before masking; (e) N/S mask image; and (f) N/S displacement field after masking.

Figure 11 .
Figure 11.Flow chart of accuracy improvement in the Landsat 8 crustal deformation field.The black polygon line denotes the area used for accuracy assessment.In Figures 11 and 12, we use image pair 10-26 September 2013 instead of the small spatial baseline pair 21 May 2013-11 May 2015 for accuracy analysis because the former is more effective for the 2013 Mw 7.7 Balochistan, Pakistan earthquake study.It has been widely used for geophysical interpretation of this earthquake and also can be used for different error source validation in our

Figure 11 .
Figure 11.Flow chart of accuracy improvement in the Landsat 8 crustal deformation field.The black polygon line denotes the area used for accuracy assessment.
comparison in Figure7has explicitly proved the advantage of the small spatial baseline method.From Figure12, both the mean and standard deviation of R1 and R2 become smaller compared to R0.It indicates that these two corresponding methods used in this paper are effective in reducing noises in coseismic deformation fields.After removing the TSA by the small radiometric baseline method, the noise level of R3 is further reduced, in which the standard deviation of E/W component reduces about 3.7% (around 0.02 m) and N/S component reduces about 15.2% (around 0.09 m).It is worth noting that, after subtracting the TSA from the N/S component of the coseismic deformation field (represented by R3), the absolute mean value increased from 0.012 m to 0.06 m (0.012 m vs. 0.06 m).As the TSA in the N/S component of the deformation field mainly have positive magnitudes, masking operation to reduce the TSA will mask amount of positive values.Thus, the mean of process result R3 shows a slight increase in the selected small area.We believe that the process result R3 would not increase the absolute mean value in large areas.As shown in Table3, the final standard deviations of E/W and N/S components are 50.4cm and 51.8 cm, respectively.From this perspective, the selected threshold (0.8 m) in the process of removing TSA is just about |μ ± 2σ| (μ-mean, σ-standard deviation, [μ − 2σ, μ + 2σ] account for 95.45% measured values), which further validates the statistical reasonability of the threshold.

Figure 12 .
Figure 12.Histograms of the process results R0, R1, R2, and R3 in the selected area.The histograms are shown as below: original COSI-Corr process result (R0, red line), result after removing orbital error (R1, black line), result after removing SA and attitude jitter distortions (R2, blue line), and result after removing TSA (R3, pink line).The vertical dotted lines represent the relevant mean (location) of histogram curves in the same color.

Figure 12 .
Figure 12.Histograms of the process results R0, R1, R2, and R3 in the selected area.The histograms are shown as below: original COSI-Corr process result (R0, red line), result after removing orbital error (R1, black line), result after removing SA and attitude jitter distortions (R2, blue line), and result after removing TSA (R3, pink line).The vertical dotted lines represent the relevant mean (location) of histogram curves in the same color.

Figure A1 .
Figure A1.Matching distribution of 16 correlation pairs using 25 Landsat 8 panchromatic images, following the uniform temporal standard of 16 days.

Figure A2 .
Figure A2.Patterns of the SA in E/W and N/S deformation fields.Clearly, patterns of these three pairs (22 April 2014-8 May 2014, 16 November 2014-2 December 2014, and 20 February 2015-8 March 2015) with much less variation are different from the ones of the rest image pairs.

Figure A2 .
Figure A2.Patterns of the SA in E/W and N/S deformation fields.Clearly, patterns of these three pairs (22 April 2014-8 May 2014, 16 November 2014-2 December 2014, and 20 February 2015-8 March 2015) with much less variation are different from the ones of the rest image pairs.

Figure A3 .
Figure A3.E/W (a) and N/S (b) deformation fields after removing the orbital error.The black ellipse lines circle the area of dominated TSA in E/W and N/S deformation fields.The symbol (+) indicates eastward and northward offsets of TSA and the symbol (−) indicates westward and southward offsets of TSA.

Figure A3 .
Figure A3.E/W (a) and N/S (b) deformation fields after removing the orbital error.The black ellipse lines circle the area of dominated TSA in E/W and N/S deformation fields.The symbol (+) indicates eastward and northward offsets of TSA and the symbol (−) indicates westward and southward offsets of TSA.

Table
22 April 2014-8 May 2014, 16 November 2014-2 December 2014 and 20 February 2015-8 March 2015) present approximate straight lines with the mean value closing to zero.From Table

Table 1 .
Spatial baseline and γ of all the correlation image pairs in this study.

Table 1 .
Spatial baseline and γ of all the correlation image pairs in this study.

Table 2 .
Radiometric baseline list of the Landsat 8 correlation image pairs.
Notes: The ∆Azi. and ∆Ele.represent the sun azimuth disparity and sun elevation disparity, respectively.

Table 2 .
Radiometric baseline list of the Landsat 8 correlation image pairs.

Table A1 .
Parameters catalog of Landsat 8 OLI Pan images used in the study.
Notes: All images are acquired in Nadir and processed into L1T level.

Table A2 .
Relationship between the TSA and Sun angles of Landsat 8 correlation images pairs.