Filtracja danych lotniczego skaningu laserowego z wykorzystaniem metody aktywnych powierzchni

ROCZNIKI GEOMATYKI 2005 m TOM III m ZESZYT 4

FILTRACJA DANYCH LOTNICZEGO SKANINGU

LASEROWEGO Z WYKORZYSTANIEM METODY

AKTYWNYCH POWIERZCHNI

1

FILTERING OF AERIAL LASER SCANNER DATA

BY ACTIVE SURFACE MODEL

Andrzej Borkowski

Katedra Geodezji i Fotogrametrii, Akademia Rolnicza we Wroc³awiu

S³owa kluczowe: skaning laserowy, filtracja, aktywne powierzchnie Keywords: laser scanning, filtering, flakes

Wstêp

Znaczn¹ czêœæ zbiorów danych lotniczego skaningu laserowego stanowi¹ odbicia od obiek-tów le¿¹cych ponad powierzchni¹ terenu, np. drzew, budynków. Punkty takie traktowane s¹ jako b³êdy grube i musz¹ byæ usuniête ze zbiorów danych, jeœli dane te maj¹ pos³u¿yæ do interpolacji numerycznego modelu powierzchni terenowej.

Pod pojêciem filtracji nale¿y rozumieæ pewn¹ formê (automatycznej) eliminacji punktów nie nale¿¹cych do modelowanej powierzchni. W ró¿nych oœrodkach naukowych opracowa-no kilka algorytmów filtracji. Ich zestawienie i ocenê mo¿liwoœci znaleŸæ mo¿na w raporcie (Sithole i Vosselman, 2003, 2004) i pracy (Borkowski, 2004). Do najbardziej skutecznych nale¿¹ algorytmy: bazuj¹ce na odpornej predykcji liniowej (Kraus, 2000; Kraus i Pfeifer, 2001, Briese et al., 2002), iteracyjnym przybli¿aniu pewnej, odpowiednio wybranej powierzch-ni startowej do danych pomiarowych (Axelsson, 1999, 2000), wykorzystapowierzch-niu operatorów matematycznej morfologii (kryterium spadków terenu) (Vosselman i Maas, 2001; Roggero, 2001; Sithole, 2001), minimalizacji energii powierzchni, zale¿nej od jej nachylenia (spadku) (Elmqvist et al., 2001; Elmqvist, 2002). W Polsce prace w zakresie filtracji danych skaningu laserowego prowadzone s¹ w Zak³adzie Fotogrametrii i Informatyki Teledetekcyjnej AGH, a bazuj¹ one na filtracji cyfrowej z wykorzystaniem FFT (Marmol, 2000, Marmol i Jachimski, 2004).

Wszystkie istniej¹ce metody filtracji maj¹ swoje zalety ale i ograniczenia; zawodz¹ w terenach o skomplikowanej strukturze, na przyk³ad w terenach zurbanizowanych zawieraj¹-cych du¿e kompleksy budowlane.

Model aktywnej powierzchni, bêd¹cy rozwi¹zaniem zadania minimalizacyjnego energii, zosta³ wprowadzony w pracy (Borkowski, 2004) jako uogólnienie modelu Snake (Kass et. al., 1987). Przez odpowiedni¹ definicjê energii zewnêtrznej mo¿na model ten uodporniæ na b³êdy grube aproksymowanych danych. W niniejszej pracy przedstawiono ideê wykorzysta-nia metody do filtracji danych skaningu laserowego.

Model aktywnych powierzchni

WyobraŸmy sobie dostatecznie ma³e, g³adkie elementy powierzchni, którym przyporz¹d-kowana jest pewna energia, która aktywnie dopasowuje powierzchniê do danych pomiaro-wych i jednoczeœnie pozwala zachowaæ po¿¹dane w³aœciwoœci geometryczne tej powierzch-ni. Dla takiego modelu zaproponowano okreœlenie flakes (Borkowski et. al., 1997; Borkow-ski, 2004).

Energia ca³kowita E_tot przypisana poszczególnym elementom powierzchni jest sum¹ energii wewnêtrznej i zewnêtrznej:

E_tot = E_int + E_ext (1)

Energia wewnêtrzna E_int opisuje w³aœciwoœci geometryczne modelowanej powierzchni i jest wa¿on¹ sum¹ jej nachylenia (membrane kernel), reprezentowanego przez kwadrat nor-my gradientu z_x + z_y i krzywizny (thin plate kernel), reprezentowanej przez kwadrat hesjanu

2 2 2

2

xx xy yy

z

+

z

+

z

_,

(2) gdzie z_x := ∂z/∂x, z_xx := ∂ 2_/∂x2_{, itd. α i β s¹ swobodnymi parametrami wagowymi}

(steruj¹cy-mi) dobieranymi zale¿nie od zastosowania i po¿¹danych w³aœciwoœci geometrycznych (g³ad-koœci) modelowanej powierzchni. Energia zewnêtrzna E_ext jest zale¿na od danych i jest defi-niowana zale¿nie od zastosowania – najczêœciej opisuje rozbie¿noœæ pomiêdzy modelowan¹ powierzchni¹ a danymi pomiarowymi. Podany w dalszej czêœci pracy model energii ze-wnêtrznej pozwala uodporniæ flakes na b³êdy grube danych pomiarowych.

Prezentowany model powierzchni aproksymuje dane pomiarowe w sposób optymalny. Oznacza to minimalizacjê energii ca³kowitej wszystkich elementów powierzchni na obszarze B⊂R2 _{obejmuj¹cym dane pomiarowe,}

(3) Warunkiem koniecznym do tego, aby funkcja z = z(x, y) by³a funkcj¹ ekstremaln¹ funk-cjona³u (3) jest spe³nienie przez tê funkcjê równania ró¿niczkowego czwartego rzêdu, tak zwanego równania Eulera (Brondstein i Semendjajew, 1962)

(4)

(

2 2

) (

2 2 2

)

int 2 , 2 x y 2 xx xy yy E =

α

z +z +

β

z + z +z

∫∫

→ = B tot x y xx xy yy dxdy z z z x z z E y x z I( ( , )) ( ; , ; , , ) min

Po obliczeniu, z wykorzystaniem podanych definicji energii, wystêpuj¹cych w (4) po-chodnych i niewielkich przekszta³ceniach otrzymuje siê równanie ró¿niczkowe czwartego rzêdu

(5) w którym parametry wagowe α i β s¹ sta³e dla danej powierzchni. Równanie to rozwi¹zywa-ne jest numerycznie.

Realizacja numeryczna

Równanie ró¿niczkowe (5) dyskretyzo-wane jest tak zwan¹ metod¹ ró¿nic skoñ-czonych, to znaczy wystêpuj¹ce pochod-ne drugiego i czwartego rzêdu aproksymo-wane s¹ przez odpowiednie skoñczone wyra¿enia. Wyra¿enia te mog¹ byæ prak-tycznie podane dla danych o okreœlonej re-gularnej strukturze. W naszym przypadku ograniczymy siê do siatki kwadratów.

Przyjmuj¹c oznaczenia poszczególnych wartoœci siatki kwadratów zgodnie z ry-sunkiem 1 otrzymamy wyra¿enia aproksy-macyjne dla poszczególnych pochodnych:

Po wstawieniu do równania (5), za³o¿eniu, ¿e ∆_x=∆_y = 1 i przekszta³ceniu otrzymujemy równanie liniowe postaci:

(6)

(

+

) (

+ +2 +

)

=0 − ∂ ∂ yyyy xyxy xxxx yy xx ext _{z z z z z} z E β α

Rys. 1. Oznaczenie danych dla siatki kwadratów

2, 2 , 2 2, 2 , 1 2, 1, , 1, 2, , 1 2, 2 , 2 2, 2 0 ext i j i j i j i j i j i j i j i j i j i j i j i j i j E _{dz cz dz bz cz bz} z az bz cz bz dz cz dz - + + + + + - -+ + - - - - + -¶ _{- + + + + + +} ¶ + + + + + + + + = , 2 1, , 1, 1 | ( 2 ) xx i j i j i j i j x z = _D

_∆

z+ - z + z -, 2 , 1 , , 1 1 | ( 2 ) yy i j i j i j i j y z = _D∆ z + - z + z -, 4 2, 1, , 1, 2, 1 | ( 4 6 4 ) xxxx i j i j i j i j i j i j x z = _∆D z+ - z + + z - z- + z -, 4 , 2 , 1 , , 1 , 2 1 | ( 4 6 4 ) yyyy i j i j i j i j i j i j y z = _D∆ z + - z + + z - z - + z -, 2 2 2, 2 , 2 2, 2 2, , 2, 2, 2 , 2 2, 2 1 | ( 2 2 4 2 2 ) 16 xyxy i j i j i j i j i j i j i j i j i j i j x y z = _{D D∆ ∆} z- + - z + + z+ + - z- + z - z+ + z- - - z - + z+

-∂

∂

–

gdzie

(7)

Równanie (6) nale¿y zestawiæ dla ka¿dego punktu danych. Mamy zatem uk³ad równañ liniowych

Az = b (8)

Jeœli przyjmiemy, ¿e modelowane dane podane s¹ na siatce kwadratów z_ij, i = 1,2,..., n ; j = 1,2,..., m to,

Wektor b zawiera pochodne energii zewnêtrznej w poszczególnych punktach danych:

gdzie E_z := ∂E_ext/∂_z. Macierz A jest macierz¹ pasmow¹ o wymiarze (n · m) × (n · m) zesta-wian¹ na podstawie wspó³czynników (7) zgodnie z równaniem (6). Mo¿na siê ³atwo przeko-naæ, ¿e jest ona macierz¹ osobliw¹. Uk³ad równañ (8) rozwi¹zywany jest iteracyjnie

z_t = (A + I)-1 _(z

t-1 + Ez, t-1) (9)

W t-tym kroku iteracji obliczane s¹ wartoœci z_t na podstawie energii zewnêtrznej z po-przedniego kroku iteracji. Pochodn¹ energii zewnêtrznej obliczamy z zale¿noœci

(10) gdzie r jest wartoœci¹ odchy³ki pomiêdzy wysokoœci¹ pomierzon¹ z_d a aproksymowan¹ z_t w t-tym kroku iteracji, r := z_d – z_t. Graficzn¹ reprezentacjê tej funkcji pokazano na rysunku 2. Energia zewnêtrzna jest proporcjonalna do odleg³oœci pomiêdzy chwilowym po³o¿eniem funkcji aproksymuj¹cej a danymi pomiarowymi i powoduje „przyci¹ganie” funkcji do danych po-miarowych. Dla dodatnich r wartoœæ ta jest t³umiona za pomoc¹ funkcji Gaussa, i dziêki temu energia zewnêtrzna oddzia³uje tylko na niewielkich odleg³oœciach tworz¹c barierê dla b³êdów grubych. Wielkoœæ tej bariery ustalana jest za pomoc¹ parametru σ. Maksimum funkcji energii zewnêtrznej znajduje siê w punkcie r_m = σ/√2. Negatywne wartoœci r nie s¹ t³umione – si³a przyci¹gania generowana przez dane jest tym wiêksza im wy¿ej ponad punk-tami pomiarowymi znajduje siê funkcja aproksymuj¹ca. Taka definicja energii zewnêtrznej zgodna jest ze struktur¹ danych skaningu laserowego, które charakteryzuj¹ siê niesyme-trycznym rozk³adem b³êdów grubych. Definicja ta mo¿e byæ wykorzystana równie¿ dla innych celów, np. aby aproksymowaæ górn¹ powierzchniê lasu wystarczy zamieniæ strona-mi poszczególne sk³adniki definicji (10).

– 25 : 4 , : ( 4 ), 2 3 1 : , : 4 8 a b c d a b a b b b = + = - + = =

Rys. 2. Energia zewnêtrzna

Przyk³ad

Prezentowan¹ metodê zastosujemy teraz do filtracji rzeczywistych danych skaningu lase-rowego (rys. 3). Dane pozyskano z OEEPE Working Group on Laser Data Acquisition; (http://www.geomatics.kth.se/~fotogram/OEEPE/oeepe_laser_main.htm) przedstawiaj¹ frag-ment zabudowy miejskiej. Nieregularnie w p³aszczyŸnie rozproszone dane skaningu lasero-wego poddano interpolacji metod¹ najbli¿szego s¹siada na regularn¹ siatê kwadratów o boku równym œredniej odleg³oœci poziomej pomiêdzy najbli¿szymi punktami w zbiorze danych pocz¹tkowych. Metoda interpolacji polega na przypisaniu do ka¿dego wêz³a wysokoœci naj-bli¿ej po³o¿onego punktu pomiarowego.

Regularyzowane dane poddano nastêpnie modelowaniu za pomoc¹ prezentowanego algo-rytmu. Jako powierzchniê startow¹ wybrano p³aszczyznê le¿¹c¹ na wysokoœci najwy¿szego punktu w zbiorze. Startow¹ powierzchnie poddano nastêpnie iteracyjnemu „dopasowaniu” do danych pomiarowych zgodnie z zale¿noœci¹ (9), przy czym energia zewnêtrzna obliczana by³a w kolejnym kroku iteracji zgodnie z modelem (10) w ka¿dym punkcie danych. Przyjêto nastê-puj¹ce parametry steruj¹ce: α = 1, β = 1 i σ2 _{= 0,1. Proces iteracyjny koñczy siê jeœli z}

i≈ zt-1.

W wyniku modelowania otrzymuje siê g³adk¹ reprezentacjê powierzchni terenu (rys. 4), woln¹ od obiektów znajduj¹cych siê ponad ni¹. Na obrze¿ach modelowanego obszaru wi-doczne s¹ pewne drobne defekty. Zjawisko to wynika z tego, ¿e do rozwi¹zania równania ró¿niczkowego potrzebne s¹ warunki brzegowe. Mo¿na je zredukowaæ poprzez do³¹czenie dodatkowych danych na obrze¿ach opracowywanego obszaru.

W ostatnim etapie eliminowane s¹ punkty zakwalifikowane jako obarczone b³êdami gru-bymi. W tym celu obliczana jest odleg³oœæ punktu pomiarowego od modelowanej powierzchni. Jeœli wielkoœæ ta przekracza dok³adnoœæ wysokoœciow¹ skaningu (lub jej wielokrotnoœæ) punkt jest eliminowany. Punkty zakwalifikowane jako terenowe pokazano na rysunku 5. S¹ to punkty le¿¹ce w odleg³oœci do 0,1 m do powierzchni modelu. Taki sam zbiór punktów otrzymano w wyniku rêcznej klasyfikacji punktów.

Rys. 3. Dane testowe

Rys. 5. Punkty zakwalifikowane jako punkty terenowe

Rys. 4. Powierzchnia terenu modelowana za pomoc¹

Zakoñczenie

Przedstawion¹ w niniejszej pracy metodê aproksymacji otrzymuje siê w wyniku rozwi¹-zania zadania wariacyjnego, w którym minimalizowane s¹ okreœlone w³aœciwoœci geome-tryczne powierzchni wyra¿one za pomoc¹ gradientu (nachylenia) i krzywizny powierzchni oraz energia zewnêtrzna opisuj¹ca rozbie¿noœæ pomiêdzy danymi pomiarowymi a modelem. Zaproponowana, asymetryczna funkcja energii zewnêtrznej pozwala uodporniæ metodê na b³êdy grube skaningu laserowego.

Przeprowadzone testy numeryczne pokaza³y wysok¹ skutecznoœæ prezentowanej meto-dy, w szczególnoœci w terenach zurbanizowanych zawieraj¹cych du¿e kompleksy budowla-ne. Metoda sterowana jest za pomoc¹ trzech swobodnych parametrów. Poprzez parametr σ

definiowany jest zakres b³êdów grubych podlegaj¹cych filtracji. Parametry α i β okreœlaj¹ rozci¹gliwoœæ i sztywnoœæ modelowanej powierzchni. Im wiêksze wartoœci tych parame-trów tym wiêksze naprê¿enie aktywnej powierzchni i wiêksze jej w³aœciwoœci filtracyjne. Dziêki temu w procesie iteracyjnym eliminowane s¹ zarówno b³êdy grube jak i filtrowane b³êdy pomiaru. W zbiorach danych skaningu laserowego zdarzaj¹ siê równie¿ odbicia le¿¹ce wyraŸnie poni¿ej powierzchni terenu. Takie, b³êdne punkty pomiarowe powstaj¹ wskutek zjawiska wielotorowoœci. Tego typu b³êdy filtrowane s¹ równie¿ dziêki odpowiedniej sztyw-noœci powierzchni.

Wysokoœci modelowanej powierzchni otrzymuje siê z rozwi¹zania uk³adu równañ linio-wych, w którym liczba równañ jest równa licznie wêz³ów regularnej siatki (n × m). Wynikaj¹ st¹d trudnoœci numeryczne. W przypadku danych skaningu laserowego, gdzie zbiory s¹ rzêdu 106 _{i wiêcej punktów, osi¹ga siê szybko kres mo¿liwoœci, nawet wydajnych}

kompute-rów. Rozwi¹zaniem tutaj jest odpowiednia segmentacja danych i stosowanie aproksymacji lokalnie, z odpowiednio dobranym obszarem wspólnym oraz zastosowanie efektywnych numerycznie algorytmów.

Oddzielnym zagadnieniem jest problem weryfikacji poprawnoœci filtracji wykonanej pre-zentowan¹ metod¹. Przeprowadzona filtracja rêczna, prezentowanego w pracy przyk³adu testowego, potwierdzi³a poprawnoœæ filtracji wykonanej metod¹ aktywnych powierzchni. W dalszej kolejnoœci wykonana bêdzie analiza porównawcza prezentowanej i innych istniej¹-cych metod filtracji. Wykorzystane zostan¹ w tym celu dane lotniczego skaningu laserowe-go, wykonanego w ostatnim czasie na obszarze Wroc³awia.

Literatura

Axelsson P., 1999: Processing of laser scanner data – algorithms and applications. ISPRS Journal of Photo-grammetry and Remote Sensing, 54(2), 138-147.

Axelsson P., 2000: DEM generation from laser scanner data using adaptive TIN models. International Borkowkski A., Burkhardt D., Meier S., 1997: Zur optimalen Approximation von Höhenprofilen. Österr.

Zeitschrift für Vermessung und Geoinformation, 182-285.

Borkowski A., 2004: Modellierung von Oberflächen mit Diskontinuitäten. Deutsche Geodätische Kommis-sion, Reihe C, Heft Nr 575.

Briese C., Pfeifer N., Dorninger P., 2002: Applications of the robust interpolation for DTM determination. Symposium ISPRS Commision III, Photogrammetric Computer Vision, Graz,. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIV / 3A, 55-61.

Bronstein I., Semendjajew K., 1962: Taschenbuch der Mathematik. B.G. Teubner, Leipzig.

Elmqvist M., Jungert E., Persson A., Soderman U., 2001. Terrain modelling and analysis using laser scanner data. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIV-3/W4, Annopolis, Maryland, 219-227.

Elmqvist M., 2002: Ground surface estimation from airborne laser scanner data using active shape models. ISPRS, Commission III, Symposium Photogrammetric Computer Vision, Graz, 114-118.

Kass M., Witkin A., Terzopoulos D., 1987: Snakes: Active contour models. Proceedings of the First Interna-tional Conference on Computer Vision, IEEE Comput. Soc. Press, 259-268.

Kraus K., 2000: Photogrammetrie. Band 3. Topographische Informatonssysteme. Dümmler, Köln. Kraus K., Pfeifer N., 2001: Advanced DTM generatin from LIDAR data. International Archives of

Photo-grammetry and Remote Sensing, Vol. XXXIV-3/W4, Annopolis, Maryland, 23-30.

Marmol U., 2003: Pozyskiwanie numerycznego modelu powierzchni topograficznej (NMPT) w oparciu o dane wysokoœciowe pochodz¹ce z lotniczego skaningu laserowego. Archiwum Fotogrametrii, Kartografii i Teledetekcji. Vol. 13B, s. 419-426.

Marmol U., Jachimski J., 2004: A FFT based method of filtering airborne laser scanner data. ISPRS Congress, Istambul, Commision3. http://www.isprs.org/commission3/wg3.

Roggero M., 2001: Airborne laser scanning: Clustering in row data. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIV-3/W4, Annopolis, Maryland, 227-232.

Sithole G., 2001: Filtering of laser altimetry data using a slope adaptive filter. International Archives of Photogrammetry and Remote Sensing, Vol. XXXIV-3/W4, Annopolis, Maryland, 203-210.

Sithole G., Vosselman G., 2003: Report ISPRS: Comparison of filters, http://www.isprs.org/commission3/ wg3.

Sithole G., Vosselman G., 2004: Experimental comparison of filter algorithms for bare-Earth extraction from airborne laser scanning point clouds. ISPRS Journal of Photogrammetry and Remote Sensing, Vo. 59, 85-101.

Vosselman G., Maas, H.G., 2001: Adjustment and filtering of raw laser altimetry data. OEEPE Workshop on Airborne Laserscanning and Interferometric SAR for Detailed Digital Elevation Models, Stockholm.

Summary

Laser scanner data sets contain a considerable number of gross errors which are reflections from objects that lie above the terrain surface, i.e. trees, houses. If Digital Terrain Model is interpolated such points need to be removed from data files. For automatic filtering of gross errors an algorithm based on active surface model (flakes) was proposed. This model constitutes a variational problem solution in which total energy comprising external and internal energy is minimized. Internal energy describes geometric properties of modelled surface and, it is the weighted sum of its inclination and curvature. External energy measures the deviation of the model from the given data. The Euler equation that is equivalent to variational problem is discretized with finite difference method and solved iteratively.

A model for external energy that enables to make flakes robust to gross errors was formulated. Filtering features of the method were presented on the laser scanner data. Hints for practical applica-tion of presented method were given.

dr hab. in¿. Andrzej Borkowski borkowski@ar.wroc.pl