Waldemar Izdebski - Wykłady z przedmiotu SIT 81
10. Numeryczny model powierzchni terenowej
Jednym z zagadnień SIP jest przedstawienie zjawisk o charakterze ciągłym jak np.
powierzchnia terenu. W ogólnym przypadku kiedy zjawisko możemy przedstawić funkcją
analityczną postaci:
z = f (x, y)
zagadnienie nie stanowi żadnych trudności, ponieważ dzięki znanej funkcji w każdym
potrzebnym punkcie P(x,y) możemy określić wartość zjawiska. W sytuacji kiedy
modelowanych zjawisk (w szczególności powierzchni terenu) nie można określić funkcją
analityczną stosujemy inne rozwiązania, oparte na wartościach zjawiska zarejestrowanych w
wybranych punktach pomiarowych. Najczęściej stosowanymi metodami przestrzennej
reprezentacji powierzchni (zjawisk) są:
" reprezentacja elementami punktowymi, dla których określono wartość zjawiska
określono w regularnej siatce kwadratów (ang grid),
" reprezentacja elementami liniowymi, dla których wartość zjawiska jest określona i
niezmienna (izolinie),
" reprezentacja w postaci elementów powierzchniowych będąca siecią
nieregularnych trójkątów TIN (ang. triangular irregular network) opartych na
punktach pomiarowych.
Schematycznie wymienione metody reprezentacji powierzchni przedstawiono na rysunku 9.1.
243
237 242
238 241
237 239 240
236
236
235
234
232 233
235 233
236 232
231
231
230
229
228
227
226
226
Rys. 10.1. Metody reprezentacji zjawisk o charakterze ciągłym
W związku z dyskretną reprezentacją powierzchni, z każdą z wymienionych wyżej
metod muszą być związane odpowiednie algorytmy interpolacyjne umożliwiające określenie
wartości zjawiska w dowolnie wybranym punkcie.
W systemach informacji przestrzennej podstawowe znaczenia ma reprezentacja
powierzchni terenu i temu zagadnieniu poświęcimy dalsze rozważania.
Numeryczny Model Terenu definiuje się jako  numeryczną reprezentację powierzchni
terenowej, utworzonej poprzez zbiór odpowiednio wybranych punktów leżących na tej
Waldemar Izdebski - Wykłady z przedmiotu SIT 82
powierzchni oraz algorytmów interpolacyjnych umożliwiających jej odtworzenie w
określonym obszarze [Gaz
dzicki 1990].
Idealne odtworzenie powierzchni terenu przez model nie jest możliwe, ponieważ ze
względów ekonomicznych, czasowych i wielkości zbiorów danych, nie da się pomierzyć ani
wyrazić całej złożoności powierzchni terenu. Podstawowymi problemami związanymi z
numerycznym modelem terenu są:
" problem odpowiedniego doboru charakterystycznych punktów powierzchni (ang.
sampling problem) w celu uzyskania jak najlepszego efektu przy minimalizacji ilości
danych,
" problem odtworzenia (przedstawienia) powierzchni na podstawie istniejących danych
(ang. representation problem).
10.1. Rodzaje modeli DTM
W praktyce podstawowe znaczenie mają dwa modele: regularny w postaci siatki
kwadratów uzupełnione charakterystycznymi punktami i liniami szkieletowymi oraz w
postaci nieregularnej siatki trójkątów (TIN). Każdy z wymienionych modeli posiada swoje
zalety i wady, które przesądzają o ich zastosowaniach. Istotą TIN jest np. przechowywanie
oryginalnych danych pomiarowych podczas gdy w modelu grid wysokości w punktach
węzłowych przeważnie są już interpolowane.
10.1.1. Model w postaci siatki kwadratów
Model oparty jest na siatce kwadratów, której punkty węzłowe posiadają określone
zij
wysokości powierzchni terenowej .
z
Pij
x
zij
i
j
P
ij
y
Rys. 10.2. Ilustracja modelu w postaci siatki kwadratów
Waldemar Izdebski - Wykłady z przedmiotu SIT 83
zij
Wysokości węzłów przechowywane w strukturze macierzy gdzie i j określają wiersz i
kolumnę macierzy. Znając interwał siatki ds , położenie jej punktu początkowego X Yo
o
możemy zawsze dla każdego węzła i j określić współrzędne terenowe X Y .
X = X + i *ds
o
Y = Yo + j *ds
Struktura taka jest wyjątkowo łatwa do przetwarzania, zabiera bardzo mało pamięci, a
algorytmy używane do modelowania terenu są stosunkowo proste. Im gęstsza siatka zostanie
zastosowana tym otrzymany model będzie dokładniejszy. Zwiększając gęstość siatki
prowadzi jednak do sytuacji, że jest ona również zwiększana w miejscach o małym
urozmaiceniu terenu, powodując tym samym znaczny wzrost nic nie wnoszących danych.
Rozwiązaniem jest uzupełnienie struktury o punkty charakterystyczne i linie szkieletowe lub
zastosowanie siatki o strukturze hierarchicznej dostosowującej gęstość do stopnia
skomplikowania rzez
by.
Wysokości w punktach węzłowych mogą pochodzić bezpośrednio z pomiaru (bezpośredniego
lub fotogrametrycznego) lub być wyznaczane z innych modeli powierzchni terenowych.
Mając model w postaci regularnej siatki kwadratów możemy interpolować wysokości w
punktach pośrednich. Wykorzystuje się do tego celu aproksymację powierzchni terenowej w
obrębie każdego kwadratu paraboloidą hiperboliczną
z = auvuv + au + avv + a
gdzie u i v są współrzędnymi układu powstałego przez równoległe przesuniecie układu x,y do
węzła i, j siatki. Wyznaczenie wysokości punktu P paraboloidy hiperbolicznej w obszarze
kwadratu może być wykonane wzorem interpolacyjnym, mającym postać ogólnej średniej
arytmetycznej:
zi, jwi, j + zi+1, jwi+1, j + zi, j+1wi, j+1 + zi+1, j+1wi+1, j+1
zp =
wi, j + wi+1, j + wi, j+1 + wi+1, j+1
gdzie zi, j , zi+1, j , zi, j+1, zi+1, j+1 są wysokościami w punktach węzłowych a wagi
wi, j , wi+1, j , wi, j+1, wi+1, j+1są polami przeciwległych prostokątów, powstających przez
Waldemar Izdebski - Wykłady z przedmiotu SIT 84
podział kwadratu liniami równoległymi do boków kwadratu i przechodzących przez rzut P
punktu P na płaszczyz
nie x,y.
u
x
ds
i+1,j+1
i+1,j
wi,j+1
wi,j
up
P
ds
wi+1,j
wi+1,j+1
v
vp
i,j
i,j+1
y
Rys. 10.3. Ilustracja interpolacji wysokości w siatce kwadratów
Cechą charakterystyczną paraboloidy hiperbolicznej jest to, że przecinając ją płaszczyznami
pionowymi u=const, lub v=const otrzymuje się linie proste.
10.1.2. Model w postaci nieregularnej sieci trójkątów
Nieregularna sieć trójkątów powstaje głownie jako efekt bezpośrednich pomiarów
terenowych, gdzie cały zakres opracowania zapełnia się trójkątami opartymi o punkty
pomiarowe. Ponieważ w tych modelach wykorzystywane są wszystkie punkty
charakterystyczne model jest stosunkowo dokładny.
Do tworzenia siatki trójkątów najczęściej wykorzystywana jest triangulacja Delaunay a.
Trójkąty tworzone są w ten sposób aby żaden z punktów nie należących do niego nie był
położony wewnątrz okręgu opisanego na trójkącie. Ilustrację zasady przedstawiono na
poniższym rysunku.
Waldemar Izdebski - Wykłady z przedmiotu SIT 85
Rys. 10.4. Ilustracja triangulacji Delaunay a
Na kolejnym rysunku przedstawiono konstrukcję diagramu (poligonu, obszaru) Voronoi,
bezpośrednio związaną z triangulacją Delaunay a. Obszar Voronoi stanowi zbiór wszystkich
punktów płaszczyzny, dla których odległość do punktu centralnego jest mniejsza od
odległości do pozostałych punktów. Jak widać ograniczenia tego obszaru stanowią odcinki
symetralnych do boków triangulacji Delaunay a.
Rys. 10.5. Ilustracja diagramu Voroni
Nazwa diagramu pochodzi od nazwiska matematyka M. G. Voronoi, który opublikował pracę
w 1908 roku. Czasami używa się również innych terminów:
Tesselacja Dirichleta od nazwiska matematyka G. L. Dirichleta (1850)
Obszar Thiessena : od klimatologa Thiessena (A. J. Thiessen, J. C. Alter  1911)
Celki Wignera-Seitza (E. Wigner, F. Seitz  1933)
Transformacja Bluma (H. Blum  1967).
Zagadnienie tworzenia nieregularnej siatki trójkątów było przedmiotem licznych badań, w
wyniku których opublikowano szereg algorytmów konstrukcji triangulacji Delaunay a:
Radial sweep  A. Mirante, N. Weingarten 1982
Recursive split  B. A. Lewis, J. S. Robinson 1978
Divide and conquer  D. T. Lee, B. J. Schachter 1980
Step by step  opublikowany m.in. przez M. J. McCullagh, C. G. Ross 1980
Hierarchical  L. De Floriani, B. Falcidieno, C. Pienovi, 1985
Incremental  opublikowany m.in. przez D. T. Lee, B. J. Schachter 1980
Waldemar Izdebski - Wykłady z przedmiotu SIT 86
Incremental delete and build  D. F. Watson 1980
Poniżej omówiono algorytm krok po kroku - step by step . Polega on na konstrukcji
triangulacji od jednego z narożników obszaru, utworzenia początkowego boku i
poszukiwaniu trzeciego punktu trójkąta. Następnie nowe dwa boki traktowane są jako boki
 początkowe . W ten sposób dodaje się kolejne punkty, aż do włączenia do triangulacji
wszystkich punktów. Oto kolejne etapy:
1. Wybiera się pierwszy bok, w ten sposób, aby żaden punkt nie znalazł się we wnętrzu
okręgu, którego średnicą jest dany bok. Można to osiągnąć wybierając dwa najbliższe
punkty w danym miejscu. Jeśli to możliwe początkowy bok powinien znajdować się
najbliżej jednej z granic obszaru.
2. Poszukuje się następnego punktu. W kole, które przechodzi przez końce odcinka
początkowego i trzeci punkt nie powinien się znalez
ć żaden inny punkt. Cel ten osiąga
się poprzez obliczanie odpowiednich kątów ą. Jako punkt tworzący trójkąt wybiera
się ten dla którego kąt ą jest największy.
3. Odcinki łączące wybrany punkt i końce odcinka początkowego stanowią drugi i trzeci
bok trójkąta.
4. Nowe boki zostają bokami początkowymi do poszukiwania następnych punktów.
Powtarzanie operacji spowoduje zabudowę trójkątami całego obszaru opracowania.
ą ą
ą
ą
Rys. 10.6. Zasada algorytmu step by step
Niekiedy przy konstruowaniu triangulacji można znalez
ć dwa lub więcej punktów, dla
których kąty są jednakowe. Wówczas przez cztery punkty przechodzi jeden okrąg (np. cztery
punkty tworzące prostokąt). Należy wybrać jeden z wariantów analizując dodatkowe
informacje.
Algorytm tworzenia nieregularnej siatki trójkątów musi mieć możliwość zmiany boków
trójkątów na podstawie znajomości spadków ze szkiców polowych oraz musi z założenia
prowadzić boki wzdłuż linii szkieletowych grzbietowych i ściekowych. Dane wygenerowane
w procesie triangulacji należy przechowywać w strukturach zapewniających szybki i
wygodny dostęp. Jednym z wariantów jest zapis danych w postaci trzech tablic:
punktów ze współrzędnymi X,Y,H,
boków ze wskaz
nikami do punktów i przyległych trójkątów,
trójkątów ze wskaz
nikami do boków.
Waldemar Izdebski - Wykłady z przedmiotu SIT 87
F
9
8
T4
E
4
B
7
T2 5
T3
3
2
T1
6 D
C
1
A
PUNKT X Y H BOK Pkt1 Pkt2 Bok1 Bok2 Bok3
"1 "2 "
A . . . 1 T1 0 A C T1 1 2 3
B . . . 2 T1 0 A B T2 3 4 5
C . . . 3 T1 T2 B C T3 5 6 7
D . . . 4 T2 T4 B E T4 4 8 9
E . . . 5 T2 T3 C E
F . . . 6 T3 0 C D
7 T3 0 D E
8 T4 0 B F
9 T4 0 E F
Rys. 10.7. Przykład zapisu danych
Mając model w postaci nieregularnej siatki trójkątów możemy interpolować wysokość punktu
P powierzchni terenowej, którego rzut P na płaszczyznę poziomą leży w trójkącie A B C .
x
B
wC
P
wB wA
A
y
C
Rys. 10.8. Ilustracja wyznaczania wysokości punktu w nieregularnej siatce trojkątów
Aproksymując powierzchnię terenową płaszczyzną przechodzącą przez punkty ABC
otrzymujemy wzór interpolacyjny w postaci:
zAwA + zBwB + zCwC
zp =
wA + wB + wC
Waldemar Izdebski - Wykłady z przedmiotu SIT 88
gdzie zA, zB , zC , są wysokościami w punktach A B C a wagi wA, wB , wC są polami
przeciwległych trójkątów, powstających przez podział trójkąta A B C odcinkami łączącymi
wierzchołki tego trójkąta z punktem P .
10.2. Tworzenie numerycznego modelu terenu
Dane do stworzenia numerycznego modelu terenu uzyskiwane są przede wszystkim z
trzech z
ródeł:
pomiarów bezpośrednich,
pomiary fotogrametryczne,
digitalizacji istniejących map.
Niekiedy wykorzystuje się również altimetrię radarową lub laserową, dla modeli
geologicznych wiercenia lub pomiary sejsmiczne.
10.2.1. Bezpośrednie pomiary terenowe
Pomiary bezpośrednie charakteryzują się wysoką dokładnością, a punkty wysokościowe
(pikiety) w łatwy sposób są wprowadzane do systemów informatycznych. Pomiary te są
jednak pracochłonne i kosztowne. Na ich podstawie otrzymuje się model nieregularny. Przy
pozyskiwaniu punktów należy zwracać uwagę na to, by dobrze charakteryzowały
powierzchnię terenu. Ilustrację przedstawiono na poniższym rysunku.
Rys. 10.9. Ilustracja wyboru punktów charakterystycznych
Punkty z lewej strony rysunku leżą na środku pochyłości lub na płaszczyzny, mogą
zostać zatem pominięte bez straty informacji o powierzchni terenu. Punkty położone po
stronie prawej leżą na szczycie wzniesienia, na punkcie przegięcia (zmiana spadku) i w
najniższym punkcie terenu. Punktów tych nie można pominąć bez utraty dokładności.
Na szkicach pomiarowych terenowych należy koniecznie zaznaczać również wszystkie
linie nieciągłości (brzegi skarp, urwiska), jak i przebieg linii szkieletowych.
Nieuwzględnienie przebiegu tych linii może prowadzić do powstania dużych błędów.
Waldemar Izdebski - Wykłady z przedmiotu SIT 89
Rys. 10.10. Ilustracja istotności linii strukturalnych
Na rysunku lewym brak jest informacji o przebiegu linii szkieletowych. Obliczenie
wysokości jako średniej arytmetycznej daje wysokość równą 10.5m. Na rysunku środkowym
znane jest położenie linii ściekowej (wysokość w zaznaczonym punkcie wynosi 10m), na
rysunku prawym znane jest położenie linii grzbietowej (wysokość 11m).
10.2.2. Pomiary fotogrametryczne
Za pomocą instrumentów fotogrametrycznych lub zaawansowanych programów
komputerowych możliwe jest automatyczne pozyskiwanie wysokości na zbudowanym
modelu. Najczęściej wysokości są pozyskiwane na siatce prostokątów lub kwadratów.
Ponieważ w ten sposób pozyskiwane wysokości nie oddają w pełni złożoności form
terenowych można zastosować automatyczne zagęszczanie siatki przy dużych zmianach
wysokości. Przy interwencji operatora możliwe jest pozyskiwanie linii strukturalnych.
Rysunki poniższe pokazują techniki zbierania wysokości.
Rys. 10.11. Przykłady rejestracji danych na instrumentach fotogrametrycznych
Próbkowanie regularne: Może być wykonywane jako profile lub w siatce kwadratów (grid).
Zaletą jest możliwość całkowitego zautomatyzowania pozyskiwania wysokości. Wadami jest
ograniczenie do terenów o małych zmianach wysokości. Liczba pozyskanych punktów jest
nieadekwatna do terenu: na terenach płaskich zbyt duża i za mała na terenach pofałdowanych.
Waldemar Izdebski - Wykłady z przedmiotu SIT 90
Metoda generuje zbyt dużą liczbę punktów, ponieważ gęstość siatki musi być mała, by
uniknąć dużych błędów.
Próbkowanie progresywne: Przy pozyskiwaniu wysokości dokonywana jest analiza i w
zależności od zmian wysokości gęstość próbkowania ulega zmianie. Zaletą jest operowanie
na mniejszej liczbie punktów przy wyższej dokładności.
Próbkowanie selektywne. Pozyskuje się dodatkowo linie strukturalne. W połączeniu z
próbkowaniem progresywnym nosi nazwę próbkowania kompozytowego. Zaletą jest wyraz
ne
poprawienie modelu terenu. Niedogodność stanowi konieczność interwencji operatora, tak
więc metoda jest jedynie częściowo automatyczna.