Fala refrakcyjna i wybrane metody interpretacyjne

Analiza cech kinematycznych fali sejsmicznej skupia się na geometrycznych aspektach przejścia fali pomiędzy źródłem a odbiornikiem (Cerveny 2001, Chapman 2004). Opis rozchodzącego się sygnału sejsmicznego można uprościć do postaci promienia sejsmicznego, posiadającego odpowiedni kierunek, zwrot i wartość długości. Podana definicja jest podstawą tzw. podejścia promieniowego w sejsmice przypowierzchniowej.

Rys. 5.1 - 1 Schemat propagacji fali refrakcyjnej oraz bezpośredniej (A) wraz przestawieniem ich hodografów (B).

88 Na drodze prostych przekształceń matematycznych można ją zastosować do opisu przejścia fali refrakcyjnej w ośrodku geologicznym (Cerveny 2001). Schemat propagacji wspomnianej fali przedstawia rysunek 5.1 - 1 - A. Warunki wystąpienia zjawiska refrakcji w ośrodku geologicznym są następujące: wzrastająca prędkość wraz z głębokością oraz odpowiedni układ warstw pozwalający na powstanie zjawiska całkowitego załamania fali. W uproszczonym modelu ośrodka, scharakteryzowanego przez wartości prędkości oraz miąższości każdej z warstw wzór hodografu fali refrakcyjnej definiuje się jako całkowity czas przejścia sygnału refrakcyjnego pomiędzy źródłem a dowolnie wybraną pozycją odbiornika (Rys. 5.1 - 1 - B):

𝑇 =^2ℎ0cos(𝑖_𝑐1) 𝑉1 +_𝑉^𝑋

2= 𝑇₀+_𝑉^𝑋

2 (5.1-1)

Gdzie 𝑉₁i 𝑉₂ są prędkościami w warstwie 1 i 2, ℎ₀ - stanowi miąższość warstwy pierwszej, x - jest odległością odbiornika od źródła oraz 𝑖_𝑐1 jest zdefiniowane wzorem:

𝑖_𝑐1 = sin(^𝑉1

𝑉₂) (5.1-2)

Przyrównując X do zera. Otrzymano:

ℎ₀ = ^𝑇0𝑉₁

2 cos(𝑖_𝑐) (5.1-3)

Gdzie 𝑇₀jest czasem cięcia. Wprowadzając do wzoru 5.1-3 wyrażenie 5.1-2. Określono zależność pomiędzy 𝑋_𝑐, 𝑋_𝑐𝑟, ℎ₀, 𝑉₁ i 𝑉₂:

2 ∗ 𝑋_𝑐~ ℎ₀ =^𝑋𝑐𝑟

2 √^𝑉2−𝑉₁

𝑉2+𝑉1 (5.1-4)

Dla n-tej ilości warstw ℎ_𝑛 jest równe:

ℎ_𝑛 = ^𝑉𝑛 cos(𝑖_𝑛)(^𝑇0𝑛+1 2 − ∑ ℎ_𝑗√𝑉¹ 𝑗²−_𝑉¹ 𝑗+1² 𝑛−1 𝑗=0 ) (5.1-5)

Przekształcając równanie (5.1-5) względem 𝑇_0𝑛 zyskuje ono postać: 𝑇_0𝑛= ^ℎ𝑛cos (𝑖𝑛) 𝑉_𝑛 −^𝑇0 2+∑ ℎ_𝑗√𝑉¹ 𝑗²−_𝑉¹ 𝑗+1² 𝑛−1 𝑗=0 (5.1-6)

W czasie analizy danych refrakcyjnych korzystano z dostępnych narzędzi interpretacyjnych. Wykorzystano metodę czasu cięcia oraz tomografię refrakcyjną. Analizy

89 przeprowadzono w dwóch programach komercyjnych. Model refrakcyjny obliczono przy użyciu systemu przetwarzania Vista Seismic Processing. Inwersję tomograficzną wykonano w programie SeisImager. Metodyki dobrano pod względem możliwości obrazowania ośrodka polodowcowego. Metoda czasu cięcia (ang. Intercept Time Method) reprezentuje uproszczoną formę analizy hodografów refrakcyjnych (Hampson 1984, Wyrobek 1956). Początkowa faza interpretacji polega na pozyskaniu hodografów fali refrakcyjnej wzdłuż profilu sejsmicznego (Rys. 5.1 - 2 - B). Każdy z wybranych zbiorów pierwszych wstąpień jest zdefiniowany powierzchniowo w kolekcji oznaczonej poligonem (Rys. 5.1 - 2 - C), nazywanym z języka angielskiego „Controll Pointem”. W każdej z kolekcji przez analizę regresji liniowej wyznacza się wartość nachylenia fali bezpośredniej i refrakcyjnej, estymując tym samym wartość prędkości refraktora i nadkładu.

Rys. 5.1 - 2 Interpretacja danych refrakcyjnych przy użyciu metody czasu cięcia . Przy użyciu wzoru 5.1-3 oraz 5.1-6 zostaje obliczona głębokość zalegania refraktora układu dwu - lub n - warstwowego. Konstrukcja modelu refrakcyjnego wymaga pozyskania wielu kolekcji pierwszych wstąpień, na bazie których estymowane są jednowymiarowe rozkłady fali P. Rezultatom zostają przyporządkowane współrzędnych centrów Controll Point’ów wzdłuż profilu sejsmicznego (Rys. 5.1 – 2 - A). Technika posiada ograniczenia.

90 Wartości prędkości pozornych refraktora są obarczone błędem w przypadku, gdy obrazowana warstwa refrakcyjna wykazuje silne zróżnicowanie morfologiczne. Bardziej złożone metody interpretacji hodografów refrakcyjnych, takie jak: metoda Plus - Minus (Hagedoorn 1959), GRM (Palmer 1981, ang. Generalized Reciprocal Method) czy metoda czasu opóźnień (Barry 1967, ang. Delay Time Method), odznaczają się podobną wrażliwością na zmiany elewacji podłoża. We wszystkich analizach materiałów refrakcyjnych - poza tomograficznymi, interpretacja przebiega przy założeniu, że refraktor nie wykazuje większych zmian elewacji niż 10 %. Jest to związane ze zmianą prędkości pozornej fali refrakcyjnej, propagującej pod upad lub z upadem warstwy. Na rysunku 5.1 - 2 „plusami” (czerwone oznaczenie) i „minusami” (niebieski marker) oznaczono dodatnie i ujemne odległości odbiorników względem pozycji źródła (offsety). Ich wzajemne zestawienie w jednej kolekcji Controll Point powoduje destabilizację gałęzi refrakcyjnych (ang. refraction branches), których liniowa postać warunkuje właściwą interpretację wartości prędkości. Ominięciem tego problemu w metodzie czasu cięcia jest zastosowanie sortowania względem absolutnych wartości offsetowych pierwszych wstąpień. Pozwala to na kontrolę estymowanej wartości prędkości mimo zróżnicowania nachylenia hodografów refrakcyjnych.

Rys. 5.1 - 3 Poglądowy rysunek przedstawiający schemat działania algorytmu tomografii refrakcyjnej.

91 Tomografia refrakcyjna jest iteracyjnym rozwinięciem metod refrakcyjnych, polegającym na dopasowaniu teoretycznych (obliczonych) czasów przyjścia do ich polowych odpowiedników (Hakan 2007, Watanabe 1999). Metoda ta odznacza się dokładnością odwzorowania modelu strefy przypowierzchniowej. Uproszczony zarys działania algorytmu tomografii przedstawia rysunek 5.1 - 3. Przedstawia on model ośrodka podzielony na strefy (cele) o jednakowej wartości prędkości. Linią szarą, ciągłą zaznaczono granicę strukturalną warstwy. Wartościami pomierzonymi w trakcie badań jest czas przyjścia fali refrakcyjnej w punkcie C. Pomiędzy punktami A, B i C zaznaczono drogę propagacji promienia sejsmicznego. Podzielono go na odcinki, tak że całkowity czas przejścia na odcinku A - B - C można obliczyć za pomocą następującej sumy:

𝑡_𝑖𝑗 = ∑ ^𝑑𝑖𝑗

𝑣𝑗

𝑁

𝑗=1 = ∑𝑁 𝑑_𝑖𝑗

𝑗=1 𝑠_𝑗 (5.1-7)

,gdzie: 𝑡_𝑖 jest czasem przejścia fali refrakcyjnej na i - tej drodze , 𝑑_𝑖𝑗 jest długością całkowitej drogi przebytej przez sygnał refrakcyjny w j - tej celi posiadającej prędkość 𝑣_𝑗.W zapisie macierzowym równanie 5.1-7 zyskuje postać:

T=D*S (5.1-8)

Równanie 5.1-8 wyjaśnia podstawową zależność pomiędzy obserwowanym czasem przyjścia fali refrakcyjnej i prędkością. Dokładne rozwiązanie tomograficzne polega znalezieniu zależności między danymi pomiarowymi a parametrami estymowanego modelu. Szerszą grupą danych pomierzonych - poza pierwsze wstąpienia fali P, są informacje o geometrii układu pomiarowego, natomiast dodatkowym parametrem - poza informacjami o geometrii warstw i ich wewnętrznym rozkładzie prędkości, jest błąd dopasowania hodografów teoretycznych do rzeczywistych. Związek ten symbolizuje następujące równanie,:

T( obs)= D*S(true) + error (5.1-9) Wektor „T(obs)” reprezentuje dane pozyskane na bazie analizy i pikowania pierwszych wstąpień, wektor „S(true)” parametry modelu, macierz D stanowi macierz układu/przejścia pomiędzy wspomnianymi wektorami, natomiast „error” oznacza błąd estymacji. Generalną ideą rozwiązania tomograficznego jest podanie zbioru parametrów estymowanego modelu, w taki sposób aby liczone na ich podstawie pierwsze teoretyczne wstąpienia, różniły się tylko w

92 nieznacznym stopniu od ich zmierzonych ekwiwalentów. Syntetyczną postać (przewidywaną) pierwszych wstąpienia można wyrazić w następujący sposób:

T(teor) = D*S(teor) (5.1-10) ,gdzie, „T(teor)” jest wektorem oznaczającym teoretyczne pierwsze wstąpienia, wyznaczone na bazie teoretycznych parametrów (np. danych tablicowych lub informacji geologicznych) ośrodka „S(teor)”, pozyskanych w procesie trasowania promienia sejsmicznego. Spełnienie założenia o satysfakcjonującym dopasowaniu między wartościami „T(obs)” i „T(teor”) należy zdefiniować w następujący sposób:

min(|T( obs)-T(teor)|)=min(|D(S(true)- S(est) )+ error|) (5.1-11) Ze względu na fakt, że wartość „error” jest różna od zera, możliwe jest tylko podanie przybliżonej informacji o ośrodku, zawartych w wartości „S(est)”. Jeśli zagadnienie jest źle postawione (np. rozwiązanie nie istnieje, lub jest niejednoznaczne), wtedy standardowym podejściem jest rozwiązanie w rozumieniu metody najmniejszych kwadratów. Oznacza to poszukiwanie minimum (proces minimalizacji) funkcjonału „F” w postaci:

F= (T(obs)- D*S(est))^T(T(obs)-D*S(est)) (5.1-12) ,gdzie:

𝑑𝐹

𝑑(𝑆(𝑒𝑠𝑡))^{= 0 (5.1-13)}

,następnie przekształcając do postaci:

D^TDS(est)- D^TDT(obs)=0 (5.1-14)

,przyrównując dalej do zera równanie zyskuje postać:

S(est)= (D^TD)^-1D^TT(obs) (5.1-15) Równanie przyjmuje postać wyrażenia pozwalającego na rozwiązanie postawionego zadania odwrotnego. Aby poprawić właściwości powyższego wyrażenia w sensie rozwiązywalności (np. wyeliminować oscylacje w rozwiązaniu, związane z błędem estymacji), należy wprowadzić czynnik określający względny udział, każdego indywidualnego błędu do błędu całkowitego predykcji. Do powyższej postaci dodaje się człon regularyzacyjny

93 M(reg)d^TM(reg)d, wyrażony przez macierz regularyzacji „M(reg)” – informacje a priori o wektorze wartości estymowanych/niewiadomych „S(est)”. Wtedy równanie otrzymuje postać:

S(est)= (D^TM(reg)d^TM(reg)d D)^-1D^TM(reg)d^TM(reg)d T(obs) (5.1-16) Jednoznaczne rozwiązanie powyższego równania w rzeczywistości może być trudne do osiągnięcia, ponieważ macierz D^TM(reg)_d^TM(reg)_dD przybiera formę osobliwą.

Uregulowanie procesu inwersji, przez wprowadzenie dodatkowych współczynników zapewniających pomocnicze informacje a priori, mogą narzucić więzy wspomnianemu procesowi, prowadząc do odwracalności macierzy D^TM(reg)d^TM(reg)dD. Poprawę jej uwarunkowania można osiągnąć poprzez wprowadzenie do powyższego wzoru wyrażenia α*I, gdzie α to tzw. współczynnik tłumienia (ang. Damping Factor), natomiast I jest macierzą jednostkową. Funkcję współczynnika tłumienia można określić jako kontrolującą postęp zmian wartości param(est) w stosunku do ich wartości referencyjnych S(est)ref

:

αj*S(est)j= αj*S(est)j^ref (5.1-17)

Podczas przeprowadzania inwersji wymagany jest początkowy zestaw oszacowanych parametrów modelu, które w przypadku powyższego równania odnoszą się do wartości czynnika S(est)j^ref. Rozkład referencyjny zostaje zastąpiony przy każdej iteracji w zakresie (j=1,2…n.). Dodatkowo, współczynnik tłumienia jest parametrem, który określa w sposób ilościowy dopuszczalną zmianę, przy każdym kroku iteracji w porównaniu z poprzednim modelem odniesienia (określa stopień regularyzacji). Parametr regularyzacyjny jest włączony przez dodanie nowego równania dla każdego parametru zgodnie z modelem:

S(est)= (D^TM(reg)_d^TM(reg)_d D+α²I)^-1(D^TM(reg)_d^TM(reg)_d T(obs) α²IS(est)^ref) (5.1-18) Podstawą metod refrakcyjnych, a także i tomograficznych jest pozyskiwanie informacji o podłożu danego ośrodka geologicznego (Oluwafemi 2009). Problem odwzorowania tomograficznego w ośrodkach – również polodowcowych, nie wynika jedynie z efektywności algorytmu inwersji, ale zależy głównie od postaci i jakości pierwszych wstąpień (Vidale 1988; Trampert 1990; Tryggvason 2002, 2006).

94