VEO Weather Routing - algorytm wyznaczania tras w oparciu o dane pogodowe, mapowe i żeglowne jachtów

Wstęp

Wyznaczanie trasy jachtu żaglowego w oparciu o dane pogodowe jest od dawna zagadnieniem istotnym w praktyce żeglowania. Dotychczas zadanie to realizowano przede wszystkim w oparciu o planowanie „w głowie” kapitana jachtu a wiec o dane pogodowe i wiedzę/intuicję żeglarza odnośnie sposobu wykorzystania wiatru i prądów morskich przez konkretny model statku. Wersja zmatematyzowana, algorytmiczna wyznaczania tras morskich była znana bardziej w świecie regatowym, gdzie każdy błąd to gorszy wynik w wyścigu. W regatach i w zastosowaniach turystycznych dominowały dotychczas algorytmy izochronowe. Algorytmy te są wykonywalne na komputerach o niskiej mocy ale posiadają szereg ograniczeń funkcjonalnych. Mianowicie nie omijają złożonej linii brzegowej - a jeśli posiadają obejścia tego problemu to tracą na jakości wyznaczonych tras lub całkowicie przestają działać.

W najbliższej przyszłości nastąpi popularyzacja wykorzystania algorytmów wyznaczania tras morskich w standardowych działaniach nawigacyjnych kapitana jachtu. Podobnie jak to miało miejsce w Google Maps i nawigacji po lądzie, podstawą matematyczną tego procesu będą algorytmy Dijkstry Algorytm Dijkstry musi jednak zostać zaadaptowany do zagadnień nawigacji morskiej. Algorytmem takim jest VEO Weather Routing prezentowany w niniejszym artykule. Logika procedury Dijkstry-ego w systemie VEO pozwala znacznie lepiej omijać lądy oraz strefy wyłączone z ruchu wodnego zyskując przy tym bardzo wysoką dokładność w zakresie rozmiaru kwadratu geograficznego, do którego zmierza pojazd. Odbywa się to kosztem mocy obliczeniowej komputerów, w tym przypadku komputerów wysokiej mocy graficznego przetwarzania danych (GPU).

Podstawy teoretyczne
Głównym narzędziem wykorzystywanym w wyznaczaniu tras jest opracowany w 1959 r. przez Edsger’a Dijkstrę tzw. algorytm Dijkstry [1]. Służy on do wyznaczania najkrótszej ścieżki w grafie obejmującej wszystkie węzły – począwszy od węzła pierwszego (źródła) lub wyznaczeniu najkrótszej ścieżki między dwoma węzłami.

Algorytm Dijkstry dokonuje podziału zbioru wierzchołków grafu na trzy zbiory: V, Q, S. Pierwszy z nich V jest zbiorem wszystkich węzłów w grafie. Q jest zbiorem wszystkich wierzchołków o wspólnej krawędzi z węzłem źródłowym. Do zbioru S trafia węzeł ze zbioru Q o najkrótszej ścieżce od środka po czym przekształca się on w węzeł środkowy, a pierwotny węzeł źródłowy przestaje być uwzględniany w dalszym działaniu algorytmu. Działanie algorytmu kończy się w momencie kiedy zbiór Q jest pusty bądź też kiedy dotrze do wcześniej zdefiniowanego węzła.

Algorytm Dijkstry dla wyznaczenia trasy między dwoma węzłami w grafie można opisać na następującym przykładzie – rysunek 1.

Rysunek 1 Procedura algorytmu Dijkstry

Wierzchołek S stanowi węzeł źródłowy natomiast wierzchołek F końcowy. Każda krawędź w grafie ma nadaną wagę (nieujemną). Działanie algorytmu Dijkstry polega na znalezieniu wierzchołka połączonego z wierzchołkiem źródłowym wspólną krawędzią o najniższej wadze, co na rysunku stanowi połączenie węzła źródłowego S z węzłem A. Następnie węzeł połączony z węzłem źródłowym przejmuje jego rolę (A). Taka iteracja zachodzi do momentu dotarcia do punktu końcowego (S->A->C->E->F). Na rysunku przedstawiona została (pogrubioną linią) ścieżka o najniższej sumie wag (4+3+3+2=12). Ich suma wyniosła 12 i przedstawia optymalną trasę (minimalny koszt dotarcia) między węzłami S i F.

Algorytm Dijkstry jest powszechnie stosowany w programach do wyznaczania tras dla jednego pojazdu. Na bazie tego algorytmu pracuje m.in. system wyznaczania tras w Google Maps.

Adaptacja morska algorytmu Dijkstry
Aby wykonać adaptację algorytmu Dijkstry do nawigacji morskiej, konieczne jest uwzględnienie następujących aspektów:

Zapewnić sposób ważenia odcinków trasy adekwatnych do specyfiki nawigacji morskiej a więc:
Uwzględnić, że pewne krawędzie będą szybsze od innych dzięki optymalnemu dopasowaniu cech żeglownych jachtu i wiatru, który wieje w danym kwadracie geograficznym w danym czasie
Uwzględnić, że docieranie jachtu do granicy kwadratu geograficznego wymaga podania nowych danych pogodowych z kolejnego okresu czasu
Zapewnić mapę zawierającą obszary po których można pływać bezpiecznie bez wpłynięcia na ląd/płyciznę
Zapewnić mapę zawierającą obszary, po których nie można pływać bezpiecznie:
żeglowanie jest ryzykowne (np. płycizny, prognozowana w danym oknie czasowym burza)
żeglowanie jest niekomfortowe (wysoka fala, silny rozkołys)
żeglowanie jest zabronione (łowiska, strefy wojskowe, inne)
Uwzględnić wpływ prądów morskich pokrywających się z pewnymi krawędziami rejsu i dodającymi jachtu prędkości
Uwzględnić znos jachtu a wiec konieczność ustawienia innego kąta żeglowania korygującą znos jachtu przez wiatr i falę
Zapewnić tak szczegółową kwadratyzację mapy aby możliwe było precyzyjne omijanie niewielkich przeszkód wodnych – np. obiektów o szerokości 100 metrów.
Zapewnić procedurę wyznaczania trasy w oparciu o podejście interacyjne – najpierw zgrubna topologia punktów brzegowych, później korekta o szczegółowe uwarunkowania mapowe i pogodowe.
Uwzględnienie części tych wytycznych opisano w dziale Implementacja. Na kolejnych rozwojach algorytmu wszystkie z nich zostaną w pełni zaimplementowane.

Implementacja
Rozwiązanie zostało zaimplementowane jako zestaw tabel oraz procedur składowanych bazy PostgreSQL z rozszerzeniem przestrzennym PostGIS oraz maszyną routingu PG_Routing.

Ze względów wydajnościowych działanie algorytmu podzielone zostało na dwie fazy:

W fazie wstępnej przygotowujemy wszystkie dane globalne oraz niezmienne w czasie, wykonywane są również bardziej czasochłonne operacje takie jak obliczanie zawartości lądu w oczku siatki, przygotowanie polar plotów jachtów, wyznaczenie krawędzi dla routingu wstępnego czy wyznaczenie wierzchołków grafu topologicznego używanego w późniejszych krokach do wyznaczania finalnej trasy
W fazie właściwej wykonywane są już tylko operacje niezbędne do wyznaczenia trasy i wyłącznie na fragmencie przestrzeni w którym trasa będzie wyznaczana. Pozwala to na znaczne przyspieszenie działania algorytmu i nie generuje niepotrzebnych operacji na fragmentach przestrzeni które nie będą odwiedzane.
W trakcie działania serwera wykonywane są również operacje cykliczne takie jak pobieranie danych pogodowych i aktualizacja ich w bazie danych.

Proces obliczeniowy składa się z dwóch faz: fazy przygotowania środowiska oraz z fazy wyznaczenia trasy.

Faza przygotowania środowiska
Przygotowanie tabeli z lądami
Tabela inicjalna z lądami tworzona jest na podstawie danych Natural Earth [3], które są importowane do bazy danych za pomocą QGIS[7].

Przygotowanie tabeli z obszarami wykluczonymi
Obszary wykluczone tworzone są na podstawie danych projektu OpenStreetMap[8] przez społeczność zgromadzoną wokół projektu OpenSeaMap[9].  Ekstrakt danych jest pobierany i filtrowany przy pomocy aplikacji osmosis, po czym wynikowy plik importowany do bazy danych procesorem osm2pgsql. Zaimportowane dane przepisywane są do tabeli zgodnych ze schematem rozwiązania.

Przygotowanie tabeli z prędkościami jachtów w stosunku do wiatru
Polar plot to wykres wiążący cechy żeglowne jachtu z siłą i kierunkiem wiatru który na niego oddziałuje, pozwalający odczytać jaką prędkość osiągnie jacht w przy danym kierunku i sile wiatru.

Przykładowy polar plot wygląda następująco:

Rysunek 2 Przykład polar plotu [4]

Tabela z prędkościami jachtów w stosunku do wiatru tworzona jest na podstawie danych dostępnych w repozytorium [5]. W ramach projektu utworzone zostało dedykowane narzędzie pozwalające na pobieranie tych danych i import ich do bazy w strukturze zgodnej ze schematem rozwiązania.

Przygotowanie tabeli z danymi pogodowymi
Na potrzeby rozwiązania utworzona została usługa, która cyklicznie pobiera aktualne dane pogodowe w formacie GRIB ze źródeł NOAA GFS[6], wyodrębnia z nich wartości konieczne do prawidłowego działania algorytmu, po czym importuje je do bazy danych.

Utworzenie tabeli z kwadratami o wielkości 0.5 stopnia
Na podstawie przeprowadzonych doświadczeń określiliśmy, że kompromisowa wielkość maksymalna kafelka to 0.5 stopnia (w przyszłości, w miarę zwiększania mocy obliczeniowych zmniejszymy ten kafelek co najmniej 10 krotnie). Wielkość ta użyta jest do utworzenia wstępnej siatki podziału.

Rysunek 3 Rozmiar kafli 0.5 stopnia

Wstawienie do tabeli fragmentów lądu przecinających się z kwadratem
Dla każdego oczka siatki wygenerowanego w poprzednim kroku wydzielany jest fragment lądu pokrywający się z nią. Geometria lądu w późniejszych krokach jest używana do określenia zawartości lądu w oczku, a podczas dalszego dzielenia siatki używana jest jako źródło lądu do mniejszych oczek siatki. To podejście znacznie skróciło czas generowania danych początkowych.

Określenie zawartości lądu w kwadracie
Na podstawie porównania pola powierzchni oczka siatki oraz pola powierzchni zawartego w niej lądu oraz obszarów wykluczonych określany jest poziom wypełnienia oczka siatki lądem w trzech wartościach:

f – full – te oczka będą pomijane przy dalszym processingu
p – partial – te oczka będą nadal dzielone w celu wyodrębnienia z nich części nie zawierających lądu
n – none – te oczka uznawane są za prawidłowe i nie są dalej dzielone

Rysunek 4 Identyfikacja kafli z lądem, z lądem i wodą, z samą wodą

Rekurencyjne dynamiczne dzielenie kwadratów o częściowej zawartości lądu
Kolejny krok ma za zadanie ekstrakcję z oczek siatki częściowo zawierających ląd ich fragmentów w których ląd nie występuje. Dzięki temu wydzielamy fragmenty przestrzeni po których jacht może się bezpiecznie poruszać. Wszystkie oczka siatki o częściowej zawartości lądu i obszarów wykluczonych dzielone są dalej, oczka o pełnej zawartości lądu są odrzucane, a oczka bez lądu pozostawiane jako prawidłowe. Proces odbywa się w rekurencyjnej pętli aż do osiągnięcia doświadczalnie określonej kompromisowej wielkości minimalnej oczka siatki. Krok ten zakłada kompromis między czasem generowania zasobu inicjalnego, czasem generowania trasy oraz dokładnością tej trasy i możliwością wyznaczenia jej na obszarach gęsto pokrytych wyspami lub w wąskich cieśninach. W kolejnych fazach rozwoju aplikacji kompromis ten może być przesunięty w kierunku dokładności w miarę uzyskiwania wyższej wydajności na etapie wyznaczania trasy.

Rysunek 5 Rekurencyjne wydzielenie kwadratów z lądem i wodą na drobniejsze elementy

Uzupełnienie relacji do tabeli zawierającej dane pogodowe
W kolejnym kroku na podstawie relacji przestrzennej tworzona jest relacja klucza obcego do tabeli zawierającej dane pogodowe. Identyfikator oczka tabeli pogodowej przepisywany jest na wszystkie elementy pochodne (krawędzie grafu topologicznego), co zapobiega wykonywaniu wielokrotnych złączeń przestrzennych i znacząco wpływa na wydajność końcowej procedury wyznaczania trasy.

Rysunek 6 Dodanie do kafli powiązania z danymi pogodowymi

Utworzenie tabeli ze zgrubnymi krawędziami
W kolejnym kroku między centroidami oczek siatki nie zawierającymi lądu wyznaczany jest graf topologiczny używany w pierwszym kroku procedury wyznaczania trasy do określenia najkrótszej trasy z pominięciem warunków pogodowych.

Utworzenie głównej siatki
W kolejnym kroku na podstawie siatki zgrubnej tworzona jest siatka główna przez podział multilinii tworzących oczka siatki na poszczególne linie, które w kolejnych krokach będą dzielone i używane do utworzenia grafu topologicznego używanego do wyznaczania docelowej trasy uwzględniającej warunki pogodowe.

Utworzenie wierzchołków topologii
Jako ostatni krok fazy przygotowania danych na utworzonej siatce tworzone są punkty używane w kolejnych krokach wyznaczania trasy jako wierzchołki grafu topologicznego po których odbywał się będzie routing. Punkty tworzone są przez podział każdej krawędzi oczka siatki na 5 części (parametr dobrany doświadczalnie).

Rysunek 7 Topologia punktów brzegowych kafli

Faza wyznaczania trasy
Wyznaczenie najkrótszej drogi po zgrubnej topologii
Pierwszym krokiem wyznaczania trasy jest wyszukanie oczek siatki przecinających się z punktem początkowym i końcowym oraz wyznaczenie zgrubnej trasy. Przeprowadzone doświadczenia pozwoliły na określenie że wzrost długości trasy ze względu na warunki pogodowe opłacalny jest jedynie w zakresie 30%, dlatego trasa końcowa musi bezpośrednio korespondować z trasą najkrótszą. Jako koszt przebycia krawędzi zastosowana została wyłącznie jej długość, a trasa wyznaczona jest za pomocą funkcji implementującej algorytm Dijkstra w rozszerzeniu PG_Routing.

Rysunek 8 Uproszczona trasa po zgrubnej topologii punktów brzegowych

Utworzenie tabeli z krawędziami
Na obecnym etapie wiemy już którędy zgrubnie prowadzi najkrótsza trasa między zadanymi punktami. W kolejnym kroku na terenie oczek siatki przebytych w ramach najkrótszej trasy wyznaczany jest gęsty graf topologiczny po którym będziemy wyznaczać docelową trasę z uwzględnieniem warunków pogodowych. Wierzchołki topologii na terenie oczek przecinających trasę zgrubną utworzone w fazie przygotowania danych łączone są krawędziami.

Rysunek 9 Szczegółowa topologia punktów brzegowych odniesiona do trasy zgrubnej

Dołączenie do topologii punktów początkowego i końcowego trasy
Jako że poprzednie wyznaczenia wykonywane były na podstawie centroidów przebywanych oczek siatki do przeprowadzenia docelowego wyszukania trasy konieczne jest podłączenie właściwych punktów początku i końca trasy do grafu topologicznego. W tym celu krawędzie na terenie oczek siatki przecinających się z punktem początkowym i końcowym są usuwane, a tworzone są nowe krawędzie które mają swój początek w punkcie początkowym trasy, a końce na uprzednio przygotowanych wierzchołkach topologii.

Rozszerzenie grafu o kafle sąsiadujące
Aby umożliwić wyznaczenie najszybszej trasy korespondującej z trasą zgrubną musimy rozszerzyć topologię o kafle sąsiadujące z tymi po których przebiega ta trasa. Procedura rozszerzająca istniejący zestaw krawędzi przygotowana jest do działania rekurencyjnego i może wielokrotnie rozszerzać topologię. Ten punkt jest kolejnym miejscem kompromisu – każde rozszerzenie grafu wydłuża czas wyznaczania końcowej trasy, ale równocześnie zwiększa prawdopodobieństwo odnalezienia trasy znacznie efektywniejszej w większej odległości od trasy zgrubnej. W kolejnych krokach rozwoju algorytmu krok ten można będzie również uzależnić od wyniku wyznaczania – jeśli trasa docelowa nie zostanie odnaleziona można do niego wrócić, rozszerzyć topologię i spróbować wyznaczyć trasę ponownie.

Rysunek 10 Szczegółowa topologia punktów brzegowych odniesiona do trasy zgrubnej

Dodanie do krawędzi informacji o pływach i wiatrach
Do wyznaczonych krawędzi dodajemy informacje o pływach i wiatrach na podstawie pobranych prognoz pogody. Informacja pogodowa przesunięta jest w czasie na podstawie szacowanego czasu dotarcia do każdej z krawędzi. W kolejnym kroku obliczany jest azymut każdej krawędzi oraz jej kąty w stosunku do wiatru zarówno w przypadku podróżowania zgodnie jak i przeciwnie do jej kierunku – parametry te będą konieczne do określenia czasu przebycia krawędzi.

Obliczenie czasu przebycia każdej krawędzi w obu kierunkach
Ostatnią operacją konieczną do przygotowania danych do wyznaczenia trasy jest obliczenie czasu przebycia każdej krawędzi w obu kierunkach. Czasy przebycia krawędzi interpolowane są na podstawie tabeli średnich prędkości jachtów dla zadanego typu jachtu, siły i kierunku wiatru, która została zaimportowana do systemu w fazie przygotowania danych.

Na przedstawionym obrazku kolor krawędzi odzwierciedla czas w jakim krawędź ta może być przebyta.

Rysunek 11 Wygenerowanie różnych tras i czasów przebycia oznaczonych kolorem

Na tym etapie najlepiej można dostrzec różnice z popularnym na rynku algorytmem izochronowym. Trasy wyznaczone algorytmem Dijkstry znajdują się w kaflach powiązanych z trasą zgrubną, na punktach brzegowych o wysokiej gęstości podczas gdy w algorytmach izochronowych  trasy te znajdują się w lejkach generowanych w pewnym kącie rozwarcia, bez powiązania z kaflami geograficznymi i ich punktami brzegowymi. Przykład lejka izochronowego w programie MAX SEA TIMEZERO:

Rysunek 12 Algorytm izochronowy dla podobnej trasy w programie MAX SEA TIMEZERO

Wyszukanie najszybszej trasy po krawędziach
Zwieńczeniem procesu przygotowania danych jest wyznaczenie trasy docelowej na przygotowanych krawędziach. Do wyznaczenia trasy – tak samo jak w przypadku trasy zgrubnej używana jest procedura implementująca algorytm Djikstry z rozszerzenia PG_Routing. W tym przypadku kosztem przebycia krawędzi jest czas jej przebycia obliczony w poprzednim kroku.

Poniższy obrazek przedstawia wyznaczoną trasę zgrubną na tle krawędzi z obliczonym czasem ich przebycia.  

Rysunek 13 Wyznaczenie trasy optymalnej (kolor jasno zielony)

Porównanie trasy zgrubnej oraz trasy docelowej na tle danych pogodowych przedstawia się następująco.

Rysunek 14 Trasa zgrubna i trasa zoptymalizowana (jasno zielona) wyznaczona automatycznie.

Przykłady obliczeń
Poniżej przedstawiono 3 kalkulacje dla przykładowych tras, dla jachtu Salona 44, dla przykładowej pogody dla wybranej daty.

Parametr „Czas rejsu na trasie uproszczonej” reprezentuje w przybliżeniu czas rejsu dla trasy wyznaczonej ręcznie, w oparciu o manualnie wyznaczone waypointy, w przeliczeniu na typową prędkość jachtu turystycznego tej klasy. Parametr „Czas rejsu na trasie zoptymalizowanej” przedstawia efekt działania algorytmu VEO Weather Routing dla kwadratyzacji map rzędu 0,5 stopnia (wraz z rekurencyjną korektą dla zadanego obszaru mapy).  Docelowo kwadratyzacja map sięgnie 0,005 stopnia przy co najmniej 20-krotnym przyśpieszeniu obliczeń dzięki zastosowaniu procesorów GPU o bardzo wysokich parametrach przetwarzania. Istnieje dalszy potencjał na zmniejszenie kwadratów mapy podlegających analizie i dalsze przyspieszenie algorytmów (zwiększania mocy obliczeniowej komputera głównego).

Parametr „Czas kalkulacji algorytmu na środowisku klasy średni PC” należy traktować jako informację porównawczą, dającą jedynie wyobrażenie o typowym obciążeniu obliczeniowym z perspektywy przeciętnych komputerów osobistych.

Trasa morska od-do

Czas rejsu na trasie uproszczonej (topologia zgrubna)

Czas rejsu na trasie zoptymalizowanej

(topologia szczegółowa)

Czas kalkulacji algorytmu na środowisku klasy średni PC

Gdynia (marina) - Łeba (port) – ok. 80 mil

16,5h

14,3h

12,8s

Ystad (Szwecja) - Rosklide (Dania) – ok. 150 mil

31,2h

26,1h

15,4s

Krk (miasto, wyspa KRK Chorwacja) - Gujak, Kornat (Chorwacja) – ok. 80 mil

16,8h

13,4h

12,8s

Tabela 1: Porównanie wyników działania algorytmu VEO Weather Routing

Podsumowanie
W artykule wykazano przydatność algorytmu Dijkstry w adaptacji morskiej zbudowanej przez autorów o nazwie VEO Weather Routing. Mechanizm ten służy do rozwiązywania problemu wyznaczania trasy morskiej pod kątem danych pogodowych, cech żeglownych jachtu – z automatycznym omijaniem lądów (i stref wyłączonych z żeglugi) o wysokiej efektywności i precyzji. Przedstawiono pełną procedurę obliczeniową w obu fazach działania algorytmu. W efekcie udowodniono istotny uzysk czasowy na trasie wyznaczonej za pomocą zbudowanego mechanizmu względem trasy wyznaczonej ręcznie. 

Literatura 

1. Dijkstra E. W., A note on two problems in connexion with graphs, nr 1(1), str. 269-271, Numerische mathematik, 1959.

2. Tan Y., Wu D., Research on Optimization of Distribution Routes for Fresh Agricultural Products Based on Dijkstra Algorithm, nr 336, str. 2500-2503, Applied Mechanics and Materials, 2013.

3. Strona internetowa: https://www.naturalearthdata.com/

4. Strona internetowa Wikipedia: https://en.wikipedia.org/wiki/Polar_diagram_(sailing)

5. Strona internetowa http://jieter.github.io/orc-data/site/

6. Strona internetowa https://www.ncei.noaa.gov/products/weather-climate-models/global-forecast

7. Strona internetowa i system QGIS: https://www.qgis.org/pl/site/

8. Strona internetowa projektu OpenStreetMaps: https://www.openstreetmap.org/

9. Strona internetowa projektu OpenSeaMaps: https://www.openseamap.org/index.php?id=openseamap&no_cache=1

Tytuł i streszczenie w języku polskim:

Tytuł

VEO Weather Routing - nowatorski algorytm wyznaczania tras w oparciu o dane pogodowe, mapowe i żeglowne jachtów 

Streszczenie

W artykule przedstawiono algorytm Dijkstry w adaptacji morskiej o nazwie VEO Weather Routing. Mechanizm ten służy do rozwiązywania problemu wyznaczania trasy morskiej pod kątem danych pogodowych, cech żeglownych jachtu – z automatycznym omijaniem lądów o wysokiej efektywności i precyzji. Przedstawiono pełną procedurę obliczeniową w obu fazach działania algorytmu. W efekcie udowodniono istotny uzysk czasowy na trasie wyznaczonej za pomocą zbudowanego mechanizmu względem trasy wyznaczonej ręcznie. 

Słowa kluczowe: Pogodowe wyznaczanie tras morskich, nawigacja morska, nawigacja jachtowa

Title

VEO Weather Routing - an innovative algorithm for calculating routes based on weather, map and navigable data of yachts

Abstract (j. ang)

The article presents the Dijkstra algorithm in a marine adaptation called VEO Weather Routing. This mechanism is used to solve the problem of determining the sea route in terms of weather data, navigable characteristics of the yacht - with automatic avoiding land with high efficiency and precision. The full calculation procedure in both phases of the algorithm's operation was presented. As a result, a significant time gain was proved on the route marked with the use of the constructed mechanism in relation to the route marked by hand.

Keywords: Weather routing, sea navigation, yacht navigation