Riccardo Campisano

Heraldo Borges

Fábio Porto

Fabio Perosi

Esther Pacitti

Florent Masseglia

Eduardo Ogasawara

Fabio Porto

Discovering Tight Space-Time Sequences

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Space and time are pervasive in our day-to-day lives. As many datasets that include both time and space data are becoming available, new opportunities to discover interesting spatiotemporal patterns arise. An event may be classified as an occurrence of a phenomenon in a given space and time. A spatiotemporal sequential pattern is a sequence of events that are constrained in space and time [START_REF] Huang | A framework for mining sequential patterns from spatio-temporal event data sets[END_REF]. Due to that, spatiotemporal sequence mining is gaining attention [START_REF] Mooney | Sequential pattern mining -Approaches and algorithms[END_REF][START_REF] Sunitha | Mining frequent patterns from spatiotemporal data sets: A survey[END_REF].

In this work, we investigate a new problem related to spatiotemporal pattern identification. We are interested in finding tight space-time sequences, i.e., sequences that are constrained in space and time that may not be frequent in the entire dataset but are frequent inside a time interval and space range (spatiotemporal blocks). The primary challenge is to discover these blocks and the frequent sequences they contain. Solving this problem has a valuable impact on many applications.

Consider, for instance, the quantified-self movement, where people wear connected bracelets giving their position and inferring their activities. By analyzing the activities of their clients' bracelets, a brand might discover some habits regarding sports and food at coarse grain (say, "people who jog in the morning step by a vegan shop once a week"). However, there might be fine-grained behaviors that cannot be extracted, using existing methods, because they concern a niche (i.e., they have extremely low support, that hinders their discovery). Such a niche could be, for instance, that "people in Manhattan who jog at 7 am and have lunch near Time Square, spend 1 to 2 minutes in front of the buildings displays". This kind of pattern is of high interest for marketing and personalized recommendations. In the meantime, the accumulation of such niches, particularly in massive data, makes their discovery a challenging and valuable target for data analytics. The challenge is to extract both the pattern (e.g., jog, lunch, buildings displays), its occurrence time (e.g., from 7 am to lunchtime), and the location where it occurs (e.g., Time Square).

This paper introduces the problem of discovering tight space-time sequences and proposes an effective and efficient solution with the Spatio-Temporal Sequence Miner (ST SM ) algorithm to address it. ST SM applies a sequential pattern mining algorithm to detect spatial ranges where sequences are frequent. Next, it composes all detected sequences inside each spatial range to discover time intervals in which these sequences are frequent.

Our group first experienced the spatiotemporal sequence mine problem while analyzing seismic surveys [START_REF]Seismic Interpretation Software & Services. Tech. rep[END_REF]. Since then we have developed the ST SM algorithm and used the seismic motivation as our use case for validating the technique. By applying constrained space-time sequence mining techniques, we managed to identify 1,500 tight space-time sequences under high support threshold, which would not have been found using general spatial-time sequence mining algorithms (GST SM ) [START_REF] Sunitha | Mining frequent patterns from spatiotemporal data sets: A survey[END_REF] for the same frequency threshold. Moreover, the algorithm was able to detect candidate areas for seismic horizons (which stand for kind of layers) and bright spots (potential areas of hydrocarbon accumulation), which are known as important characteristics for the domain experts [START_REF] Zhou | Practical Seismic Data Analysis[END_REF].

In addition to this introduction, this paper is organized into five more sections. Section 2 discusses related works. Sections 3 and 4 formalize the problem and present the definitions needed to describe our algorithm. Section 5 gives the ST SM algorithm to identify sequential pattern constrained in space and time. Section 6 details the experimental evaluation. Finally, Section 7 concludes the paper.

Related works

Sequential pattern mining is a broad field of research embracing several approaches [START_REF] Mooney | Sequential pattern mining -Approaches and algorithms[END_REF][START_REF] Han | Frequent pattern mining: Current status and future directions[END_REF]. Over the last decades, some works have been developed for sequential mining of spatiotemporal datasets [START_REF] Sunitha | Mining frequent patterns from spatiotemporal data sets: A survey[END_REF]. We have conducted a systematic map study querying Scopus using the keywords ("sequence mining ∨ "sequential pattern ) ∧ ("space-time ∨ "spatiotemporal ). We identified 98 articles that were related to spatiotemporal sequential pattern mining, which are grouped according to different types of datasets: (i) trajectory datasets, (ii) event-based datasets, and (iii) emerging patterns datasets.

Trajectories datasets can be classified as a collection of events of the same moving object at different spatial locations and times. The goal is to retrieve similar trajectories that might reveal underlying traveling patterns of moving objects in the dataset [START_REF] Huang | A framework for mining sequential patterns from spatio-temporal event data sets[END_REF]. It includes works that search for routes frequently used [START_REF] Giannotti | Trajectory pattern mining[END_REF]. Since commonly the position of moving objects may be inaccurate, some approaches address the position uncertainty to better recognize common trajectories [START_REF] Li | Mining probabilistic frequent spatio-temporal sequential patterns with gap constraints from uncertain databases[END_REF]. Although trajectory problems are relevant, they differ from our problem. We are not interested in moving objects that have the same behavior. Instead, we are interested in regions and time intervals when some events are related (constrained) and relevant (with high support).

Event-based datasets correspond to the majority of related work. The goal is to find sequences constrained by space and time. For that, data is partitioned according to spatial or temporal dimensions, and events are related whenever they preserve certain proximity [START_REF] Tsoukatos | Efficient mining of spatiotemporal patterns[END_REF]. Huang et al. [START_REF] Huang | A framework for mining sequential patterns from spatio-temporal event data sets[END_REF] define sequence index as the significance measure for evaluating identified spatiotemporal sequential patterns. Such approach improves the mining algorithm to identify constrained sequences. Julea et al. [START_REF] Julea | Unsupervised spatiotemporal mining of satellite image time series using grouped frequent sequential patterns[END_REF] studied satellite images and target patterns with a high connectivity measure. Alatrista-Salas et al. [START_REF] Alatrista-Salas | Spatio-sequential patterns mining: Beyond the boundaries[END_REF] expand the search-space for identifying sequences by looking at neighbor partitions. Li et al. [START_REF] Li | Prediction of human activity by discovering temporal sequence patterns[END_REF] address the problem of event prediction from identified spatiotemporal sequence patterns. Finally, Batu et al. [START_REF] Batu | A Non-Parametric Algorithm for Discovering Triggering Patterns of Spatio-Temporal Event Types[END_REF] avoid partitioning space and time by computing the most relevant sequences in the dataset. Although approaches for event-based datasets are in the same category of our paper, they differ since all identified sequences are frequent in the entire dataset. In our work, sequences may only be frequent inside spatiotemporal blocks (i.e., a time interval and space range).

Finally, emerging pattern datasets correspond to solve the problem where data are continuously added to the database. Consider the patterns identified in the initial dataset and the patterns that are identified in an updated dataset. Previously identified patterns may become irrelevant, and new patterns may emerge. Some initiatives in emerging spatiotemporal datasets have been developed so far [START_REF] Chen | Constraint-based sequential pattern mining: The consideration of recency and compactness[END_REF][START_REF] Tsai | A change detection method for sequential patterns[END_REF][START_REF] Aydin | Spatiotemporal event sequence mining from evolving regions[END_REF]. These problems are complementary to ours since all identified patterns in both datasets (initial or updated) have high support. The problem tackled in this paper may find tight space-time sequence in both datasets, which may emerge or not.

Considering the aforementioned related works, to the best of our knowledge, we are not aware of similar works that studied tight space-time sequences from spatiotemporal datasets.

Formalization

Let t = <v 1 , v 2 , • • • , v n >
be a time-stamped sequence (TS), i.e., a sequence of items, where |t| = n is the number of items in t. A time index j is an integer value between 1 and n that is related to item v j . A sequence

s = <w 1 , w 2 , • • • , w k > is included in another sequence t = <v 1 , v 2 , • • • , v n > if there exist integers i 1 < i 2 < • • • < i k such that w 1 = v i1 , w 2 = v i2 , • • • , w k = v i k .
Let P = {p 1 , p 2 , ..., p m } be a set of positions, a spatial time-stamped sequence (STS) d is a couple (p, t) where p ∈ P is a position and t is the associated TS. A STS dataset D is a set of STS.

An STS d = (p, t) is said to support a sequence s if s is included in t. The support of a sequence s in D is the number of STS in D in which s is included. The frequency of a sequence s in D is the fraction of STS in D that supports s: f req(s, D) = sup(s,D)

|D|

. Given a user's minimum threshold γ ∈]0..1], a sequence s is said to be frequent if f req(s, D) ≥ γ. Let I be the set of all possible sequences s included in D, the set F of frequent sequences in D concerning γ is denoted by F (D, γ) = {s ∈ I : f req(s, D) ≥ γ}.

A spatial range (or simply range) r = (p s , p e ) is defined by a start position p s and an end position p e . We define the set of all potential ranges over D as P R. The set of STS that belong to range r is defined as T r(r) = {d : d ⊆ D, p s ≤ d.p ≤ p e }. The frequency of s over T r(r) in STS of a range r is denoted by f req(s, r) whenever it is clear from the context. The same approach is valid for support of s over T r(r), which is denoted by sup(s, r). A ranged sequence sr is a triple (s, r, f r) where s is a sequence, r is a range, and f r is the frequency of the sequence s over the range r: f r = f req(s, r). Let k be the size of s; then s is called a k-sequence.

A time interval (or simply interval) i = (i s , i e ) is defined by a start time i s and an end time i e . The length of an interval i is given by:

|i| = i e -i s + 1. Given a interval i, a sequence s = <w 1 , w 2 , • • • , w k > is a subsequence of another sequence t = <v 1 , v 2 , • • • , v n >: s = subseq(t, i) iff i s ≥ 1 ∧ i e ≤ n, |i| = k and ∀j ∈ [1..k], w j = v is+j-1
. By definition, s is also included in t. We define the set of all possible intervals over D as P I. A block b is a couple (r, i) where r is a range (r ∈ P R) and i is an interval (i ∈ P I). The size of a block b is the product of range size with interval length: |b| = |b.r|•|b.i|. We define the set of all possible blocks over a range r as P B(r). Moreover, a bounding box between two blocks b i , b j : bb(b i , b j ) is the minimum block b that contains them.

The occurrences of a sequence s in a block b refer to the positions in which all instances of s is present in b.

Thus, occur(s, b) is a set of blocks o, such that ∀b j ∈ o, b j ⊆ b | |T r(b j .r)| = 1 ∧ ∀d ∈ T r(b j .r), subseq(d, b j .i) = s.
: if req(s, b) = isup(s,b) |b| . Given a user's minimum threshold δ ∈]0..1], a sequence s is said to be item-frequent in a block b if if req(s, b) ≥ δ. A blocked sequence sb is a triple (s, b, if r) where s is a sequence, b is a block, and if r is the item-frequency of s over b: if r = if req(s, b).
Considering an STS dataset D, the problem we address is to find sequences in D that are frequent in constrained spatial range and time interval. The goal is to discover these ranges and intervals and the frequent sequences they contain.

In the following definitions, we introduce the notions of solid-ranged sequence and solid-blocked sequence that are fundamental for STSM algorithm. Their intuition is to respectively support the identification of spatial range and timespace blocks where a pattern is frequent.

STSM

This section is devoted to the presentation of ST SM (Spatio-Temporal Sequence Miner), our algorithm designed for the identification of solid-blocked sequences in the spatiotemporal dataset. The notion of kernels (range-kernels and blockkernels) introduced in this section allows for extracting solid-blocked sequences efficiently. First, we give an overview of our main algorithm in Section 5.2. Later, each relevant function is described in the following sections.

STS definitions

Given γ and δ, user's minimum thresholds, respectively, for frequency and itemfrequency, we introduce the characteristics of solid-ranged sequence in Definition 1 and solid-blocked sequence in Definition 2.

Definition 1 Let sr be a ranged sequence of range r, sequence s, and frequency f r. Then, sr is called a solid-ranged sequence iff the following conditions hold:

1) f r ≥ γ 2) ∀r 2 ∈ P R | r ⊆ r 2 , we have either a) or b) or both: a) f req(s, r 2 ) < γ b) sup(s, r 2 ) = sup(s, r) 3) ∀r 2 ∈ P R such that r 2 ⊂ r, sup(s, r 2 ) < sup(s, r)
The first condition in Definition 1 ensures that sr represents a sequence that is frequent over its associated range. The second condition ensures that the size of r is maximal. If a larger range exists, then, in this range, s is not frequent, or the support of s is the same (i.e., it is not worth extending the range from r to r 2 , since the extension is not going to contribute to the support of s). Finally, the third condition ensures that the size of r is minimal. In fact, s is supported by the first and last STS in r, so if a smaller range exists where s is frequent, the support is going to be lower anyway (i.e., relevant STS supporting s would have been dropped from the range and should be kept).

Let k be the size of s then sr is a k-solid-ranged sequence. SR k is the set of all k-solid-ranged sequences.

Definition 2 Let sb be a blocked sequence with a block b, sequence s, and itemfrequency if r. Then, sb is called a solid-blocked sequence iff the following conditions hold:

1) ∃ sr ∈ SR |s| | b.r ⊆ sr.r and s = sr.s 2) if r ≥ δ 3) ∀b 2 ∈ P B(r) | b ⊆ b 2 , we have either a) or b) or both: a) if req(s, b 2 ) < δ b) isup(s, b 2 ) = isup(s, b) 4) ∀b 2 ∈ P B(r) | b 2 ⊂ b, isup(s, b 2 ) < isup(s, b) 5) sup(s, b) > 1
The first condition ensures that the range of sb is within the range of a solidranged sequence sr. In other words, candidate blocks for solid-blocked sequences are identified from computed solid-ranged sequences. The second condition ensures that s corresponds to a sequence that is item-frequent in b. The third condition ensures that the size of b is maximal. If a larger block exists that contains b, s is not item-frequent, or the item-support of s is the same (i.e., it is not worth extending the block from b to b 2 , since the extension is not going to contribute to the item-support of s). The fourth condition ensures that the size of b is minimal. In fact, an item of s is present in the first and last range in b.r, and in the first and last interval in b.i, so if a smaller range or interval exists where s is solid-item-frequent, the support is going to be lower anyway. Finally, the fifth condition avoids trivial solid-blocked sequences that contain only a single occurrence of s in b.

Let k be the size of s; then sb is a k-solid-blocked sequence. SB k is the set of all k-solid-blocked sequences.

General principle

Algorithm 1 aims to generate solid-blocked sequences (SB) from size 1 to k. It receives as input an STS dataset D, a solid range threshold γ and a solid block threshold δ. It has three main functions: (i) candidate generation, (ii) identification of solid-ranged sequences, (iii) identification of solid-blocked sequences. The algorithm starts from candidates ranged sequences of size 1 (detailed in Section 5.3). They are built from all distinct items presented in D considering its entire range.

It starts a repeat-until loop that computes all solid-ranged sequences SR k with a frequency greater than or equal to γ (detailed in Section 5.4). Then, candidate ranged sequences of size k + 1 are computed from solid ranges SR k . Once we get empty candidates ranged sequences of size k + 1, the loop stops.

Finally, all solid-blocked sequences SB k with item-frequency greater than or equal to δ are computed (detailed in Section 5.5) from identified solid-ranged sequences.

Algorithm 1 Spatio-Temporal Sequence Miner

1: function ST SM (D, γ, δ) 2: C1 ← generateCandidates(D, nil) 3: k ← 0 4: repeat 5: k ← k + 1 6: SR k ← solidRangedSequences(D, C k , γ) 7: C k+1 ← generateCandidates(D, SR k ) 8:
until

C k+1 = ∅ 9:
for

(i ∈ {1 • • • k}) do 10: SBi ← solidBlockedSequences(D, SRi, δ) 11: return {SB1, • • • , SB k }

Generation of candidates

Frequent pattern mining algorithms aim at providing efficient algorithms on larger datasets. The generation of candidates should accomplish the mission: start with solid-ranged sequences of size 1 (SR 1 ) and explore the support of larger solid-ranged sequences with a limited number of scans over the dataset.

To this end, we compute range frequencies for candidate solid-ranged sequences of size k from computed solid-ranged sequences of size k -1 in only one scan. At the end of each scan, solid-ranged sequences of size k are identified.

ST SM introduces a spatial aware counting step for generated candidates. Let us consider xr(s x , r x , f r x ) to be a ranged sequence that is not a solid-ranged sequence. Any ranged sequence superset yr = (s y , r y , f r y ), such that s x ⊆ s y ∧ r x ⊆ r y cannot be a solid-ranged sequence (i.e., since f r x < γ then f r y < γ). Thus, the generate candidates follows an apriori-like principle with an additional filter on the possible intersection of the candidates (i.e., if two solid-ranged sequences of size k have a common subsequence, but their ranges do not intersect, they are not considered for generating a new candidate).

Selection of solid-ranged sequences

In this section, we define a range kernel, which is the basis of solid-ranged sequences described in Algorithm 2. A range kernel for a sequence s is a range where its frequency is greater or equal to γ. Let K(s, r, γ) be the set of range kernels for the sequence s over a range r concerning the minimum threshold γ. K(s, r, γ) is defined as follows: Let k ⊆ r be a subrange such that s is included in T r(k), i.e., sup(s, k) > 0. Let k s be the first position in which s is present in r. If k s does not exist, then K = ∅. If k s exists, then let P be the set of positions such that ∀p ∈ P, p ∈ r ∧ p > k s ∧ f requency(s, [k s ..p]) < γ (in other words, P is the set of positions in r such that extending the range k up to any of those positions leads to a frequency less than γ for s). If P is empty, then k e is defined as the last position in which s is present in r, and K(x, r, γ) = {k}. Otherwise (i.e. P = ∅), let q ∈ P such that ∀p ∈ P, p > q (q is the first position such that frequency of s is lost on [k s ..q]). Then, k e is defined as the last occurrence of s in [k s ..q] and K(s, r, γ) = {k} ∪ K(s, r -[k s ..k e ], γ).

The mechanics of computing range kernels is encapsulated by function rangeK of Algorithm 2. Intuitively, range kernels of a sequence are the longest ranges such that: (i) The first and last records support the sequence; (ii) The frequency of the sequence is always greater than or equal to γ when it is counted from the beginning until the end of the range.

Algorithm 2 solid-ranged Sequences After extracting range kernels for all candidate sequences C k , whenever possible, Algorithm 2 merges them. Algorithm mergeK finds solid-ranged sequences of s over r concerning threshold γ. Given two range kernels t and u, Algorithm mergeK tries to merge them into v using function mRS(t, u) as long as f req(s, v) ≥ γ is satisfied, which simply combines their ranges: t, u. In this case, v is included in the sker list, and both t and u are marked for removal. Such approach is repeated until no more pairs of kernels can be merged. Finally, in lines 6-8 of Algorithm 2, all merged kernels are used to produce solid-ranged sequences SR k .

Selection of solid-blocked sequences

In this section, we define a block kernel, which is the basis of solid-blocked sequences described in Algorithm 3. A block kernel is a block b in a solid range sr that contains a set of occurrences of s in b such that the item-frequency for s in b is greater or equal to δ. Let B(sr, δ) be the set of all possible block kernels in sr that has item-frequency not less than δ: B(sr, b)

: b k ∈ P B(sr)∧if req(s, b k ) ≥ δ.
Given a block kernel b ∈ B(sr, δ), let o be the set of occurrences of sr in b: o = occur(s, b). If |o| = 1, then according to this definition, if req(s, o) ≥ δ. We name this case as trivial block kernel. However, if |o| > 1, then there exists two block kernels b i and b

j (b i ⊂ b, b j ⊂ b), such that if req(s, b i ) ≥ δ ∧ if req(s, b j ) ≥ δ ∧ o = occur(s, b i )∪occur(s, b j ) ∧ occur(s, b i )∩occur(s, b j ) = ∅. Let o = occur(s, b) of s in a block kernel b. Whenever |o| > 1, (s, b, if req(s, b)) is a solid-blocked sequence.
From the previous definition, Algorithm 3 computes all solid-blocked sequences. It starts computing all trivial block kernels. Then, Algorithm mergeK finds solid-blocked sequences of s over sr concerning threshold δ. Given two block kernels t and u, Algorithm mergeK tries to merge them into v using function bb(t, u), which produces a bounding box for t and u. If v has an if req for s greater than or equal to δ, v is included in the stkernel list, and both t and u are marked for removal. Such approach is repeated until no more pairs of kernels can be merged. Finally, in lines 6-8 of Algorithm 3, all merged block kernels are used to produce solid-blocked sequences SB k .

Algorithm 3 solid-blocked Sequences

Merging kernels

Let us consider that we are provided with a sequence s and L, the set of kernels of s concerning a threshold σ. Algorithm mergeK can be used to compute both solid-ranged sequences of s over r concerning γ and solid-blocked sequences of s over sr concerning δ. To do so, Algorithm mergeK receives as input mf unc, which is a function that merges kernels, and thresf unc, which is a function that computes the frequency used to evaluate if a merged kernel respects threshold σ.

Proof for merge kernels Lemma: Let K be the set of range kernels of s on r concerning γ. Algorithm mergeK makes it possible to find all the solid-ranged sequences sr = (s, r, f r) on r.

Proof Let k ∈ K, be a range kernel of s after Algorithm mergeK (i.e. k cannot be merged with any other range kernel in K), then:

1. f req(s, k) ≥ γ. According to solid range definition, s is frequent in each kernel. Furthermore, if k is the result of a merging process, then Algorithm mergeK checks the frequency of s on the resulting range.

2. ∀q such that k ⊆ q, we have one of the following cases:

Algorithm 4 Algorithm Merge Kernels if thresf unc(s, v) ≥ σ then 8:

L ← L + v 9: K ← K + t + u 10:
mergeable ← true 11:

for k ∈ K do 12:

L ← L -k 13:

K ← ∅ sup(s, q -k) > 0, then s is not frequent in q (otherwise, let us consider k the kernel to which s is present in q, then k and k would have been merged). sup(s, q -k) = 0, then sup(x, q)=sup(x, k) (in this case, x may remain frequent in q or not, depending on the size of q). 3. According to Definition 2.1, s is supported by the first and the last STS in k. Then, s has lower support on any sub-range of k.

Based on the three observations above, let SR x = {(s, r, f r), ∀r ∈ K)} be the set of ranged sequences corresponding to all the merged kernels of s on P R concerning γ, then SR x is the set of all solid-ranged sequences sr = (s, r, f r) on P R concerning γ.

A similar proof is analogous to mergeK applied in solid-blocked sequences.

Toy Example

Figure 1 shows an example of spatiotemporal dataset D to better illustrate the definitions and mechanics of ST SM algorithm. Each position p i is associated with an STS d i in dataset D. Let us consider a frequency threshold γ = 1 2 given by the user. In this context, the only frequent sequence for D is <a> with a frequency of 1 2 using GST SM . This dataset presents other sequences that can be identified by ST SM according to our definitions. Consider the ranges r

1 ([d 1 ..d 2 ]), r 2 ([d 8 ..d 10 ]), r 3 ([d 2 ..d 5 ]
), and r 4 ([d 7 ..d 8 ]). Both r 3 and r 4 are, respectively part of solid-ranged sequences sr 3 (<e, o>, r 3 , 3 4 ) and sr 4 (<i, u>, r 4 , 1). In fact, considering our generate candidate algorithm, solid-ranged sequence sr 3 is obtained from (<e>, r 3 , 3 4 ) and (<o>, r 3 , 3 4 ). Additional, r 1 and r 2 are not solid-ranged sequences. They correspond to range kernels for <a> identified by algorithm solidRangedSequences. Later, they are merged to build a solid-ranged sequence (<a>, [d 1 ..d 10 ], 1 2 ). As it is, the sequence <a> may not appear to be interesting, but if frequency threshold is increased to γ = 2 3 , both r 1 and r 4 become part of a solid-ranged sequence for <a>, whereas no sequence is identified by general spatial-time sequence mining algorithms (GST SM ).

(
Taking into account the relationship between time and space, intuitively, sequence <a> still does not appear to be interesting. In fact, the goal of ST SM is to identify tight space-time sequences. Consider blocks b

1 = ([d 2 ..d 5 ], [v 3 ..v 5 ]) and b 2 = ([d 7 ..d 8 ], [v 1 ..v 5 ]). The occurrences of <i, u> in b 2 are {([d 7 ], [v 1 ..v 2 ]), ([d 8 ], [v 4 ..v 5 ])}.
Adopting an item-frequency threshold δ = 2 5 , both blocks are part of solid-blocked sequences sb 1 (<e, o>, b 1 , 1

2 ) and sb 2 (<i, u>, b 2 , 2 5 ). Additionally, sequence a is not produced by ST SM algorithm as a solid-blocked sequence for such threshold.

We can observe that using solid-ranged sequences and solid-blocked sequences definitions together with ST SM algorithm, it is now possible to find sequences that would not be found using higher frequency thresholds. According to the example, the sequence <e, o> is supported by the entire dataset with a frequency of 3 10 occurs on a particular range (i.e. [d 2 ..d 5 ]) and interval (i.e. [v 3 ..v 5 ]) with a frequency of 3 4 and item-frequency of 1 2 . Also, it filters unwanted sequences that are not tight in space-time, such as <i, u> in block ([d 1 ], [v 5 ..v 6 ]). Sequence <i, u> i is only solid-blocked sequence in b 2 .

Experimental evaluation

This section presents the experimental evaluation. Section 6.1 starts presenting the dataset used as a proof of concept. Section 6.2 presents the exploratory analysis of ST SM . Finally, Section 6.3 discusses the behavior of the algorithm ST SM when compared to GST SM .

Dataset

The dataset chosen as input data is the inline 401 of the Netherlands Offshore F3 Block reflection seismic survey [START_REF]Seismic Interpretation Software & Services. Tech. rep[END_REF]. It is composed of 651 inlines. The inline contains 951 spatiotemporal series of timely-ordered observations collected in different positions on the surface. Each of spatiotemporal series contains 462 observations, a sample of every four milliseconds, representing the different depths of the subsurface.

Seismic surveys are commonly collected by producing pertubations (shots) on the ground or in the ocean. Such shots generate elastic waves that are propagated in the Earth interior. They are reflected and refracted in each material layer and are returned to the surface in a given position and time. A set of receivers positioned on the surface (onshore or offshore) at fixed aligned intervals records the waveform and wave arrival times. The collected wave amplitudes and frequencies are related to the reflection location, and the wave traveled distance. Thus, the early received signals correspond to upper subsurface layers, and viceversa. As it can be depicted in Figure 2, for example, each spatiotemporal series can be interpreted as vertical columns of pixels, where negative amplitude values correspond to lighter pixels whereas positive amplitude values correspond to darker pixels. The inline was made available6 using SAX with an alphabet of size 7 [START_REF] Lin | A symbolic representation of time series, with implications for streaming algorithms[END_REF].

Exploratory Analysis

The algorithm proposed in this work allows for users to set solid range threshold γ and solid block threshold δ constraints. Finding adequate values for these thresholds depends on the characteristics of input dataset/application. Lower values for such constraints can lead to the identification of a large number of non-useful frequent sequences. Conversely, higher values for these thresholds can result in the detection of a small number of frequent sequences that may become too small to be interesting. All code, experimental data, and results are made available 6 .

We explored the combination of solid range threshold γ (60%, 70%, 80%, 90%, and 100%) with solid block threshold δ (10%, 20%, 30%, 40%, and 50%). For sake of comparison, consider the identified sequence <a, a, g, g>. In this scenario, when γ is under 80%, a single solid range for this sequence is detected covering the entire dataset. Values between 80% and 100% divide the dataset into more fine-grained solid ranges. A consequence of having better solid range division is that the computation required for identifying solid-blocked sequences is reduced. Relaxing γ from 80% to 60% led to an increase 18% in computing time without improving the detection of a solid-blocked sequence. Adopting the extreme value of 100% may, however, limit too much the detection of solid block sequences. Figure 2.a depicts sequence occurrences for γ = 80%.

Similarly, examining the solid block threshold δ, when relaxing it to 10%, many larger solid-blocked sequences are identified with many overlapping among them. Conversely, when increasing the δ threshold to higher values (30%), it is possible to observe a significant decrease in the number of identified solid-blocked sequences, many of them lose their overlapping. In Figure 2.a, δ is set to 20%.

Discussions

The goal of this paper is to provide a way to discover sequences that are frequent in tight time and space. By applying ST SM , we managed to identify sequences that would not have been found using traditional spatiotemporal techniques using the same support. We ranked the identified sequences by density. The density of a sequence s was computed by the mean block size of all solid-blocked sequences for s. A representative high-ranked sequence (<e, e, f >) is depicted in Figure 2.b.

Despite the simplicity of used ranking and the absence of any significant seismic concept in its logic, it was good enough to prioritize interesting areas. In the seismic domain, it is of particular importance to identify horizons, i.e., zones of major nonconformities between geologic boundary materials, usually associated to different lithologies, which produces variations in the collected reflection values [START_REF] Zhou | Practical Seismic Data Analysis[END_REF]. ST SM algorithm was able to identify known horizons. The majority of high-ranked sequences made available in our site match known seismic horizons previously known for this dataset [START_REF]Seismic Interpretation Software & Services. Tech. rep[END_REF]. Additionally, ST SM was also able to find potential areas of accumulation of hydrocarbon, known as bright spots. Bright spots are rare patterns that occur when there is an inversion of the wave phase. Such pattern is important not only because gas can be interesting for exploration, but also since it can be near to rocks containing oil. For example, in Figure 2.a, the discovered patterns, evidenced in red, are plotted on the top of yellow marks obtained from the previously known bright spots for this dataset [START_REF]Seismic Interpretation Software & Services. Tech. rep[END_REF]. Figure 2.a shows that ST SM algorithm has discovered patterns that cover large parts of known bright spots.

Figure 2.b shows a comparison between GST SM and ST SM . The sequence <e, e, f > is identified by both algorithms. In GST SM it has high support, whereas in ST SM it is mid-term ranked results obtained the chosen configuration. It can be observed that although many solid-blocked sequences (depicted in red) are covering the majority of known horizons, many isolated occurrences identified by GST SM (marked as black) correspond to noise. Considering Figure 2, it is possible to conclude two important differences between GST SM and ST SM . Not only ST SM can identify sequences that GST SM is not capable of identifying without lowering too much the frequency threshold, but it also can filter-up sequence occurrences regarding space and time. Such approach allows for researchers to focus on tight sequences.

Conclusion

This paper presented a challenging problem with high potential impact: spatiotemporal sequential mining that focuses on discovering frequent tight patterns, along with their constraints of frequency in both space and time. Our formal background provides us with the necessary conditions of such discovery. We proposed ST SM , a novel algorithm that performs an efficient extraction of these patterns and their spatiotemporal constraints. To the best of our knowledge, this is the first study on this topic. Our experiments, using a real-world seismic dataset, highlighted significant insights according to the feedback of specialists in the domain. ST SM allows for extracting patterns that follow areas of major nonconformities between geologic boundary materials of the subsurface, such as horizons and bright spots. The existence of these discovered patterns was confirmed comparing them with a previously known survey on this dataset, while their meaning was assessed by domain experts. Owing to its novelty and the characteristics of the extracted patterns, ST SM is opening new research for tight spatial-time sequences discovery.

Although we have evaluated ST SM algorithm using seismic data, it has been conceived to be sufficiently generic for different spatiotemporal datasets (such as IoT domain). The plot of detected solid-blocked sequences according to space and time, as shown in this paper, can aid specialists to analyze identified sequences. Additionally, in scenarios in which a large number of frequent constrained sequences are produced, we can explore specialized ranking functions to prioritize the most interesting ones.

  The support of a sequence s in a block b, sup(s, b) is the number of occurrences of s in b: |occur(s, b)|. The item-support of a sequence s in a block b, is the product of support of a sequence s in a block b with the length of the sequence s: isup(s, b) = sup(s, b) • |s|. The item-frequency of a sequence s in a block b is the fraction of item-support over the size of b

Fig. 2 .

 2 Fig. 2. Identified solid-blocked sequence <a, a, g, g> for inline 401, alphabet size 7, solid range threshold γ = 80% and solid block threshold δ = 20%. Its density was 206. Solid-blocked sequences are marked in red. The results follow the yellow pattern produced using the previously known bright spots for this dataset [5]. (a) Comparison of quality between GST SM and ST SM for sequence <e, e, d, d> in inline 401 using alphabet size 7, with support of 80% for GST SM and with solid range threshold (γ) of 80% and solid block threshold (δ) of 20% for ST SM . Identified occurrences are marked as red when identified by ST SM and as black in GST SM . Although occurrences from ST SM correspond to seismic horizons, many occurrences from GST SM correspond to noise. (b)

  Fig. 1. Dataset D for ten spatial-time series with six observations in each.

	t	D	d 1	d 2	d 3	d 4	d 5	d 6	d 7	d 8	d 9 d 10
	v 1	a	b	c	d	t	q	i	g	a	h
	v 2	k	l	m	n	p	q	u	s	t	v
	v 3	w	e	e	x	y	m	a	r	d	a
	v 4	h	o	o	g	e	i	e	i	c	b
	v 5	i	j	k	l	o	z	n	u	z	p
	v 6	u	a	r	s	t	µ	c	d	f	a
			(sr 1 )								(sr 2 )
						sr 3 )			(sr 4 )

https://github.com/eogasawara/stsm

Acknowledges

The authors would like to thank CAPES, CNPq, and FAPERJ for partially funding this paper.


Zhiqi Tang
email: zhiqitang@tecnico.ulisboa.pt

Rita Cunha

Tarek Hamel

Carlos Silvestre

Going through a window and landing a quadrotor using optical flow

This paper considers the problem of controlling a quadrotor to go through a window and land on a planar target, using an image-based controller, with only a camera and an Inertial Measurement Unit (IMU) as sensors. The maneuver is divided into two stages: crossing the window and landing on the target plane. For the first stage, a control law is proposed that guarantees that the vehicle will not collide with the wall containing the window and will go through the window with non-zero velocity along the direction orthogonal to the window, keeping at all times a safety distance with respect to the window edges. For the landing stage, the proposed control law ensures that the vehicle achieves a smooth touchdown, keeping at all time a positive height above the target plane. For control purposes, the centroids of the images of a collection of landmarks (corners) for both the window and the target are used as position measurement. The translational optical flow relative to the wall, window edges, and target plane is used as velocity cue. To achieve the proposed objective, no direct measurements of position or velocity are used and no explicit estimate of the height above the target plane or of the distance to the wall is required. Simulation results are provided to illustrate the performance of the proposed controller.

I. INTRODUCTION

.

The use of vision for obstacle avoidance has been a active topic of research. The work in [3] uses methods with pushbroom stereo to automatically avoid obstacles. In [1], the authors address the problem of state estimation, control, and planning for aggressive flight through a narrow window using only a single camera and an IMU.

 [START_REF] Ho | Optical-flow based self-supervised learning of obstacle appearance applied to mav landing[END_REF]introduce a novel setup of self-supervised learning, in which optical flow cues serve as a scaffold to learn the visual appearance of obstacles in the environment. In [START_REF] Mcguire | Efficient optical flow and stereo vision for velocity estimation and obstacle avoidance on an autonomous pocket drone[END_REF], stereo vision and optical flow are used for obstacle avoidance during the landing process, requiring estimation of velocity and depth. A method for ego-motion estimation using optical flow is proposed for UAV flying in cluttered environments [START_REF] Briod | A method for ego-motion estimation in micro-hovering platforms flying in very cluttered environments[END_REF].

In this paper, we consider the problem of landing a quadrotor on a planar target that is placed inside a room and requires the vehicle to first go through a window and only then land on the target. This scenario occurs in many practical applications, such as search and rescue in an earthquakedamaged building [START_REF] Michael | Collaborative mapping of an earthquake-damaged building via ground and aerial robots[END_REF], recharging or battery change during long missions, package delivery using UAVs, etc. This work is builds on previous work reported in [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF], which used optical flow to land a UAV in cluttered environments. Compared with [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF], where obstacles were represented by spherical balls, the present work considers the more realistic scenario of crossing a window, while avoiding the window edges and the wall that contains the window. Compared with other vision-based methods, the proposed controller requires no estimation of 3D quantities, using directly image centroids as position measurements and optical flow as a velocity cue.

The proposed control law comprises two stages. The first accomplishes the goal of going through the window and the second the goal of landing on the target. For control purposes, collections of landmarks are used as features, which are given by the edges and corners of the window for the first stage and the corners of the landing area for the second stage. The corresponding spherical image measurements are combined into centroid vectors, which provide information about position. For velocity, it is assumed that both the target plane and the window plane have enough texture to provide reliable optical flow measurements, which are then used to compute the translational optical flow relative to the window plane, window edges, and target plane. With no direct measurement of position or velocity, the control law proposed for going through the window ensures that no collision with the wall or window edges will occur and the vehicle will align with the center line orthogonal to the window, crossing it with non-zero velocity. The control law proposed for landing on the target is an improvement with respect to the one used in [START_REF] Serra | Landing of a quadrotor on a moving target using dynamic imagebased visual servo control[END_REF], with the centroid vector now directly given by the image centroid. The proof of convergence is also simplified by relying on an alternative storage function that considers all directions of motion in conjunction, requiring no ultimate bound on the height above the target and keeping the same property that the height will remain positive for all time.

1 The body of the paper consists of six parts. Section II presents the kinematics and dynamics of a quadrotor vehicle. Section III introduces the environment and presents the image features that are used in the control law. Section IV proposes two control laws, one for landing in an obstacle-free environment and the other for flying through the window. The combination of these two control laws in a practical case is also presented in this section. Section V shows simulation results obtained with the proposed controller. Conclusions are presented in Section VI.

= v v = F (1)
Ṙ = RS(Ω)

I Ω = -S(Ω)IΩ + Γ (2)
where F ∈ IR 3 expressed in {I} denotes the force applied to the vehicle, Γ ∈ IR 3 expressed in {B} is the torque and I is the tensor of inertia. The matrix S(Ω) denotes a skewsymmetric matrix. The force F combines thrust and gravity components and can be described as

F = ge 3 - 1 m T Re 3 ( 3
)
where g is the gravitational acceleration, m is the mass of the vehicle, the variable T represents the thrust magnitude, and e 3 = [1 0 0] .

III. ENVIRONMENT AND IMAGE FEATURES

Using an image-based controller to perform a landing maneuver and to avoid the wall and the window edges requires the definition of adequate images features. Some assumptions regarding the environment and the setup are established.

Assumption 1: The camera is attached to the center of mass of the vehicle so that the camera reference frame coincides with the body-fixed frame.

Assumption 2: The landing target lies on a textured plane called the target plane. The normal to the target plane, with coordinates in the inertial frame denoted by η t is aligned with the gravitational direction (η t ≡ e 3 ).

Assumption 3: The target window has the form of a square and lies on a textured wall called the window plane, with normal in the inertial frame denoted by η w . Assumption 4: The orientation matrix R is known in order to represent all image information as well as the dynamics of the system in the inertial frame.

Assumption 5: Assume that the control input is F or more specifically the thrust vector (T Re 3 ). This implicitly means that any desired orientation in the direction e 3 , R d e 3 = mge 3 -F d T

(with T = mge 3 -F d ), can be obtained after a short transient. In the language of automatic control there is a high gain inner-loop controller that ensures that Re 3 ≈ R d e 3 [START_REF] Hérissé | Landing a vtol unmanned aerial vehicle on a moving platform using optical flow[END_REF]. Consequently, the control input of the translational dynamics described by ( 1) is directly given by F = F d .

Both target plane and window plane are placed in the environment, as shown in Figure 1. It is assumed that the vehicle is able to recognize the target and the window by landmarks on the target and the edges of the window respectively. The background texture on both target plane and window plane can be used to obtain information about velocity relative to the background and also to avoid collision with both planes. With the initial position outside the room containing the target, the main objective of the controller is to drive the vehicle through the window and land it safely on the target, avoiding collisions along the way.

A. Image features of the target plane

The target on the target plane is shown in Figure 1. The axes of {I} are given by (u t , ρ t , η t ), where ρ t = η t × u t , and the origin of {I} is placed at the center of the target. As shown in Figure 1,s t i denotes the position of ith corner of the target relative to the inertial frame. Note that η t s t i = 0 and ∑ i s t i = 0. Define the position vector of ith corner of the target relative to {B} as

P t i = s t i -ξ . (4)
Using the spherical projection model for a calibrated camera, the spherical image points of target's corners can be expressed as

p t i = P t i P t i .
A visual feature that encodes information about the position of the vehicle relative to the target plane is the centroid vector

q t := - 1 n t n t ∑ i=1 p t i ( 5
)
where n t is the number of the corners.

B. Image features of the window plane

A rectangular window with r of width is placed on a textured wall as shown in the Figure 1. The axes representing the window are given by (η w , ρ w , u w ), where ρ w = u w × η w . As shown in Figure 1,s w i denotes the position of ith corner of the target relative to the inertial frame. Similar to Section A, the spherical image points of the corners of the window can be expressed as

p w i = P w i P w i ( 6
)
where P w i = s w iξ . Let l i denotes the ith edge of window as shown in Figure 2. Using the image of l i and by exploiting the fact that the window has a rectangular shape, it is straightforward to get the directions u w and ρ w and consequently η w . As described in [START_REF] Mahony | Image-based visual servo control of aerial robotic systems using linear image features[END_REF], in binormalized Euclidean Plucker coordinates, l i can be represented by its unit direction u w (resp. ρ w ) and the unit direction h i , normal to the plane defined by the origin of the camera/body-fixed frame and the line l i . The unit vector h i can be obtained directly from the images [START_REF] Mahony | Image-based visual servo control of aerial robotic systems using linear image features[END_REF]. Using the fact that l 1 and l 3 (resp. l 2 and l 4 ) are parallel in the inertial frame, we get the measure of the direction u w (resp. ρ w ) from the following relationships:

ρ w = ± h 1 × h 3 h 1 × h 3 u w = ± h 2 × h 4 h 2 × h 4 . (7)
Define a visual centroid vector q := 1 n w ∑ n w i=1 p w i , then the normal vector to the window plane can be obtained by

η w = ± u w × ρ w u w × ρ w . ( 8
)
The direction of η w is decided such that η w q(0) > 0. Let

d o := η w P i w (9)
denote the distance along the direction η w between the vehicle and the window plane.

We can also extract the direction to the nearest point to each edge

p i e = ±(h i × ρ w ), i = {1, 3}, p i e = ±(h i × u w ), i = {2,
W := {ξ : q(ξ ) ≤ ε}, (12)
where ε > 0. The bound ε must be such that for ξ ∈ W , the condition 1 n w ∑ n w i=1 P i w < r 2ε also holds, meaning that the region defined by W does not contain the window edges.

The visual feature that encodes information about the position of the vehicle relative to the window and is used for feedback is given by

q w =: -1 n w ∑ n w i=1 p w i (1 -α w (t)) 1 η w p w i + α w (t) η w p e η w p w i , (13)
where n w (n w = 4) is the number of corners of the window and the weight α w (t) is given by

α w (t) =      0 , if q > ε 1 δ (ε -q ), , if ε -δ ≤ q ≤ ε 1 , if q < ε -δ , (14)
which ensures that there is continuity in the visual feature q w (t). Noting that η w p e = d o d e , and substituting (6) into (13), q w can be rewritten as

q w (t) = (1 -α w (t)) χ(t) d o (t) + α w (t) χ(t) d e (t) (15)
where χ denotes the position of the vehicle relative to the window's center, such that

χ = - 1 n w n w ∑ i=1 P w i = ξ - n w ∑ i=1 s w i , (16)
and ∑ n w i=1 s w i is constant.

C. Image Kinematics and Translational Optical Flow

The kinematics of any observed points on the target plane can be written in the inertial frame as:

Ṗt = -v (17)
where P t denotes a point on the textured ground of the target plane, not to be confused with P i t , which refers to the ith corner of the target. So the kinematics of the corresponding image point p t = P t P t can be written as

ṗt = -π p t v P t ( 18
)
For the point on the textured ground of target plane, define the height of the vehicle above the target plane d t := η t P t and, equation (18) can be rewritten as

ṗt = -cos θ t π p t v d t (19)
where cos θ t = η t p t . It can be shown that the translational optical flow used for landing is

φ t (t) = v(t) d t (t) , (20)
which is obtained from the integral of ṗt in the normal direction to the observed surface over a solid angle [START_REF] Serra | Landing of a quadrotor on a moving target using dynamic imagebased visual servo control[END_REF]. Similarly, we can exploit the kinematics of observed image points from the textured surface of the window plane, to obtain the following optical flow measurement

φ w = (1 -α w (t)) v(t) d o (t) + α w (t) v(t) d e (t) (21)
which will be used as velocity cue for going through the window.

IV. CONTROLLER DESIGN

A. Landing in obstacle free environment Theorem 1: Consider the system (1) subject to the control input

F t = -k t 1 q t -k t 2 φ t . (22)
Then for any initial condition such that d t (0) > 0 and for some positive gains k t 1 and k t 2 the following assertions hold 1) the height d t (t) is positive and bounded for all the time, 2) d t (t) and ḋt (t) converge to zero asymptotically, 3) the states (ξ (t), v(t)) are bounded and converge to zero asymptotically.

Proof: See Appendix A.

B. Going through the center of the window

To accomplish the goal of going through the window, while avoiding the wall and window edges, we propose to use the following control law

F w = σ (t)(-k w 1 π η w q w -k w 2 π η w φ w -k w 3 η w (η w φ w + φ * )), (23)
with

σ (t) = 0 , if η w q ≤ 0 1 , if η w q > 0, (24)
yielding F w = 0 for d o ≤ 0, which indicates vehicle already crossed the window. Note that for η w q > 0, the resulting closed-loop system can be written as

   χ = v v = -k w 1 π ηw χ d w -k w 2 π ηw v d w -k w 3 η w (η w v d w + φ * ), (25)
where d w is such that

1 d w = ((1 -α w ) 1 d o + α w 1 d e ) ( 26
)
which according to the definition of α w given in (14) yields

d w =      d o , if q > ε d o d e (1-α w )d e +α w d o , if ε -δ ≤ q ≤ ε d e , if q < ε -δ (27)
Theorem 2: Consider the system (1) with the control input given by (23). Then for any initial condition such that d w (0) > 0, desired optical flow φ * > 0, and positive gains k w 1 , k w 2 and k w 3 such that

k w 2 2 k w 1 > r 2
, the following assertions hold: 1) the height d w (t) never crosses zero as long as q ≥ ε -δ ; 2) the states (π η w χ(t), π η w v(t)) are bounded and converge to zero exponentially as long as q ≥ εδ .

3) there is a finite time T > 0 when the system crosses the window d o (T ) = 0, with strictly negative velocity ḋo (T ) and inside the region where q(T ) < εδ .

Proof: See Appendix B.

C. Landing on the target plane after crossing the window

The double goal of crossing the window and landing on the target can be achieved by simply applying the control laws F w and F t in sequence, with an adequate trigger to switch from F w to F t . Under the assumption that the initial position of the vehicle is such that it must cross the window to reach the landing target, the trigger to start using F t is given by the condition η w q(t) ≤ 0, which is equivalent of having d o (t) ≤ 0.

By Theorem 2, we know that d o (t) = 0 can only be reached when the vehicle position is inside W and consequently inside the window area. Moreover, when the crossing occurs the velocity ḋo (t) is strictly negative, which implies that these are adequate initial conditions to start the landing maneuver using F t . To avoid any chattering behavior due to noise or external disturbances, the switching from F w to F t is triggered only once by the first occurrence of η w q(t) ≤ 0.

V. SIMULATION

In this section, simulation results are shown to illustrate the behavior of the closed-loop system using the designed controller. A high-gain inner-loop controller based on backstepping is used to control the attitude dynamics [START_REF] Tang | Homing on a moving dock for a quadrotor vehicle[END_REF]. It generates the torque inputs in order for the vehicle to follow a desired thrust direction R d e 3 , which is provided by the outerloop controller. In the simulation, the position of the center of the window is [-5 -

t 1 = 3, k t 2 = 3, k w 1 = 3, k w 2 = 3, k w 3 = 2
, and φ * = 0.1 s -1 . As shown in Fig. 3, with different initial positions and initial velocities, the vehicle avoids the wall successfully and aligns with the center line of the window for crossing. Finally the quadrotor lands on the center of the target. VI. CONCLUSION This paper considers the problem of controlling a quadrotor to go through a window and land on planar target, using an image-based controller that relies on optical flow measurements. For control purposes, the centroids of the images of a collection of landmarks (corners and edges) for the window and target are used as position measurement. The translational optical flow relative to the wall, window edges, and target plane is used as velocity measurement, meaning that no direct measurements or explicit estimates of 3-D position or velocity are used. With the initial position outside the room containing the target, the proposed control law guarantees that the quadrotor aligns with the center line orthogonal to the window, crosses it with non-zero velocity and finally lands on the planar target successfully avoiding collisions with the wall and edges of the window. A proof of convergence of the overall control scheme is provided and simulation results show the effectiveness of the proposed controller. Future work will focus on implementing and experimentally testing the proposed solution, taking into consideration the presence of disturbances and measurement noise.

APPENDIX

A. Proof of Theorem 1

Proof: The proof follows a reasoning very similar to that of Theorem 1 in [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF]. Recalling (1) and applying the control input (22), we can write the closed-loop system as

   ξ = v v = -k t 1 q t (ξ ) -k t 2 v d t . ( 28
)
Before proceeding to the proof of item 1), we define a positive definite storage function L 2 (ξ , v) and show that if d t (t) remains positive, L2 is negative semi-definite, which implies that the solutions remain bounded for all t ≥ 0. Let L 2 be given by

L 2 (ξ , v) = k t 1 L 1 (ξ ) + 1 2 v 2 (29)
where L 1 (ξ ) is the radially unbounded function

L 1 (ξ ) = 1 n t n t ∑ i=1 ( P t i (ξ ) -P t i (0) ). ( 30
)
To show that L 1 (ξ ) is positive definite, we note that

∂ L 1 ∂ ξ = q t ∂ 2 L 1 ∂ ξ 2 = Q (31)
where Q = 1 n t ∑ n t i=1 1 P t i π p t i is positive definite, as long as at least two of the vectors p t i are non-collinear. It follows that L 1 is a convex function of ξ , with a global minimum attained when ∂ L 1 ∂ ξ = q T t = 0, or equivalently, when ξ = 0. Since L 1 (0) = 0 is the global minimum of the function, we can conclude that L 1 (ξ ) is positive definite. Noting that, L1 = k t 1 q t v, it follows that

L2 = - k t 2 d t v 2 (32)
which is negative semi-definite as long as d t remains positive and implies that the states ξ (t) and v(t) remain bounded for all t ≥ 0.

Proof of Item 1: Using (28) and the fact that d t (t) = -η t ξ and ḋt = -η t v yields dt = -k t 2 ( ḋt

d t + α) (33)
with

α(t) = - k t 1 k t 2 η t q t = k t 1 k t 2 1 n t ∑ n t i=1 d t P i t , α(t) ∈ (0, 1), ∀t (34)
Taking the integral of (33) and applying a similar analysis as in [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF] (Th. 1, Item 1), we can conclude that d t remains positive for all time.

Proof of Item 2:

Since the time evolution of d t (t) is described by (33) and α(t) is bounded, we can conclude that d t (t) and ḋt (t) converge to zero asymptotically by direct application of the result in [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF] (Th. 1, Item 2).

Proof of Item 3:

To show that v(t) converges to zero, note that

v = -k t 2 v d t -k t 1 q t . (35)
Since d t (t) is converging to zero and q t (t) is bounded, we can immediately conclude that v(t) converges to zero by direct application of [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF] (Lemma 3).

To show that ξ (t) is also converging to zero, we can use an analysis similar to that in [START_REF] Rosa | Optical-flow based strategies for landing vtol uavs in cluttered environments[END_REF] (Th. 1, Item 4) applying it to q t (t) instead of ξ (t) .

B. Proof of Theorem 2

Proof: We consider only the case where σ = 1, that is, the situation in which the vehicle is going through the window without collision.

Proof of Item 1: We first recall that while q(t) > ε (and σ = 1), the dynamics of d o can be written as By the choice of gains k w 1 and k w 2 such that

do = -k w 3 ( ḋo d o + φ * ). (
(k w 2 ) 2 k w 2
> r 2 and recalling that χ(t) < r 2ε when q(t) ≤ ε, we have that the condition on d w holds when q(t) ≤ ε. From the proof of item 1), it is already proved that as long as q(t) > ε, d o (and therefore d w ) is decreasing and converging to zero which in turn implies that there is a finite time after which d w (t) = d o (t) < is converging to zero. This implies that there exists a finite time T after which d 0 (t) = 0 with ḋo (t) < 0.

  II. MODELING Consider a quadrotor vehicle equipped with an IMU and a camera. To describe the motion of the vehicle, two reference frames are introduced: the inertial reference frame {I} fixed to the earth surface and the body-fixed frame {B} attached to the quadrotor's centre of mass. Let R = I B R ∈ SO(3) denote the rotation matrix that transforms vectors from the frame {B} to {I} and ξ ∈ IR 3 the position vector of the vehicle expressed in {I}. Let v ∈ IR 3 denote the translational velocity expressed in {I} and Ω ∈ IR 3 the angular velocity expressed in {B}. The kinematics and dynamics of the quadrotor vehicle are then described as ξ

Fig. 1 .

 1 Fig. 1. Target plane and window plane.

  Fig. 2. Window.

  4}.Define the closest point from the vehicle to the edges of the window as P e . Then, its direction p e = P e P e is given byp e = arg max = P e = d 2 o + π η w P e 2(11)denote the closest distance from the vehicle to the edges of the window, where π x denotes the projection operator defined for any x ∈ S 2 := {x ∈ R 3 : x x = 1} as π x = I 3xx . Then define the region W such that

  4 20] m and the position of the center of the target is [0 0 0]. Different initial conditions are set as ξ (0) = [-5.5 -2.5 21] m with v(0) = [0.5 0 0] m/s, ξ (0) = [-4 -5.5 18] m with v(0) = [2 0 0] m/s and ξ (0) = [-4.5 -3.5 17] m with v(0) = [2 1 2] m/s, which means the quadrotor starts outside of the room containing the target. And gains are chosen following values k

Fig. 3 .

 3 Fig. 3. 3-D plot of the quadrotor trajectories.

F 1 :

 1 = π η w F w , χ 1 := π η w χ and v 1 := π η w v denote the components orthogonal to η w . From (1) and (23) L 4 is decreasing provided that 0 < d w <

w 3 d

 3 L 4 becomes negative definite. This implies that the state (χ 1 , v 1 ) is bounded and converges to zero exponentially.Proof of Item 3: Consider now the case for which q(t) < εδ while σ = 1. In this case the dynamics of d o takes the formdo = -k w 3 ( ḋo d e + φ * ),(41)Recalling that in this case the vehicle is inside the window region and hence one can ensure that d e is lower bounded, d e ≥ c > 0, implying that do ≤ -k e ( ḋo + cφ * ). Defining a new variable z := ḋo + cφ * , we obtain ż τ) dτ)z(0), and ḋo (t) ≤ -cφ * + o(t) where o(t) = exp(-k w dτ)( ḋo (0) + cφ * ) (43)

  36)Proceeding analogously to the Proof of Item 1 in Theorem 1, it is straightforward to verify that, while q(t) > ε, d o (t) remains positive and bounded and never crosses zero, and d o (t) and ḋo (t) are converging to zero.When εδ ≤ q(t) ≤ ε, the dynamics of d o (t) becomes

		do = -k w 3 ((1 -α w )	ḋo d o	+ α w	ḋo d e	+ φ * )	(37)
	with α w a continuous valued function in (0, 1) defined by
	(14). The above equation can be rewritten as
			do (t) = -k w 3 b(t)	ḋo (t) d o (t)	+ φ (t)	(38)
	with b(t) = (1-α w (t))d e (t)+α w (t)d o (t) d e (t) following storage function:	and φ (t) = φ * b(t) . Define the
	S =	1 2	ḋ2 o , with Ṡ = -k w 3 b(t) ḋo	ḋo (t) d o (t)	+ φ (t)
	that implies that ḋo is bounded as long as d

o (t) is positive (i.e. σ = 1). Using the fact that for σ = 1 one has d o (t) is positive, and therefore b(t) and φ (t) are uniformly continuous and upper and lower bounded by positive constants. Now by a direct application of [10]-(Th. 5.1), we can prove that d o (t) and ḋo (t) are still converging to zero and d o (t) never crosses zero. Proof of Item 2: Let

T. Hamel is with I3S-CNRS, Institut Universitaire de France, Nice-Sophia Antipolis, France. thamel@i3s.unice.fr

[START_REF] Ho | Optical-flow based self-supervised learning of obstacle appearance applied to mav landing[END_REF] C. Silvestre is with the Department of Electrical and Computer Engineering of the Faculty of Science and Technology of the University of Macau, Macao, China, on leave from the Instituto Superior Técnico, Universidade de Lisboa, Portugal. csilvestre@umac.mo

The authors in *This work was partially supported by the Project MYRG2015-00126-FST of the University of Macau; by the Macao Science and Technology, Development Fund under Grant FDCT/026/2017/A; by Fundac ¸ão para a Ciência e a Tecnologia (FCT) through Project UID/EEA/50009/2013 and Project PTDC/EEI-AUT/5048/2014; and by the EQUIPEX project Robotex. The work of Z. Tang and R. Cunha was supported by FCT through Ph.D. Fellowship PD/BD/114431/2016 under the FCT-IST NetSys Doctoral Program and the FCT Investigator Program (IF/00921/2013), respectively. 1 Z. Tang and 2 R. Cunha are with the Institute for