STATYSTYCZNE I GEOSTATYCZNE METODY BADANIA DEFORMACJI POWIERZCHNI

Marcin GIL¹, Paweł FRĄCKIEWICZ², Katarina PUKANSKÁ³, Peter BLIŠTÁN⁴

STATYSTYCZNE I GEOSTATYCZNE METODY BADANIA DEFORMACJI POWIERZCHNI

1. WPROWADZENIE

Obecne technologie pomiarowe są efektywne zarówno w przypadku, gdy ob- serwacje obejmują tylko subiektywnie wybrane punkty, jak i chmury pikiet, gdy pozyskiwane są tysiące, a nawet miliony danych. Przy pomiarach, których celem jest budowa modelu powierzchni, znaczenie ma zarówno dokładność wyznaczo- nych współrzędnych, jak i liczba oraz rozmieszczenie pikiet. Modele generowane z danych pozyskanych innymi metodami różnią się od siebie. Różnice są znaczące, gdy porównywane modele powierzchni są generowane na podstawie danych pozy- skanych metodami klasycznymi z efektami analizy chmury punktów.

Na etapie wizualizacji powierzchni dane punktowe zamieniane są na pola ciągłe w postaci rastrowej. Jakkolwiek metody interpolacji automatycznej różnią się w szczegółach, to zasada ich działania jest podobna. W każdym przypadku o jako- ści modelu decyduje rozmieszczenie punktów pomiarowych i ich wpływ na oto- czenie. Punkty tworzą siatki regularnych węzłów, ciągi punktów w liniach pomia- rowych lub nieregularne rozkłady punktów pomiarowych. Wartość zmierzona w danym punkcie powierzchni jest informacją opisującą konkretnie miejsce, a ma- tematyczna funkcja, zastosowana przy budowie modelu tylko aproksymuje otocze- nie tego punktu.

Wybór modelu interpolacji ma znaczenie praktyczne. Dokładność i wiarygod- ność analiz przestrzennych wykonywanych na wygenerowanych modelach zależy od tego, jak dokładnie model reprezentuje rzeczywistą powierzchnię. Znalezienie najlepszego modelu jest problemem badawczym.

1 mgr inż. Marcin Gil – Politechnika Świętokrzyska, aleja Tysiąclecia Państwa Polskiego 7, 25-314 Kielce, PL.

2 mgr inż. Paweł Frąckiewicz – Politechnika Świętokrzyska, aleja Tysiąclecia Państwa Polskiego 7, 25-314 Kielce, PL.

3 dr hab. inż. Katarina Pukanská, prof. PŚk – Politechnika Świętokrzyska, aleja Tysiąclecia Państwa Polskiego 7, 25-314 Kielce, PL.

4 dr hab. inż. Peter Blištán, prof. PŚk – Politechnika Świętokrzyska, aleja Tysiąclecia Państwa Pol- skiego 7, 25-314 Kielce, PL.

1. CHARAKTERYSTYKA BADANEGO TERENU

Jako obiekt badań wybrano teren położony w miejscowości Mozgawa w gminie Pińczów w województwie świętokrzyskim. Lokalizację obszaru badań zaprezen- towano na mapie topograficznej (rys. 1).

Rys. 1. Lokalizacja obszaru badań

Morfologicznie jest to niecka ablacyjna powstała wskutek spłukiwania materia- łu w kierunku strumienia znajdującego się w kierunku południowym. Zarówno szczyt, jak i dno niecki ze względu na swój kształt (teren płaski o niewielkim, regu- larnym spadku) wykorzystywane są na cele rolnicze. Zbocza niecki zostały zabez- pieczone roślinnością trawiastą oraz niskim zakrzewieniem w celu zapobiegnięcia możliwości dalszego spływu materiału. Południowo-zachodni fragment zbocza został dodatkowo obsadzony drzewami iglastymi (rys. 1).

3. METODY OBSERWACJI POWIERZCHNI TOPOGRAFICZNEJ

Pomiary powierzchni topograficznej można przeprowadzać różnymi metodami, tj. klasycznymi, obserwacjami satelitarnymi GNSS, skanowaniem lub metodami fotogrametrycznymi. Wynikiem ostatnich dwóch są chmury punktów, które przy odpowiednio dobranej rozdzielczości dają efekt bliski rzeczywistemu modelowi mierzonej powierzchni. Jakkolwiek klasyczne metody pomiaru tej zalety nie mają, to ich zastosowanie nie wyklucza stworzenia modelu o wysokiej dokładności. W takich przypadkach duże znaczenie ma dobór odpowiednich algorytmów interpolacyjnych.

Obserwacje można prowadzić metodami punktowymi, obszarowymi lub meto- dami mieszanymi. Metody punktowe polegają na pomiarze wcześniej zastablizo- wanych punktów w charakterystycznych przekrojach lub punktach przewidywa- nych zmian, a sama liczba punktów jest niewielka. Metody obszarowe polegają na pomiarze regularnych lub nieregularnych siatek punktów zlokalizowanych na całej powierzchni, a subiektywna ocena doboru lokalizacji punktów wpływa na jakość wygenerowanego modelu. Metody mieszane wykorzystują elementy obu metod.

Połączenie umożliwia kontrolowanie obszaru we wcześniej sygnalizowanych punktach oraz pozwala dobrać pozostałe w taki sposób, aby powstały model był jak najwierniejszym odzwierciedleniem rzeczywistej powierzchni.

Wybrany obszar pomierzono dwukrotnie dwiema metodami:

1. Pomiar tachimetryczny wykonano w listopadzie 2017 roku za pomocą ta- chimetru Stonex R1 Plus (rys. 2a), mierzącego kąty z dwusekundową do- kładnością (tabela 1). Zbiór obserwacji tworzy 435 punktów.

2. Pomiar RTN GNSS (jeden zbiór obserwacji) wykonano w listopadzie 2018 roku za pomocą nowoczesnego odbiornika Sokkia GRX1 (rys. 2b) o parame- trach przedstawionych w tabeli 2. Zbiór obserwacji składa się z 857 punktów.

Rys. 2. Tachimetr Stonex R1 Plus (a), odbiornik Sokkia GRX1 (b) a)

b)

Tabela 1. Parametry techniczne tachimetru Stonex R1 Plus

Pomiar kąta Dokładność wyświetlania 1” / 0,0002^gon / 0,005 mil

Dokładność 2”

Luneta Powiększenie lunety 30×

Pole widzenia lunety 1°30’

Pomiar odległości

Bezlustrowy 1,0 ~500 m

Na pryzmat 1,0 ~3000² 2000³ 1500⁴ m klasa 1 ponad 5000 m klasa 3

Na folię (60 × 60 mm) 1,0 ~ 800 m

Na mini pryzmat 1,0 ~ 800 m

Dokładność pomiaru odległości

Bezlustrowy ±3 mm + 2 ppm

Na pryzmat ±2 mm + 2 ppm

Na folię (60 × 60 mm) ±3 mm + 2 ppm

Tabela 2. Parametry techniczne odbiornika Sokkia GRX1

Liczba kanałów 72

Śledzone sygnały GPS L1, L2, L2C; GLONASS L1, L2

Częstotliwość 20 Hz (stand. 5Hz)

Antena Zintegrowana

Dokładności RTK H : +/–10 mm + 1 ppm; V : +/–20 mm + 1 ppm Dokładność STATIC H : +/–3 mm + 1 ppm; V : +/–10 mm + 1 ppm

Modem GSM/GPRS Wbudowany

Modem radiowy UHF Wbudowany

Porty Bluetooth × 2 , RS-232

Pamięć SD/SDHC

Format danych TPS, NMEA, RTMC, CMR

Kontroler SHC250

Zasilanie BDC58 Li-ion 4300 mAh/7.2 VDC

Czas pracy 7,5 h

Temperatura pracy –40°C + 65°C

Wodoszczelność IP 67

Na rysunku 3 przedstawiono widok na badany obszar z jego północno- wschodniego zbocza podczas wykonywania pomiaru w listopadzie w 2017 roku (rys. 3a) oraz w listopadzie 2018 roku (rys. 3b).

a) b)

Rys. 3. Obiekt badań: a) listopad 2017 rok, b) listopad 2018 rok

4. METODYKA BADAŃ

W przypadku pomiarów powierzchniowych obserwacje poszczególnych punk- tów dają tylko przybliżenie rzeczywistego kształtu obserwowanej powierzchni.

Jakość modelu z nich wygenerowanego zależy od wyboru charakterystycznych punktów oraz liczebności zbiorów danych. Obszarowe badania powierzchni wy- magają dużej liczby obserwacji, ale z uwagi na duże objętości danych nie wymaga- ją dużej dokładności.

Analizę metod interpolacji zrealizowano przy użyciu programu komputerowego Surfer. Program oferuje 12 metod interpolacji opartych na modelu GRID. Metody bazują na różnych algorytmach, skutkując różnym rozmieszczeniem prognozowa- nych wartości poza punktami pomiarowymi. Do analizy wykorzystano pięć najczę- ściej stosowanych metod interpolacji, takich jak [2, 3, 5, 9–13]:

1. Metoda odwrotnej odległości IDP (Inverse Distance to a Power), która pole- ga na obliczeniu średniej ważonej obserwacji w otoczeniu. Wagi obserwacji są odwrotnie proporcjonalne do odległości pomiędzy punktami pomierzo- nymi a punktem interpolowanym.

2. Zmodyfikowana metoda Sheparda MSM (Modified Shepard’s Method), która jest uogólnieniem metody odwrotnych odległości. Algorytm wykorzystuje dwa rodzaje funkcji interpolacyjnych: wierne, w których parametr funkcji jest zgodny z pomierzonym parametrem, oraz wygładzające, w których wartość wejściowa nie jest precyzyjnie zlokalizowana na generowanej powierzchni.

3. Metoda minimalnej krzywizny MC (Minimum Curvature), która tworzy ela- styczną powierzchnię z minimalną liczbą wygięć. Algorytm generuje możli- wie jak najgładszą powierzchnię, tworząc siatkę, poprzez wielokrotne zasto- sowanie równania algorytmu na punkcie pomiarowym w celu wygładzenia.

4. Metoda radialnych funkcji bazowych RBF (Radial Basis Function), w której wartości rzeczywiste zależą od odległości od punktu początkowego. Metoda pozwala na wygenerowanie gładkiej powierzchni z dużych zbiorów danych.

5. Metoda krigningu (K), czyli metoda interpolacji oparta na geostatystyce, w której wyznacza się błąd interpolacji zwany wariancją krigingu. Algorytm krigingu jest skuteczny, ponieważ może zrekompensować dane w zbiorze, nadając tym obszarom mniejszą wagę w ogólnej prognozie. Pozwala rów- nież na ekstrapolację poza obszarem danych.

Błąd interpolacji v (1) powstaje w wyniku zastąpienia wyników obserwacji siat- ką interpolacji. Jest to różnica pomiędzy wartością wyinterpolowaną a wielkością fizycznie pomierzoną.

v = Zinterpolacja – Z_pomiaru (1) Dobrą miarą dokładności interpolacji jest suma kwadratowych różnic v. Przy analizie różnych wariantów metod interpolacji czy danych pomierzonych suma ta stanowi obiektywną miarę dopasowania wygenerowanej powierzchni względem pomierzonych punktów (2). Ograniczeniem metody jest brak odporności na punkty odstające w zbiorze danych. Dlatego ocena wielkości błędów interpolacji musi zostać poprzedzona oceną jakościową (wizualną).

2 1

min

n i







 (2) Wielkość błędów wyznaczano przy użyciu funkcji Grid/Residuals w progra- mie Surfer [10]. Celem wykonanych analiz było określenie wiarygodności uzy- skanych modeli powierzchni terenu przy zastosowaniu metod statystycznych i geostatystycznych. Następnie porównano modele przy różnej liczbie punktów obserwacji wynikających z zastosowania różnych metod pomiarowych przy za- stosowaniu tej samej geometrii siatki GRID oraz porównano błędy interpolacji modeli [1, 4, 6–8].

Analizy sumy kwadratów błędów różnic przeprowadzono dla modeli interpola- cyjnych przy rozdzielczości siatki GRID dobranej przez program Surfer w oparciu o algorytm metody i najmniejszą odległość pomiędzy parą punktów. Dla pomiaru tachimetrycznego o mniejszym zagęszczeniu oraz pomiaru GNSS o większym zagęszczeniu program dobrał siatkę interpolacyjną 100 rzędów na 33 kolumny, co daje siatkę o oczku ~2,5 × 2,5 m. Analizy pozwalają sprawdzić wpływ stopnia zagęszczenia na jakość modelu i wielkość błędu interpolacji.

Testowane modele algorytmu odwrotnej odległości generowano w oparciu o za- łożenie współczynnika anizotropii (anisotropy ratio) równego 1 i kąta anizotropii (anisotropy angle) równego 0. W metodzie Sheparda wybrano zasięg południkowy (Range 1) i równoleżnikowy (Range 2) równy w obu kierunkach celem osiągnięcia współczynnika anizotropii równego 1. Dla algorytmu minimalnej krzywizny zasto-

sowano współczynnik anizotropii równy 1, maksymalny błąd interpolacyjny (ma- ximum residual) 0,022 oraz dla maksymalnej liczby iteracji (maximum iteration) 100 000. W algorytmie radialnych funkcji bazowych wybrano współczynnik anizo- tropii równy 1, kąt anizotropii 0 oraz funkcję wielokwadratową (multiquadric).

Modele według algorytmu krigingu analizowano w oparciu o współczynnik anizo- tropii równy 1oraz kąt anizotropii 0. Wybrano liniowy model wariogramu oraz nie analizowano wpływu trendu (typedrift) [10].

Pierwszym etapem kontroli i analizy poprawności modelu była ocena wyglądu mapy. Przykładowe modele interpolacyjne na podstawie pomiaru tachimetryczne- go i pomiaru GNNS zaprezentowano na rysunku 4.

Rys. 4. Model powierzchni z pomiaru tachimetrycznego według algorytmu krigingu (a), model powierzchni z pomiaru GNSS według algorytmu krigingu (b)

W celu oceny jakościowej porównano mapy warstwicowe dla badanych algo- rytmów. Na rysunku 5 przedstawiono mapy z pomiaru tachimetrycznego, a na rysunku 6 zwizualizowano mapy z pomiaru GNSS. Wielkości błędów modelu za- prezentowano pod mapami warstwicowymi

Rys. 5. Mapy warstwicowe z tachimetrii według: a) algorytmu IDP, b) algorytmu MSM, c) algorytmu MC, d) algorytmu RBF, e) algorytmu K

Rys. 6. Mapy warstwicowe z GNSS według: a) algorytmu IDP, b) algorytmu MSM, c) algorytmu MC, d) algorytmu RBF, e) algorytmu K

5. ANALIZA ALGORYTMÓW INTERPOLACJI

Analizę jakościową zastosowanych algorytmów interpolacji przeprowadzono pod kątem czytelności i ogólnej zgodności z rzeczywistym ukształtowaniem tere- nu. Jej efektem są następujące spostrzeżenia:

 algorytm odwrotnych odległości (IDP) spowodował nagromadzenie izolinii wokół pomierzonych pikiet, szczególnie na dnie niecki. Powierzchnia jest niezgodna z rzeczywistym ukształtowaniem. Przebieg warstwic jest zabu- rzony ze względu na algorytm, w którym pojedynczy punkt o lokalnym eks- tremum wpływa na wynik na punkty wokół, a nie tylko na najbliższy punkt;

 algorytm zmodyfikowanej metody Sheparda (MSM) powoduje na brzegach nierzeczywisty przebieg formy terenu. Największe niedoskonałości modelu algorytm wygenerował na górnych krawędziach niecki. Algorytm łączy eks- tremalne punkty na brzegach, nie nawiązując do ogólnego trendu zmienności danych;

 algorytm najmniejszej krzywizny (MC) wygenerował warstwice o reali- stycznym przebiegu z minimalną ilością wygięć przy maksymalnym jej wy- gładzeniu;

 algorytm radialnych funkcji bazowych (RBF) spowodował czytelność i zgodność mapy z rzeczywistym przebiegiem terenu. Zastosowana funkcja wielokwadratowa pozwoliła na uniknięcie osobliwości macierzy lub jej złe- go uwarunkowania w przypadku niewielkich odległości pomiędzy pikietami;

 algorytm krigingu dzięki nieobciążonemu liniowo oszacowaniu wygenero- wał warstwice zgodne z rzeczywistym ukształtowaniem terenu. Zastosowane przybliżenie powierzchni na podstawie metody polegającej na minimalizacji wariancji krigingu jest wiarygodne i czytelne.

Analiza graficzna czytelności i wiarygodności map warstwicowych skłania do uznania większej wiarygodności modeli opartych na metodzie najmniejszej krzy- wizny, radialnych funkcji bazowych oraz krigingu. Najtrafniej ukształtowanie tere- nu przedstawia model oparty na geostatytstycznej estymacji. Modele oparte na różnej liczbie punktów sugerują podobne oceny. Algorytmy RBF i K przy większej liczbie obserwacji dają zbliżone do siebie przebiegi warstwic.

Wykonano analizę ilościową sumy kwadratów błędów interpolacji w oparciu o ich minimalizację (tabela 3). Najmniejszą wartość Σv² uzyskano dla zmodyfiko- wanej metody Sheparda. Jest mniejsza od 2 m² przy siatce interpolacyjnej GRID

~2,5 × 2,5 m. Maksymalna różnica pomiędzy wartością pomierzoną a modelową nie przekracza 0,40 m. Drugim pod względem dokładności algorytmem jest RBF, dla którego Σv² z pomiaru tachimetrycznego jest większe od 2 m², a dla pomiaru GNSS większe od 5 m². Skrajna różnica pomiędzy modelem a pikietą nie przekra- cza 0,30 m. Algorytm krigingu uzyskał błędy na poziomie poniżej 10 m² z pomiaru tachimetrycznego oraz ponad 11,5 m² dla pomiaru GNSS. Różnice w modelu z tachimetrii nie przekraczają 46 cm, a z GNSS 40 cm. Metoda najmniejszej krzy- wizny daje błędy odpowiednio ponad 15 m² i 26 m². Maksymalne różnice nie prze-

kraczają wartości 0,90 m. Największe wartości Σv² przekraczające 100 m² uzyska- no dla metody odwrotnych odległości. Największe różnice pomiędzy wartością pomierzoną a wymodelowaną wynoszą ponad 1,6 m dla pomiaru tachimetrycznego oraz 1,15 m dla GNSS.

Tabela 3. Sumy kwadratów błędów z zastosowanych algorytmów interpolacji Algorytm Σv², m² z tachimetrii Σv², m² z GNSS

IDP 106,06 104,29

MSM 1,79 1,47

MC 15,23 26,20

RBF 2,25 5,61

K 9,93 11,61

Dla metod, których ocena graficzna była najsłabsza (IDP oraz MSM), nieznacz- nie zmniejszył się błąd interpolacji przy większej ilości punktów z metody GNSS.

Dla pozostałych badanych metod wartości błędów uległy zwiększeniu. Przy różnej liczbie punktów błąd interpolacji zwiększył się o 149% dla algorytmu RBF, o 72%

dla MC oraz 16,9% dla krigingu.

Analiza błędów interpolacji wskazuje na małą wiarygodność metody odwrotnej odległości. Pozostałe algorytmy dają błędy na poziomie akceptowalnym w prakty- ce. Najlepsze rezultaty uzyskano przy zastosowaniu metody radialnych funkcji bazowych oraz metody krigingu. Obie mogą być stosowane przy analizach po- wierzchni niecek oraz wyrobisk poeksploatacyjnych.

6. WNIOSKI

Nowoczesne metody pomiarowe dają możliwość szybkiego pomiaru gęstych chmur punktów, które odzwierciedlają rzeczywiste ukształtowanie terenu. Wyma- gają jednak kosztownego sprzętu i specjalistycznego oprogramowania oraz jedno- stek komputerowych o dużej mocy obliczeniowej. Post-processing nowoczesnych metod staje się żmudny i dodatkowo nie zawsze kończy się powodzeniem. Alterna- tywą są klasyczne metody pomiarowe, do których należą np. tachimetria czy GNSS. Czas trwania pomiaru klasycznymi metodami jest dłuższy, ale opracowanie wyników jest dużo łatwiejsze, wymaga jedynie zastosowania odpowiednich algo- rytmów obliczeń i modelowania danych.

Algorytmy interpolacji w różnym stopniu uwzględniają wpływ trendu ukształ- towania terenu. Największe błędy interpolacji związane są z punktami skrajnymi modelowanej powierzchni. Na krawędzi obszarów występują wygładzenia po- wierzchni najmniej zbliżone do jej rzeczywistego przebiegu.

Algorytm odwrotnych odległości daje nierealistyczny przebieg mapy warstwi- cowej. Pojedyncze, lokalne, ekstremalne punkty powodują zaburzenie trendu zmienności danych. Metoda efektywna dla powierzchni o małej zmienności danych sąsiadujących ze sobą punktów. Zmodyfikowana metoda Sheparda daje najmniej- sze błędy interpolacji, jednak algorytm nie uwzględnia ogólnego trendu otrzyma- nych danych. Przebieg uzyskanego modelu znacznie odbiega od rzeczywistości i na etapie oceny graficznej wymaga odrzucenia. Jeżeli dane mają duże zmiany wartości powierzchni na niewielkich odległościach, algorytmy IDP oraz MSM stają się nieefektywne.

Algorytm najmniejszej krzywizny daje błędy interpolacji do zaakceptowania w praktyce, odzwierciedla rzeczywiste ukształtowanie badanego obszaru przy jego maksymalnym wygładzeniu. Najmniejsze błędy interpolacji uzyskuje się, stosując algorytm radialnych funkcji bazowych oraz krigingu, które prezentują powierzch- nię najbardziej zbliżoną dla rzeczywistej.

Wykonane analizy potwierdziły efektywność metody radialnych funkcji bazo- wych oraz krigingu przy zastosowaniu nieregularnej siatki pomiarowej ze względu na niehomogeniczność zastosowanych algorytmów. W przypadku mniejszej liczby punktów bardziej efektywny jest algorytm krigingu, przy większej ilości danych modele dają bardzo zbliżone graficzne wyniki. Błąd interpolacji jest mniejszy dla metody radialnych funkcji bazowych.

Zastosowanie algorytmu RFB albo K w przypadku badania powierzchni terenu jest najbardziej wydajne. Wygenerowane powierzchnie przebiegają najbliżej danych pomiarowych przy realistycznym i wiarygodnym odwzorowaniu powierzchni. Ze względu na stosowanie nowoczesnych metod pomiarowych w postaci skanerów laserowych lub dronów należy prowadzić dodatkowe analizy w celu możliwości porównania pomiarów pochodzących z klasycznych metod z chmurami punktów.

LITERATURA

[1] Bartuś T., Słomka T., Geostatystyczna estymacja parametrów jakości węgla brunatnego w polu Bełchatów wykorzystująca znajomość zmienności lokalnej, „Gospodarka Surow- cami Mineralnymi” 2009, t. 25, z. 2, s. 5–22.

[2] Cressie N.A.C., Statistics for spatial data, revised edition, Wiley & Sons, New York 1993.

[3] Davis J.C., Statistics and Data Analysis in Geology, Wiley & Sons, New York 1986.

[4] Goldsztejn P., Skrzypek G., Wykorzystanie metod interpolacji do numerycznego kreśle- nia map powierzchni geologicznych na podstawie nieregularnie rozmieszczonych da- nych, „Przegląd Geologiczny” 2004, Vol. 52, nr 3, s. 233–236.

[5] Hengl T., Minasny B., Gould M., A geostatistical analysis of geostatistics,

“Scientometrics” 2009, Vol. 80, s. 491–514.

[6] Mucha J., Wasilewska M., Dokładność interpolacji zawartości siarki i popiołu w wy- branych pokładach węgla kamiennego GZW, „Gospodarka Surowcami Mineralnym”

2005, t. 21, s. 5–21.

[7] Mucha J., Metody geostatystyczne w dokumentowaniu złóż, Wydawnictwo AGH, Kra- ków 1994.

[8] Mitka B., Piech I., Modelowanie powierzchni terenu, „Infrastruktura i Ekologia Tere- nów Wiejskich” 2012, nr 3, PAN, Oddział w Krakowie, s. 167–180.

[9] Namysłowska-Wilczyńska B., Geostatystyka. Teoria i zastosowania, Oficyna Wy- dawnicza Politechniki Wrocławskiej, Wrocław 2006.

[10] Surfer – User’s Guide, Golden Software, 2002.

[11] Zawadzki J., Metody geostatystyczne dla kierunków przyrodniczych i technicznych, Oficyna Wydawnicza Politechniki Warszawskiej, Warszawa 2015.

[12] Zawadzki J., Cieszewski C.J., Zasada M., Lowe R.C., Applying geostatistics for investigations of forest ecosystems using remote sensing imagery, “Silva Fennica”, t. 39, s. 599–617.

[13] Zhou F., Guo H.-C., Ho Y.-S., Wu C.-Z., Scientometric analysis of geostatistics using multivariate methods, “Scientometrics” 2007, Vol. 73, s. 265–279.