remote sensing An Objective Method with a Continuity Constraint for Improving Surface Velocity Estimates from the Geostationary Ocean Color Imager

: Mapping surface currents with high spatiotemporal resolution over a wide coverage is crucial for understanding ocean dynamics and associated biogeochemical processes. The most widely used algorithm for estimating surface velocities from sequential satellite observations is the maximum cross-correlation (MCC) method. However, many unrealistic vectors still exist, despite the utilization of various ﬁltering techniques. In this study, an objective method has been developed through the combination of MCC and multivariate optimum interpolation (MOI) analysis under a continuity constraint. The MCC method, with and without MOI, is applied to sequences of simulated sea surface temperature (SST) ﬁelds with a 1/48 ◦ spatial resolution over the East China Sea continental shelf. Integration of MOI into MCC reduces the average absolute differences between the model’s ‘actual’ velocity and the SST-derived velocity by 19% in relative magnitude and 22% in direction, respectively. Application of the proposed method to Geostationary Ocean Color Imager (GOCI) satellite observations produces good agreement between derived surface velocities and the Oregon State University (OSU) regional tidal model outputs. Our results demonstrate that the incorporation of MOI into MCC can provide a signiﬁcant improvement in the reliability and accuracy of satellite-derived velocity ﬁelds.


Introduction
Ocean surface current observations have critical importance in understanding ocean dynamics, air-sea interactions, and biogeochemical processes. Particularly in coastal waters, real-time current observations with a high spatiotemporal resolution take much more significance through their impacts on search and rescue missions, environmental monitoring, and optimal ship routing. Traditional direct measurements of surface currents are both time-and cost-consuming, and it is impossible to obtain two-dimensional, continuous surface current fields with wide coverage. The prospect of observing nearly instantaneous surface velocity fields without temporal aliasing over a large area has prompted numerous attempts to infer advective surface currents directly from sequential satellite imagery [1][2][3][4][5][6][7].
A well-known method called the maximum cross-correlation (MCC) method is by far the most widely used to obtain high spatial resolution surface currents in coastal areas from thermal infrared and ocean color satellite observations since it is less sensitive to errors and does not need precise values [2,[8][9][10]. For example, using sequential sea surface satellite (SST) imagery in the British Columbia coastal ocean, Emery et al. [2] succeeded in retrieving the surface currents by tracking the displacements in SST pattern. By combing chlorophyll a with SST images off the California coast, Crocker et al. [8] demonstrated that spatial and temporal coverage of satellite-derived surface currents could be increased significantly. Note that polar-orbiting satellites acquire only one or two images a day, which poses challenges for application in coastal waters with strong tidal currents. The Geostationary Ocean Color Imager (GOCI), on the other hand, could provide up to eight hourly snapshots per day with a spatial resolution of 500 m [11]. The unprecedented GOCI capability is adequate for deriving changing surface currents in the China seas (here referring to the Bohai Sea, Yellow Sea (YS), and East China Sea (ECS)) where the tidal currents dominate [7,10,[12][13][14][15][16][17][18]. Substantial efforts have been carried out to retrieve the hourly variability of surface currents using various GOCI-derived hourly ocean color products, such as total suspended matter (TSM) [13,15,17], chlorophyll a [4,6], normalized water-leaving radiances, and a diffuse attenuation coefficient [10,12]. Based on GOCIderived hourly surface currents, Hu et al. [14][15][16] extracted semi-diurnal M2 tidal currents and mapped snapshots of the Changjiang (Yangtze River) plume, the YS warm current, and the seasonal surface circulations of the Taiwan Strait. The satellite-derived surface currents have been validated with extensively independent in situ current measurements, primarily from an acoustic Doppler current profiler (ADCP), surface drifters in drogue at 15 m, and high-frequency radars [4,10,12,[14][15][16]18].
While the great advantage of the MCC method in quantitatively mapping nearly instantaneous surface currents over large spatial and temporal scales has been well-documented, obviously, erroneous incoherent vectors still exist due to mismatching of surface features under the influence of clouds, image registration, and/or ocean dynamical processes. To identify and remove the unrealistic vectors, several techniques are proposed, such as a minimum cross-correlation coefficient limit [2], a reasonable velocity filter [9], a nearest neighbor filter [19], and a reciprocal filter [20]. However, their performance still needs to be improved, and some reliable vectors are often discarded.
The main objective of this study is to develop a robust objective method for improving satellite-derived surface velocity estimates without data loss. The paper is organized as follows: Section 2 describes the data and methodology; Section 3 evaluates the validity and accuracy of the proposed method by comparing the 'ground truth' with the MCC-derived velocities from high-resolution simulated SST; Section 4 applies the newly developed method to the GOCI image over the ECS shelf, and the resulting velocities are compared with a well-validated regional tidal simulation; discussions and conclusions are given in the last section.

GOCI Data
GOCI has a spatial coverage of around 2500 km × 2500 km, centered at (130 • E, 36 • N) in northeast Asia, and eight spectral bands (six visible and two near-infrared (NIR)). It can acquire eight hourly snapshots between 00:15 and 07:45 UTC [11]. Level-1B data with rare cloud cover for 5 April 2011, 13 May 2013, and 16 September 2013 are downloaded from the Korea Institute of Ocean Science and Technology (KIOST; http://kosc.kiost.ac.kr, accsess on 1 July 2021), and processed using its official GOCI Data Process Software to generate TSM. To remove noise, a median filter with a 5 × 5 pixel window is applied to the TSM images, following Hu et al. [18]. Figure 1 shows the composite of eight hourly GOCI-derived TSM images on 13 May 2013 in the Changjiang estuary and its adjacent waters. The distinction in TSM between the nearshore and offshore waters is very remarkable. High TSM is mainly restricted to shallow waters (<30 m depth).

Ocean Model Data
The global Massachusetts Institute of Technology general circulation model (MITgcm) on a latitude-longitude polar cap grid (named as LLC4320) has a horizontal resolution of 1/48 • and provides hourly hindcasts of temperature and ocean currents from September 2011 to November 2012 (https://data.nas.nasa.gov/ecco/data.php, accessed on 18 December 2021). It is jointly driven by both six-hourly atmospheric fields from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System and hourly tidal forcing from the 16 most significant tidal constituents [21,22]. In this study, we randomly selected MITgcm-simulated surface temperatures and currents at hourly intervals during spring and neap tides and for different seasons. Details regarding the evaluation of the utility of the proposed objective method can be seen in Section 2.3.
The Oregon State University (OSU) regional barotropic tidal model with data inversion (China Seas 2010) [23,24] has a spatial resolution of 1/30 • . It assimilates altimetry sea level observations (https://www.tpxo.net, accessed on 18 December 2021). The simulated tidal currents in the China seas have been validated with in situ observations from historical moored ADCP and HF radar. The average absolute differences between OSU outputs and total surface currents (combination of both tidal and wind-driven components) are about 0.5 cm/s in magnitude and 10 • in direction, with a correlation coefficient of 0.93 [12]. These findings suggested that wind-driven velocities play a very limited role in the total surface currents in the China seas, where tidal currents significantly dominate [16,25]. In this study, the OSU model tidal currents on 13 May 2013 are used for comparison with GOCI-derived surface velocities.

MCC Method
The MCC method is commonly used to estimate hourly total surface currents from the sequential tracer images [4,6,7,10,12,13]. Briefly, a sub-source box in an initial image is cross-correlated with the same size box in a subsequent image, searching for the most likely location in the second image by the maximum cross-correlation coefficient. The water parcel displacement corresponding to the first sub-source box is assumed to be the distance between its initial location and the sub-source center in the second image with the maximum cross-correlation. Following Hu et al. [15], we used a 20 × 20 pixel template (first image) window and a 36 × 36 pixel search (second image) window. To eliminate erroneous vectors, three subsequent steps are performed ( Figure 2). Vectors with the maximum cross-correlations < 0.9 are discarded. A reciprocal filter that combines the information of the forward and backward MCC tracking between two sequential images is employed [26]. Vectors with speeds larger than realistic values in this study region (i.e., 2.5 m/s) are ignored. More details are given in Hu et al. [18]. The performance of the MCC method is well-demonstrated with various types of surface current measurements in previous studies [4,10,12,13,15,18].

Multivariate Optimum Interpolation Analysis
Using multivariate optimum interpolation (MOI) analysis based on the Gauss-Markov Theorem [27,28], the MCC-derived velocity components, u and v, are fitted simultaneously to obtain a stream function, assuming horizontal non-divergence. The stream function (Ψ) can be formulated for multivariate functional fitting, where Ψ(r i ) is the optimal estimate of the stream function, N is the number of the MCCderived velocities, r ik is the distance (|r i −r k |) between analysis grid (r i ) and observation point (r k ), W u r ik and W v r ik are the weights for u and v components, and u(r k ) and v(r k ) are the MCC-derived velocities. Once the weights are obtained, the optimal stream function can be calculated from (1) and the surface velocities from the following equation, where u and v are the zonal and meridional velocity components. By minimizing the error variance of Equation (1) with respect to the weights. The resulting equations are where denotes the sum over all the analysis points, and l denotes the observation point. A Gaussian autocorrelation function for the stream function is assumed, where L is the 0.5-degree radius of influence, k and l denote the analysis and observation points. The velocity autocorrelation functions Ψ,u , Ψ,v , and u,v can be easily obtained from Equation (2). The weights are then estimated from Equations (3) and (4) using the least-squares method.

Synthetic Tests
To assess the MOI effectiveness, surface velocities based on the MCC method with and without MOI are estimated from the MITgcm simulated sequential SST and compared directly with the model velocity vectors as a benchmark. Since the model velocity data are 'ground truth' with no sources of errors, the errors of the SST-derived velocities are attributed to the methods for estimating the surface velocities. In order to quantitatively evaluate the MOI performance, angular and relative magnitude errors for motion estimations are used following Chen [29]. Assuming w = 0 in the vertical direction, the average angular and relative magnitude errors (∆θ = AAE and ∆V V = AME) between the model velocityv and an estimate v are defined as: where the relative magnitude error is a dimensionless quantity. We define The bracket indicates averaging over all locations. Moreover, a skill (S) as a more strict measure of the goodness-of-fit is used following Hu, Wang, Pan, He, Miyazawa, Bai, Wang and Gong [16], where u o and u m are the simulated and derived velocities, respectively. A perfect agreement has S = 1. Figure 3 shows the comparisons between the benchmark velocity vectors from the numerical model and the MCC-derived velocity fields with and without MOI for the two time intervals on 13 September 2011, namely 12:00-14:00 UTC and 19:00-21:00 UTC. The spatial patterns of the three velocity vector fields are mutually similar, though there are few vectors within the highly turbid coastal zone. Surface velocity fields for different tidal phases can be well-captured with both MCC and MCC+MOI. Strong surface velocities with speeds > 1 m/s are mostly restricted to shallow waters (<50 m depth), whereas weak surface currents occur in relatively deep waters. On the other hand, the overall agreement with the benchmark vectors is much better when incorporating MOI into the MCC method. Those evident erroneous vectors can be effectively corrected from the MCC-derived velocity fields. Relative to the benchmark vectors, AAEs are greatly reduced from 10 • to 7 • for 12:00-14:00 UTC and from 8 • to 7 • for 19:00-21:00 UTC, respectively. AMEs are also decreased from 0.29 to 0.25 and from 0.29 to 0.22 for the two time intervals, respectively. The skills are significantly improved from 0.90 to 0.94 for 12:00-14:00 UTC and from 0.89 to 0.91 for 19:00-21:00 UTC, respectively.
To further assess the accuracy of the MCC method with and without MOI, Figure 4 shows histograms of misfits in magnitude and direction of the MCC-derived surface velocity fields with and without MOI for 12:00-14:00 UTC and 19:00-21:00 UTC on 13 September 2011. There exist significant differences between MCC alone and MCC+MOI. Relative to the MCC vectors, the MCC+MOI vectors agree better with the benchmark velocities. Overall, both ∆V and ∆θ for MCC+MOI vectors have narrower ranges and larger peaks. The histograms for MCC+MOI are closer to zero (∆V = ∆θ = 0). These results indicate that the incorporation of MOI into MCC can greatly improve the accuracy of surface velocity estimates, together with no data loss.  To avoid the possibility that the above good results are a mere coincidence, more cases are tested for different background conditions, such as in spring and neap tides and four seasons (See Supplementary Materials). The qualitative and quantitative comparisons consistently show that ∆V and ∆θ can be decreased considerably through the combination of MOI with MCC, indicating that the effectiveness of MOI is realistic and robust.

Application to GOCI Images
The intended application of MOI is to improve estimates of sea surface velocities from GOCI image sequences. To examine the performance of the combined solution (MCC+MOI) for GOCI satellite observations, hourly surface currents over the ECS shelf are reconstructed from the velocity vectors retrieved from GOCI-derived TSM images and directly compared with the OSU model currents both qualitatively and quantitatively. Figure 5 shows the comparisons of velocity vectors from the OSU model with those estimated by the MCC method without and with MOI on 13 May 2013. Note that we only plot the OSU model velocities on the grids where there exist MCC velocities. The three velocity fields have very similar features. The velocity vector transitions from the flood tide to the ebb tide are very apparent. However, it is obvious that there exist many more unrealistic vectors from MCC alone than from MCC+MOI ( Figure 5, middle panel). The combined solutions have smoother results. The MCC+MOI can produce great improvement in GOCI-derived surface velocities. Statistical measures for the comparison between the model simulations and velocity estimates from MCC and MCC+MOI on 13 May 2013 is shown in Figure 6. The absolute differences in the magnitude and direction from MCC+MOI are smaller than those from MCC alone. Averaged over four time intervals, AAEs decrease greatly from 31 • to 18 • while AMEs are from 0.7 to 0.5. The skill of 0.76 from MCC+MOI is also much larger than that from MCC (0.37). The statistical analysis indicates that the combined solution (MCC+MOI) has better performance than MCC alone.
To confirm the robustness of MCC+MOI, another two cases from relatively cloud-free GOCI satellite observations are examined. The results are compared with the corresponding OSU model outputs (See Supplementary Materials). Similar findings to those in Figures 5  and 6 are obtained. This further corroborates that the velocity estimates can be greatly improved through the combination of MCC with MOI. Although an evaluation of the realism of the smoother velocity fields is limited by the lack of comparison with in situ data, the solutions estimated by MCC+MOI are realistic in comparison with the OSU model.

Discussions and Conclusions
Previous studies suggested the existence of many unrealistic velocity vectors in MCCderived surface currents [2,7]. Great efforts are made to eliminate erroneous vectors using various techniques [6,9,30]. The problem, however, remains unconquered. In this study, a robust technique for estimating surface velocity fields in our study has been developed through the incorporation of MOI into the MCC method. Under the constraint of two-dimensional continuity requirement, a stream function is fitted by MOI based on MCC-derived surface velocities and then used to reconstruct the MCC velocity estimates in order to improve reliability and accuracy of the MCC method without data loss.
The high-resolution MITgcm solutions are used as a benchmark to examine the utility of both MOI and MCC combinations. Three independent measures of goodness-of-fit (i.e., angular and magnitude errors and skill) identically indicate that the MCC+MOI combination significantly outperforms the MCC alone for estimating surface velocities. Application of the MCC+MOI to actual GOCI image sequences with both magnitude and angular errors exhibits a significant decrease, as shown in Figure 6. These results suggest that the combination of MOI with MCC could not only effectively correct many erroneous velocity vectors but improve the validity and accuracy of velocity estimates as well.
Both velocity fields from the model-simulated SST and satellite-derived TSM images show that the MCC method is not capable of effectively detecting displacements of the tracer (SST/TSM) patterns in highly turbid coastal waters (<20 m depth) (Figures 1 and 3). This suggests that processes other than horizontal advection dominate the heat and sediment budget and make major contributions to the changing SST and TSM patterns. Note that the basic assumption in the MCC approach is that horizontal advection is assumed to be the dominant mechanism contributing to the linear displacement of tracer features in the time interval between the images. The non-advective mechanisms include surface heating and cooling for SST (sinking and resuspension for TSM), vertical convection (upwelling and downwelling), mixing, and diffusion [2,31]. The time interval between the two models SST/TSM is short. It seems unlikely for diffusion and surface heating (cooling) to change the SST pattern, as reported by Wahl and Simpson [32], based on idealized model experiments. They found that horizontal diffusion has only a minor effect at time intervals of 24 h and less. In this shallow zone (<20 m depth), where tidal currents are extremely strong, the tide can easily induce intense mixing, which may be primarily responsible for affecting the MCC estimates through altering the SST and TSM patterns.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/rs14010014/s1, Figure S1. The hourly tidal heights at (122 • E, 31 • N) from 10 September to 15 October 2011. The start time for each case during spring and neap tides is marked by a bold gray line. The time is in Greenwich. Figure S2. The averaged velocity fields from model (left), MCC (middle), and MCC+MOI (right) for 12:00-14:00 UTC and 19:00-21:00 UTC on 30 September 2011 (during spring tide), overlying with the MITgcm model SST for 12:00 UTC (a), 14:00 UTC (b,c), 19:00 UTC (d), and 21:00 UTC (e,f). Figure S3. Same as Figure S2, but for 7 October 2011 (during neap tide). Figure S4. Same as Figure S2, but for 1 January 2012 (in winter). Figure S5. Same as Figure S2, but for 1 April 2012 (in spring). Figure S6. Same as Figure S2, but for 1 July 2012 (in summer). Figure S7. Same as Figure S2, but for 1 October 2011 (in autumn). Figure S8. Histograms of magnitude and angle differences between the MITgcm model velocities and those from either MCC or MCC+MOI on 30 September 2011 (during spring tide). For this histogram comparison, the total number of model vectors equals those estimated using MCC and MCC+MOI. Figure S9. Same as Figure S8, but for 7 October 2011 (during neap tide). Figure S10. Same as Figure S9, but for 1 January 2012 (in winter). Figure S11. Same as Figure S8, but for 1 April 2012 (in spring). Figure S12. Same as Figure S8, but for 1 July 2012 (in summer). Figure S13. Same as Figure S8, but for 1 October 2011 (in autumn).  Figure S15. Same as Figure S14, but for 16 September 2013. Figure S16. Comparison of AME, AAE, and skill for velocities from the OSU model and from GOCI imagery using MCC and MCC+MOI for four different time intervals on 5 April 2011. Figure S17. Same as Figure S16