Wprowadzenie do metodyki interpretacji badañ georadarowych
przy u¿yciu procedury modelowania numerycznego
Tomis³aw Go³êbiowski*
Introduction to GPR surveys interpretation with procedure of numerical modelling. Prz. Geol., 52: 563–568.
S u m m a r y . The paper presents results of numerical modelling of the electromagnetic wave field in the ground based on the Finite Dif-ference Time Domain Method (FDTD). Theoretical background of the electromagnetic wave field in the ground is presented. Tech-niques of building of the discrete model and the numerical simulation are described. Computer program REFLEX-W, used for modelling, is also briefly characterized. A sample measurement echogram for shallow georadar (GPR, ground penetrating radar) sur-veys is presented together with the synthetic one. Usefulness of the FDTD procedure in the interpretation of echograms is analysed and evaluated. The main goal of the paper is to introduce a reader into the procedure of numerical modelling, a modern method of the GPR surveys interpretation.
Key words: GPR measurements, numerical modelling
W artykule przedstawiono zastosowanie metody FDTD (Finite Difference Time Domain Method) — metoda ró¿nic skoñczonych w domenie czasu, do modelowania elektroma-gnetycznego pola falowego w gruncie. Na œwiecie metoda ta znajduje zastosowanie m.in. w procesie interpretacji wyników pomiarów georadarowych (Bergmann i in., 1998; Carcione, 1996; Holliger, 2002; Roberts & Daniels, 1997).
Metoda modelowania numerycznego dla celów geora-darowych nie by³a dotychczas stosowana w naszym kraju, wiêc dok³adny opis procedury modelowania zawarty w arty-kule pozwoli na wprowadzenie nowego narzêdzia do inter-pretacji pomiarów georadarowych. Na etapie interinter-pretacji wyników badañ geofizycznych wspó³pracuj¹ ze sob¹ m.in. geofizycy, geolodzy, hydrogeolodzy, z tego te¿ powodu g³ównym celem artyku³u jest wprowadzenie w sposób przy-stêpny czytelnika (nie zawsze geofizyka) do nowoczesnej metody wspomagaj¹cej interpretacje badañ geofizycznych jak¹ jest procedura modelowania numerycznego.
Odpowiednio du¿o uwagi i miejsca w artykule poœwiê-cono opisowi technologii budowy modelu numerycznego i sposobowi prowadzenia obliczeñ, co z jednej strony w spo-sób doœæ przystêpny pozwala zrozumieæ sam¹ technikê modelowania, a z drugiej pozwala na odpowiednie odniesie-nie wyników modelowania (które s¹ zawsze pewn¹ idealiza-cj¹ rzeczywistoœci) i ich korelacjê z wynikami pomiarów.
W artykule pominiêto dok³adny opis aparatury pomiaro-wej, podstaw fizyczno-matematycznych metody rowej oraz techniki prowadzenia pomiarów georada-rowych, a tak¿e sposobu ich przetwarzania, poniewa¿ zagad-nienia te wykraczaj¹ poza g³ówne za³o¿enie artyku³u. Bli¿sze informacje na wy¿ej wymieniony tematy mo¿na znaleŸæ m.in. w pracach: Karczewskiego (1997) i Annana (2001).
Przyk³adow¹ analizê numeryczn¹ przeprowadzono dla wybranego profilu georadarowego, który po przeprowa-dzeniu standardowego przetwarzania i interpretacji dawa³ wieloznaczne wyniki. Przy dotychczas stosowanych meto-dach interpretacyjnych, echogram ten zosta³ odrzucony, ze wzglêdu na zbyt du¿¹ niejednoznacznoœæ interpretacji. Artyku³ pokazuje, ¿e nawet dla tak „s³abych” wyników pomiarowych mo¿emy próbowaæ dokonywaæ ich
inter-pretacji z zastosowaniem procedury modelowania numerycznego.
Wybrany jako przyk³ad do modelowania profil pomiaro-wy jest pomiaro-wynikiem badañ miejskiej infrastruktury podziem-nej (rury, kable, kana³y itp.) we Wroc³awiu.
Podstawy teoretyczne numerycznego modelowania pola georadarowego
Matematyczny opis pól elektromagnetycznych jest oparty na funkcjach ci¹g³ych. W celu rozwi¹zania zadania w sposób numeryczny nale¿y przekszta³ciæ model ci¹g³y w model dyskretny o parametrach skupionych. Wprowadza-my nastêpuj¹c¹ notacje: rozmiar oczek siatki)x, )y, )z; krok czasowy)t; numeracja wêz³ów siatki i, j, k; rozmiar modelu i× )x, j × )y, k × )z.
Aby rozwi¹zaæ zagadnienie propagacji fal w modelu o skoñczonych rozmiarach i w skoñczonym czasie, dokonu-jemy dyskretyzacji przestrzennej z krokiem )x, )y, )z oraz dyskretyzacji czasowej z krokiem)t. W ten sposób otrzymujemy siatkê obliczeniow¹ (ryc. 1), w wêz³ach któ-rej obliczamy wartoœci sk³adowych pola elektromagne-tycznego.
Analiza numeryczna elektromagnetycznego (e.m.) pola falowego wymaga rozwi¹zania równañ Maxwella przy okreœlonych warunkach brzegowych. Równania te
y z x Ez (i, j, k) Ez Ez Ez Ex Ex Ex Ex Ey Ey Ey Ey Hy Hy Hx Hx H z Hz (i 1, j, k)+ (i, j 1, k)+ (i, j+1 -1, k ) (i+1, j+1 -1, k ) (i+1, j, k )-1
Ryc. 1. Oczko trójwymiarowej siatki obliczeniowej Fig. 1. One cell of 3-dimensional grid
*Wydzia³ Geologii, Geofizyki i Ochrony Œrodowiska, Akademia Górniczo-Hutnicza, al. Mickiewicza 30, 30-059 Kraków; goleb@2-0.pl
opisuj¹ rozk³ad sk³adowych pola falowego oraz ich zale-¿noœæ od rodzaju Ÿród³a. Oœrodek, w których przewodnoœæ F > 0 jest oœrodkiem stratnym (Morawski & Gwarek, 1998). W oœrodkach stratnych i anizotropowych, jakimi s¹ utwory geologiczne, równania Maxwella przybieraj¹ nastêpuj¹c¹ postaæ (Bergmann i in., 1998; Carcione, 1996):
Ñ x E = -¶ + ¶B Mt S Ñ x H = ¶ ¶D Jt + je¿eli: D = g*¶ ¶ E t B= :* ¶ ¶ H t J = F*¶ ¶E Jt + S to: Ñ x E = –:*¶ ¶ 2 2 H M t + S Ñ x H = F*¶ ¶ e ¶¶ 2 E E J t + * t + S 2 gdzie:
Ñ — operator nabla, E — sk³adowa elektryczna e.m. pola falowego,H —– sk³adowa magnetyczna e.m. pola falowego, D — indukcja elektryczna; B — indukcja magnetyczna,J — gêstoœæ pr¹du w oœrodku; JS— gêstoœæ pr¹du w Ÿródle,MS— gêstoœæ pr¹du magnetycznego,: — przenikalnoœæ magnetyczna oœrodka, g — sta³a dielek-tryczna oœrodka,F — przewodnoœæ elektryczna oœrodka, t — czas.
Obliczenia numeryczne e.m. pola falowego mog¹ byæ wykonywane albo w domenie czasu lub w domenie czêsto-tliwoœci. Najczêœciej obecnie stosowan¹ metod¹ modelo-wania pola e.m. jest metoda FDTD (metoda ró¿nic skoñczonych w domenie czasu), a algorytm numeryczny dla tej metody zosta³ opracowany przez Yee w 1966 r. (Yee, 1966). W metodzie tej ró¿niczkowanie po zmiennych prze-strzennych (x, y, z) i po czasie (t) w równaniach Maxwell’a aproksymowane jest ró¿nicami skoñczonymi zgodnie z ogólnym schematem (dla dowolnej funkcji F):
¶ ¶ F x dF dx F F x i dF dx F i j i j i j i j i j ® -1 2/ , » , – -1, +1 2/ , » +1, – D F x i j, D
( )
¶ ¶ 2 2 1 1 2 2 F x F F F x i j i j i j ® -, – , + + , DAnalogicznie dokonujemy zamiany ró¿niczek na ró¿-nice skoñczone dla pozosta³ych zmiennych tj. y, z i t. Zgod-nie z powy¿szym, mo¿emy obliczyæ odpowiednie sk³adowe pola e.m. w uk³adzie 3D tj. Ex, Ey, Ez, Hx, Hy, Hz, dla dyskretnych wartoœci po³o¿enia i czasu. Pole elek-tryczne jest obliczane w krokach czasowych n+1, nato-miast magnetyczne w n+1/2. Ze wzglêdu na obszernoœæ opisu poni¿ej podano przyk³ad wyliczenia sk³adowej Ex pola e.m.:
(
)
(
(
)
)
(
)
E i j k i j k t i j k i j k t x n x x x x + = × × + × × 1 0 5 0 5 , , , , – , , , , , , e s e s D D(
) (
)
(
)
(
)
i j k E i j k t i j k t i j k H i x n x x z n , , , , , , , , , , , × + + + × × × + D D e 0 5 s 0 5(
)
(
)
j k H i j k y z n + é ë ê ê + 1, – 0 5 , , – , D(
)
(
)
– , , – , , , , H i j k H i j k z y n y n + + + ù û ú ú 0 5 0 5 1 D gdzie: Ex n+1(i,j,k) — wartoœæ sk³adowej elektrycznej Ex e.m. pola falowego w wêŸle siatki obliczeniowej o wspó³rzêd-nych i,j,k w n+1 kroku czasowym.
Analogicznie postêpujemy dla wszystkich pozosta³ych sk³adowych wektorów E i H.
Podczas obliczeñ w domenie czasu wielkoœci)x, )y, )z i )t s¹ ze sob¹ powi¹zane. Rozmiar oczka siatki obli-czeniowej ()l) musi byæ min. 10-krotnie mniejszy ni¿ naj-mniejszej d³ugoœci fali (8) wystêpuj¹cej w sygnale Ÿród³owym. Najmniejsz¹ d³ugoœæ fali (a co za tym idzie najwiêksz¹ czêstotliwoœæ — fmax) sygna³u nale¿y odczytaæ z widma czêstotliwoœciowego tego sygna³u. Czêsto dla uproszczenia u¿ywa siê wzoru:
fmax= 3× fs 8 = c fmax× e Dl = l 10 gdzie:
fs— czêstotliwoœæ g³ówna w sygnale Ÿród³owym,
c — prêdkoœæ œwiat³a.
Aby zachowaæ zbie¿noœæ i stabilnoœæ rozwi¹zania numerycznego nale¿y uwzglêdniæ tzw. warunek Courant-Freidrich-Lewy (CFL):
( ) ( ) ( )
Dx 2+ Dy 2+ D 2 > ×c Dt=Dt× 1 × zm e
Na brzegu modelu obliczeniowego ustala siê warunki brzegowe. W przypadku modelowania pól falowych przyj-muje siê warunki typu ABC (Absorbing Boundary Condi-tion) — powoduj¹ one t³umienie fali padaj¹cej na brzeg modelu (Kosloff & Kosloff, 1986; Mur, 1981). Wzrost przewodnoœci brzegu modelu powoduje t³umienie fali elektromagnetycznej, zgodnie ze wzorem:
E(x,y,z,t) = E0×ej(w bt- ×k r)×e-ak r× a msn= 2 a
[
]
s e »1 7 3, e × dB m/ gdzie:E — wartoœæ amplitudy sk³adowej elektrycznej e.m. pola falowego w punkcie x, y, z modelu w chwili t; E0— wartoœæ amplitudy w Ÿródle; $ — wspó³czynnik fazy; " — wspó³czynnik t³umienia;k — wektor kierunku ruchu fali; r — wektor od Ÿród³a do punktu x, y, z; < — prêdkoœæ fali elektromagnetycznej w gruncie wyznaczana ze wzoru:
< = c
r
gdzie:
gr— wzglêdna przenikalnoœæ elektryczna oœrodka.
W programie REFLEX-W, który u¿yto do modelowa-nia, brzeg modelu jest okreœlony jako obszar wzrostu prze-wodnoœci od wartoœci zadanej w modelu do wartoœci maksymalnej Fmax. Wartoœæ Fmax oraz szerokoœæ obszaru
t³umi¹cego s¹ wyznaczane z poni¿szych zale¿noœci: Fmax= A× × ×w e e0 s1 =s +s s× max B – i lub Fi=F× De ln i B max × æèç öø÷ æ èç ö ø÷ s s gdzie:
B — wspó³czynnik okreœlaj¹cy szerokoœæ brzegu modelu, o wartoœciach 50 lub 100 oczek siatki; A — wspó³czynnik okreœlaj¹cy charakter zmian wartoœciF w obszarze brzegowym, o wartoœciach: B/50 dla liniowej zmiany lub B/20 dla zmiany wyk³adniczej;T — czêstoœæ drgañ,Fi— wartoœæ przewodnoœci elektrycznej w i-tym wêŸle wewn¹trz obszaru t³umienia na brzegu modelu.
W modelu dyskretnym budowê geologiczn¹, niejedno-rodnoœæ i anizotropowoœæ uzyskujemy poprzez przypisy-wanie ro¿nych parametrów materia³owych:
gr[-] — wzglêdna przenikalnoœæ elektryczna oœrodka,
:r[-] — wzglêdna przenikalnoœæ magnetyczna oœrodka,
F [S/m] — przewodnoœæ elektryczna oœrodka,
do odpowiednich oczek siatki obliczeniowej. Na gór-nym brzegu modelu wprowadza siê warstwê o parametrach powietrza. Dla pró¿ni wielkoœci powy¿sze maj¹ okreœlon¹ wartoœæ: g0=8,85 · 10 -12 [F/m] :0=1,25 · 10 -6 [H/m] F =0 [S/m] W oœrodku materialnym jakim jest górotwór wartoœci g, :, F, okreœlamy nastêpuj¹co:
g = gr× g0 : = :r× :0 F > 0
ród³o wymuszaj¹ce propagacjê fali (Ÿród³o punktowe lub fala p³aska), modelowane jest poprzez zaczepienie wektorów pr¹dowych w okreœlonych oczkach siatki, a ich zachowanie sterowane jest przez odpowiedni¹ funkcjê Ÿród³a np. sygna³ sinusoidalny, sygna³ Gaussa, sygna³ Ric-kera, sygna³ Blackmana-Harrisa. Najbardziej zaawanso-wan¹ technik¹ budowy Ÿród³a wymuszaj¹cego propagacjê fali jest budowa dodatkowej gêstej siatki obliczeniowej
tzw. sub-grid i skonstruowanie w niej odpowiedniego Ÿród³a wielowektorowego symuluj¹cego rzeczywist¹ radiacjê energii z anteny — tzw. radiation pattern (Carcio-ne, 1998). W programie REFLEX-W dysponujemy kilko-ma typami Ÿróde³ fal e.m.: Ÿród³o punktowe, fala p³aska, exploding reflector. Rozwi¹zanie numeryczne dla Ÿród³a typu exploding reflector jest oparte na zasadzie Huygensa i w uproszczeniu polega na przeniesieniu odpowiednio przeliczonego Ÿród³a na granicê odbijaj¹c¹ i modelowanie propagacji fali od momentu jej odbicia. Dok³adny opis tego typu Ÿród³a mo¿na znaleŸæ m.in. w pracy Carcione i in. (2002).
Przyk³adowy echogram pomiarowy
Jak wspomniano we wstêpie do przetestowania zaproponowanej technologii modelowania numerycznego wybrano „s³aby” echogram pomiarowy, który zosta³ odrzu-cony w procesie interpretacji badañ georadarowych. Bada-nia te przeprowadzono we Wroc³awiu w celu inwentaryzacji podziemnych kabli, rur i kana³ów pod chodnikiem i jezdni¹.
Pomiary wykonano aparatur¹ RAMAC szwedzkiej fir-my MALA GeoScience. W badaniach u¿ywano anteny ekranowanej 500MHz. Bior¹c pod uwagê wysok¹ czêsto-tliwoœci sygna³ów wysy³anych z tej anteny oraz œredni¹ wartoœæ prêdkoœci dla badanego oœrodka (vsr=0,08 m/ns)
mo¿na próbowaæ lokalizowaæ cia³a anomalne o rozmiarach kilkunastu centymetrów do g³êbokoœci kilku metrów. Ze wzglêdu na wystêpowanie utworów ilasto-gliniastych w badanym rejonie nale¿y spodziewaæ siê du¿ego t³umienia w oœrodku.
Na echogramie odleg³oœæ miêdzy trasami wynosi 0,04 m, sygna³ zosta³ spróbkowany z krokiem 0,2ns, i zastoso-wano 16-krotne sk³adanie. Po przeprowadzeniu wstêpnego przetwarzania z u¿yciem procedur (REFLEX-W Manual, 2003):
usuniêcie przesuniêæ tras usuniêcie ringingu usuniêcie dudnieñ wyg³adzenie medianowe dobór wzmocnieñ u¿ycie funkcji kroskorelacji przesuniêcie czasowe 0 korekcja faz
wyznaczono z echogramu prawdopodobne po³o¿enie kana³ów, kabli i rur. Wstêpna analiza echogramu pozwoli³a na budowê modelu numerycznego, który nastêpnie u¿yto w procesie interpretacji.
Wyniki przeprowadzonych analiz numerycznych
Wyniki modelowania numerycznego prezentowane w artykule uzyskano przy u¿yciu programu REFLEX-W nie-mieckiej firmy Sandmeier-Geo.
Pakiet oprogramowania REFLEX-W, pozwala na: przetwarzania danych georadarowych i sejsmicznych, pro-wadzenie analiz CMP, interpretacjê 3D danych georadaro-wych i sejsmicznych, modelowanie numeryczne przy pomocy metody FDTD oraz na analizê danych tomogra-ficznych i refrakcyjnych (REFLEX-W Manual, 2003). Modu³ do symulacji numerycznych pozwala na modelowa-nie pola w uk³adzie dwu- i trójwymiarowym, z uwzglêd-nieniem w analizach numerycznych anizotropii i niejednorodnoœci modelowanego oœrodka. Program daje mo¿liwoœæ modelowania ró¿nych pól falowych: elektro-magnetycznego, sejsmicznego (wraz z polem fali SH) oraz
t[ns]60 00 10 2 3 4 5 20 1 30 50 40 6 7 0 1 2 3 4 90 80 70 100 x [m] h [m]
Ryc. 2. Echogram pomiarowy po przetwarzaniu Fig. 2. Real echogram after processing
akustycznego. W analizowanym modelu mo¿na uwzglêd-niæ topografiê, informacje otworowe oraz mo¿na umieœciæ w nim dowoln¹ strukturê geologiczn¹ i cia³o anomalne o dowolnym kszta³cie. Program pozwala na obliczenie echo-gramów syntetycznych jak równie¿ obrazu propagacji fali w oœrodku (REFLEX-W Manual, 2003).
Przeprowadzone symulacje mia³y za zadanie pokazaæ mo¿liwoœci narzêdzi jakim jest modelowanie numeryczne w procesie interpretacji wyników georadarowych. Na przyk³adowym echogramie pomiarowym (ryc. 2) próbo-wano z pomoc¹ modelowania numerycznego identyfiko-waæ rodzaj obiektów infrastruktury podziemnej oraz okreœliæ ich po³o¿enie i rozmiary.
W czasie analiz numerycznych przeprowadzono kilkana-œcie symulacji dla wybranego profilu, dobieraj¹c odpowied-nio geometriê modelu i sta³e materia³owe. W symulacjach numerycznych badano wp³yw ró¿nych cia³ anomalnych tj. rur, kabli i kana³ów o ró¿nych przekrojach i po³o¿eniu, aby okreœliæ ich obraz na echogramie syntetycznym.
Modelowanie by³o prowadzone w dwuwymiarowej tarczy obliczeniowej w uk³adzie x–y (uk³ad 2D). Do anali-zy pranali-zyjêto model o wymiarach x = 8 m (odpowiednio do d³ugoœci profilu) oraz y = 4 m (odpowiednio do g³êbokoœci dla wyznaczonej prêdkoœci œredniej). Model zosta³ zdyskre-tyzowany siatk¹ obliczeniow¹ o oczkach kwadratowych, o wymiarze)x = )y = 0,01 m, zgodnie z zale¿noœciami opisa-nymi w czêœci teoretycznej. Górny brzeg modelu opisano jako warstwê powietrza o parame-trach: gr=1, :r=1, F=0[S/m]. Na
dolnym, lewym i prawym brzegu modelu przyjêto t³umi¹ce warun-ki brzegowe — stref szerokoœci 50 oczek z liniowym wzrostem przewodnoœci od wartoœci zada-nej w modelu doFmax.
W czasie symulacji testowych analizowano ró¿ne dostêpne w pro-gramie typy Ÿróde³ fal elektroma-gnetycznych, opisane w czêœci teoretycznej. W obliczeniach koñcowych Ÿród³o fali e.m. (symu-lacja anteny nadawczej) by³o modelowane poprzez procedurê exploding reflector, która daje naj-lepsze wyniki dla zastosowanego uk³adu pomiarowego. Parametry Ÿród³a dobrano tak, aby odpowia-da³y sygna³om wysy³anym przez antenê 500 MHz.
Opis w³aœciwoœci gruntu oraz cia³ anomalnych (rury, kable, kana³y) przeprowadzono poprzez przypisanie do wêz³ów siatki odpowiednio dobranych sta³ych materia³owychgr,F i :r.
Dla wszystkich materia³ów w modelu przyjêto sta³¹ wartoœæ :r=1. Parametry materia³owe dla
elementów stalowych — rury oraz kable elektryczne — przyjê-to jako: gr=1, F=10[S/m]. Na
bazie informacji geologicznych okreœlono, ¿e utwory buduj¹ce przypowierzchniow¹ warstwê gruntu to mieszanina i³u, piasku gruboziarnistego, gliny i ¿wiru. Dla uwzglêdnienia wyraŸnego przechodzenia jednych utworów w inne np. gliny w ¿wir, testowa-no w analizach numerycznych zmienne parametry materia³owe w odpowiednich czêœciach mode-lu mode-lub stosowano gradientow¹ zmianê tych wartoœci.
4,0m 0 1 2 3 x [m]4 5 6 7 t[ns] 60 0 10 20 30 50 40 90 80 70 100 rura z PCV pipe of PCVf=0,1 m, =3, =1, =0,0001 [S/m]e m s rura metalowa metal pipe f=0,1 m, =1, =1, =10 [S/m]e m s pustka z powietrzem
void with air e=1, =1, =0 [S/m]m s
betonowy fundament foundation of concrete e=5, =1, =0,005 [S/m]m s wilgotny piasek wet sand e=9, =1, =0,01 [S/m]m s 0 0,5 0,4 1,0 3,0 0 0,9 1,9 2,1 3,4 5,4 6,0 8,0m A B
Ryc. 4. A — model testowy o parametrach gruntu piaszczystego, z cia³ami anomalnymi; B —
echogram syntetyczny dla modelu (A)
Fig. 4. A — numerical model, wet sand ground with anomaly bodies, B — synthetic
echo-gram for the model (A)
4,0m 0 1 2 3 x [m]4 5 6 7 t[ns] 60 0 10 20 30 50 40 90 80 70 100 rura z PCV pipe of PCVf=0,2 m, =3, =1, =0,0001 [S/m]e m s rura metalowa metal pipe f=0,2 m, =1, =1, =10 [S/m]e m s kana³ z powietrzem
sewer with air e=1, =1, =0 [S/m]m s
cia³o metalowe metal body e=1, =1, =10 [S/m]m s i³+glina+piasek+¿wir loam+clay+sand+gravele=15, =1, =0,05 [S/m]m s 0 0,5 0,4 0,7 1,0 1,5 0 0,9 1,9 2,1 3,4 5,4 6,0 8,0m A B
Ryc. 5. A — model testowy o parametrach gruntu bêd¹cego mieszanin¹ i³u, gliny, piasku,
¿wi-ru, z cia³ami anomalnymi, B — echogram syntetyczny otrzymany dla modelu (A)
Fig. 5. A — numerical model: loam, clays, and gravel ground with anomaly bodies, B —
syn-thetic echogram for the model (A)
4,0m 0 1 2 3 4 5 6 7 x [m] t[ns] 60 0 10 20 30 50 40 90 80 70 100 rura z PCV pipe of PCVf=0,5 m, =3, =1, =0,0001 [S/m]e m s rura metalowa metal pipe f=0,5 m, =1, =1, =10 [S/m]e m s kana³ z powietrzem
sewer with air e=1, =1, =0 [S/m]m s
kana³ wype³niony wod¹ sewer with water e=80, =1, =0,5 [S/m]m s i³+glina loam+claye=20, =1, =0,1 [S/m]m s 0 0,5 0,4 0,7 1,0 1,5 0 0,9 1,9 2,1 3,4 5,4 6,0 8,0m A B
Ryc. 3. A — model testowy o parametrach gruntu ilasto-gliniastego, z cia³ami anomalnymi, B
— echogram syntetyczny dla modelu (A)
Fig. 3. A — numerical model of a loam and clay ground with anomaly bodies, B — synthetic
Wyznaczono krok czasowy do obliczeñ)t=0,01ns — odpowiedni do czêstotliwoœci, zgodnie z wczeœniejszym opisem (warunek CFL), a obliczenia prowadzono w oknie czasowym T=100ns.
Na rycinach 3–5 przedstawiono kilka wybranych modeli testowych oraz model koñcowy (ryc. 6) dla przedstawione-go na ryc. 2 echogramu pomiaroweprzedstawione-go. Na koñcu rozdzia³u zamieszczono wyniki symulacji dla modelu koñcowego w postaci echogramu syntetycznego dla sk³adowej Ey elektro-magnetycznego pola falowego.
Analiza wyników symulacji numerycznych
W analizach numerycznych pomijamy modelowanie fal bezpoœrednich tj. fali powietrznej i gruntowej, ponie-wa¿ fale te nie nios¹ informacji o obiektach podziemnych. Dodatkowo fale te skutecznie utrudniaj¹ interpretacjê p³ytko po³o¿onych cia³ anomalnych.
Zastosowanie w badaniach anteny 500 MHz daje nie-wielk¹ mo¿liwoœæ rejestracji niejednorodnoœci w pasie przypowierzchniowym (ok. 0,5 m) — w asfalcie i podsyp-ce ¿wirowej. Dla potwierdzenia zarejestrowanych w tym
pasie niejednorodnoœci nale¿y przeprowadziæ pomiary antenami wysokoczêstotliwoœciowymi np. 1GHz, co da mo¿liwoœæ interpretacji iloœciowej i jakoœciowej.
Aby najlepiej dopasowaæ parametry materia³owe modelu numerycznego do sk¹pych informacji geologicz-nych, testowano ró¿ne modele górotworu: wilgotny, zailo-ny piasek (ryc. 4), i³ i glina (ryc. 3), model o parametrach uœrednionych (ryc. 5), model z gradientow¹ zmian¹ para-metrów od zailonego piasku do mieszaniny i³u i gliny (ryc. 6). Badano tak¿e efekt od modeli warstwowanych. Dobie-rano odpowiednio wartoœci materia³owe oœrodka w prze-dziale: gr = 6–20, F = 0,001–0,1[S/m] aby nachylenie
ramion hiperbol na echogramie syntetycznym by³o zgodne z echogramem pomiarowym. Wszystkie dodatkowe reflek-sy o ma³ych amplitudach na echogramie pomiarowym s¹ wynikiem niejednorodnoœci oœrodka (np. wk³adki ¿wirowe w piasku czy glinie), ale nie posiadaj¹c dok³adnych infor-macji geologicznych nie mo¿na ich zamodelowaæ.
Na echogramie pomiarowym widaæ silne t³umienie oœrodka ze wzglêdu na jego budowê (i³y, gliny), a tak¿e zauwa¿a siê niewielki zasiêg g³êbokoœciowy badañ (do ok. 1 m) co wynika z silnego t³umienia fal wysokoczêstotliwo-œciowych. Efekty te potwierdzaj¹ siê na echo-gramach syntetycznych, co wskazuje na poprawny dobór parametrów materia³owych modelu.
Aby dopasowaæ kszta³t hiperboli w modelu do rzeczywistego kszta³tu z pomiarów zmienia-no œrednicê rur i testowazmienia-no dwa materia³y: stal i PCV. W wyniku okaza³o siê, ¿e nie ma znacze-nia czy s¹ to rury z wod¹ czy gazem lub czy jest to plik metalowych przewodów. Na kszta³t hiperboli i intensywnoœæ refleksów ma wp³yw œrednica i materia³ (ryc. 3–5); w wyniku ekspe-rymentów numerycznych ustalono, ¿e s¹ to dwie rury stalowe o œrednicach rzêdu 10–20 cm, po³o¿one w maksimum hiperbol tj. na x = 2,1 m i x = 3,4 m.
Aby okreœliæ od jakiego typu cia³a nast¹pi³o odbicie w lewej (x = 0,9–1,9 m) i prawej (x = 5,4–6,0 m) czêœci echogramu pomiarowego sprawdzano w symulacjach efekty od: zakopa-nych obiektów stalowych (ryc. 5), fundamen-tów betonowych (ryc. 4), kana³ów lub pustek wype³nionych powietrzem lub wod¹ (ryc. 3). Dobierano tak¿e optymalny kszta³t samych obiektów, aby generowany przez nie efekt by³ zbli¿ony do rzeczywistego. Na echogramie pomiarowym nie zauwa¿a siê dyfrakcji na kan-tach kana³ów (pustek), wiêc prawdopodobnie maj¹ one ³aman¹ liniê stropu. Nie s¹ to tak¿e kana³y z pó³okr¹g³ym sklepieniem, poniewa¿ efekt na echogramie by³by podobny jak od rur. Ani na echogramie pomiarowym ani na synte-tycznym nie mo¿na uchwyciæ dolnej granicy pustki, wiêc do modelowania przyjêto kszta³t kwadratowy o wysokoœci równej szerokoœci (znanej z echogramu pomiarowego). Aby nie pojawia³a siê dyfrakcja na rogach kwadrato-wych pustek, zast¹pione je lini¹ ³aman¹ (ryc. 6). Na granicy pustki (kana³u z powietrzem) nastê-puje kilkukrotny wzrost prêdkoœci fali, a na
metalowe rury instalacji wodnej lub gazowej =0,15 m, =1, =1, =10 [S/m]
gas or water metal pipes
f e m s
kana³y lub pustki wype³nione wod¹ =80, =1, =0,5 [S/m] sewers or voids with water
e m s
¿wir(gravel e) =10, =1, =0,001 [S/m]m s
warunek brzegowy typu ABC Absorbing Boundary Condition powietrze air e=1, =1, =0 [S/m]m s asfalt (asphalt)e=5, =1, =0,002 [S/m]m s wilgotny piasek+¿wir =9, =1, =0,05 [S/m] wet sand+gravel e m s i³+glina =20, =1, =0,1 [S/m] loam+clay e m s 0 0,9 1,9 2,1 3,4 5,4 6,0 8,0 m 0 0,6 0,4 0,7 1,0 1,7 4,0 m D Dx= y=0,01 m smax s 50 oczek siatki grid points
Ryc. 6. Model koñcowy do obliczeñ numerycznych dla echogramu
pomiaro-wego (ryc.2)
Fig. 6. Final model for numerical simulation for a real echogram (Fig. 2)
t[ns] 60 0 0 10 2 3 4 5 20 1 30 50 40 6 7 90 80 70 100 x [m]
Ryc. 7. Syntetyczny echogram dla modelu koñcowego (ryc. 6) — sk³adowa Ey
elektromagnetycznego pola falowego
Fig. 7. Synthetic echogram for final model (Fig. 6) — Ey component of
echogramie zmniejsza siê w pustce czêstotliwoœæ impul-sów i pojawia siê efekt interferencji fal — widoczny na ryc. 3–5. Bior¹c to pod uwagê przy analizie echogramu pomia-rowego mo¿na domniemywaæ, ¿e kana³y (lub pustki) wype³nione s¹ wod¹.
Prawdopodobnie pomiêdzy rurami jest jeszcze jeden obiekt, na co wskazuje poziomy pas silnych refleksów, lecz echogram pomiarowy daje zbyt ma³o informacji, aby taki obiekt zamodelowaæ.
Zaproponowany model (ryc. 6) wydaje siê najlepiej odzwierciedlaæ budowê oœrodka gruntowego oraz geome-triê, po³o¿enie i rodzaj cia³ anomalnych, zwa¿aj¹c na bar-dzo ubog¹ informacjê geologiczn¹ i s³aby materia³ pomiarowy.
Podsumowanie
W artykule zaprezentowano technologiê wspomagania interpretacji wyników georadarowych przy u¿yciu mode-lowania numerycznego. Jak widaæ na podstawie zamieszczo-nego przyk³adu, zastosowanie modelowania numeryczego mo¿e przyczyniæ siê do ³atwiejszej interpretacji wyników pomiarowych, zw³aszcza w przypadku „s³abego” mate-ria³u pomiarowego. W przypadku, gdy po przeprowadzeniu przetwarzania echogramów nie uzyskamy zadawalaj¹cych wyników, mo¿na próbowaæ zastosowaæ opisywan¹ metodê w procesie interpretacji. Metody symulacyjne daj¹ najlepsze wyniki, gdy jesteœmy w stanie dok³adnie opisaæ geometriê i parametry materia³owe modelu. W przypadku zastosowania modelowania numerycznego w interpretacji geofizycznej znacz¹cym utrudnieniem w stosowaniu modelowania jest niewielka iloœæ informacji o analizowanym oœrodku. Infor-macje te (zazwyczaj punktowe) pochodz¹ z otworów lub ods³oniêæ w miejscu prowadzenia pomiarów. W przypadku prowadzenia badañ w rejonie o s³abo rozpoznanej geologii i parametrach materia³owych oœrodka, modelowanie numeryczne sprowadza siê do ca³ej masy testów minimali-zuj¹cych rozbie¿noœæ pomiêdzy uzyskanym wynikiem pomiaru, a otrzymanym echogramem syntetycznym. W takim przypadku zazwyczaj przyjmuje siê parametry mate-ria³owe przybli¿one (tabelaryczne), a wyniki takiej symula-cji uwzglêdnia siê w procesie interpretasymula-cji z du¿¹ ostro¿noœci¹ i jedynie w skali jakoœciowej. Zupe³nie inaczej wygl¹da sytuacja w momencie posiadania bogatego mate-ria³u geologicznego oraz badañ laboratoryjnych okreœlaj¹cych elektromagnetyczne parametry utworów buduj¹cych analizo-wany oœrodek. W takim przypadku symulacje numeryczne pozwalaj¹ na interpretacjê iloœciow¹, a tak¿e mo¿na pokusiæ siê o bardzo daleko posuniête wnioski interpretacje np. dok³adne kszta³ty i rozmiary cia³ anomalnych, okreœlenie
materia³u z jakiego wykonany jest dany obiekt (stal, beton, tworzywo syntetyczne, drewno itp.).
Jako, ¿e jest to pierwszy artyku³ opisuj¹cy metodê modelowania numerycznego z zastosowaniem w procesie interpretacji danych georadarowych, autor skupi³ swoj¹ uwagê na technologii modelowania i jednym przyk³adzie pomiarowym. Autor pracuje obecnie nad zastosowaniem modelowania numerycznego dla badañ georadarowych w uk³adzie 3D w oœrodkach niejednorodnych, anizotropo-wych z ograniczon¹ informacj¹ geologiczn¹ i stochastycz-nym rozk³adem parametrów modelowanego oœrodka. W nastêpnych artyku³ach zamieszczane bêd¹ wyniki zastoso-wania opisywanej metody w procesie interpretacji danych georadarowych dla coraz bardziej skomplikowanych sytu-acji geologiczno-in¿ynierskich.
Autor pragnie podziêkowaæ mgr in¿. J. Ziêtkowi z Zak³adu Geofizyki AGH za udostêpnienie wyników pomiarów georada-rowych, przeprowadzonych przez niego we Wroc³awiu.
Praca zosta³a zrealizowana w ramach grantu KBN nr: 5 T12B 036 23
Literatura
ANNAN A.P. 2001 — Ground Penetrating Radar. Workshop Notes, Sensor and Software Inc., Canada.
BERGMANN T., ROBERTSSON J. & HOLLIGER K. 1998 — Finite difference modeling of electromagnetic wave propagation in dispersive and attenuating media. Geophysics, 63: 856–867.
CARCIONE J.M. 1996 — Ground Penetrating Radar: wave theory and numerical simulation in lossy, anisotropic media. Geophysics, 61: 1664–1677.
CARCIONE J.M. 1998 — Radiation patterns for 2D GPR forward modelling. Geophysics, 63: 424–430.
CARCIONE J.M., PINERO F.L. & ZAMPARO M. 2002 — Exploding reflector concept for Ground Penetrating Radar modeling. Ann. Geoph., 45: 473–478.
HOLLIGER K. 2002 — Finite difference modeling of Ground Penetrating Radar data. 9-th International Conference on Ground Penetrating Radar, Santa Barbara, California, USA.
KARCZEWSKI J. 1997 — Zastosowanie metody georadarowej do wykrywania zaburzeñ naturalnych i antropogenicznych w warstwach przypowierzchniowych. Arch. Wydz. Geol. Geof. Och. Œrod., AGH, Kraków.
KOSLOFF R. & KOSLOFF D. 1986 — Absorbing boundaries for wave propagation problems. Jour. Comput. Phys., 63: 363–376. MORAWSKI T. & GWAREK W. 1998 — Pola i fale elektromagne-tyczne. WNT, Warszawa.
MUR G. 1981 — Absorbing boundary condition for finite difference approximation of time domain electromagnetic field equations. IEEE Trans. Electromag. Comput., 23: 1073–1077.
REFLEX-W Manual 2003 — K.J. Sandmeier, Karlsrue, Germany. ROBERTS R.L. & DANIELS J.J. 1997 — Modeling near-field GPR in three dimensions using FDTD method. Geophysics, 62: 1114–1126. YEE K.S. 1966 — Numerical solution of initial boundary value pro-blems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propagat., 14: 302–307.
