*Article* **Extremely Robust Remote-Target Detection Based on Carbon Dioxide-Double Spikes in Midwave Spectral Imaging**

#### **Sungho Kim 1,\*, Jungsub Shin 2, Joonmo Ahn <sup>2</sup> and Sunho Kim <sup>2</sup>**


Received: 6 April 2020; Accepted: 15 May 2020; Published: 20 May 2020

**Abstract:** Infrared ship-target detection for sea surveillance from the coast is very challenging because of strong background clutter, such as cloud and sea glint. Conventional approaches utilize either spatial or temporal information to reduce false positives. This paper proposes a completely different approach, called carbon dioxide-double spike (CO2-DS) detection in midwave spectral imaging. The proposed CO2-DS is based on the spectral feature where a hot CO2 emission band is broader than that which is absorbed by normal atmospheric CO2, which generates CO2-double spikes. A directional-mean subtraction filter (D-MSF) detects each CO2 spike, and final targets are detected by joint analysis of both types of detection. The most important property of CO2-DS detection is that it generates an extremely low number of false positive caused by background clutter. Only the hot CO2 spike of a ship plume can penetrate atmosphere, and furthermore, there are only ship CO2 plume signatures in the double spikes of different spectral bands. Experimental results using midwave Fourier transform infrared (FTIR) in a remote sea environment validate the extreme robustness of the proposed ship-target detection.

**Keywords:** ship detection; false alarm; farbon dioxide peaks; midwave infrared; FTIR

#### **1. Introduction**

Remote-ship detection is important in various applications, such as maritime navigation [1], coast guard searches [2], and homeland security procedures [3]. Infrared images are frequently adopted, especially due to their operational capability, day or night.

The key issue is how to reduce the high number of false detections caused by cloud clutter and sea surface glint while maintaining an acceptable detection rate [4,5]. Cloud edges produce false detections, and sea glint is similar to small infrared targets, both of which degrade detection performance. Since the 1990s, a lot of methods have been proposed to reduce the false positives by using either spatial information or temporal information from infrared images.

Spatial information processing, such as background subtraction, can be a feasible approach to reducing background clutter. Background images are estimated by spatial filters, such as the mean filter [6], the least mean square filter [7,8], the median filter [9], and the morphological filter [10]. These filters use neighboring pixels to estimate background pixel values, but with different strategies. In particular, mean subtraction filter (MSF)-based target detection is the most simple, but is weak when it comes to cloud edges. Non-linear filters, such as max-median, morphology, and data-fitting, show better cloud clutter suppression capabilities [11,12]. Local directional Laplacian-of-Gaussian filtering can remove false positives from cloud edges [13]. In the spatial

information approach, classification with a spatial shape feature can discriminate clutter. Hysteresis thresholding [14], statistics-based adaptive thresholding [15], the Bayesian classifier [16], and support vector machines [17] are well-known methods. Voting by various classifiers can enhance dim-target detection rates [18]. Spatial frequency information can be used to remove low-frequency clutter, such as the three-dimensional fast Fourier transform (3D-FFT) spectrum [19]. The wavelet transform is effective against sea glint [20], and low pass filter can also deal with sea glint [21]. An adaptive high-pass filter can reduce cloud and sea glint clutter [22], while spatial target-background contextual information enhances the target signature and reduces background clutter, which is effective for sea glint reduction [23]. Spatial multi-feature fusion can increase the clutter discrimination capability [24].

Temporal information processing, such as track-before-detection, can remove slowly moving cloud clutter and fast-blinking sea glint [25]. Temporal profiles are useful to discriminate between a fast-moving target and slow-moving cloud [26–28]. The temporal contrast filter is suitable for detecting small, supersonic, infrared targets [29], and a 3D matched filter for wide-to-exact searches can remove cloud clutter very fast [30]. A power-law detector for image frames is effective for target detection in dense clutter [31]. Since a previous frame can be regarded as background, a weighted autocorrelation matrix update using the recursive technique can be useful in order to eliminate sea glint [32]. An advanced adaptive spatial-temporal filter can achieve tremendous gain [33], and principal component analysis in multi-frames can alleviate false positives from sun glint [34].

Synthetic aperture radar (SAR) with deep learning can be a useful approach for remote ship detection application [35–37]. Faster R-CNN [35], DRBox [36], and Contextual CNN [37] can provide improved detection performance but these approaches require huge number of training dataset to reduce false positives.

As mentioned above, the conventional approaches in the spatial- or temporal-image domains have their own pros and cons, depending on the situations and conditions. The spatial filter based approach relies on target shapes and intensity distributions, which is ambiguous if targets are far away. The temporal filter based approach assumes fast target motion with stationary sensors. If this assumption is not satisfied, it will generate many false positives in the maritime environment.

This paper proposes a novel spectral–spatial signature analysis-based ship detection method in midwave hyperspectral images, instead of the conventional spatial or temporal methods. Based on the combustion process, CO2 emissions show a double-spike signature at around 4.16 μm (2400 cm−1) and 4.34 μm (2300 cm−1). The hot CO2 emission band is broader than that which is absorbed by normal atmospheric CO2, which generates CO2-double spikes. The directional-mean subtraction filter (D-MSF) detects each CO2 spike, and final targets are detected by joint analysis of the detection of both spikes. So, we call the proposed method carbon dioxide-double spike (CO2-DS) detection.

The contributions from this paper can be summarized as follows. First, the phenomenon of the CO2-double spike emission in the midwave band is analyzed systematically. Second, a novel spectral–spatial analysis is proposed to detect remote ships in the maritime environment based on this analysis. Third, CO2-DS detection can suppress background clutter (cloud, sea glints) completely, which leads to an extremely low false positive rate in remote-ship detection. Finally, the performance of CO2-DS detection is demonstrated in a real sea environment with real ships.

The remainder of this paper is organized as follows. Section 2 explains the proposed CO2-DS method, including analysis of carbon dioxide-double spike radiation. Section 3 evaluates the clutter suppression performance of the CO2-DS method in the maritime environment. The paper concludes in Section 4.

#### **2. Proposed Ship CO**<sup>2</sup> **Plume Detection**

#### *2.1. Signature Analysis of Carbon Dioxide-Double Spikes*

Before explaining the details of the proposed CO2-DS, the double-spike phenomenon in midwave spectral bands should be introduced [38,39]. Figure 1 summarizes the overall mechanism of CO2-double spikes in received spectral radiance from ship plumes. Hot CO2 emits thermal energy in wider spectral band than that of normal air CO2, as shown in Figure 1a. Atmospheric spectral transmittance shows strong absorption at 2300–2380 cm−<sup>1</sup> by normal CO2 in the atmosphere, as shown in Figure 1b. Therefore, the received spectral radiance in a TELOPS Fourier transform infrared (FTIR) camera [40] shows a double-spike signature: the first spike is around 2276 cm−<sup>1</sup> (4.39 μm), and the second spike is around 2393 cm−<sup>1</sup> (4.18 μm).

**Figure 1.** Operational concept of remote ship detection and the double spike-signature generation mechanism of hot carbon dioxide: (**a**) spectrum of a hot CO2 emission by a ship, (**b**) spectral atmospheric transmittance, and (**c**) received spectral radiance of hot CO2 from a ship's plume.

The measurement of CO2-double spikes can be derived from radiative transfer calculated with Equation (1). We adopted the radiative-transfer equation used in the Moderate-Resolution Atmospheric Radiance and Transmittance (MODTRAN) [41]. In general, at-sensor received radiance in the midwave infrared (MWIR) region consists of CO2-emitted radiance, transmitted background radiance, and total atmospheric path radiance (thermal+solar components).

$$L\_{\text{target}}(\lambda) = \pi(\lambda) \left[ \varepsilon(\lambda) B\_{\text{CO}\_2}(\lambda, T\_{\text{CO}\_2}) + (1 - \varepsilon(\lambda)) L\_{\text{bg}}(\lambda) \right] + L\_{\text{s}}^{\uparrow}(\lambda) + L\_{\text{atm}}^{\uparrow}(\lambda) \tag{1}$$

*Ltarget*(*λ*) is the observed at-sensor target radiance; *λ* is wavelenth; *ε*(*λ*) is spectral CO2 gas emissivity; *BCO*<sup>2</sup> (*λ*, *TCO*<sup>2</sup> ) is the spectral radiance of the hot CO2 gas, assuming a blackbody in the Planck function with the combusted CO2 gas temperature (*TCO*<sup>2</sup> ); *τ*(*λ*) is the spectral atmospheric transmittance, and *L*<sup>↑</sup> *<sup>s</sup>* (*λ*) and *L*<sup>↑</sup> *atm*(*λ*) are the spectral upwelling solar and thermal path radiance, respectively, reaching the sensor.

According to the spectral emissivity of combusted CO2 gas [38], *ε*(*λ*) is relatively high. Therefore, the term *Lbg* can be removed. According to MWIR radiometric characteristics [42], the contribution of solar radiance, *L*<sup>↑</sup> *<sup>s</sup>* (*λCO*<sup>2</sup> ), from air scattering is very small, even for very dry conditions (less than 2% at 5 μm) [42]. Ignoring the solar term, we can simplify Equation (1) to Equation(2):

$$L\_{\text{target}}(\lambda) = \tau(\lambda) \left[ \varepsilon(\lambda) B\_{\text{CO}\_2}(\lambda, T\_{\text{CO}\_2}) \right] + (1 - \tau(\lambda)) B\_{\text{atm}}(\lambda, T\_{\text{atm}}) \tag{2}$$

where the definition of thermal upwelling, *L*<sup>↑</sup> *atm*(*λ*), is replaced by the multiplication of blackbody radiation of the atmosphere, *Batm*(*λ*, *Tatm*) and 1 − *τ*(*λ*). The observed target spectral radiance, *Ltarget*, directly depends on both emissivity *ε*(*λ*) of hot CO2, and atmospheric transmittance, *τ*(*λ*), that should be analyzed specifically.

In addition, neighboring pixels to target hot CO2 has only atmospheric spectral path radiance at wavelenth *λ*. Therefore, the received background signal can be expressed with Equation (3):

$$L\_{m\text{-}high}(\lambda) = (1 - \tau(\lambda))B\_{atm}(\lambda, T\_{atm})\tag{3}$$

If we apply center-surround difference by subtracting Equation (3) from Equation (2), which is the same as a mean subtraction filter (MSF) operation for a specific band image, the target-only signal can be extracted with Equation (4). This physical property is used in the D-MSF to extract the target (hot CO2) signature:

$$L\_{\text{targetOnly}}(\lambda) = L\_{\text{target}}(\lambda) - L\_{\text{neighbor}}(\lambda) = \pi(\lambda) \left[ \varepsilon(\lambda) B\_{\text{CO2}}(\lambda, T\_{\text{CO2}}) \right] \tag{4}$$

Conventional ships are powered by diesel engines, and the corresponding hydrocarbon combustion produces water vapor, H2O) and CO2 as expressed in Equation (5) [39,43]:

$$\rm C\_nH\_{2n+2} + (3n+1)O\_2 \to (n+1)H\_2O + nCO\_2 \tag{5}$$

The first important parameter is spectral emissivity *ε*(*λ*) in Equation (2). The normal temperature of combusted CO2 gas from ships is approximately 600 K by infrared signature suppression (IRSS) [44]. The emission band-width of CO2 gas is proportional to the temperature of the CO2 gas [39]. Figure 2 demonstrates this phenomenon by changing CO2 temperature using the high-resolution transmission molecular absorption database (HITRAN) simulation on the web (http://hitran.iao.ru/molecule/). As the CO2 temperature increases, the absorption band range is broader. The absorption band range of normal atmospheric temperature (300 K) is 2300–2380 cm−<sup>1</sup> (or 4.20–4.35 μm), where the inverse acts as atmospheric transmittance. The spectral unit can be either the wavelength (*λ*) [μm] or wavenumber (*k*) [cm<sup>−</sup>1] by unit selection (*k* = 10,000/*λ*). For hot CO2, the spectral absorption coefficient is the same as the spectral emission coefficient , *ε*(*λ*), in Equation (2).

**Figure 2.** Spectral absorption coefficient variations of carbon dioxide according to different gas temperatures.

The second important parameter is spectral atmospheric transmittance *τ*(*λ*) in Equation (2). In the MODTRAN simulation in the MWIR band, the spectral transmittance of the CO2 band (2320–2375 cm−<sup>1</sup> or 4.21–4.31 μm) decreases abruptly with distance, as shown in Figure 3. The average transmittance in the CO2 band is 0.13, 0.005, and 0 at 5, 20, and 100 m, respectively. Note that the atmospheric transmittance band of zero value coincides with the spectral absorption coefficient of CO2 at 300 K in the top part of Figure 2. In addition, it remains only the atmospheric path radiance in Equation (2) in the CO2 absorption band, *τ*(*λ*) = 0), where there is no target or background signals.

**Figure 3.** Spectral transmittance according to ship-camera distance in the MWIR band: (**a**) 5, (**b**) 20, and (**c**) 100 m.

Multiplication of the spectral emissivity of hot CO2 and atmospheric transmittance produces double spikes, as shown in the bottom portion of Figure 1. The spectral feature where a hot CO2 emission band is broader than that which is absorbed by normal atmospheric CO2 generates CO2 double spikes. The midwave spectral information is acquired via TELOPS MWIR hyperspectral camera [40]. It can provide calibrated spectral radiance images with a high spatial and spectral resolution from a Michelson interferometer in the short-wave to midwave bands (1.5–5.6 μm).

#### *2.2. CO*2*-DS-Based Ship Plume Detection*

Motivated by the double-spike phenomenon of hot CO2, a novel ship-detection method is proposed, as shown in Figure 4. The proposed method is called CO2-DS because it is based on the carbon dioxide-double spike feature. The proposed CO2-DS approach consists of three steps: spectral band selection, spatial small-target detection, and the final ship detection. The spectral band selection sets the search range in the spectral domain focusing on the CO2 absorption band. Then, each CO2 spike is detected by a carefully designed spatial filter (the directional-mean subtraction filter). Finally, the target ship is confirmed using a joint detection operation.

**Figure 4.** Overall processing flow of remote-ship detection using spectral–spatial information of hot CO2 plumes.

In Step 1, the specific spectral band range should be defined to attain successful remote-ship detection with extremely few false positives. The proposed CO2-DS method is based on the atmospheric CO2 absorption band (2320–2375 cm−<sup>1</sup> or 4.21–4.31 μm) that should be included. In addition, double-spike bands emitted by hot CO2 should be included. The first spike range is below 2300 cm<sup>−</sup>1, and the second spike range is above 2380 cm<sup>−</sup>1.

The details of Step 2 and Step 3 in Figure 4 are described in Figure 5. Each spike signal is detected by a parallel search algorithm: the first spike detection and the second spike detection.

**Figure 5.** Details of the proposed CO2-DS: first spike detection, second spike detection, and joint detection.

The starting wavenumbers are 2300 cm−<sup>1</sup> and 2380 cm−<sup>1</sup> because there is no signal at all except for atmospheric path radiance. Given a wavenumber, a test image of this band is probed by D-MSF-based small infrared target detection [4]. D-MSF is adopted in this paper because it is a well-proven method and robust against thermal noise in the maritime environment. Figure 6 summarizes the detailed flow of the D-MSF graphically. Given a hypothesized spectral band image, *Iλ*(*x*, *y*), a mean subtraction filter, *h*(*x*, *y*), is used to enhance the signal-to-clutter ratio (1). Then, directional local median *I<sup>D</sup> <sup>λ</sup>* (*x*, *y*) is used to estimate row directional background (2) from the MSF result, *IMSF <sup>λ</sup>* (*x*, *y*). Subtracting (2) from (1) in Figure 6 produces a clearer signature image. The region of interest (ROI) is generated by an initial thresholding (*th*1) and eight-nearest neighbor-based clustering. Final detection is achieved by applying adaptive thresholding *th*2 in constant false alarm rate (CFAR) detector to the subtracted image with the given ROI. It can provide a signal-to-clutter ratio (SCR) as well as the ROI of the target.

If there is no detection, the wavenumber (*k*1) of the first spike is decreased by Δ, and that of the second spike is increased by Δ, the value of which is 13 cm−<sup>1</sup> in this paper. If there is a detection, each neighboring band image is probed additionally to obtain a maximum signal score.

**Figure 6.** Operational flow of the directional-mean subtraction filter (D-MSF) for small-target detection in a band image.

Finally, in Step 3, the joint detection process is performed from the first and second spike detections: [*x*1, *y*1, *w*1, *h*1, *SCR*1], [*x*2, *y*2, *w*2, *h*2, *SCR*2], where [·] denotes [column, row, width, height, SCR] of the detected target. If two detections exist, then the final target information is merged with Equation (6). If the merged SCR is larger than a predefined threshold (*th*3), then it is declared a final target.

$$Final\,\,\text{target}\,\,\text{info} = \left[\frac{\mathbf{x}\_1 + \mathbf{x}\_2}{2}, \frac{y\_1 + y\_2}{2}, \frac{w\_1 + w\_2}{2}, \frac{h\_1 + h\_2}{2}, (SCR\_1 + SCR\_2)\right] \tag{6}$$

Figure 7 shows the spectral profile of a remote ship's plume and an enlarged view focusing on the atmospheric CO2 absorption band. The numbers (1), (2), (3), and (4) represent the probing band locations of the proposed CO2-DS method; (1) and (2) belong to the first spike band, and (3) and (4) belong to the second spike band. Note that the first spike is not clear in remote cases but it does not matter in our CO2-DS method since D-MSF is used. Figure 8 shows the corresponding band images, and the max SCR is used to select a band image for each spike: (2) and (3) are selected. The bottom image in Figure 8 represents the result of joint detection from the detection results of the first spike and the second spike.

**Figure 7.** Spectral profile of a remote ship's plume with the probed band positions: first spike: (1) and (2); second spike: (3) and(4).

**Figure 8.** Probed band images of first and second spikes corresponding to (1) and (2), and (3) and (4), respectively. The final target is detected via joint detection from the selected max signature for each spike.

#### **3. Experimental Results**

#### *3.1. Experiment 1: Analysis of Signature Variation*

MWIR hyperspectral images of remote ships were acquired with the open path TELOPS Hyper-Cam Mid-Wave Extended (MWE) model [40]. It could provide calibrated spectral radiance images with a high spatial and spectral resolution from a Michelson interferometer in the shortwave to midwave band (1.5–5.6 μm). The spatial image resolution was 320 ×256, and the spectral resolution was up to 0.25 cm−1. The noise-equivalent spectral radiance (NESR) was 7 [nW/(cm<sup>2</sup> · sr · cm−1)], and the radiometric accuracy was approximately 2 K. The field of view was 6.5 × 5.1 degrees.

**Figure 9.** Comparison of CO2 signatures in terms of distance: (left) near CO2 (78 m), (right) remote CO2 (2 km).

In the first experiment, the effect of signature variation was compared, depending on target distance. Figure 9 shows the spectral profiles of hot CO2 plumes, both near (78 m) and remote (2 km). A strong double-spike phenomenon can be observed for near CO2, as shown in Figure 9 (left side) due to relatively weak signal attenuation. However, a remote target does not show such a distinctive double-spike pattern as seen in Figure 9 (right side), which makes it clear that only the spectral profile-based approach was not effective at detecting targets. This conclusion was confirmed through a real experiment applying the well-known spectral angle mapper-based detector [45] to the spectral profiles. Figure 10 shows the spectral profile-based detection results. The left column is the training region; the middle column is a representative spectral profile after training (just the mean of spectral profiles); and the right column is SAM-based detection results indicated by red dots. For the near target, the hot CO2 plume region could be detected correctly with a similarity threshold of 0.8, as shown in Figure 10a. On the other hand, the remote target signature was unclear, which led to false positive detections with a similarity threshold of 0.98, as shown in Figure 10b.

#### *3.2. Experiment 2: Detection Parameter Analysis*

In the second experiment, key parameters of the proposed CO2-DS method were analyzed and evaluated. There are four parameters; spectral band range for probing, initial threshold (th1) and CFAR threshold (th2) in the D-MSF, and final threshold (th3) from joint detection. The spectral band range for robust ship detection should be determined analytically based on Figure 11. Basically, the atmospheric CO2 absorption band should be included. Since the air temperature is usually 300 K, this band is fixed at [2320–2380 cm−1] as shown in the bottom part of Figure 11. The first-spike band starts at 2300 cm−<sup>1</sup> and ends at 2180 cm−<sup>1</sup> considering band offset. The second-spike band starts at 2380 cm−<sup>1</sup> and ends at 2450 cm−<sup>1</sup> considering band offset. The band images below 2180 cm−<sup>1</sup> or above 2450 cm−<sup>1</sup> show background clutter with target information that should be excluded to ensure extremely few false positives.

**Figure 11.** Analysis of spectral band ranges for ship detection by visualizing corresponding band images.

The selection of threshold1 (*th*1) is used to determine an ROI in order to investigate whether it is a target or not, as shown in Figure 6. Figure 12 shows the image for *IMSF <sup>λ</sup>* (*x*, *y*), where the bright blob represents the hot CO2 target, and others are background pixels. The 3D surf views of the target patch and background patch show the signal levels, and the lower left graph of Figure 6 displays the signal intensity profile corresponding to the dotted red line in the image. Note that the max signal is <sup>20</sup> <sup>×</sup> <sup>10</sup>−5[W/m2 ,sr,cm−1] and the max background signal is 5 <sup>×</sup> <sup>10</sup>−5[W/m<sup>2</sup> , sr,cm−1]. Considering the safety margin of factor 2, *th*1 was set at 1.0 <sup>×</sup> <sup>10</sup>−4[W/m<sup>2</sup> , sr,cm−1]. The threshold 2 parameter (*th*2) can be set by analyzing the SCR values for different targets, as shown in Figure 13. Most SCRs of a spike band image are larger than 12.8, so, *th*2 can be set at 10.0 for a safe detection capability. The last threshold (*th*3) was varied to calculate the receiver operating characteristic (ROC) curve in order to compare different target-detection methods.

**Figure 12.** Parameter analysis of threshold 1 (*th*1) from directional mean subtraction image.

**Figure 13.** Parameter analysis of threshold 2 (*th*2) from statistics of SCR values of different targets.

#### *3.3. Experimental 3: Performance Evaluation*

Based on these parameter analyses, the final detection performance was evaluated by comparing the proposed CO2-DS, D-MSF [4], and the high-boost multi-scale contrast measure (MLCM) [46]. MLCM is based on the human visual system, and greatly improves the detection rate for small infrared targets. The test images consist of 53 hypercubes acquired by a TELOPS MWIR hyperspectral camera. Three different ships with distance ranges between 1.5 km and 4 km were considered. For the baseline methods, test images were prepared by making broad band images in the spectral range of 3.5–5.6 μm. The ROC metric was used for a fair detection performance based on false positives per image (FPPI). Figure 14 summarizes detection performance. The proposed method (the red solid line) showed

an ideal detection curve. MLCM ranked second, and D-MSF showed the worst result. Visual inspection was conducted at the indicated FPPI (0.65/image). With a normal homogeneous background, the three methods showed reasonable working performance, as seen in Figure 15. Magenta circles represent ground truth, cyan rectangles are detection results from the proposed CO2-DS, and yellow rectangles represent baseline detection results. If there was strong cloud clutter, D-MSF generated many false alarms, as shown in Figure 16. In a strong sea-clutter environment, MLCM detected a lot of sea-glints, as shown in Figure 17. Note that the proposed CO2-DS showed stable detection results without any false positives from sky or sea clutter. In addition, the ship was located at 4 km (near the horizon) and was detected successfully by our proposed method.

**Figure 14.** Receiver operating characteristic (ROC)-based comparison of detection performance.

**Figure 15.** Comparative detection examples for a homogeneous background: (**a**–**c**) detection results, (**d**–**f**) filtered maps.

**Figure 16.** Comparative detection examples for cloud-clutter backgrounds: (**a**–**c**) detection results, (**d**–**f**) filtered maps.

**Figure 17.** Comparative detection examples for sea-clutter backgrounds: (**a**–**c**) detection results, (**d**–**f**) filtered maps.

#### *3.4. Experiment 4: Detection Performance Analysis*

It can be useful to analyze the effect of weather condition and target distance to the detection performance. The proposed CO2-DS method is based on Equation (4), especially atmospheric transmittance (*τ*(*λ*)). Therefore, if we use the MODTRAN-based atmospheric transmittance calculation and measured background noise, we can conduct the Monte Carlo simulation. Figure 18 represents the MODTRAN-based atmospheric transmittance calculation by changing distance at a specific wavenumber (2390.16 cm−1) or wavelength (4.184 μm). The atmosphere model was selected either "MidLatitude Summer" or "MidLatitude Winter" to consider the effect of humidity and temperature. The distance length changed from 0.1 km to 6.0 km with 0.1 km interval.

**Figure 18.** Moderate-Resolution Atmospheric Radiance and Transmittance (MODTRAN) environment for atmospheric transmittance calculation.

The proposed CO2-MS detector uses target-background difference method (mean subtraction filter). Therefore, the path radiance term can be removed. Figure 19 shows a filtered map for the spectral band image at 2393 cm−1. The maximum target signal-background at 2 km is 6.1742 <sup>×</sup> <sup>10</sup>−4[W/m<sup>2</sup> · sr · cm−1]. If one uses transmittance of 0.18 at 2 km, the original signal-background is 3.4 <sup>×</sup> <sup>10</sup>−3[W/m2 · sr · cm−1]. Background noise can be modeled as the rectified Gaussian noise model with 0 mean and standard deviation of 1.3964 <sup>×</sup> <sup>10</sup>−5[W/m2 · sr · cm−1]. Th2 (SCR) is set as 10 in this simulation. Figure 20 shows the atmospheric transmittance and Monte Carlo simulation-based detection performance according to the target distance and weather conditions. Normally, the transmittance in mid-latitude summer is lower than that in mid-latitude winter because of higher water vapor contents (higher humidity), which leads to shorten the detection range as Figure 20b. If we fix target distance at 5 km, the detection rate is reduced from 95.4% to 5.2%.

The proposed CO2-MS detection method shows extremely low false positives in ship detection. In addition, CO2-MS is based on the thresholding and clustering. So, it can detect multiple targets. Because the CO2-MS can use both the spectral and spatial information, the size of the ship (fishing boat, military ship, or cargo ship) can be identified if each ship has different type of fuel, engine, gas volume, and temperature. However, it has several limitations such as expensive spectral measurement device, weak to weather condition such as high humidity.

**Figure 19.** MODTRAN environment for atmospheric transmittance calculation.

**Figure 20.** (**a**) Atmospheric transmittance vs. distance, (**b**) Monte Carlo simulation-based detection results.

#### **4. Conclusions**

Remote-ship detection is important in various applications, such as navigation and surveillance. Midwave infrared image-based ship detection is deployed due to its long-range detection capability in the maritime environment. However, MWIR images show strong response to cloud clutter and sea glint, which generate many false positives, even in state-of-the-art methods. This paper proposed a completely different approach by fusing spectral and spatial information in the MWIR carbon dioxide absorption band with the double-spike phenomenon. The atmospheric CO2 absorption band (2320–2375 cm−1) can block all background radiance, and remote hot CO2 can penetrate through the spectral double-spike band, which is detected by the proposed carbon dioxide-double spike algorithm. Experimental results in a real cluttered sea environment with remote ships showed excellent detection performance with extremely few false positives if the ships' diesel engines are working.

**Author Contributions:** The contributions were distributed between authors as follows: S.K. (Sungho Kim) wrote the text of the manuscript and programmed the hyperspectral detection method. J.S., J.A. and S.K. (Sunho Kim) provided the midwave infrared hyperspectral database, operational scenario, performed the in-depth discussion of the related literature, and confirmed the accuracy experiments that are exclusive to this paper. All authors have read and agree to the published version of the manuscript.

**Funding:** This research was funded by ADD grant number UE191095FD, 2020 Yeungnam University Research Grants, NRF (NRF-2018R1D1A3B07049069) and the APC was funded by NRF.

**Acknowledgments:** This study was supported by the Agency for Defense Development (UE191095FD). This work was supported by the 2020 Yeungnam University Research Grants. This research was also supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of Education (NRF-2018R1D1A3B07049069).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Aircraft and Ship Velocity Determination in Sentinel-2 Multispectral Images**

#### **Henning Heiselberg**

National Space Institute, Technical University of Denmark, 2800 Kongens Lyngby, Denmark; hh@space.dtu.dk; Tel.: +45-4525-9760

Received: 31 May 2019; Accepted: 26 June 2019; Published: 28 June 2019

**Abstract:** The Sentinel-2 satellites in the Copernicus program provide high resolution multispectral images, which are recorded with temporal offsets up to 2.6 s. Moving aircrafts and ships are therefore observed at different positions due to the multispectral band offsets, from which velocities can be determined. We describe an algorithm for detecting aircrafts and ships, and determining their speed, heading, position, length, etc. Aircraft velocities are also affected by the parallax effect and jet streams, and we show how the altitude and the jet stream speed can be determined from the geometry of the aircraft and/or contrail heading. Ship speeds are more difficult to determine as wakes affect the average ship positions differently in the various multispectral bands, and more advanced corrections methods are shown to improve the velocity determination.

**Keywords:** Sentinel-2; multispectral; temporal offsets; ship; aircraft; velocity; altitude; parallax; jet stream

#### **1. Introduction**

Surveillance for marine and air space situation awareness is essential for monitoring and controlling traffic safety, piracy, smuggling, fishing, irregular migration, trespassing, spying, icebergs, shipwrecks, the environment (oil spill or pollution), etc. Dark ships and aircrafts are non-cooperative vessels with non-functioning transponder systems such as the automatic identification system (AIS) for ships or automatic dependent surveillance (ADS-B) for aircrafts. Their transmission may be jammed, spoofed, sometimes experience erroneous returns, or simply turned off deliberately or by accident. Furthermore, AIS and ADS-B land based and satellite coverage is sparse at sea and at high latitudes. Therefore, other non-cooperative surveillance systems as satellite or airborne systems are required.

The Sentinel-2 satellites under the Copernicus program [1–3] carry multispectral imaging (MSI) instruments that provide excellent and freely available imagery with pixel resolutions down to 10 m. The orbital period is 5 days between the Sentinel-2 (S2) satellites A + B, but as the swaths from different satellite orbits overlap at higher latitudes, the typical revisit period for each satellite is two or three days in Europe and almost daily in the Arctic. S2 MSI has the potential to greatly improve the marine and airspace situational awareness, especially for non-cooperative ships and aircrafts.

Ship detection, recognition, and identification in optical satellite imagery has been studied in a number of papers with good results [4–10]. The resolution and sensitivity are generally better and the number of multispectral bands is larger, but clouds reduce the continuous coverage. Ship positions, and their length, breadth, form and heading can be determined accurately and the multispectral reflections can be fingerprints for ID. Ship speeds have only been determined from satellite imagery in a few cases where Kelvin wakes are observed [7].

Detection, tracking and speed determination of vehicles on ground has been studied in video sequences and images recorded with time intervals fx. by change detection. Heights of tall buildings or altitudes of clouds [11–13] and other static or slow moving objects have been determined by shadow lengths or parallax methods, when the images are recorded from a flying platform as an aircraft, drone or satellite with time delay imaging or from two platforms.

Recently, parallax effects were observed for aircrafts and their condensation trails (contrails) in Sentinel-2 color images, referred to as "plainbows" [14]. The work presented here is novel when it comes to exploiting the temporal offsets in Sentinel-2 MSI, and specifically for determining aircraft velocities and altitudes, and ship speeds. We consider this work as the first analysis of such effects due to temporal offsets, as we have not been able to find any studies with scientific analyses or applications. In this work, we outline the basic physics for moving objects in satellite multispectral images with temporal offsets, the parallax effects and influence of jet streams. We primarily consider aircrafts and ships, but the analysis also applies in principle to all kinds of moving vehicles including cars and helicopters, and also clouds, auroras, etc. The basic formulas are derived, and as proof of principle, we show a number of representative examples for both aircrafts and ships. From Sentinel-2 multispectral images with known temporal offsets we calculate the resulting speed, heading, altitudes, etc. Subsequently, we test our results by comparing to data from the navigation systems ADS-B or AIS.

In Section 2 we analyze how the temporal offset affects moving objects in S2 multispectral images and include the parallax effect for aircrafts as well as the effect of jet streams. After a description of the S2 MSI offset times and resulting apparent velocities, we calculate aircraft velocities and altitudes from S2 multispectral images in Section 3. In Section 4 we turn to ship velocities, which are not affected by parallax but are slower and more difficult to determine accurately due to long wakes. We describe a simple but effective correction method, which improves the speed calculations considerably.

In the summary and outlook we suggest ways to improve the model by better position determination and understanding of spectral dependence of object reflections [15]. In addition, calculations for a large number of ships and aircrafts are required for better statistics and improving the model by fine-tuning the parameters. This could also lead to an annotated database useful for machine learning methods.

#### **2. Satellite Images and Method of Analysis**

The S2 multispectral images were analyzed using dedicated software developed specifically to detect ships and aircrafts in large images with different pixel resolutions and determine their precise position and orientation as described in [7,8]. As detection is not the focus of this paper, we mainly describe the MSI temporal offsets and how they affect the multispectral images for velocity calculations. A flow chart for the algorithm is shown in Figure 1.

**Figure 1.** Flow chart of algorithm: For a given scene the m = 0, ... , 9 high resolution S2 multispectral images are selected. Search and detection of objects above a threshold T is done for the red band m = 3. For each object, a region of interest (RoI) around the object center <sup>→</sup> *c* <sup>3</sup> is chosen for all ten images, where m = 4–9, are resampled to 10 m pixel size. The center <sup>→</sup> *<sup>c</sup> <sup>m</sup>* and length <sup>→</sup> *Lm* are calculated for m = 0, ..., 9. From 2D linear regression (3) of the variance <sup>σ</sup>, the apparent velocity <sup>→</sup> *V* is found and inserted in (9–15), whereby the aircraft velocity and altitude or ship velocity is found. By comparing to ADS-B or AIS, the velocities, altitudes and positions are validated.

#### *2.1. Sentinel-2 Multispectral Images*

S2 carries the Multispectral Sensor Imager [1–3] that records images in 13 multispectral bands (see Table 1) with different resolutions and time offsets. As we are interested in small object detection and tracking, we focused on analyzing the high resolution images, the 4 bands with 10 m and the 6 bands with 20 m pixel resolution. These are mega- to giga-pixel images with 16 bit grey levels.


**Table 1.** Spectral bands for the Sentinel-2B Multispectral Imager. The 10 high resolution bands *m* = 0, ..., 9 are ordered according to temporal offset [1].

We analyzed several S2 level 2A processed images from 2019 covering Copenhagen airport in Denmark, see Figure 2, and the Channel between Amsterdam and London. These images are convenient because several aircrafts are usually present after takeoff or before landing. In addition, there are a larger number of ships present in the strait of Øresund surrounding the airport and in the Channel.

**Figure 2.** East part of Denmark in the red band B4 (m = 3) from 5 April 2019, 10:30 a.m. UTC. The red box is a ROI in the strait of Øresund surrounding Copenhagen airport.

In the S2 multispectral images *Im*(*i*, *j*), the spatial coordinates <sup>→</sup> *r* = (x, y) are the pixel coordinates (*i*,*j*) multiplied by the pixel resolution *l* = 10 m, 20 m or 60 m as given in Table 1 for the 13 bands. The 10 high resolution multispectral images with pixel size 10 m or 20 m are ordered according to temporal offset *tm*, *m* = 0, 1, ..., 9. As shown in Table 1, they range from 0 s to 2.085 s in temporal offset. Due to the odd and even detector array in MSI, the offsets are either delayed or advanced, respectively. The imaging sequence is such that the offsets are reversed in stripes along track within the image [13].

#### *2.2. Object Detection and Position Determination*

To detect an object, its reflection must deviate from the background. For proving the principle of velocity determination, we chose for simplicity a region of interest with sea background, which is usually darker than the objects and therefore makes detection easier. The multispectral variant background over land will require a more elaborate detection algorithm, but has the potential of determining velocities of driving vehicles as well.

When the sea covers more than half of the image after land removal, the median reflection value provides an accurate and robust value for the average background. For detection, we chose the red band *m* = 3 (see Figure 2) because it has high resolution and average time offset such that temporal offset objects will appear around the red center (see Figures 3–5). In addition, solar reflections from ships and aircrafts generally have high contrast in red with respect to the sea background. For each object a small region, e.g., 100 × 100 pixels or smaller, is extracted around the central object coordinate, such that it covers the object extent including movement, wakes or contrails. The same 1 km × 1 km region is then extracted for the 10 high resolution bands *m* = 0, 1, 2, 3 with spatial resolution 10 m and *m* = 4, ..., 9 with 20 m. The latter are corrected for the different resolutions. For each band the median value is chosen as the background.

The pan-sharpening technique [7,16] for increasing the resolution in lower resolution images can only be applied to static images. As moving objects change pixel position in the multispectral images due to the temporal offset, we could not apply pan-sharpening in this analysis.

For each multispectral image, the object is defined spatially by the pixels with reflections above the background value plus a threshold T, which depends on target type as will be discussed below. The central object position <sup>→</sup> *c* <sup>m</sup> = (x, y), length *Lm*, breadth *Bm* and orientation/heading angle θ*<sup>m</sup>* are calculated in each band *m* = 0, ..., 9 by weighting the object pixels with their reflection *Im*(*i*, *j*) and calculating the first moments in *i* and *j*, as described in detail in [7].

Unfortunately, ship wakes and aircraft contrails can corrupt the position determination considerably. Both generally move the central position backwards with respect to vessel direction by an amount that varies with the band. At the same time, the object length is extended. We corrected for this effect to the first order by adding the distance from the average center position to the ship front, which is half of the object length *Lm*, in the vessel heading direction, i.e., the vector <sup>→</sup> *Lm* = *Lm*(cos θ*m*, sin θ*m*).

$$
\overrightarrow{\dot{r}}\_{\text{m}} = \overrightarrow{\dot{c}}\_{\text{m}} + \frac{1}{2}\overrightarrow{L}\_{\text{m}} \tag{1}
$$

This position is now at the front of the object in each band *m* as shown in the images below.

#### *2.3. Multispectral Temporal O*ff*sets and Velocity Determination*

We define the apparent velocity as the change in position as observed in the multispectral images divided by the band dependent time delay. An object moving with apparent velocity <sup>→</sup> *V* will ideally be recorded in band *m* at position

$$
\overrightarrow{r}\_m \sim \overrightarrow{r}\_V + \overrightarrow{V} \cdot \mathbf{t}\_m \tag{2}
$$

here, <sup>→</sup> *r <sup>V</sup>* is the vessel position at zero temporal offset—ideally the blue band *m* = 0 for which *t*<sup>0</sup> = 0. For ships, the apparent velocity <sup>→</sup> *V* is simply the vessel speed and direction, whereas for aircrafts,

the parallax effect due to satellite motion must be included as will described below. Currents and jet streams also influence V.

For example, the aircrafts shown in Figures 4 and 5 fly with apparent speeds V ~ 200 m/s, and move a distance of ~400 m or 40 pixels in the time interval of 2.085 s. Consequently, the aircrafts show up as ten "pearls on a string" when the high resolution multispectral bands are plotted all together in a (false) color image.

The object positions <sup>→</sup> *r* m that are calculated for each multispectral image *m* = 0, ..., 9 will generally scatter around the linear prediction of Equation (2). We define the variance as the mean square average of distance deviations

$$
\sigma^2 = \frac{1}{10} \sum\_{m=0}^{9} \left( \overrightarrow{r}\_m - \overrightarrow{r}\_V - \overrightarrow{V} \cdot t\_m \right)^2 \tag{3}
$$

By minimizing this variance, which is equivalent to a two-dimensional linear regression, we obtain the best fit values for vessel position <sup>→</sup> *<sup>r</sup> <sup>V</sup>* and the apparent velocity vector <sup>→</sup> *V* = *V*(cos θ*V*, sin θ*V*). An estimate for the uncertainty in apparent velocity is the lowest standard deviation in distance divided by the temporal interval

$$
\sigma\_V = \sigma / \mathbf{t}\_\theta \tag{4}
$$

Typically, the positions are accurate up to a pixel size *l* or less in each multispectral image. The two-dimensional linear regression fit of <sup>→</sup> *r <sup>m</sup>* vs. *tm* therefore has a standard deviation less than σ ∼ *l*/ √ 10. Dividing by the temporal offset interval *t*<sup>9</sup> = 2.085 s, we obtain the uncertainty for the apparent velocity of σ*<sup>V</sup>* = *l*/*t*<sup>9</sup> √ 10 ~5–10 kph. This is comparable to speeds of slow ships, whereas typical aircraft cruise velocities are 800–1000 kph, and gives a relatively accurate aircraft velocity determination.

#### *2.4. Comparison to AIS and ADS-B*

By matching positions of aircrafts from ADS-B and ships from AIS at the same time and positions, we can identify the vessels and compare size, heading, velocities and altitudes. Unfortunately, the positioning systems sometimes lag or the updating is delayed or is infrequent. Therefore, we find that ships and aircraft do not always match precisely at the correct position and time. In addtion, the S2 overflight time included in the file name is not the local image recording time. In Figure 2, the Copenhagen regions are recorded five minutes later than the time 10 h 30 m 29 s given in the file name, and the local recording time is 14 s later from north to south for the descending track.

#### **3. Aircraft Velocities**

For low flying aircrafts and ships the parallax effect is negligible and the apparent velocity <sup>→</sup> *V* is just the aircraft velocity <sup>→</sup> *VAC*. For high altitude aircraft, we need to consider the satellite orbit, velocity and viewing angle in order to correct for the parallax effect. In addition, the jet stream must be considered as it affects the contrails.

**Figure 3.** (**a**) Illustration of the satellite-aircraft-ground parallax effect. (**b**) Earth map with projected satellite orbits. The orbits cross the vertical lines (latitudes φ) at angles θ*S*2.

#### *3.1. Satellite Direction w.r.t. Ground*

We used standard mathematical and celestial convention where angles are measured from the equator counter-clockwise. In navigation, angles are measured from the North Pole clockwise and thus differ by 90◦ and angular direction.

The S2 satellites fly in a sun-synchronous orbit at mean altitude HS2 = 786 km with speed VS2 = 7.44 km/s. Their polar orbit is slightly retrograde, descending with inclination angle i = −98.62◦ on the dayside. Due to Earth's curvature, the orientation of the satellite track θ*S*<sup>2</sup> with respect to latitude φ is (see Figure 3b)

$$\cos\theta\_{S2} = \cos i / \cos \phi \tag{5}$$

At the equator, φ = 0◦ and θ*S*<sup>2</sup> = *i*, but at maximum polar S2 latitude φ = 180◦ − *i* = 81.38◦, the S2 satellite flies straight west, i.e., θ*S*<sup>2</sup> = ±180◦. For the images around Copenhagen, φ 55◦ and we find <sup>θ</sup>*S*<sup>2</sup> −105◦. At Amsterdam, <sup>φ</sup> 52, 5◦ and <sup>θ</sup>*S*<sup>2</sup> −104◦. The resulting satellite velocity is <sup>→</sup> *VS*<sup>2</sup> = *VS*2(cos θ*S*2, sin θ*S*2), w.r.t. ground.

#### *3.2. Parallax E*ff*ect*

The satellite movement during the multispectral temporal delays causes a parallax effect (see Figure 3a) that moves objects at an altitude H northeast-wards in direction −θ*S*<sup>2</sup> with respect to the ground. Stationary objects such as tall buildings, clouds, balloons, stalling aircrafts, etc. will be moved by an apparent velocity <sup>→</sup> *V* = − → *VS*<sup>2</sup> · *H*/*HS*<sup>2</sup> with respect to ground due to their parallax. Determining *V* from a linear regression of Equation (3), we find the object altitude

$$H = \frac{V}{V\_{S2}} \cdot H\_{S2} \tag{6}$$

The parallax has recently been exploited for determining altitudes and movement of, for example, volcanic plumes [11,12].

The parallax effect moves the flight paths opposite to the satellite direction, and separates each multispectral band such that an aircraft (with its contrails) appears as a multispectral rainbow, when plotted in a false color image as shown in Figures 4 and 5. Previously [14], the RGB contrails were named "plainbows". We named our ten multispectral contrails as a "planebows". Contrails are usually observed at high altitudes 7.5–12 km.

#### *3.3. Moving Objects*

When the object moves, its velocity must be added to the apparent velocity. In the absence of wind, a moving object, such as an aircraft at altitude *HAC* with velocity <sup>→</sup> *VAC* = *VAC*(cos θ*AC*, sin θ*AC*), will appear to have velocity

$$
\overrightarrow{V} = \overrightarrow{V}\_{\text{AC}} - \overrightarrow{V}\_{\text{S2}} \cdot \frac{H\_{\text{AC}}}{H\_{\text{S2}}} \tag{7}
$$

or 
$$\begin{pmatrix} \cos \theta\_{AC} \\ \sin \theta\_{AC} \end{pmatrix} V\_{AC} = \begin{pmatrix} \cos \theta\_V \\ \sin \theta\_V \end{pmatrix} V + \begin{pmatrix} \cos \theta\_{S2} \\ \sin \theta\_{S2} \end{pmatrix} V\_{S2} \cdot \frac{H\_{AC}}{H\_{S2}} \tag{8}$$

The heading of the aircraft θ*AC* is given by the aircraft orientation angles θ*m*, *m* = 0, ..., 9 as described in Section 2.2 and [7]. The aircraft heading angle θ*AC* must be determined independently for the ten multispectral bands θ*m*. The four high resolution bands generally provide a consistent and robust average aircraft heading angle θ*AC*. When contrails are visible (see Figure 3b), the contrail and aircraft directions are the same θ*CT* = θ*AC*, and they provide a more accurate heading.

From the two equations in (8), we obtain the aircraft velocity

$$V\_{AC} = \frac{\sin(\theta\_V - \theta\_{S2})}{\sin(\theta\_{AC} - \theta\_{S2})} \cdot V \tag{9}$$

and the aircraft altitude

$$H\_{\rm AC} = \frac{\sin(\theta\_V - \theta\_{\rm AC})}{\sin(\theta\_{\rm AC} - \theta\_{\rm S2})} \cdot \frac{V}{V\_{S2}} \cdot H\_{\rm S2} \tag{10}$$

These relations can also be obtained from simple sine relations for the triangles in Figures 4 and 5. When the aircraft and satellite directions are parallel, θ*S*<sup>2</sup> = θ*AC*, the aircraft velocity and altitude cannot be determined separately. At zero altitude (as will be discussed below for ships) there is no parallax effect so that θ*AC* = θ*<sup>V</sup>* and *VAC* = *V*.

Note that Equations (9) and (10) are invariant to the orientation of the coordinate system as only relative angles appear. They are also clock- vs. counter clockwise invariant.

**Figure 4.** Planebows. Scale in pixels. (**a**) A slow and low flying aircraft near Amsterdam airport after takeoff recorded in 10 multispectral images with offsets. Red, green and blue are true colors whereas the seven remaining colors are overlayed with false colors, e.g., m = 1 is yellow. Numbers 0, ..., 9 indicate the central aircraft position in each band (but moved 10 pixels down). The triangle shows the apparent (<sup>→</sup> *V*), aircraft (<sup>→</sup> *VAC* ) and satellite (<sup>→</sup> *V G <sup>S</sup>*<sup>2</sup> <sup>=</sup> <sup>→</sup> *VS*2*HAC*/*HS*<sup>2</sup> ) velocity vectors. (**b**) Planebow of a fast flying aircraft at high altitude near Amsterdam with strong contrails from each wing motor.

In Figure 4a, a slow and low flying aircraft is shown after take-off from Amsterdam. From linear regression of Equation (3), we find its apparent velocity V = 586 ± 2 kph and direction θ*<sup>V</sup>* = 128◦. The aircraft orientation angle θ*AC* = 133◦ is determined by the object orientation. From Equations (9) and (10), we find an aircraft speed of *VAC* = 552 kph, at altitude *HAC* = 1.764 m. According to live flight tracking, ADS-B, the aircraft at that time and position, has a speed of 507 kph at altitude 1.875 m. Considering delays in the flight tracking updates and neglecting wind speeds, the agreement is fair.

In Figure 4b another aircraft with strong contrails is captured near Amsterdam, where we expect little jet stream. Thus the long contrails show the aircraft heading θ*CT* = θ*AC* = −134◦. The contrails do, however, corrupt the determination of the aircraft central positions, and we therefore remove them in the object images by setting the threshold T sufficiently high above the contrail but below the aircraft reflections. The central object positions in each band can then be used for determining the aircraft positions <sup>→</sup> *r <sup>m</sup>* as shown in Figure 4b. The resulting apparent velocity is *V* = 557 ± 8 kph, and the apparent direction θ*<sup>V</sup>* = −155◦. From Equations (9) and (10), we find an aircraft speed of *VAC* = 855 kph and altitude *HAC* = 11.231 m. According to the flight tracking system, an aircraft at that time and position has speed *VAC* = 833 kph and altitude *HAC* = 11.582 m, both in good agreement with our calculations.

The large apparent velocities *V* can in turn be exploited for a search for aircrafts in S2 images as they are the only fast moving objects.

#### *3.4. Jet Stream and Contrails*

The polar jet stream circulates eastward in a meandering way as illustrated in Figure 5a. It lies between latitudes 50–60◦ at altitudes 9–12 km, and are a few hundred km wide. The jet typically has a speed of ~100 kph but can exceed 400 kph. Flight information systems show that aircrafts benefitting from the jet stream typically fly a hundred kph faster east- than westwards.

**Figure 5.** (**a**) Illustration of Earth's meandering polar and subtropical jet streams. (**b**) An aircraft over Copenhagen where the contrails are affected by the polar jet stream (see text).

The jet stream (and winds in general) with velocity <sup>→</sup> *VJS* will sweep an aircraft and its contrails along. The aircraft orientations θ*<sup>m</sup>* are the same as the contrail direction θ*CT*. The aircraft orientation in the image no longer matches the true aircraft heading θ*AC* with respect to ground as shown in Figure 5b. If the contrail angle is used in Equations (9) and (10) for the aircraft angle, it can erroneously lead to supersonic aircrafts flying well above the altitude of most commercial airliners. Such extreme

values only apply to the Concorde and a few other aircrafts. Fighter jets can be excluded in Figure 5b as the size of the aircraft is too large.

In order to correct for wind speeds, we introduce two auxillary quantities (see Figure 5b), namely the velocity in the contrail direction

$$V\_{\rm CT} = \frac{\sin(\theta\_V - \theta\_{S2})}{\sin(\theta\_{\rm CT} - \theta\_{S2})} \cdot V \tag{11}$$

and the jet stream correction to the parallax effect

$$
\Delta = \frac{\sin(\theta\_{IS} - \theta\_{CT})}{\sin(\theta\_{CT} - \theta\_{S2})} \cdot V\_{JS} \tag{12}
$$

The aircraft altitude is then (see Equation (10))

$$H\_{AC} = \frac{H\_{S2}}{V\_{S2}} \left(\frac{\sin(\theta\_V - \theta\_{CT})}{\sin(\theta\_{CT} - \theta\_{S2})} \cdot V - \Delta\right) \tag{13}$$

The aircraft velocity can be determined from

$$V\_{AC}^2 = V\_{CT}^2 + \Delta^2 - 2V\_{CT} \cdot \Delta \cdot \cos(\Theta\_{CT} - \Theta\_{S2}) \tag{14}$$

Finally, the aircraft heading angle can be determined from Equation (9).

In Figure 5b the aircraft engines under each wing create two long contrails with angle θ*CT* = −24◦. The apparent velocity is V = 1070 ± 2 kph with direction θ*<sup>V</sup>* = 2◦. From Equation (11) we find *VCT* = 1035 kph. Assuming that the jet stream heads east θ*JS* = 0◦ with speed *VJS* = 200 kph, we obtain Δ = 82 kph, *HAC* = 11.514 m, *VAC* = 1.025 kph and θ*AC* = −19◦. The aircraft speed relative to the jet stream is thus 825 kph. These numbers are compatible with normal aircraft cruising altitude and speed, however, the unknown jet stream speed was fitted.

#### **4. Ship Velocities**

Low objects on the surface of the Earth, such as ships, have no parallax. In addition, ocean currents are usually slow and can be neglected. The apparent velocity is then simply the ship velocity

$$V\_{Slip} = \text{ } V \tag{15}$$

Likewise, the vessel heading angle is the apparent direction θ*<sup>V</sup>* and is the same for all bands. By overlaying all multispectral images as shown in Figures 6 and 7, the heading angle θ*<sup>V</sup>* can be determined accurately, especially for large ships and/or when ship wakes are long.

Ships sail much slower than aircrafts, typically ~10 m/s (19.4 knots or 36 kph). Therefore, the ships in Figures 5 and 6 move less than two pixels in the temporal offset interval *t*<sup>9</sup> = 2.085 s, which requires accurate determination of ship positions.

#### *4.1. Short Wakes*

When ship wakes are short, they have less effect on the estimated central positions and both the central and the corrected positions of Equation (1) can be used for determining the ship velocity. An example is shown in Figure 6 where <sup>→</sup> *c* <sup>m</sup> and <sup>→</sup> *r* m are plotted with black and red numbers, respectively. Both are temporally ordered correctly and yield the ship velocity *V* = 15 ± 1 kph. According to AIS, a ship at that time and position is sailing at a speed of 15 kph in the same direction.

**Figure 6.** (**a**) Ship with short wake. As all ten multispectral ship images almost overlap, they appear white. The central positions shown with black numbers *m* = 0, ..., 9 are expanded in (**b**). The temporal ordering appears approximately correct.

#### *4.2. Wake Corrections*

Unfortunately, ship wakes can be longer than the ship and we find that they can corrupt the position determination considerably. Ship wakes move the apparent central position backwards with respect to vessel direction by an amount that typically is larger for the lower wavelengths. This is observed in Figure 7 where the ordering is not temporal but rather spectral, i.e., according to band wavelength as: *m* = 0, 2, 3, 4, 6, 7, 1, 8, 5, 9 (see Table 1). Wake reflection seems to decrease gradually towards the infrared [9,10,15]. If one simply uses the central positions for determining V, one obtains a ship speed of 69 ± 9 kph, which is much too large. Using the front positions of Equation (1), their ordering is closer to temporal and the ship speed is only 30 ± 4 kph. According to AIS, a ship at that time and position has a speed of 25 kph in the same direction. The length correction to central positions not only improves the accuracy of the velocity determination but also reduces the uncertainty.

The uncertainty in velocity is typically 5–10 kph, which sets the limit on how slow ship speeds can be determined. To improve the position accuracy, the more advanced calculation of <sup>→</sup> *r* m and study of the spectral dependence of wakes is required.

The two examples in Figures 6 and 7 may give the wrong impression, that faster ships create longer wakes, which is not always the case. Wake length may depend on e.g., ship type and size, surface winds and background. The temporal offsets are therefore useful complementary but are also partly correlated information.

**Figure 7.** (**a**) As Figure 6 but for a fast ship creating a long wake, which corrupts the temporal ordering. Adding ship lengths Lm to central position as shown with red numbers in (**b**) improves the temporal ordering and velocity determination.

#### *4.3. Kelvin Waves*

Kelvin waves from large and fast ships are occasionally observed in S2 images [7]. A sailing ship creates Kelvin waves bounded by cusp-lines separated by an angle of ±arcsin(1/3) = ±19.47◦ on each side of the ship and its wake. The Kelvin wavelength is related to the ship speed *V* as [17]

$$
\lambda = \frac{2\pi V^2}{g} \tag{16}
$$

where *g* = 9.8 m/s2 is the gravitational acceleration at the surface of the Earth. The wavelength can be determined by a Fourier analysis of the image, and the ship speed follows from Equation (16).

In Figure 8, 10 Kelvin waves are observed within ca. 50 pixels, i.e., λ = 50 m, and we obtain the ship speed V = 8.8 m/s or 32 kph. The apparent ship speed from temporal delays is V = 29 ± 3 kph, when corrected for wakes. According to AIS, a ship at that position and time is sailing north with 30 kph in good agreement with both satellite results.

The Kelvin waves are stationary relative to the ship, i.e., travel with ship velocity as seen from the satellite. The time delay of 2 s between the first and last multispectral corresponds to ~20 m in this case or almost half a wavelength. The Kelvin wave in the last image is therefore in anti-phase with the first and tends to interfere destructively if the multispectral images are added or plotted on top of each other. The Fourier transform of the individual band images is therefore optimal. It also improves the transform if the ship is masked, which is straight forward to do as the ship central position, width, breadth and heading are known.

The time offsets also cause a Kelvin wave phase shift in the multispectral bands

$$
\phi\_m = \frac{2\pi V t\_m}{\lambda} = \frac{g t\_m}{V} \tag{17}
$$

Slow ships result in large phase shifts. The multispectral phases can also be found in the Fourier analysis.

**Figure 8.** Ship with up to a dozen observable Kelvin waves in its wake.

#### **5. Discussion**

The above results can be taken as proof of principle for our novel method utilizing the temporal offsets in the Sentinel-2 multispectral imager for determining velocities and altitudes of moving objects. For aircrafts in particular, the accuracy is excellent; a few kph compared to cruise velocities of the order of 1000 kph, and a few hundred meters uncertainty in altitude compared to the standard 10 km cruise altitude of commercial airliners. Contrails improve heading direction determination but can also debase the positions unless the threshold is set correctly between aircraft and contrail reflections. Jet streams when present were shown to affect velocities and altitudes significantly. Unfortunately, the jet stream velocity cannot be determined from the S2 images but requires separate atmospheric information or alternatively, the jet stream can be estimated by requiring the aircraft to fly with standard cruise speed or altitude of commercial airliners.

The method was also shown to apply to ships with similar uncertainty of a few kph, which is sufficient for fast ships but comparable to slow ships. However, wakes cause a serious systematic error as they corrupt the temporal towards spectral ordering. Utilizing measurements of ship length and adding them can correct for part of this error. Yet, a much better understanding of the correlations between wakes, speed and spectral reflections is required and needs further investigation.

Our method is limited in the sense that it only utilizes the central position and length of objects for each band. The image contains much more spectral information on the object extent and form that is not utilized. Yet, it works surprisingly well even for small and slow objects such as ships with a variety of complex wakes.

Moving objects over land was not considered in this work. The more complex background will require better algorithms for removing the background in each spectral image, which is outside the scope of this work.

#### **6. Summary and Outlook**

We have described the basic physics behind moving objects in satellite multispectral images with temporal offsets, parallax effects and influence of jet streams. The basic formulas were derived and as proof of principle, a number of representative examples were shown for aircrafts and ships. The analysis serves as a proof of principle and provides a working model.

From apparent velocities the resulting aircraft speed, heading, and altitudes were calculated accurately and compared to data from the navigation systems ADS-B with good agreement. Jet streams can influence aircraft speeds and altitudes and the jet stream velocity must be determined independently or fitted.

Ship velocities are not affected by parallax but difficult to determine accurately for slow ship speeds or when long wakes are present. We described a simple but effective correction method, which improves the calculation of ship speeds considerably when compared to AIS. The detailed influence of thresholds, backgrounds, object lengths and contrails and wake reflections in the different multispectral bands should be studied in more detail in order to further improve the position determination. In addition, wake lengths may depend on e.g., ship type and size, surface winds and background. The temporal offsets are useful complementary information and the correlation to wake and ship lengths and widths provides additional information. Comparison to wake detection and velocity determination in Synthetic Aperture Radar radar images [18] should also be studied.

For better statistics, a large number of ships and aircrafts is required, where velocities and altitudes of aircrafts and ships are calculated and compared to AIS and ADS-B data with improved trajectory prediction [19]. This would also allow for improving the model by fine-tuning and optimizing parameters such as thresholds, better wake corrections and possibly introduce non-equal weights in the linear regression analysis of Equation (3). The large set of ships and aircrafts would also build an annotated database necessary for training machine learning algorithms [20–23]. Convolutional neural networks could be trained on this database so as to attempt to extract many more parameters and possibly refine the estimation of altitudes and velocities.

**Funding:** This work received no external funding.

**Acknowledgments:** We acknowledge the Arctic Command Denmark for support and interest. Images contain modified Copernicus Sentinel data from 2019 [1].

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
