Image velocimetry for clouds with relaxation labeling based on deformation consistency

Cloud tracking has been used to observe atmospheric motion. It has been widely used for operational weather forecasting (e.g. Schmetz et al (1993)), and also it is the most used technique to estimate atmospheric motion on the planets covered with thick clouds such as Venus, Jupiter, and Saturn since the 1980s (e.g. Limaye and Suomi (1981), Ingersoll et al (1981) and Sromovsky et al (1983)). The task of cloud tracking is basically the same as that of particle image velocimetry (PIV), so its techniques have been adapted in many ways (e.g. Asay-Davis et al (2009)). Cloud tracking is usually used to estimate horizontal wind, since clouds are mainly advected horizontally. This treatment is especially valid when features having large horizontal scales are treated, since the vertical extent of atmospheric motion is generally limited geometrically and by stratification.

Despite the long tradition, cloud tracking is still a demanding and often difficult task, so manual tracking by human operators' eyes is still widely used especially for Venus (e.g. Hueso et al (2012) and Khatuntsev et al (2013)). Unlike PIV in a laboratory, frame rate is limited for satellite cloud observation. It is especially the case for planetary satellites like Akatsuki because of the long distances and the limitation of satellite payload. Factors that make planetary cloud tracking difficult include instrumental noise or error, positioning error, cloud morphology itself (sometimes it is featureless or only has blurred large-scale features), deformation, formation, diminishment, and fluctuation by traveling atmospheric waves. As for the Earth, sparseness of clouds is also a factor, which makes clouds at different altitudes often observed nearby or overlappingly.

The present paper reports a cloud-tracking algorithm developed for the Venus orbiter Akatsuki, which reached Venus and started observation on Dec 7, 2015 (Nakamura et al 2016). Our tracking is based on a cross-correlation method by (Ikegawa and Horinouchi (2016), hereinafter IH16), in which multiple cross-correlation surfaces (CCSs) are overlaid to reduce erroneous template matching and increase the signal-to-noise (S/N) ratio. In this study, the IH16 method is extended by combining a new formulation of the relaxation labeling process to correct erroneous match. We further propose post processes to eliminate or to detect velocities that are not associated with the flow.

Figure 1(a) shows one of the images used in this study. It was taken by the ultraviolet imager (UVI) onboard Akatsuki with the filter for 365 nm. Although the camera shots by $1024\times 1024$ pixels, shown there is a sub-image of $200\times 200$ pixels, since the orbiter is at relatively long distance from Venus. Akatsuki has four imaging cameras, and three of them have multiple filters, enabling to provide images of clouds at different altitudes on day-side and/or night-side. See Nakamura et al (2011) for more on the observation by Akatsuki.

Zoom In Zoom Out Reset image size

The relaxation labeling, also called the relaxation method (RM), is a class of mechanisms to assign labels to objects or scenes, which has been explored extensively since the study of Rosenfeld et al (1976). The RM has been applied to correct erroneous template matching in optical flow detection and particle tracking velocimetry (PTV) (Wu 1995, Baek and Lee 1996, Ohmi and Li 2000, Evans 2000), and it is also used in our previous studies (Ogohara et al 2012, Kouyama et al 2012). In the RM for template matching, the candidates of displacement vectors, which are associated with local peaks in CCSs, are first computed, and selection is made by evaluating the closeness of the vectors over neighboring grid points (here, the expression grid points refers to the centers of template regions at which motion vectors are defined). Like these studies, we use the simple cross-correlation method that does not treat image deformation. In this study, we introduce to the RM a measure to evaluate the consistency with the neglect of deformation, which we shall call the deformation consistency. It gives a quantitative basis for the scoring of 'no-match' cases; unless the no-match possibility (i.e. the possibility that none of the candidates associated with a template represent the actual displacement) is considered, its impact is not large. Although we have not tried, the present method is considered suitable for the first stage of the two-stage PIV in which the affine transformation is considered at the second stage (e.g. Scarano (2002)). An important study to note in this respect is Tokumaru and Dimotakis (1995). They proposed a method to estimate motion by using two successive images, where deformation is explicitly treated. The method is obtained by modifying a variational method pioneered by Horn and Schunck (1981), and it can optimize the mapping between two images through a global constraint. On the other hand, our method utilizes local consistency to correct (if possible) or reject erroneous template matching.

Our method secures the overall consistency, and it is sufficient in some cases. However, it does not necessarily eliminate the inconsistency among all of the nearby grid-point pairs. Therefore, we propose an optional post process to eliminate the inconsistency entirely. We also propose another post process to detect velocities that are not associated with the flow. The two post processes are closely related. The former classifies the motion vectors into groups, within each of which the consistency is perfect, and the latter expand the groups. Their important application is to find waves in fluids such as atmospheric gravity waves and to quantify their phase velocities, as shown in what follows. Also, to separate flows at different altitudes is another possible application.

In the relaxation labeling, 'labels' are assumed to be independent. It implies that the candidates of displacement vectors should represent isolated or well-separated peaks of CCSs. However, this requirement and actual methods to select independent candidates are not found in previous studies on the RM for template matching. We discuss this issue and present a method to find well-separated peaks of a surface. Limiting candidates to the highest points of well-separated peaks is a prerequisite for the post processes stated above, and it is also required in the precision estimation proposed by IH16, which is also used to process data from Akatsuki.

The rest of the paper is organized as follows. Section 2 briefly introduces the particular cross-correlation method used in this study. The new RM formulation is presented in section 3. Section 4 presents the necessity and an algorithm to find well-separated peaks of CCS. Section 5 presents our post processes. Results using Akatsuki images are shown in section 6, and conclusions are drawn in section 7.

In this section, we briefly review the cross-correlation method by IH16 and the error estimation methods therein, on which the present study is based. Minor updates made are also noted here. IH16 used ultraviolet images from the Venus Express orbiter, which was operated from 2006 to 2014. IH16 did not use RM.

In their method, multiple CCSs are superposed before deriving cloud motion vectors (CMVs). The superposition is conducted by averaging CCSs in two ways. One is by superposing CCSs over nearby grid points, which is similar to the one introduced by Hart (2000) but for using summation (to average) rather than multiplication. The particular superposition used in IH16 and in this study as well is to use five CCSs, one at the middle and the adjacent four. This superposition is referred to as the spatial smoothing of CCSs. In addition, they used multiple (more than three) images obtained subsequently, for which interpolation is needed before averaging. CCSs are calculated over the image combinations as in Sciacchitano et al (2012) but for considering the movement of template centers expanding with time. In the implementation in IH16, an approximation to consider uniform motion was used (equation (10) of IH16), but in this study, the full treatment based on equation (9) of IH16 is used.

IH16 proposed two methods to estimate uncertainties. One of them utilizes the acuteness of CCS peaks, where the region of a peak above the lower bound of confidence associated with its summit is regarded as statistically 'indistinguishable'. Many of clouds on Venus do not exhibit sharp boundaries, so the precision of tracking depends on cloud size, or cloud feature scale. The greater the scale is, the more difficult it is to track precisely, and at the same time, the flatter the CCS peak gets. The present method exploits this relation. As in IH16, a parameter termed to indicate the precision is derived for each CMV. For an idealized CCS that have a single-peaked Gaussian shape, is inversely proportional to the 'peak to correlation energy ratio' (PCE) used by Xue et al (2014).

As in figures 1(b) and (c), clouds on Venus often have elongated streaky features especially at mid- to high-latitude. The CCS peaks associated with the movement of such a feature is also elongated. The formulation by IH16 explicitly treats it by fitting CCS peaks by elliptic paraboloids (the parameter mentioned above is derived by combining multiple parameters that provide directional information; see IH16).

Our method was implemented by enhancing the cloud tracking system developed by Ogohara et al (2012), and we use the pre-process therein as follows. The UV brightness is first interpolated onto longitude-latitude grid points with a resolution of 0.125°. It is then divided by cosine of the solar incident angle at each grid point, and a two-dimensional Gaussian band-pass filter is applied. The filter uses the convolution with Gaussian functions in the form of $\exp \left[\frac{\left(\lambda -{{\lambda}_{i}}\right){{}^{2}}{{\cos}^{2}}{{\phi}_{i}}}{2{{\sigma}^{2}}}+\frac{\left(\phi -{{\phi}_{i}}\right){{}^{2}}}{2{{\sigma}^{2}}}\right]$, where ${{\lambda}_{i}}$ and ${{\phi}_{i}}$ are longitude and latitude of a grid point and λ and ϕ are longitudinal and latitudinal variables, respectively. The band-pass filter is realized by low-pass filtering with $\sigma ={{0.25}^{\circ}}$ and high-pass filtering with $\sigma ={{3}^{\circ}}$.

Following a convention, the RM is formulated as follows:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{k}}(i):=\frac{{{p}_{k}}(i){{q}_{k}}(i)}{\underset{j=0}{\overset{{{I}_{k}}}{\sum}}\,{{p}_{k}}\left(\,j\right){{q}_{k}}\left(\,j\right)}, \end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{q}_{k}}(i)\equiv \underset{l\in {{G}_{k}}}{\sum}\,\underset{{{i}^{\prime}}=0}{\overset{{{I}_{l}}}{\sum}}\,{{r}_{kl}}\left(i,{{i}^{\prime}}\right){{p}_{l}}\left({{i}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 2 }$$

Here, the symbol $:=$ means iterative updates, and p_k(i) is the score of ith candidate at the grid point k (the center of a template region). We assign i  =  0 to the no-match (no solution) case and i  =  1, 2,.., I_k to the actual displacement vector candidates associated with the peaks of the CCS for k; I_k is the number of the candidates. The score is calculated in terms of compatibility with the candidates at neighboring grid points, denoted symbolically as a set of grid points G_k. The term ${{r}_{kl}}\left(i,{{i}^{\prime}}\right)$ represents the compatibility between the ith candidate at the grid point k and the ${{i}^{\prime}}$th candidate at the grid point l. The formulation of ${{r}_{kl}}\left(i,{{i}^{\prime}}\right)$, which we introduce below, characterizes the RM. The initial value of p_k(i) is not very important, and it can be set equally to $\frac{1}{{{I}_{k}}+1}$. In our actual implementation, it is set by using cross-correlation values as in Evans (2000).

Generally in the RM, the score of one of the candidates at each grid point approaches to 1 through iteration, and the scores of the other candidates approach to 0. Therefore, as underscored by Hummel and Zucker (1983), p_k(i) cannot be regarded as probability, even though it is normalized so that the summation over i is 1. This nature implies that the results of the RM is not very sensitive to the details in the formulation of ${{r}_{kl}}\left(i,{{i}^{\prime}}\right)$ or its parameter values as long as it is based on the closeness between the displacement vectors.

However, care is needed when the no-match possibility is considered, since there is competition between match and no-match, which must be evaluated somehow. In our formulation, no-match cases are incorporated as i  =  0. Therefore, the task is to define the compatibility between match and no-match cases (i  =  0 and ${{i}^{\prime}}\ne 0$) in a comparable manner to the compatibility between the regular match and match cases ($i\ne 0$ and ${{i}^{\prime}}\ne 0$).

We now introduce our formulation of compatibility based on the deformation consistency. For convenience, ${{r}_{kl}}\left(i,{{i}^{\prime}}\right)$ is first separated as

$$\begin{eqnarray}&&{{r}_{kl}}\left(i,{{i}^{\prime}}\right)\equiv {{w}_{kl}}\,{{c}_{kl}}\left(i,{{i}^{\prime}}\right),\end{eqnarray} \tag{ 3 }$$

where w_kl is a weight that is optional, and ${{c}_{kl}}\left(i,{{i}^{\prime}}\right)$ is the actual compatibility dependent on i and ${{i}^{\prime}}$. It is not needed to introduce the weight, but following the web site www.skycoyote.com/FITSFlow/ and Kouyama et al (2012), it is introduced as ${{w}_{kl}}=\frac{s}{{{D}_{kl}}}$, where D_kl is the distance between the grid points k and l in terms of the number of CMV grid points (D_kl  =  1 if the points k and l are next to each other), and s is the normalization factor derived to make ${\sum}_{l\in {{G}_{k}}}{{w}_{kl}}=1$ just for cleanliness. The weight w_kl is introduced in order to weaken the contribution to equation (1) from grid points that are relatively far from the point k. The weight is irrelevant to the following argument, since it is constant for a given combination of k and l.

The compatibility among match and match cases (i  >  0 and ${{i}^{\prime}}>0$) is formulated as follows:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{c}_{kl}}\left(i,{{i}^{\prime}}\right)\equiv \exp \left(-a\frac{\parallel \boldsymbol{d}_{{k}}(i)-\boldsymbol{d}_{{l}}\left({{i}^{\prime}}\right){{\parallel}^{2}}}{L_{kl}^{2}}\right). \end{array}\end{eqnarray} \tag{ 4 }$$

Here, $a\equiv \ln 2$, $\boldsymbol{d}_{{k}}(i)$ represents the ith displacement vector at k (measured in the same way as D_kl), $\parallel \parallel$ expresses the L₂ norm of vector (thus, if $\boldsymbol{d}=\left({{d}_{x}},{{d}_{y}}\right)$, $\parallel \boldsymbol{d}\parallel =\sqrt{d_{x}^{2}+d_{y}^{2}}$), and L_kl is the half-length half-maximum specified as

$$\begin{eqnarray}&&{{L}_{kl}}=\alpha {{D}_{kl}},\end{eqnarray} \tag{ 5 }$$

where α is a tunable parameter, which is set to 0.5 in our implementation.

Equation (5) is regarded as a simple realization of the deformation consistency. When the candidates i and ${{i}^{\prime}}$ are selected at the grid points k and l, respectively, non-zero $\boldsymbol{d}_{{k}}(i)-\boldsymbol{d}_{{l}}\left({{i}^{\prime}}\right)$ indicates expansion-contraction and/or rotation of the line segment between these points, indicating deformation including rotation, over the time between the initial and end times of tracking. Since we use the simple cross-correlation method with no affine transformation, we cannot track features if deformation is too large, so the combination of displacement vectors that indicates a too large deformation should be rejected as being inconsistent. Equations (4) and (5) mean that ${{c}_{kl}}\left(i,{{i}^{\prime}}\right)$ is 0.5 when the length of $\boldsymbol{d}_{{k}}(i)-\boldsymbol{d}_{{l}}\left({{i}^{\prime}}\right)$ is equal to $\alpha {{D}_{kl}}$, α times the initial distance between the points. Since $0\leqslant {{c}_{kl}}\left(i,{{i}^{\prime}}\right)\leqslant 1$, we express that the compatibility is neutral if it is 0.5.

Given the compatibility formulated as above, it is reasonable to give a neutral compatibility when no-match cases are considered, so naïvely it can be set as follows:

(Type A)

$$\begin{eqnarray}&&{{c}_{kl}}\left(i,{{i}^{\prime}}\right)=0.5\quad \text{if}~i=0~\text{or}~{{i}^{\prime}}=0.\end{eqnarray} \tag{ 6 }$$

In this case, q_k(0) is always 0.5 throughout the iteration. On the other hand, q_k(i) for i  >  0 tends to be smaller than that at the beginning of iteration, since ${{p}_{l}}\left({{i}^{\prime}}\right)$ is more or less uniform over ${{i}^{\prime}}$, and positive support (meaning ${{c}_{kl}}\left(i,{{i}^{\prime}}\right)>0.5$) is expected to occur at limited ${{i}^{\prime}}$s. Therefore, the initial support for the no-match possibility can be too much. Hence, we devised another formulation:

(Type B)

$$\begin{eqnarray}{{c}_{kl}}\left(i,{{i}^{\prime}}\right)=\left\{\begin{array}{*{35}{l}} 0.5 & \text{if}~i=0~\text{and}~{{p}_{l}}\left({{i}^{\prime}}\right)~\text{is}~\text{maximum}~\text{over}~{{i}^{\prime}}=0,1,...{{I}_{l}}, \\ 0 & \text{otherwise}~\text{and}~\text{if}~i=0~\text{or}~{{i}^{\prime}}=0. \end{array}\right. \end{eqnarray} \tag{ 7 }$$

In some cases, it is conceivable that this formulation results in the reduction of no-match scores too fast, so it is recommended to modify the iterative updates as follows:

(Amendment of equation (1) for Type B)

$$\begin{eqnarray}&&{{p}_{k}}(0):=\max \left[\frac{{{p}_{k}}(0){{q}_{k}}(0)}{\underset{j=0}{\overset{{{I}_{k}}}{\sum}}\,{{p}_{k}}\left(\,j\right){{q}_{k}}\left(\,j\right)},\,\text{Second}\,\text{greatest}\,{{p}_{k}}(i)\,\text{over}\,i=1,..,{{I}_{k}}\right].\end{eqnarray} \tag{ 8 }$$

This treatment ensures that the decrease of the no-match score is not faster than that of the second-likely match score over the course of iteration.

As mentioned earlier, the final state of the iteration equation (1) is normally p_k(i)  =  1 for one of i in 0, 1,.., I_k, so equations (6) and (7) are equivalent at that state, meaning that a solution of another is a solution of the other. However, there is a hysteresis with the initial value, and in general the number of no-match results becomes greater with equation (6) than with equation (7).

As stated in section 2, we use the superposition of CCSs obtained from multiple images. In this case, the time intervals of image pairs are not constant. However, CCSs are interpolated onto the lags for the greatest time interval. It indicates that the deformation consistency is measured uniformly in terms of the greatest time interval.

The relation between the present formulation and those in earlier studies can be summarized as follows. The RM by Wu (1995) is the closest among the studies introduced in section 1. The study consolidated the no-match treatment as in this study and used equation (6). However, deformation is not considered there, and the compatibility parameter becomes weaker than neutral only when the directions of velocities differ by more than 90°. Therefore, it is not Galilean invariant; adding a uniform flow changes the result. The other studies do not make the consolidation. However, since there are tunable parameters in Ohmi and Li (2000) (e.g. R_c), it might be possible to realize a similar effect to ours. However, it is not clear how we can achieve it. The no-match possibility is not considered in Evans (2000), Ogohara et al (2012) and Kouyama et al (2012).

The candidates in the RM for template matching are associated local peaks of CCSs. However, to the authors' knowledge, how to find them has not been described in previous papers. In practice, a simple nine-point comparison is used in our earlier implementation (Ogohara et al 2012, Kouyama et al 2012), although it is not written in these papers. There, if a cross-correlation value at a CCS grid point is greater than those at the neighboring eight points, it is classified as a local peak. It appears that the same or similar algorithms are used in many studies.

This simple method would be practically sufficient if CCS peaks are peaky enough. However, it may capture local bumps on a ridge as in figures 2(a) and (b). The defect of this method is especially evident in figure 2(b). In this example, the number of candidates is limited to nine, so no candidate is assigned to the secondly peak of the CCS. Note that this kind of limitation is employed in many implementations.

Zoom In Zoom Out Reset image size

The elongated shape of the primary peak in figure 2(b) is due to the streaky feature of high-latitudinal clouds. In this case, cloud motion is well-defined only to the direction perpendicular to the streak. Therefore, to set candidates along the ridge itself is reasonable. The problem, however, is that it is done in an ad hoc way; the distribution and the number of local peaks varies significantly case by case and are not adequately constrained. The proper way to treat this situation should be either (1) to treat the peak as a continuous region and to develop a theoretically consistent treatment or (2) to set only one candidate for each well-separated peak and evaluate the directional uncertainty. Here, we choose (2) and leave (1) out of scope. Note that a method to evaluate the directional uncertainty is proposed by IH16. Note also that the option (1) means that the candidates resulting from its discretization are not independent, which must be taken into account in RM.

4.1. Formulation and algorithm

We define a peak as well-separated from others if it has a deep-enough skirt that is not connected to any of higher peaks. Here, the criterion for being deep-enough is that the depth from the local maximum is greater than a specified parameter d. It is illustrated in figure 3 for one-dimensional cases. The points shown by filled circles are the highest points, called as 'summits' here, of the well-separated peaks, and those shown by cross marks are local bumps because their skirts are shallower than d, so they are regarded as a part of peaks. In other words, d is the minimum gap depth required to separate peaks.

Zoom In Zoom Out Reset image size

An algorithm to implement it is described in terms of programming as follows. The procedure takes a two-dimensional array holding 'height' data, which is the cross-correlation values in the context of the present study. In the following description, the term grid points refers to array indexes.

• 1.  
  Make a list of the grid points sorted in the descending order by height values.
• 2.  
  Repeat the following until a sufficient number of peaks are found or all the grid points are classified (or, if a threshold is provided, all the unclassified grid points have heights lower than it).

  • (a)  
    Find all the grid points that have not been classified and are connected to the list's first (i.e. highest) unclassified point (referred to as P) and have height values greater than or equal to h  −  d, where h is the height at P. Here, a group of points are regarded as connected if they are contiguous through vertically and/or horizontally adjacent (next) points (the criterion of connection can be altered).
  • (b)  
    If the grid points are connected to an existing peak, they are classified as its skirt. Otherwise, they are classified as a new peak, and the point P is recorded as its summit location.

We provide a source code in C language as a supplementary material of this paper (stacks.iop.org/MST/28/085301/mmedia). It is efficient, where connection finding is implemented with a recursive call to avoid unnecessary array scan; note that efficiency is important, since this procedure is conducted for each of all CCSs. The program returns not only the summit points but also the map of peaks, which is used for uncertainty evaluation to limit the region of surface fitting.

In our implementation, d is set to a constant (0.05), since the results were insensitive to its choice in our case. However, one can elaborate it, for instance, by using the confidence interval, which is a function of the peak correlation value.

The results of the RM can be further improved and/or used to draw additional information. Here, we describe two methods that exploit the deformation consistency. The both can be implemented as post processes.

5.1. Elimination of inconsistent CMVs / grouping consistent CMVs

5.2. Wave detection / flow detection at multiple altitudes

The proposed approach described in previous sections is applied to Akatsuki UVI images obtained on Dec 9, 2015. The size of template images are set to ${{6}^{\circ}}\times {{6}^{\circ}}$. The set of neighboring grid points used for the RM, G_k, consists of the eight surrounding points. The CCS peaks treated are limited to those having values greater than 0.5. When the no-match possibility is considered, the type B compatibility equation (7) and the associated modification equation (8) are applied.

Figure 4 shows the results obtained using the three 365 nm images taken at 14, 16, 18 UTC on Dec 9, 2015, so three combinations of images (14–16, 14–18, 16–18 h) are used to superpose CCSs. The spatial smoothing of CCSs is also used. Therefore, each CMV is derived from a synthetic CCS obtained from 15 CCSs (3 image pairs  ×  5 nearby CCSs by the spatial smoothing).

Zoom In Zoom Out Reset image size

Figure 4(a) shows CMVs obtained without considering the no-match possibility. It suggests that the flow is predominantly westward, which is known as the super-rotation of the Venusian atmosphere (see e.g. Hunten et al (1983), for review). To ease comparison of CMVs, a constant value of 100 m s⁻¹ is added and shown in figure 4(b). The CMVs generally exhibit gradual spatial variation, but at some places, changes are abrupt (e.g. vectors around $\left(\lambda,\phi \right)\simeq \left({{105}^{\circ}},-{{10}^{\circ}}\right)$, $\left({{100}^{\circ}},{{28}^{\circ}}\right)$, and $\left({{117}^{\circ}},-{{20}^{\circ}}\right)$, where λ and ϕ are longitude and latitude, respectively). The RM incorporating the no-match possibility (section 3) eliminates these vectors as shown in figure 4(c).

The elimination is easily understood by using figure 5, which visualizes the deformation indicated by the CMVs at around $\left({{105}^{\circ}},-{{10}^{\circ}}\right)$. It suggests that a significant deformation occurs in 4 h, if all the CMVs shown in the figure are admitted. This deformation violates the premise of the cross-correlation method in which template regions are simply slid with no deformation (note that the first and the last images used are separated by 4 h). The proposed RM eliminates the four CMVs shown in the red color in figure 5(a). It is obvious that if they are removed from the grid, it does not indicate a significant (inconsistent) deformation. The situation is similar for the other CMVs in figure 5(b) that do not exist in figure 5(c).

Zoom In Zoom Out Reset image size

The above results show that the IH16 method (CCS superposition) provides reasonable results without post processes. For the sake of demonstration, results obtained with the conventional cross-correlation method in which each CMV is derived from a single CCS are shown in figure 6.

Zoom In Zoom Out Reset image size

Consideration of the no-match possibility eliminates the four erroneous vectors at around $\left({{115}^{\circ}},-{{20}^{\circ}}\right)$, as seen in the difference between figures 6(a) and (b). However, the three erroneous vectors at around $\left({{97}^{\circ}},{{30}^{\circ}}\right)$, which have significant southward velocity, still remain. This is because the three support with each other and not many vectors are obtained nearby. It leaves the deformation inconsistency between the CMV at $\left({{100}^{\circ}},{{28}^{\circ}}\right)$ (southwestward) and that at $\left({{100}^{\circ}},{{25}^{\circ}}\right)$ (westward). As seen in figure 6(c), the three CMVs are separated by the post-process proposed in section 5.1.

Akatsuki UVI has two filters: one is the 365 nm filter, with which the images used above are obtained, and the other is the 283 nm filter. Fukuhara et al (2017) found that stationary gravity waves whose phase velocity is close to zero are observed at this wavelength as well as in the images obtained by the long infrared camera (LIR), another imager on Akatsuki. Figure 7 shows the CMVs obtained with the proposed RM and the two post processes in section 5. The black arrows show the largest group, representing atmospheric flow, which is consistent with the 365 nm results. The red arrows show the CMVs of the second group. They are small, indicating nearly stationary movement. Actually, the red arrows are obtained at where a stationary gravity wave packet is observed. This result demonstrates the ability of the post-process in section 5.2 to detect waves. Note that both of the two kinds of velocities are obtained at one of the grid points.

Zoom In Zoom Out Reset image size

Cloud tracking has been widely used to observe atmospheric winds on the Earth and other planets having clouds. PIV techniques have been adapted, but cloud tracking is still demanding and often difficult especially for clouds on Venus.

We have been developing a cloud tracking system for the Venus orbiter Akatsuki by using the cross-correlation method by IH16, in which multiple CCSs are overlaid to reduce erroneous template matching and to increase the S/N ratio. IH16 also provides methods to estimate uncertainty. We also use the relaxation labeling technique, also called the RM, to select candidate displacement vectors from the peaks of CCSs.

In this study, we developed methods to improve and enhance the correlation-based cloud tracking. It was made by exploiting consistency inherent in the type of cross-correlation method where template sub-images are slid without deformation, which we call the deformation consistency. Namely, if the resultant displacement vectors indicate a too-large deformation, it is contradictory to the assumption of the method.

The compatibility parameter in the RM was formulated in terms of the deformation consistency. Its value is decreased if the fractional change of grid-point distances is increased, and a neutral value is given to the compatibility between match and no-match cases. With this formulation, the RM tends to select no-match (no data) rather than to select displacement vectors where candidates over neighboring grid points do not have a deformation consistent solution.

We also developed two new post processes based on the deformation consistency; one clusters the displacement vectors into groups within each of which the consistency is perfect, and the other extends the groups by using the original candidate lists. The group having the largest number of vectors represents the flow, unless the tracking is severely unsuccessful or the cloud morphology is dominated by waves. Erroneous vectors, if any, are separated out. If clouds are sparse and the images used in tracking capture clouds at multiple altitudes with significantly different velocities, the clustering can separate them into different groups.

When images capture clouds modulated by a wave such as atmospheric gravity wave, its phase velocity may create peaks in CCSs, and the post-process may group them. The second post-process is useful to depict the extent of the area exhibiting wave signature, and both flow and phase velocities may be identified simultaneously.

It should be stressed that the clustering is based entirely on the deformation consistency. To physically interpret the results would require one to inspect the original images. Also, care is needed for possible degradation by interference.

In the RM and post processes, it is important to limit the candidates corresponding to the highest points of well-separated peaks in CCSs. It is also required in a uncertainty estimation by IH16. An algorithm to realize it was described in this paper.

All the methods described in this paper were implemented, and their effectiveness was demonstrated with images obtained on Dec 9, 2015 by the ultraviolet imager UVI onboard Akatsuki.

We thank Drs Masahiro Takagi, Hiroki Kashimura, Takeshi Imamura, Kensuke Nakajima, Naoki Sato, Takehiko Satoh, Javier Peralta, and Yoshi-Yuki Hayashi for helpful discussion. We also thank the two anonymous reviewers for their comments that helped us improve the paper. This study is partially supported by the JSPS Grant-in-Aid 16H02231 and 16H02225.