4. Badania otworowe i wielokanałowa analiza fal powierzchniowych
4.3 Adaptacja schematu przetwarzania MASW do obrazowania ośrodka polodowcowego
polodowcowego
Do analizy zapisów fali powierzchniowej użyto komercyjnego oprogramowania SeisImager oraz GeopsyPack. Przetwarzanie danych rozpoczęto od transformacji pozyskanych rekordów sejsmicznych do dziedziny częstotliwość - prędkość fazowa. Procedurę przeprowadzono zgodnie z definicją transformacji (4.2-7). Wprowadzono przedziały prędkości fazowej: 50 -1000 m/s i częstotliwości 2- 60 Hz. Trendy dyspersyjne pikowano na mapach znormalizowanej amplitudy określonych wzorem (4.2-9). Dane charakteryzowały się wielomodalnym trendem. Szczególną ich właściwością było współwystępowanie mód wyższych przy nagłym wygaszeniu mody podstawowej. Powodowało to oddzielenie trendów modalnych i trudności w oznaczeniu ostatecznego wzorca krzywej dyspersji. W niektórych przypadkach mody wyższe posiadały charakter szczątkowy. W początkowej fazie przetwarzana założono, że są to naturalne zniekształcenia wynikające ze zniekształcenia rejestrowanej fali powierzchniowej. Zjawisko to utożsamiono z efektem występowania silnych kontrastów impedancji akustycznych (ew. obecność eratyków w materiale polodowcowym) i rozpraszania fali powierzchniowej. Z tego względu przyjęto kryterium, definiujące modę wyższą jako podniesienie prędkości fazowej względem mody
51 podstawowej o zakresie 8 Hz. Warunek przyjęto na podstawie obserwacji. Modę podstawową zdefiniowano, śledząc najwyraźniejszy trend krzywej dyspersji. Nie zawsze oznaczało to, że poddano analizie modę podstawową, ponieważ rozstrzygniecie tej kwestii w ośrodku polodowcowym nie jest jednoznaczne (Park 1995b, 2001). Po zakończeniu pozyskiwania krzywych polowych, rozpoczęto proces ich dopasowywania do krzywych teoretycznych. Adaptacja procesu przetwarzania polegała na wprowadzeniu informacji „a priori” do początkowych rozkładów fizycznych m określających ośrodek polodowcowy. Wartość prędkości poprzecznej dla każdego z rekordów definiowano jako:
𝛽𝑝𝑜𝑐𝑧= 0.86*𝐶𝑅, (4.3-1)
,gdzie 𝐶𝑅 pochodziły z wcześniej wypikowanych krzywych dyspersji. Pobierano je z interwałem 2 Hz (Rys. 4.3 - 1).
Rys. 4.3 - 1 Przykład krzywej dyspersji. Definicja zasięgu warstw przebiegała przy wykorzystaniu wzoru (4.3-1):
𝜆𝐶𝑅 = 𝐶𝑅
𝑓𝐶𝑅 (4.3-2)
,gdzie 𝑓𝐶𝑅 to częstotliwość fali propagującej z prędkością fazową 𝐶𝑅, natomiast 𝜆𝐶𝑅 to długość fali. Miąższość falowodu h można zdefiniować jako wartość w przedziale < 0.3 𝜆𝐶𝑅
52 ,0.5 𝜆𝐶𝑅>. Dla uproszczenia obliczeń przyjęto jego górną wartość. Rozkład gęstości pozornej obliczono na podstawie wzoru (Meyne 1999):
𝜌 = 1 +0.614+58.7(log(𝑧)+1.095)/𝑉1
𝑠 (4.3-3)
Gdzie z oznacza głębokość, 𝑉𝑠- prędkość fali poprzecznej. Dla potrzeb parametryzacji procesu inwersji przekształcono go do postaci wyrażenia 4.3-4:
𝜌𝐶𝑅 = 1 + 1
0.614+58.7(log(𝜆𝐶𝑅/2)+1.095)/𝐶𝑅 (4.3-4) Wykorzystanie wzoru 4.3-4 pozwoliło na uzależnienie gęstości od postaci krzywej dyspersji.
Początkowy rozkład fali P starano się powiązać ze znanymi wartościami 𝐶𝑅, przez korelację wartości Vp i Vs pozyskanych na podstawie studium literaturowego (Załącznik A, Tab. 1 - 4) oraz wyników profilowań otworowych (Rys.4.1-2 i 4.1-3).
Zależność pomiędzy prędkością fali poprzecznej i podłużnej w utworach przypowierzchniowych jest niejednoznaczna i lokalnie uwarunkowana. Nie wynika to jedynie z właściwości litologicznych czy geomechanicznych ośrodka, ale przede wszystkim jest związane z występowaniem stref zawodnionych, które odgrywają szczególną rolę w opisie i estymacji przypowierzchniowych pól prędkości. W przypadku utworów strefy przypowierzchniowej - w tym utworów polodowcowych, określenie relacji pomiędzy Vp i Vs należy rozpocząć od nakreślenia cech pól prędkości na bazie, których oblicza się wspomniane współczynniki. Postępująca konsolidacja materiału przypowierzchniowego powoduje wzrost prędkości fali podłużnej i poprzecznej wraz z głębokością. Przyrost wartości prędkości odbywa się przy zmianach we właściwościach mechanicznych ośrodka, wyrażonych ilościowo za pomocą modułów sprężystości. Rozchodzenie się fali poprzecznej w jednorodnej izotropowej półprzestrzeni zależy wprost proporcjonalnie od zmian wartości modułu sprężystości postaciowej, a odwrotnie proporcjonalnie do zmian gęstości w ośrodku (Aki i Richard 2002). Wyjaśnia to brak jej reakcji na obecność stref saturacji w materiale nieskonsolidowanym (West 2001). Odwrotnie reaguje fala podłużna, której rozchodzenie zależy od modułu ścinania a także od modułu sprężystości objętościowej, co skutkuje reakcją na wypełnienie porów. Wspomnianą zależność podaje teoria Biota - Gassman’a (Biot 1956). Wrażliwość prędkości fali P na zmiany saturacji utożsamia się ze zmianą w module sprężystości objętościowej. Równocześnie neguje ona wpływ wypełnienia porowego na zmiany w wartości prędkości fali poprzecznej. Opisana zależność w skałach zwięzłych
53 znalazła zastosowanie głównie w kontekście złożowym. Charakter sedymentów powierzchniowych trudno wiązać z określeniem zwięzłości, jednak generalne sformułowanie Biota - Gassman’a o braku związku pomiędzy wartością prędkości Vs i saturacją w sejsmice płytkiej pozostaje aktualne. Ośrodek porowaty jakim są utwory polodowcowe powoduje powstawanie dużych kontrastów pomiędzy Vp i Vs, związanych z obecnością stref saturacji. Mieszanina piasków, pospółek i materiału ilastego oraz pośrednich form osadowych jest szczególnie podatna na infiltrację wody. Zawodnienie ma znaczący wpływ na określenie zależności pomiędzy Vp i Vs, jednak nie jest to jedyny warunek określający ich wzajemną relację (Zimmer 2003). W sejsmice inżynierskiej pole prędkości fali poprzecznej i podłużnej uzyskuje się na podstawie kilku wariantów powierzchniowych (Gabriels 1987, Foti 2000, 2003, 2004, Hunter 2002, Paillet 1995) i profilowań otworowych (Hunter 2002, Stephenson 2005, Paillet 1995, Wair 2012, Xia 2002). Bazują one na konkretnych zjawiskach falowych, często różniących się sposobem propagacji, rozdzielczością i zasięgiem penetracji strefy przypowierzchniowej. Związek pomiędzy wspomnianymi stałymi zależy od ilościowej zbieżności rozkładów prędkości fali P i S (Bergamo 2014; Bodet 2014; Castagana 1985; Gibbs 1999, 2000, 2001; Huffman 2001; Zimmer 2003). Trudności w znalezieniu odpowiedniej formuły pozwalającej na ilościowe przeliczenie wartości prędkości Vs i Vp wynikają z niedoskonałości postaci otrzymanych pól prędkości oraz ich różnic strukturalnych. Literatura poświęcona obrazowaniu sedymentów polodowcowych przedstawia rozkłady stosunku Vp/Vs w postaci dwuwymiarowych map głębokościowych. Nieliczne publikacje skupiają się określeniu związku ilościowego pomiędzy Vp i Vs w postaci funkcji liniowych. Szczególnie cenna pod tym względem okazała się publikacja autorstwa Salem’a i in. 2000, który na zasadzie eksperymentu sejsmicznego podjął próbę znalezienia zależności pomiędzy Vp i Vs w utworach polodowcowych. Ich postać jest następująca:
Vp/Vs= 1.948+0.00205*Vs (4.3-5)
Vp/Vs= 1.6113+0.0007*Vp (4.3-6)
Funkcje liniowe podane przez Salem’a (4.3 - 5) i (4.3 - 6) dotyczą zarówno refraktora jak i nadkładu. Można je traktować jako uśrednienie wyrażające związek pomiędzy Vp/Vs(Vs) i Vp/Vs(Vp) w całym zakresie penetracji fali refrakcyjnej. Rysunek 4.3 - 2 przedstawia wykres zbiorczy zależności pomiędzy wartością prędkości podłużnej i poprzecznej w utworach polodowcowych. Czerwonymi okręgami o zmiennym promieniu określono ilość wartości Vp i Vs na podstawie publikacji (Cox 2009). Niebieski wykres powstał na bazie wyników profilowań otworowych fali podłużnej (P-P, Rys.4.1 - 2) i
54 poprzecznej (SV, Rys. 4.1 - 3), na którym średnice okręgów wzrastają wraz z wartością prędkości.
Rys. 4.3 - 2 Zbiorczy wykres zależności pomiędzy wartościami prędkości fali podłużnej i poprzecznej.
Wykres 4.3 - 2 podzielono na trzy strefy A, B i C. Zakres A odpowiada zakresowi prędkości fali podłużnej i poprzecznej w strefie aeracji. Obszar B stanowi strefę przejściową pomiędzy strefą aeracji i saturacji. Poligon C określa strefę saturacji. W każdym z obszarów zaznaczono barwnymi krzyżami punkty reprezentatywne. Na ich podstawie dokonano analizy regresji liniowej. Wyprowadzono trzy wzory określające zależność pomiędzy Vp i Vs Rys. 4.3 – 2 - (1),Rys. 4.3-2- (2) oraz Rys.4.3-2 - (3):
55 Vp= 1.98*Vs+60 (4.3-7)
Vp=Vs*1.93+650 (4.3-8)
Vp=Vs*2.09+1191 (4.3-9)
Stosunek prędkości fali podłużnej do prędkości fali poprzecznej wiąże w sposób bezpośredni wartości prędkości Vp i Vs. Jest wskaźnikiem zawodnienia w nieskonsolidowanych utworach przypowierzchniowych. Pozwala na określnie strefy silnego przyrostu wartości Vp. Rozbieżności w wartościach prędkości Vp i Vs, przekładają się na zakres zmienności współczynnika Vp/Vs. W gruntach słabo zwięzłych oscyluje on w przedziale 1,4 - 8,9, co podkreśla fakt złożoności relacji pomiędzy wartościami Vp i Vs w utworach nieskonsolidowanych. W celu potwierdzenia określonych wzorów regresyjnych na rysunku 4.3 - 3 przeprowadzono konstrukcje wykresu zmienności wartości Vp/Vs w przedziale wartości prędkości 100 - 2500 m/s. Za krzywe referencyjne obrano funkcje podane przez Salem i in. (Rys. 4.3 - 3). Na jego podstawie analizowano przyrosty wartości Vp/Vs w odniesieniu do dwóch pól prędkości: Vp i Vs. Barwnymi wykresami kołowymi (czerwonym i zielonym) oznaczono typowe przedziały wartości Vp/Vs(Vp) oraz Vp/Vs(Vs). Na wykres naniesiono wartości funkcji Salem’a oznaczone przez „c” i „d”. Poniżej wartości prędkości 500 m/s (w strefie aeracji) zależność wartości fali poprzecznej i podłużnej można przybliżyć przez stałą wartość Vp/Vs =2. Iloraz prędkości fali podłużnej do poprzecznej wobec zmiany wartości prędkości Vs (wykres zielony) przyrasta dwukrotnie szybciej w porównaniu do rozkładu Vp/Vs(Vp). Jest to efekt spadku wartości prędkości fali S przy równoczesnym wzroście prędkości fali P. Poniżej zwierciadła wód gruntowych (Vp> 1500 m/s) przeliczenie wartości Vs na Vp należy przeprowadzić przy Vp/Vs ≥ 3. Przyjęcie wspomnianych stałych przeliczeniowych jest właściwe w przypadku niezaburzonego wzrostu prędkości Vp i Vs wraz z głębokością. W sedymentach nieskonsolidowanym prędkość fali S przyrasta znacznie wolniej w porównaniu do prędkości fali P wraz z głębokością. Wymaga to wprowadzenia poprawki na skok wartości prędkości Vp względem Vs w strefie zawodnionej. W tym celu wprowadzono dwie proste „a” i „b” (Rys. 4.3 - 3). Wskazują one przesunięcie ostałą wartość przyrostu Vp/Vs(Vp) względem Vp/Vs(Vs), równą 1320 m/s . Wynika stąd, że dla szerokiego przedziału stosunku Vp/Vs, istnieje możliwość oszacowania wartości :
56 Rys. 4.3 - 3 Wykres zależności pomiędzy stosunkiem Vp/Vs względem wzrastającej wartości
prędkości.
Typowa wartość prędkości fali powierzchniowej w utworach nieskonsolidowanych wynosi ok. 300 m/s (środek strefy oznaczonej kwadratem, Rys. 4.3 - 3). Na jej przykładzie przedstawiono drogę transformacji stałej wartości prędkości fali S do wartości fali P w strefie zawodnionej. Na wykres 4.3 - 3 wprowadzono dane otworowe (Rys. 4.1 - 2 i 4.1 - 3). Obszar estymacji zawężono do obszaru ograniczonego przez funkcje 4.3 - 5 i 4.3 - 6. Wykazano, że wyniki profilowań otworowych są zbieżne z literaturową populacją punktów w ograniczonym
57 zakresie prędkości. Różnią się one znacznie w zakresie niskich prędkości od rozkładów zaczerpniętych z literatury. Wskazuje to na szczególnie złożoną charakterystykę prędkościową strefy aeracji. Część ich wartości uwzględniono w analizie regresji. Wyprowadzenie wzoru przejścia pomiędzy Vs i Vp przedstawiono w postaci odcinka A - B, który zlokalizowano na krzywej regresji liniowej „w”. Powstał ona na bazie dwóch chmur punktów związanych ze strefa aeracji (przedział Vp: 100-500 m/s) i zawodnioną ( Vp > 1500 m/s). Charakteryzuje go wzór regresyjny:
Vp/Vs= 0.00083*Vs+ 2.05 (4.3-11)
Który można przekształcić do postaci:
Vp/Vs=1204.1/Vs+2.05 (4.3-12)
Na podstawie którego relację pomiędzy Vp i Vs można określić jako:
Vp=2.05*Vs+ 1204.1 (4.3-13)
Do oszacowania ostatecznej postaci zależności pomiędzy Vp i Vs posłużono się zależnościami 4.3 - 7, 4.3 - 9 oraz 7.3 - 13. Pomięto wzór 4.3 - 8 ze względu na brak pokrycia jego wartości z danymi otworowymi. Uśredniając wartości współczynnika kierunkowego oraz wyrazu wolnego wspomnianych funkcji liniowych, otrzymano ostateczna formułę przeliczenia wartości Vs na Vp:
Vp=2.07*Vs+ 1197.4 (4.3-14)
Uwzględniając wyrażenie 4.3-1, otrzymano:
𝛼𝑝𝑜𝑐𝑧=2.07*0.86*𝐶𝑅,+ 1197.4 (4.3-15) Ze względu na ilość posiadanych rekordów, kontrolowano proces inwersji w wariancie MASW przez zachowanie jednolitego schematu parametryzacji. Zdefiniowano go iteracyjnie. Jako model referencyjny wykorzystano otworowy rozkład prędkości SV (Rys. 4.1 - 3). Rekordy znajdujące się w bliskim otoczeniu poddano kilkukrotnej analizie, każdorazowo zmieniając podejście do pikowania i ilościowej definicji stałych kontrolujących proces inwersji. Wynik najbardziej zbliżony do modelu otworowego uznano za końcowy, a schemat powielono w przetwarzaniu pozostałych rekordów sejsmicznych. Wszystkie wyniki inwersyjne 1 - D zostały uzyskane przy zastosowaniu następujących założeń:
58 1. Miąższość oraz ilość warstw uzależniono od zakresu spektralnego krzywej dyspersji
(wzór 4.3-2).
2. Ilość iteracji w procesie inwersji dla krzywych gładkich (mono- modalnych) określono jako 5 , natomiast w przypadku krzywych wielo - modalnych była równa 10.
3. Wartości startowe prędkości fali poprzecznej zdefiniowano jako 𝛽𝑝𝑜𝑐𝑧= 0.86*𝐶𝑅 wartości prędkości fazowych, pozyskanych w czasie pikowania krzywych dyspersji.
4. Rozkład początkowy wartości prędkości fali podłużnej w strefie aeracji ustalono na 𝛼𝑝𝑜𝑐𝑧=2*0.86*𝐶𝑅, natomiast w strefie saturacji 𝛼𝑝𝑜𝑐𝑧= 2.07*0.86*𝐶𝑅,+ 1197.4.
5. Głębokość poziomu wodonośnego określono na podstawie profilowań otworowych prędkości fali P - punkt A na rysunku 4.1 - 2 (podrozdział 4.1).
6. Pozorną gęstość ośrodka zdefiniowano wykorzystując wartości Vp i wzór podany przez Meyne 1999 (4.3-3).
7. Nie wprowadzono elementów stałych inwersji w wartościach rozkładu początkowego fali P, dopuszczono możliwą inwersje prędkości fali S oraz w celu podkreślenia detali rozkładów wykorzystano surowy (ang. robust’owy) typ inwersji.
8. Wartości miąższości warstw h i gęstości pozorne ρ oznaczono jako stałe.
9. Rezultaty eksportowano do postaci plików tekstowych (ASCII). Na ich podstawie założono siatkę interpolacyjną - grid, którego cele w poziomie zachowały wielkość interwału odbiorników, pionowo były równe miąższości warstw po inwersji.
59 Rys. 4.3 - 4 Rozkład średniego błędu kwadratowego względem długości profilu badawczego. Przyjęty schemat przetwarzania pozwolił na dopasowanie krzywych teoretycznych do polowych z wartościami średnich błędów kwadratowych (wzór 4.1-22) przedstawionymi na rysunku 4.3 - 4. Każda z wielkości jest średnią dwudziestu krzywych dyspersji, odnoszących się do 20 metrowego odcinka profilu. Zgodnie z założoną geometrią pomiarową, obszary stanowiły samodzielne obszary estymacyjne. Jego zmienność zawiera się w przedziale 3-22%. Względna zmienność pomiędzy na mniejszym i największym błędem dopasowania wynosi 19 %. Poszukiwano przyczyn zróżnicowania w rozkładzie błędu RMS. W tym celu wykonano wykres zależności pomiędzy długością profilu i zakresem spektralnym krzywej dyspersji (Rys. 4.3 - 5). Odniesiono go do wartości błędu RMS. Wykres podzielono na dwie strefy: określając dominację mody podstawowej (< 30 Hz) oraz obecność mód wyższych w zakresie (30 - 45 Hz).
60 Rys. 4.3 - 5 Wykres zależności pomiędzy długością profilu i zakresem spektralnym krzywych
dyspersji.
Największe wartości błędów dopasowania znajdują się w strefie dominacji mody podstawowej. Po przeglądzie krzywych dyspersji, zjawisko to utożsamiono ze fragmentaryczną formą trendu dyspersyjnego mody podstawowej. Wykres wskazuje, że błąd RMS może osiągać niewielkie wartości w czasie inwersji wielomodalnej. Warunkiem jest czytelność mód wyższych w dziedzinie prędkość fazowa - częstotliwość. Modalna postać fali powierzchniowej dominuje w utworach polodowcowych, ponieważ większość wartości błędu została przypisana do przedziału 30 - 45 Hz (65% wyników inwersyjnych). Na rysunku 4.3 - 5 linią przerywaną zaznaczono strefę kumulacji wartości błędu RMS.
61 Rys. 4.3 - 6 Wykres rozdzielczości metody MASW wyrażony w metrach, zestawiony z
krzywą zasięgu oraz dyspersji.
Na potrzeby kolejnych etapów przetwarzania obliczono rozdzielczość wariantu wielokanałowej analizy fal powierzchniowych. Przedstawiono ją w postaci wykresu 4.3 - 6. Linią niebieską oznaczono zasięg metody, natomiast czerwoną rozdzielczość wyrażona w metrach. Krzywa czarna przedstawia postać krzywej dyspersji, którą uznano za wzorcową. W strefie do 10 metrów głębokości MASW posiada wysoką rozdzielczość. Średnio w zakresie od 1 do 10 metrów p. p. g., metoda MASW jest zdolna wykryć anomalie prędkości fali poprzecznej o rozmiarze 1.5 metra, poniżej 15 metrów głębokości dochodzi do siedmiokrotnego spadku możliwości detekcyjnych wariantu MASW.
62