Index of /rozprawy2/10107

Pełen tekst

(1)AKADEMIA GÓRNICZO-HUTNICZA IM. STANISŁAWA STASZICA WYDZIAŁ GEOLOGII, GEOFIZYKI I OCHRONY ŚRODOWISKA KATEDRA GEOFIZYKI. ROZPRAWA DOKTORSKA. ANALIZA ZMIAN ZAPISU SEJSMICZNEGO Z OFFSETEM JAKO NARZĘDZIE DO IDENTYFIKACJI AKUMULACJI GAZU ZIEMNEGO NA PRZYKŁADACH Z ZAPADLISKA PRZEDKARPACKIEGO. mgr inż. Artur Tatarata. PROMOTOR prof. dr hab. inż. Kaja Pietsch. KRAKÓW 2009.

(2) PODZIĘKOWANIA Chciałbym przede wszystkim serdecznie podziękować Pani Profesor Kai Pietsch, pod której opieką moja praca doktorska powstawała, za okazaną pomoc oraz cenne wskazówki udzielone podczas przygotowywania pracy. Dziękuję Jej również za wieloletnią opiekę naukową. Dziękuję także moim kolegom z Pracowni Interpretacji Danych Sejsmicznych, Katedry Geofizyki WGGiOŚ AGH – dr inż. Marcinowi Kobylarskiemu, mgr inż. Pawłowi Marcowi, mgr inż. Ryszardowi Hodiakowi oraz dr inż. Edycie Frankowicz za współpracę i owocne dyskusje pomocne przy realizacji moich badań.. II.

(3) Spis treści 1 Wstęp..........................................................................................................................................................................................................2 2 Podstawy teoretyczne zjawiska zmian amplitud z offsetem.............................................................................................4 2.1 Rozchodzenie się fal sprężystych na granicy izotropowych ośrodków...............................................................4 2.2 Petrofizyczne uwarunkowania zróżnicowania prędkości fal podłużnych i poprzecznych oraz gęstości w ośrodku skalnym.....................................................................................................................................................................6 2.3 Klasyfikacja zmian amplitud z offsetem ..........................................................................................................................16 2.4 Ogólna charakterystyka metod wykorzystania rejestracji sejsmicznych w analizie AVO ........................22 3 Metody analiz zmienności amplitud z offsetem w kontekście identyfikacji akumulacji gazu.........................28 3.1 Atrybuty intercept, gradient, curvature i ich kombinacje.........................................................................................29 3.1.1 Ekstrakcja atrybutów A, B i C z kolekcji kątowej WPG ..................................................................................29 3.1.2 Interpretacja atrybutów A, B i C...............................................................................................................................42 3.2 Fluid Factor ...................................................................................................................................................................................47 3.3 Ocena wiarygodności obliczeń atrybutów AVO ...........................................................................................................52 3.4 Impedancja elastyczna.............................................................................................................................................................57 3.5 Stałe materiałowe z AVO ........................................................................................................................................................66 3.6 Impedancja Poissona................................................................................................................................................................71 3.7 Atrybuty polaryzacyjne...........................................................................................................................................................73 3.8 Podsumowanie............................................................................................................................................................................79 4 Interpretacja złożowa danych sejsmicznych przed składaniem z obszaru zapadliska przedkarpackiego.....................................................................................................................................82 4.1 Szkic geologiczny wybranych obszarów zapadliska przedkarpackiego ...........................................................83 4.2 Wykorzystane dane geofizyczne i ich przetwarzanie................................................................................................86 4.2.1 Zdjęcie sejsmiczne Lubliniec-Aleksandrów-Cieszanów 2D ..........................................................................86 4.2.2 Zdjęcie sejmiczne Sokołów Małopolski-Smolarzyny 3D (obiekt Pogwizdów).....................................90 4.2.3 Zdjęcie sejsmiczne Ujkowice-Batycze 3D .............................................................................................................94 4.3 Interpretacja wybranych fragmentów zdjęć sejsmicznych.....................................................................................99 4.3.1 Lubliniec-Aleksandrów-Cieszanów 2D ..................................................................................................................99 4.3.2 Sokołów Małopolski – Smolarzyny 3D................................................................................................................ 121 4.3.3 Ujkowice–Batycze 3D.................................................................................................................................................. 134 4.4 Interpretacja złożowa – podsumowanie ...................................................................................................................... 144 5 Wnioski................................................................................................................................................................................................ 147 Literatura................................................................................................................................................................................................ 149.

(4) 1 Wstęp Sejsmika jest główną metodą geofizyczną wykorzystywaną w poszukiwaniach złóż ropy naftowej i gazu ziemnego. Pozwala ona na obrazowanie budowy geologicznej, co ułatwia pośrednie wnioskowanie o obecności akumulacji węglowodorów. Od początku lat siedemdziesiątych, na skutek wzrostu jakości rejestracji, możliwe stały się bezpośrednie poszukiwania węglowodorów. Jednym z mierzalnych efektów wskazujących na obecność węglowodorów są zmiany wielkości amplitud zapisu sejsmicznego w funkcji odległości źródło-odbiornik (offset). Na obszarze zapadliska przedkarpackiego akumulacje gazu ziemnego nawiercane są w coraz głębszych strefach, bardzo często w utworach cienkowarstwowych, co powoduje, że dużego znaczenia nabiera poprawna lokalizacja i charakterystyka stref anomalnych zapisu sejsmicznego, wywołanych nasyceniem gazem ziemnym. Duża zmienność parametrów petrofizycznych utworów miocenu, w których zlokalizowane są pułapki powoduje, że do złożowej interpretacji zapisu sejsmicznego powinny być wykorzystane wszelkie dostępne dane sejsmiczne, nie tylko tradycyjnie interpretowane przekroje sejsmiczne, ale również rzadziej wykorzystywane dane sejsmiczne przed składaniem. Celem pracy jest zaprezentowanie możliwości interpretacji danych sejsmicznych przed składaniem w kontekście bezpośrednich poszukiwań gazu, opierając się na przykładach z zapadliska przedkarpackiego. Praca składa się z trzech głównych części. W rozdziale 2 przedstawiono zagadnienia teoretyczne odnośnie wpływu parametrów petrofizycznych charakteryzujących ośrodek geologiczny na rejestrowane zmiany amplitud zapisu sejsmicznego w funkcji offsetu. Przedstawiono również w zarysie wpływ akwizycji i przetwarzania danych sejsmicznych na wynikowy zapis amplitudowy. W rozdziale 3 zebrano dostępne techniki interpretacji danych sejsmicznych przed składaniem. Opierając się na założonych modelach sejsmogeologicznych, sparametryzowanych na warunki z zapadliska przedkarpackiego, przeprowadzono symulacje ekstrakcji informacji o ośrodku geologicznym. Przeanalizowano uwarunkowania testowanych metod interpretacyjnych wynikające z charakteru analizowanego zapisu sejsmicznego oraz parametrów ośrodka. Ważnym aspektem badań było przedstawienie wpływu sposobu wykonywania analiz na finalny wynik. Rozdział 4 przedstawia interpretację wybranych przykładów danych sejsmicznych z obszaru zapadliska przedkarpackiego. Przeanalizowano przydatność testowanych atrybutów bazujących na zmienności amplitudy z offsetem w kontekście bezpośrednich poszukiwań akumulacji gazu ziemnego. W zakończeniu przedstawiono wnioski dotyczące zarówno metodyki wyznaczania atrybutów jak i ich przydatności do interpretacji złożowej w warunkach zapadliska przedkarpackiego.. 2.

(5) Do realizacji pracy wykorzystano system Hampson-Russell oraz język obliczeń technicznych MATLAB. Wykorzystane dane sejsmiczne i geofizyki otworowej uzyskano dzięki uprzejmości PGNiG S.A. oraz spółek: Geofizyka Kraków Sp. z o.o. i Geofizyka Toruń Sp. z o.o.. 3.

(6) 2 Podstawy teoretyczne zjawiska zmian amplitud z offsetem 2.1 Rozchodzenie się fal sprężystych na granicy izotropowych ośrodków. W sejsmice poszukiwawczej do badania ośrodka skalnego wykorzystujemy zjawisko rozchodzenia się dwóch głównych typów fal sprężystych: fal podłużnych – P i fal poprzecznych – S. Ruch fal P polega na rozchodzeniu się zagęszczeń i rozrzedzeń (zmian objętości) ośrodka w kierunku równoległym do kierunku rozchodzenia się czoła fali. Fale P rozchodzą się we wszystkich typach ośrodków tj. ciałach stałych, cieczach i gazach. Ruch fal S polega na rozchodzeniu się zmian kształtu ośrodka, dla którego ruch cząstek odbywa się w płaszczyznach prostopadłych do kierunku rozchodzenia się fali. Powoduje to, że fale S nie rozchodzą się w ciałach pozbawionych sprężystości postaci: cieczach i gazach. W przypadku jednorodnego i izotropowego ośrodka do opisu jego własności sprężystych i określenia prędkości propagacji fal P i S wystarczy znajomość gęstości oraz dwóch modułów sprężystości tj. modułu odkształcenia objętości – K lub stałej Lamégo – λ oraz modułu odkształcenia postaci – µ (Plewa, 1992). Prędkości fal P i S wyrażają się wtedy wzorami:. VP =. λ + 2µ = ρ VS =. K+. 4 µ 3. ρ. µ ρ. (2.1.1) (2.1.2). Do opisu własności sprężystych często również używa się współczynnika Poissona – σ, będącego stosunkiem względnego odkształcenia poprzecznego do względnego odkształcenia wzdłużnego. Współczynnik Poissona wiąże relację prędkości fal P i S i wyraża się wzorem:. σ=. 1 (VP2 - 2VS2 ) 2 (VP2 - VS2 ). (2.1.3). Fundamentalnymi zjawiskami będącymi podstawą metod sejsmicznych są odbicie i załamanie padających fal na granicy rozdziału dwóch ośrodków skalnych o różnych własnościach sprężystych. Na granicy dochodzi również do konwersji części energii padającej fali podłużnej na energię fal poprzecznych: załamanej i odbitej (rys. 2.1.1). Wzajemne relacje pomiędzy padaniem kątami padania i załamania poszczególnych typów fal opisuje prawo Snelliusa: sin θ1 sin θ 2 sin Φ1 sin Φ 2 = = = = pp (2.1.4) VP1 VP 2 VS 1 VS 2 gdzie: p p – parametr promienia sejsmicznego, objaśnienia – patrz rys. 2.1.1.. 4.

(7) Rys. 2.1.1 Typy fal powstających przy padaniu fali podłużnej (pP) na granicę rozdziału dwóch ośrodków. Fale odbite: podłużna (oPP) i poprzeczna (oPS), fale załamane: podłużna (zPP) i poprzeczna (zPS). Kąt padania fali pP (θ1), kąt załamania fali zPP (θ2), kąt odbicia fali oPS (Φ1), kąt załamania fali zPS (Φ2).. Zjawisko konwersji energii na granicy pociąga za sobą zmiany amplitud poszczególnych typów fal, co stanowi podstawę analiz amplitud z offsetem (amplitude versus offset – AVO). W przypadku ośrodka jednorodnego, izotropowego i padającej płaskiej fali podłużnej amplitudy fal odbitych: podłużnej – AoPP i poprzecznej – AoPS oraz załamanych: podłużnej – AzPP i poprzecznej – AzPS opisywane są równaniami Knotta– Zoeppritza (Knot (1899) i Zoeppritz (1919) za Sheriff i Geldhart (1995)). Równania te wiążą wspomniane amplitudy z własnościami ośrodka górnego (VP1, VS1, ρ1) i dolnego (VP2, VS2, ρ2) oraz kątami: padania fali P – θ1, załamania fali P – θ2, załamania fali S – Φ2 i odbicia fali S – Φ1. Wygodną macierzową formę równań przedstawili Aki i Richards (1980):  sin θ1  − cos θ 1   sin 2θ 1   cos 2Φ1 . cos Φ1 sin Φ1 VP1 cos 2Φ1 VS 1 V − S 1 sin 2Φ1 VP1. − sin θ 2 − cos θ 2 ρ 2VS22VP1 sin 2θ 2 ρ!VS21VP 2 ρV − 2 P 2 cos 2Φ 2 ρ1VP1. cos Φ 2    AoPP   − sin θ1  − sin Φ 2  AoPS   − cos θ1  ρ 2VS 2VP1 (2.1.5) = − cos 2Φ 2    AzPP   sin 2θ1  ρ1VS 1     ρV − 2 S 2 sin 2Φ 2   AzPS  − cos 2Φ1   ρ1VP1. 5.

(8) 2.2 Petrofizyczne uwarunkowania zróżnicowania prędkości fal podłużnych i poprzecznych oraz gęstości w ośrodku skalnym. Jednoznaczne powiązanie prędkości fal podłużnych VP i poprzecznych VS oraz gęstości objętościowej ρ z daną facją geologiczną jest bardzo trudne z powodu wpływu na nie wielu czynników natury petrofizycznej i litologicznej. Poniżej przedstawiono w zarysie wpływ na VP, VS i ρ poszczególnych czynników, prezentując także spotykane w literaturze modele teoretyczne i empiryczne łączące prędkości fal oraz gęstości w ośrodku skalnym z jego parametrami takimi jak nasycenie gazem, porowatość i zailenie. Prędkości rozchodzenia się fal sprężystych oraz gęstości matrycy skalnej Proporcje poszczególnych minerałów w matrycy skalnej warunkują prędkości fal sprężystych. Wartość średnią prędkości obliczamy z relacji: n. Vma = ∑ AiVi. (2.2.1). i =1. gdzie: Ai – objętościowa zawartość danego minerału, Vi – prędkość dla danego minerału. Uśrednione prędkości i gęstości poszczególnych matryc skalnych zebrane z różnych publikacji zestawiono w tabeli 2.2.1. Tab. 2.2.1 Prędkości fal: P i S, stosunek VP/VS, współczynnik Poissona σ, moduły sprężystości: K, λ, μ oraz gęstość objętościowa ρ dla wybranych minerałów skałotwórczych – uśrednione (zestawione za Bała, 2006). VP. VS. K. λ. μ. ρ. m/s. m/s. GPa. GPa. GPa. g/cm3. kwarc. 5780. 3876. 1,50. 0,09. 35,4. 8,9. 39,8. 2,65. ił*. 3175. 1645. 1,93. 0,32. 17,2. 12,4. 7,2. 2,65. anhydryt. 6096. 3280. 1,86. 0,30. 67,8. 46,5. 32,0. 2,97. klacyt. 6410. 3322. 1,93. 0,32. 71,5. 51,5. 29,9. 2,71. dolomit. 6993. 3831. 1,82. 0,29. 84,2. 56,1. 42,1. 2,87. halit. 4608. 2628. 1,75. 0,26. 25,6. 15,8. 14,7. 2,13. minerał. VP/VS. σśr. * – zailenie często uwzględnia się oddzielnie i wprowadza do obliczeń poza matrycą skalną. Prędkości rozchodzenia się fal sprężystych oraz gęstości mediów Wartości prędkości i gęstości mediów nasycających pory skalne, zależą od składu chemicznego w przypadku węglowodorów, warunków termobarycznych, zasolenia wody złożowej, sposobu wypełnienia przestrzeni porowej przez różne media. 6.

(9) (m.in. Wang i Nur, 1986, Batzle i Wang, 1992 za Bała, 2006). Uśrednione prędkości i gęstości poszczególnych mediów zebrane z różnych publikacji zestawiono w tabeli 2.2.2. Tab. 2.2.2 Prędkość fali P i gęstość objętościowa dla wybranych mediów – uśrednione (zestawione za Bała, 2006). VP ρ medium m/s g/cm3 Solanka. 1600. 1,03. Ropa. 1350. 0,85. Gaz. 707. 0,10. Modele teoretyczne łączące parametry sprężyste i zbiornikowe skał porowatych zestawione za Bała i Cichy (2006) Najprostszym modelem jest model Wylliego (Wyllie et al., 1956), który uwzględnia w opisie ośrodka dwie fazy: matrycę i medium nasycające pory. Prędkość fali P w skale VP wyraża się wzorem: 1 Φ 1- Φ = + (2.2.2) V V f Vma gdzie: Vma – prędkość fali P w matrycy, Vf – prędkość fali P w medium porowym, Φ – porowatość. Przy jednorodnym rozkładzie porów i braku zailenia może model ten być stosowany do określenia porowatości w zakresie 0% do 30% z pomiarów geofizyki otworowej. Stosuje się również poprawki na zailenie (zailenia nie wlicza się do matrycy skalnej), obecność węglowodorów oraz zwięzłość skały (Jarzyna et al., 1999). Model empiryczny Raymera-Hunta-Gardnera (1980) wyraża prędkość fali w skale wzorem:. V = (1 − Φ ) 2 Vma + ΦV f. (2.2.3). gdzie: oznaczenia jak we wzorze Williego (2.2.2) Może być stosowany dla skał o porowatości międzyziarnowej z przedziału 0% do 37%. W obliczeniach uwzględnia się również obecność węglowodorów i skład mineralny skały. Połączenie modeli zaproponowanych przez Biota (1956) i Gassmana (1951) zaproponował Krief et al. (1989). Model Biota-Gassmana wiąże moduły Ksat i μsat (formacji skalnej) z porowatością ośrodka Φ oraz modułami szkieletu KS i μS, matrycy skalnej w stanie suchym Kdr i μdr oraz medium złożowego Kf. Moduły sprężystości formacji skalnej wyrażają się wzorami: K sat = K S (1 − β ) + β 2 M. (2.2.4). 7.

(10) µ sat = µ S (1 − β ). (2.2.5). 1 β −Φ Φ = + M KS Kf. (2.2.6). gdzie: β – współczynnik Biota (zależny od porowatości) Prędkości fali P i S wyrażają się wzorami: VP2 =. 1. ρ. [ ρ sVPs2 (1 − β ) + Mβ 2 ]. VS2 =. ρs 2 VS (1 − β ) ρ. (2.2.7) (2.2.8). Model BG nie uwzględnia sposobu rozłożenia wolnych przestrzeni w skale. Kształt porów uwzględnia natomiast model Kustera-Toksöza (1974), w którym niejednorodność ośrodka została opisana przy użyciu rozkładu statystycznego współczynnika kształtu porów (αm). W modelu ośrodek geologiczny jest aproksymowany poprzez dwie jednorodne fazy: pierwszą stałą odpowiadającą matrycy skalnej (Kma, µma, ρma) i drugą fazę medium wypełniające pory (Kf, µf, ρf). Rozkład fazy porowatej jest opisany za pomocą rozkładu statystycznego poprzez współczynnik kształtu porów. Moduły sprężystości dla efektywnego ośrodka porowatego K* i µ* są dość skomplikowanymi funkcjami modułów matrycy i medium oraz współczynnika kształtu porów (szerzej patrz: Wu (1966), Toksöz et al. (1976) za Bała i Cichy (2006)). Prędkości fal P i S można wyliczyć ze wzorów 2.1.1 i 2.1.2 wstawiając do nich K* i µ*. Powyższa teoria zakłada brak połączeń między porami. Istnieją jednak modyfikacje teorii KT pozwalające na oddziaływanie porów ze sobą (Cheng i Toksöz (1979) za Bała i Cichy (2006)). Prace Bały (1989, 1994) uzupełniają model o możliwość uwzględniania rozłożenia w przestrzeni porowej dwóch faz: wody i gazu. Powyższe 4 modele tj. W, RHG, BG i KT zostały zaimplementowane w programie ESTYMACJA (opracowanym w Katedrze Geofizyki AGH), który posłużył do wykonania obliczeń przy realizacji pracy. Wpływ nasycenia przestrzeni porowej gazem na prędkość rozchodzenia się fal sprężystych i gęstość objętościową W bezpośrednich poszukiwaniach sejsmicznych bardzo często odwołujemy się nie do bezwzględnych wartości prędkości fal, gęstości objętościowej czy modułów sprężystości, ale do wzajemnych różnic tych parametrów pomiędzy skałą nasyconą i nienasyconą. Względną zmianę wartości Δx/x danego parametru X przyjęło się określać jako stosunek różnicy jego wartości w warstwie dolnej x2 i górnej x1 do średniej arytmetycznej obu wartości: x2 − x1 ∆x = (2.2.9) x 0,5( x1 + x2 ). 8.

(11) Przyjmuje się, że zmiana wartości danego parametru ma znak dodatni, gdy jego wartość rośnie przy przejściu z warstwy górnej do warstwy dolnej. Przedstawiony sposób opisywania względnych zmian jest używany w pracy. Przedstawiono (za Bała, Cichy 2006) porównanie prędkości fali podłużnej obliczonych dla matrycy „kwarcowej” o porowatości 25% dla wyżej omówionych modeli w zależności od nasycenia przestrzeni porowej gazem (2.2.1). Prędkości obliczone z modelu W wyraźnie maleją wraz ze wzrostem nasycenia, natomiast spadek względna zmiana prędkości ΔVP/VP przy wzroście nasycenia od 0 do 1 wynosi prawie -1,1, co jest ekstremalnie dużą wartością. Model RHG również zakłada spadek prędkości przy wzroście nasycenia jest on jednak największy przy wzroście nasycenia od 0 do 0,1 i zmiana wynosi -0,11, aby przy rosnących nasyceniach od 0,1 do 1 poziom prędkości zmieniał się już nieznacznie. Model BG zakłada spadek prędkości do wartości minimalnej dla nasycenia około 0,2 (zmiana -0,05) następnie przy zwiększającym się nasyceniu następuje niewielki jednostajny wzrost prędkości. Model KT zakłada spadek prędkości fali przy przejściu od nasycenia wodą do pełnego nasycenia gazem (zmiana wynosi -0,24), przy czym zmiana -0,12 zaznacza się już przy nasyceniu gazem rzędu 0,1. W modelu KT założono wariant, w którym gaz i woda nie mieszają się ze sobą. Generalizując, powyższe modelowania wyraźnie pokazują, że prędkość fali podłużnej nie zmienia się jednostajnie ze wzrostem nasycenia gazem. Największy spadek pojawia się już przy nasyceniu rzędu 0,1-0,2.. Rys. 2.2.1 Porównanie prędkości fali podłużnej – VP od nasycenia umownego piaskowca Sg (patrz tekst) gazem obliczonych według różnych modeli (zmodyfikowane za Bała, Cichy 2006). Prędkość fali S (rys. 2.2.2) wyliczoną zgodnie z modelem BG dla wybranego piaskowca o porowatości 0,25 przy zamianie nasycenia z wody na gaz niewiele rośnie od 2360 do 2510 m/s. Względna zmiana prędkości wynosi 0,06 (Bała, 2007).. 9.

(12) Rys.2.2.2 Zależność VS2–VP2 dla piaskowca o porowatości Φ=0,25 w funkcji nasycenia wodą – SW (linia czarna). Zaznaczono linie odpowiadające całkowitym nasyceniom wodą (niebieska) i gazem (czerwona) dla różnych porowatości (zaczerpnięte z Bała, 2007).. Gęstość objętościowa ρ skały zmniejsza się liniowo wraz ze wzrostem nasycenia, co wynika z relacji (rys. 2.2.3): ρ = ρ ma (1 − Φ ) + ρ f Φ (2.2.10) gdzie: ρma – gęstość matrycy, ρf – gęstość medium, Φ - porowatość Względna wielkość zmiany gęstości ∆ρ/ρ jest zależna od objętości wolnych przestrzeni w skale (porów, szczelin) i jest rzędu -0,05 przy porowatości 5% i –0,4 przy porowatości 40%. Znajomość gęstości objętościowej mogłaby być potencjalnym wskaźnikiem nasycenia skały gazem. Rho (5%) Rho (20%) Rho (35%). 2.9 2.7. Rho (10%) Rho (25%). Rho (15%) Rho (30%). ρ [g/cm3]. 2.5 2.3 2.1 1.9 1.7 1.5 0. 0.2. 0.4. Sg. 0.6. 0.8. 1. Rys. 2.2.3 Zależność gęstości objętościowej skały – ρ od nasycenia wolnych przestrzeni gazem – Sg. Gęstość matrycy skalnej założono na 2,6 g/cm3, wody złożowej na 1,03 g/cm3 i gazu na 0,1 g/cm3. Parametr krzywych - objętość wolnych przestrzeni w skale.. 10.

(13) Wpływ rozłożenia mediów złożowych w przestrzeni porowej na prędkość rozchodzenia się fal Prędkość rozchodzenia się fal w skale zależy także od sposobu ułożenia mediów w przestrzeni porowej skały. Na rys. 2.2.4 przedstawiono trzy warianty sposobu wzajemnego rozłożenia wody i gazu w piaskowcu o porowatości 0,2 i współczynnikach kształtów porów α=0,2 według model KT (za Bała, 2007): • mieszanina - wariant 1 •. media niezmieszane (np. pęcherzyki gazu w wodzie (Han et al., 2002) – wariant 2. •. media występują w oddzielnych porach – wariant 3. W przypadku występowania porów o jednakowym współczynniku kształtu dla wariantów 1 i 3 obserwujemy liniowy spadek prędkości fali P ze wzrostem nasycenia gazem. Dla wariantu 2 prędkość maleje bardzo szybko i osiąga minimum przy nasyceniu około 0,1, następnie zaczyna jednostajnie wzrastać wraz ze wzrostem nasycenia gazem. W przypadku wariantu 2 moduł odkształcenia objętościowego dla medium złożonego z wody i gazu Kf można obliczyć z równania Wooda (1957): S 1− Sw 1 = w − (2.2.11) K f Kw Kg gdzie: Kw i Kg – moduły sprężystości objętościowej wody i gazu, Sw – objętość wody w mieszaninie Z równania wynika, że już niewielka zawartość gazu powoduje znaczny spadek modułu Kf, co pociąga za sobą spadek prędkości fali P. Wariant 2 odpowiada warunkom złożowym.. Rys. 2.2.4 Zależność prędkości fali P od nasycenia gazem obliczone według modelu KT dla poszczególnych wariantów sposobu rozmieszczenia wody i gazu (patrz tekst) dla a) Kg = 0,00015 GPa; b) Kg = 0,01089 GPa; c) Kg = 0,05 GPa (zaczerpnięte z Bała, 2007).. 11.

(14) Wpływ obecności minerałów ilastych na prędkość rozchodzenia się fal (wg Bała, 2007) Według modelu BG i RHG wartości prędkości fal P i S spadają dla wzrastającej zawartości minerałów ilastych (rys. 2.2.5). Przy zmianie zawartości minerałów ilastych od 0 do 0,8 i przy porowatości 0,2 względna zmiana VP wynosi -0,5 natomiast VS -0,77. Porowatość nie wpływa znacząco na wielkość spadku obu prędkości. Stosunek VP/VS w przypadku modelu BG rośnie wraz ze wzrostem zailenia. Względna zmiana stosunku przy wspomnianym wyżej przedziale zmienności iłów wynosi 0,18 przy porowatości 0,05 i osiąga wartość 0,5 przy porowatości 0,25. Model KT pozwala uwzględnić sposób rozłożenia substancji ilastej w skale (rys. 2.2.6). Przyjmując, że frakcja ilasta występuje w szkielecie skalnym wraz ze wzrostem zailenia (od 0 do 0,8) spada prędkość fali P i S (względne zmiany wynoszą odpowiednio -0,38 i -0,30). Obserwowany jest niewielki wzrost stosunku VP/VS wraz ze wzrostem zailenia (względna zmiana 0,18). W przypadku występowania zailenia w porach wzrost zawartości frakcji ilastej powoduje znacznie większy spadek prędkości rozchodzenia się fal P i S. Względne zmiany wynoszą w przypadku fali P -1,32, natomiast w przypadku fali S obliczone prędkości przy dużych zaileniach wydają się być już nierealistyczne. Minerały ilaste wykazują silne własności anizotropowe (Katahara, 1996 za Bała, 1997) co powoduje, że stosunek VP/VS minerałów ilastych może zmieniać się w przedziale 1,795 do 2,5.. Rys. 2.2.5 Zależność prędkości fal P i S oraz stosunku VP/VS w zależności od zawartości substancji ilastej dla wybranego modelu piaskowca. Model BG – krzywe ciągłe, model RHG – krzywe przerywane (zaczerpnięte z Bała, 2007).. 12.

(15) Rys. 2.2.6 Zależność prędkości fal P i S oraz stosunku VP/VS w zależności od zawartości substancji ilastej dla wybranego modelu piaskowca (z Bała, 2007). Model KT uwzględniający zailenie w szkielecie – krzywe czarne oraz zailenie w porach – krzywe fioletowe).. Relacje empiryczne łączące prędkości fali podłużnej i poprzecznej W przypadku braku pomiaru w otworze czasu interwałowego fali S istnieje możliwość uzyskania informacji o prędkości fali S poprzez zastosowanie zależności VS=f(VP) uzyskanej na drodze empirycznej, wykorzystując dostępne korelacje VS–VP. Najbardziej znaną relację, przy założeniu prędkości fal w km/s, zaprezentował Castagna (1985): VP = 1,16VS + 1,36 (2.2.12) Innymi powszechnie stosowanymi zależnościami są zależności Greenberga i Castagni (1992), dla których VP jest powiązane z VS przy użyciu relacji dla różnych litologii (tab. 2.2.3). Tab. 2.2.3 Zależności Greenberga i Castagni dla różnych matryc skalnych. matryca relacja VS–VP [km/s] piaskowiec. VS=0,80416VP-0,85588. iłowiec. VS=0,766969 VP-0,86735. wapień. VS=-0,05508VP2-1,01677VP-1,03049. dolomit. VS=0,858321VP-0,07775. 13.

(16) Wartości VS są średnią dwóch wagowanych średnich: arytmetycznej i harmonicznej: n. Varytm = ∑ VSi wi. (2.2.13). i =1. Vharm =. 1. (2.2.14). n. ∑V i =1. Si. wi. Inną spotykaną relacją wiążącą prędkości fal sejsmicznych jest równanie Kriefa (1990) o ogólnej postaci: VP2 = aVS2 + b. (2.2.15). Wartości a i b zależą od litologii i rodzaju medium złożowego (tab. 2.2.4). Otrzymane wartości prędkości wyrażone są w km/s. Tab. 2.2.4 Współczynniki równania Kriefa dla różnych litologii litologia a. b. piaskowiec (solanka). 2,213. 3,857. piaskowiec (gaz). 2,282. 0,902. piaskowiec (zailony). 2,033. 4,894. Relacje empiryczne łączące prędkości fali podłużnej i gęstość objętościową W przypadku braku gęstości objętościowej możemy dokonać jej estymacji z prędkości fali podłużnej (i vice versa) stosując zależności Gardnera (1974) wyrażoną ogólnym wzorem: ρ = aVPb. (2.2.16). gdzie: a=0,23 i b=0,25 odpowiadają ogólnie skałom osadowym (VP w ft/s i ρ w g/cm3). Zależność 2.2.16 została również wykorzystana przez Castagna (1993), który zaproponował wartości a i b dla różnych litologii (tab. 2.2.5). Tab. 2.2.5 Współczynniki a i b do równania Gardnera zaproponowane przez Castagna dla różnych litologii (VP w km/s i ρ w g/cm3) litologia a b piaskowiec/iłowiec (średnia). 1,741. 0,25. iłowiec. 1,750. 0,265. piaskowiec. 1,660. 0,261. wapień. 1,500. 0,225. dolomit. 1,740. 0,252. anhydryt. 2,190. 0,160. 14.

(17) Podsumowując powyższe rozważania dokonano zestawienia względnych zmian parametrów takich jak: ΔVP/VP, ΔVS/VS, Δρ/ρ, Δλ/λ, Δμ/μ oraz Δν/ν na granicy dwóch umownych ośrodków geologicznych: skała ilasta – piaskowiec nasycony gazem, skała ilasta – piaskowiec nasycony wodą (tab. 2.2.6). Wartości VP, VS i ρ parametryzujące ośrodek założono za Pietsch et al., (2007). Z punktu widzenia poszukiwań naftowych przedstawiona analiza pokazuje, że nasycenie gazem powoduje największą zmianę stałej Lamégo – λ i oraz współczynnika Poissona – σ.. 1660. 2.45. 1.81. 8.5. 6.8. 0.27. PC/WOD. 3470. 2010. 2.35. 1.73. 9.3. 9.5. 0.24. IŁOWIEC. 3000. 1660. 2.45. 1.81. 8.5. 6.8. 0.27. PC/GAZ. 2670. 1850. 2.25. 1.44. 0.6. 7.7. 0.04. Δλ/λ. Δμ/μ. Δν/ν. VP/VS. Δρ/ρ. 0.15. 0.19. -0.04. 0.10. 0.09. 0.34. -0.12. -0.12. 0.11. -0.09. -0.20. -1.72. 0.13. -1.52. 0.26. 0.08. 0.05. 0.30. 1.81. 0.21. 1.40. ΔZ/Z. ΔVS/VS. IŁOWIEC. ΔVP/VP. MODEL. σ. ρ[g/cm3]. 3000. μ [GPa]. VS [m/s]. IŁOWIEC. MODEL. λ [GPa]. VP [m/s]. Tab. 2.2.6 Zestawienie zmian wybranych parametrów na granicy dla dwóch modeli: iłowiec – piaskowiec nasycony wodą, iłowiec – piaskowiec nasycony gazem.. PC/WOD IŁOWIEC PC/GAZ Δ. 15.

(18) 2.3 Klasyfikacja zmian amplitud z offsetem. Wprowadzenie w latach 60 i 70 do sejsmicznych poszukiwań naftowych metody profilowań wielokrotnych, polegającej na rejestracji fal odbitych od tego samego fragmentu reflektora i następnie ich sumowaniu, umożliwiło polepszenie jakości sejsmiki poprzez podniesienie stosunku sygnału do szumu. Dodatkowo dało początek analizie zmian wielkości rejestrowanych amplitud z odległością źródło odbiornik - AVO. Najbardziej znaną pierwszą pracą w obszarze analizy AVO przedstawił Ostrander (1984), który zaobserwował wzrost wielkości amplitud wraz z offsetem dla piaskowca nasyconego gazem. Od tego czasu obserwuje się wzrost zainteresowania metodą AVO w przemysłowych poszukiwaniach naftowych. Skomplikowany kształt wykresów zmienności amplitud z kątem padania fali opisywanych przez równania Knotta-Zoeppritza wymagał wprowadzenia uproszczonej klasyfikacji. Biorąc pod uwagę, że głównym celem analizy AVO są bezpośrednie poszukiwania węglowodorów, wprowadzono podział anomalii w odniesieniu do charakteru zmian amplitud zarejestrowanych w stropie piaskowca nasyconego węglowodorami, w nadkładzie którego znajduje się nieprzepuszczalna skała ilasta. Podstawą klasyfikacji stanowiły dwa fakty: znak współczynnika odbicia przy pionowym padaniu fali na granicę oraz kierunek zmian amplitud, czyli fakt czy amplitudy zarejestrowane dla dalszych offsetów maleją czy rosną. Pierwotnie dokonano podziału na trzy klasy. Następnie klasyfikacja została rozszerzona i obecnie w literaturze spotyka się z pięcioma klasami AVO oznaczanymi zwyczajowo dużymi cyframi rzymskimi (rys. 2.3.1).. Rys. 2.3.1 Klasyfikacja zmian amplitud z offsetem – klasy AVO. Zestawione i zmodyfikowane za Rutherford i Williams (1989), Ross i Kinman (1995) oraz Castagna i Swan (1993).. W przypadku występowania anomalii klasy I, przy pionowym padaniu fali na strop piaskowca nasyconego gazem mamy do czynienia z dużym dodatnim 16.

(19) współczynnikiem odbicia oraz spadkiem amplitudy dla większych kątów. Klasa I jest charakterystyczna dla starych, skonsolidowanych piaskowców o niewielkiej porowatości. Dla klas II i IIp, przy pionowym padaniu fali na strop piaskowca nasyconego gazem znak współczynnika może być dodatni (IIp) lub ujemny (II), a moduł amplitudy zwykle jest niewielki. Wraz ze wzrostem kąta padania w przypadku klasy IIp amplitudy maleją i następnie stają się ujemne. W przypadku klasy II, wraz ze wzrostem kąta padania, ujemne niewielkie amplitudy zwiększają swoją wartość bezwzględna. Klasy IIp i II są charakterystyczne dla piaskowców, w których wzajemne relacje porowatości i stopnia nasycenia gazem powodują, że wartości prędkości fali P i S i/lub gęstości są porównywalne z otoczeniem (Ross i Kinman, 1995). Klasa III jest najczęściej spotykana, w przypadku której przy pionowym padaniu fali na strop piaskowca nasyconego gazem mamy wyraźny współczynnik ujemy, zwiększający swój moduł wraz ze wzrostem kąta padania. Klasa III jest charakterystyczna dla piaskowców średnio skonsolidowanych, położonych relatywnie płytko. W klasie IV przy pionowym padaniu fali na strop piaskowca nasyconego gazem mamy wyraźny współczynnik ujemy, którego moduł nieznacznie maleje dla rosnących kątów padania. Klasa IV występuje w przypadku osadów nieskonsolidowanych, w nadkładzie których występuje warstwa twarda (Castagna i Swan, 1997, Foster et al., 1997). Powyższą klasyfikację niektórzy autorzy odnoszą do wszystkich piaskowców, nie koniecznie do tych wypełnionych gazem . Chcąc sformalizować matematycznie powyższą klasyfikację można odnieść charakter zmian amplitud do usytuowania danego odbicia na płaszczyźnie x: RPP(0), y: G (rys. 2.3.2). G – gradient jest nachyleniem linii aproksymującej wielkości amplitud w układzie: x: sin2θ, y: RPP(0). Powyższa interpretacja wynika z aproksymacji równań Knotta-Zoeppritza przyjętego przez Shueya (1985) – szerzej patrz rozdz. 3: RPP (θ ) ≈ RPP (0) + G sin 2 θ. (2.3.1). Rys. 2.3.2 Klasy AVO na płaszczyźnie RPP(0)–G (Castagna et al., 1997).. 17.

(20) Klasa I wiąże się z występowaniem anomali typu dim spot. Wzrost nasycenia gazem powoduje zmniejszenie się współczynników odbicia w stropie i spągu piaskowca, co odpowiada na złożonym przekroju czasowym niewielkiem amplitudom w strefie nasyconej. Klasa III i IV wiąże się z występowaniem anomalii typu bright spot. Nasycenie gazem powoduje występowanie silnie ujemnej amplitudy w stropie i silnie dodatniej w spągu piaskowca. Anomalia typu bright spot może wykazywać zmianę fazy, gdy impedancja akustyczna w piaskowcu nasyconym wodą jest większa od impedancji nadkładu, natomiast ta jest większa niż impedancja piaskowca nasyconego gazem. Klasy II i IIp występują kiedy impedancje piaskowca nasyconego wodą, gazem oraz nadkładu są porównywalne. Występowanie poszczególnych klas AVO uwarunkowane jest wzajemnymi relacjami współczynnika Poissona – σ i impedancji akustycznej, odzwierciedlając generalne zmiany wartości wspomnianych parametrów z głębokością (rys. 2.3.3).. Rys. 2.3.3 Wzajemne relacje współczynnika Poissona i impedancji akustycznej determinujące występowanie klas AVO z głębokością (Bacon et al., 2007).. Na wykresie korelacji krzyżowej x: RPP(0), y: G można wyróżnić kilka trendów związanych z parametrami petrofizycznymi ośrodka (Simm et al., 2000). Na rys. 2.3.4 przedstawiono przykład odbicia od granicy iłowiec-piaskowiec nienasycony odpowiadający I klasie AVO (niebieski punkt). Wspomniany punkt (jak i wszystkie inne) otoczone są elipsami, opisującymi szum obecny na kolekcjach WPG o wyraźnym trendzie (linia czarna). Zwiększanie porowatości powoduje przesuwanie się punktu w stronę środka układu współrzędnych (trend porowatości). Wprowadzenie gazu powoduje przesunięcie trendu w lewo (punkty z żółtymi elipasami). Odbicia pochodzące od spągu piaskowca znajdują się w lewym górnym rogu rysunku. W przypadku analizy. 18.

(21) dużej liczby odbić punkty odpowiadające poszczególnym kontaktom znajdują się wewnątrz elipsy z wyraźnie zaznaczonym generalnym trendem tła dla granic iłowiec – piaskowiec nienasycony (wet trend). Castagna i Swan (1997) dostrzegli zależność nachylenia linii trendu od średniego stosunku VP/VS (rys. 2.3.5).. Rys. 2.3.4 Różne trendy pojawiające się na wykresie korelacji krzyżowej x: RPP(0) y: G (zmodyfikowane za Simm et al., 2000).. Rys. 2.3.5 Położenie linii trendu aproksymującej odbicia dla kontaktów skał nienasyconych od średniego stosunku VP/VS (Castagna i Swan, 1997). Gęstości przeliczone zostały z równania Gardnera.. Czynniki komplikujące odpowiedź AVO Pomijając wpływy akwizycji i przetwarzania danych sejsmicznych na odpowiedź AVO (rozdz. 2.4) jest ona warunkowana również wpływami anizotropii i cienkowarstwowości ośrodka skalnego.. 19.

(22) Thomsen (1986) uznał, że dla większości ośrodków geologicznych mamy do czynienia ze słabą anizotropią. Efekty anizotropowe wpływają na odpowiedź AVO dla ośrodków silikoklastycznych (Wright, 1987, Kim et al., 1993). Słabą anizotropię można opisać stałymi ε, γ, δ zdefiniowanymi jako: C − C 33 C − C 44 ε = 11 γ = 66 (2.3.2, 2.3.3) 2C 33 2C 44. δ=. (C13 + C 44 ) 2 − (C 33 − C 44 ) 2 2C 33 (C 33 − C 44 ). (2.3.4). gdzie: Cxx – składowe tensora sztywności C zdefiniowanego dla medium ośrodka anizotropowego o pionowej osi symetrii (transversely isotropic) jako:. C11 2C11 − C12 C13 C= 0 0 0. 2C11 − C12 C11 C13 0 0 0. C13 C13 C 33 0 0 0. 0 0 0 C 44 0 0. 0 0 0 0 C 55 0. 0 0 0 0 0. (2.3.5). 0,5(C11 − C12 ). Stałe Thomsena można wyrazić poprzez prędkości fal P i S rejestrowane w różnych kierunkach: V (90°) − VP (0°) ε= P (2.3.6) VP (0°). ε=. VSH (90°) − VSV (0°) VSH (90°) − VSH (0°) = VSV (90°) VSH (0°). (2.3.7). Stała ε opisuje anizotropie fali P, stała γ najlepiej opisuje różnicę pomiędzy pionową i poziomą polaryzacją fali S propagującej w kierunku poziomym, natomiast stała δ jest trudna do jednoznacznej definicji. Dla ośrodka anizotropowego o pionowej osi symetrii wartość współczynnika odbicia można rozłożyć na część izotropową i anizotropową oraz zakładając słabą anizotropię i niewielkie offsety (Banik, 1987) wyrazić następująco: RPP (θ ) ≈ RPP _ IZO. sin 2 θ + ∆δ 2. (2.3.8). Iłowiec wykazujący własności anizotropowe w nadkładzie piaskowca o własnościach izotropowych wywołuje wpływ na wielkości amplitud odbić od granicy ośrodków. Charakter zmian współczynnika odbicia zależy od impedancji akustycznej charakteryzującej piaskowiec (Blangy, 1994). Obecność cienkich warstw powoduje interferencję sygnałów odbitych od stropu i spągu analizowanej cienkiej warstwy, co wywołuje zakłócenie relacji amplitud (zjawisko tuningu). W przypadku wzrostu offsetu i czasu przebiegu fali separacja odbić od stropu i spągu warstwy staje się trudniejsza (zakłócenia amplitud są większe 20.

(23) dla większych offsetów). Teoretyczną odpowiedź AVO - Rt(θ) pojedynczej cienkiej warstwy można wyrazić wzorem (Lin i Phair, 1993): Rt (θ ) = ω0 ∆T (0)(cosθ ) RPP (θ ) (2.3.9) gdzie: ω0 – częstotliwość dominująca, ∆T(0) – czas podwójny dla przypadku pionowego padania na granicę ,θ – kąt padania fali na granicę. W przypadku występowania dużych kontrastów we własnościach sprężystych pomiędzy cienką warstwą i otoczeniem, wpływ na wielkość amplitud mają także odbicia wewnątrz warstwowe, a także nakładanie się różnych mod fal konwertowanych (Simmons i Backus, 1994). Opis wpływu powyższych czynników na odpowiedź AVO jest niezwykle trudny i wymaga modelowań elastycznych pełnego pola falowego.. 21.

(24) 2.4 Ogólna charakterystyka metod wykorzystania rejestracji sejsmicznych w analizie AVO. Znakomita większość badań sejsmicznych wykorzystuje wzbudzanie i rejestrację fal podłużnych. Znając geometrię pomiaru i rejestrując czasy dojścia fal odbitych od granic sejsmicznych możemy wyznaczyć prędkość rozchodzenia się fali P w ośrodku skalnym i pośrednio wnioskować o rodzaju skały. W celu poszerzenia uzyskanych informacji o prędkość fali poprzecznej i gęstości można odczytać informację (najczęściej pośrednią) o obu parametrach, zawartą w rejestrowanych amplitudach fali podłużnej (rys. 2.4.1). Na etapie akwizycji danych sejsmicznych w celu pozyskania rejestracji pod kątem analizy AVO zwraca się szczególną uwagę na wielkości pokrycia (ilość tras w kolekcji danego wspólnego punktu głębokościowego - WPG), rozkład offsetów i ewentualnie azymuty par źródło-odbiornik (w celu uwzględnienia efektów anizotropowych). Dla analiz AVO najkorzystniejsze byłoby posiadanie rejestracji odpowiadających danemu punktowi odbicia o dobrej jakości (brak szumów koherentnych i niekoherentnych) o możliwie jak najszerszym spektrum offsetów. Dodatkowo zwraca się uwagę na wpływ długości grup geofonowych (arrays), od stosowania których obecnie się odchodzi, bowiem ich długość wywiera wpływ na rejestracje na poszczególnych offsetach (Allen i Peddy, 1993). Efekt ma znaczenie szczególnie dla płytkich warstw, zwiększając tłumienie sygnału dla dalekich offsetów w miarę wzrostu długości grupy. Szerzej problemy wpływu relacji źródło-odbiornik na odpowiedź AVO przedstawiono w pracy Martineza (1993).. Rys. 2.4.1 Idea rejestracji danych sejsmicznych dla celów analizy AVO.. Procedury przetwarzania rejestracji sejsmicznych pod kątem analizy AVO, uwzględniają wszystkie główne etapy przetwarzania charakterystyczne dla powszechnie stosowanej metody pokryć wielokrotnych. Istnieją jednakże pewne etapy, na których przetwarzanie musi być przeprowadzane z uwzględnieniem szczególnej 22.

(25) parametryzacji procedur, a także ostrożności. Poniżej przedstawiono zarys tych procedur, które są istotne z punktu widzenia analizy AVO. Dodatkowo z uwagi na fakt, że do interpretatora trafiają już dane po wykonaniu pierwszych etapów przetwarzania, na które nie ma on zwykle wpływu, przedstawiono poniżej bardziej szczegółową analizę procedur przetwarzania, dostępnych w systemie Hampson-Russell CE8R3.1. wykorzystanym do analiz AVO przy realizacji pracy. Niejednokrotnie wspomniane procedury określa się jako post-processing do celów analizy AVO, którego celem jest korekta i/lub uzupełnienie przetwarzania o niezbędne etapy, które wydają się dla interpretatora nieodzowne. Przetwarzanie danych pod kątem AVO ma na celu osiągnięcie dwóch generalnych celów: • odtworzenie względnych relacji amplitud na kolekcjach WPG • poprawnej geometrycznej lokalizacji rejestrowanych tras. Odtworzenie i zachowanie względnych relacji amplitud stanowi fundamentalny aspekt przetwarzania danych sejsmicznych do AVO. Korekcja dywergencji sferycznej może wpływać na niezerooffsetowe odbicia (Ursin, 1990). Dekonwolucja powinna uwzględniać ideę procedur powierzchniowo-zgodnych (Taner i Koehler, 1981). Średnie zmiany amplitud z offsetem wynikające np. z grupowania geofonów powinny być usuwane. Refleksy wielokrotne należy usuwać stosując filtrację Radona - liniową lub paraboliczną (Hampson, 1986) unikając filtracji f-k, która usuwa amplitudy z tras odpowiadających bliskim offsetom i w niewielkim stopniu zmieniając dla dalszych offsetów, co powoduje powstanie artefaktów AVO. Poprawną geometryczną lokalizację refleksów zapewnia czasowa migracja przed składaniem (prestack time migration), która powinna być stosowana w metodyce przetwarzania do analiz AVO (Resnik et al., 1987). Ma to znaczenie szczególnie przy obecności upadów. Migracja zmniejsza promień strefy Fresnela, zwiększając rozdzielczość poziomą (Berkhout, 1985). Dodatkowo migracja usuwa dyfrakcje powstałe od uskoków i ostrych wyklinowań stratygraficznych, które to interferują z analizowanymi refleksami. Równie istotne może być zastosowanie poprawki DMO, szczególnie przy obecności dużych upadów. Technikę ekstrakcji informacji AVA w obecności upadu przedstawił Shang et al. (1993). Procedury przetwarzania dostępne w systemie Hampson-Russell CE8R3.1 i wykorzystywane przy analizach AVO (na podstawie Hampson-Russell (2008)). Chcąc odczytać informację o własnościach ośrodka geologicznego z rejestrowanych amplitud, należy dokonać transformacji kolekcji wspólnego punktu głębokościowego z domeny offsetu (AVO) do domeny wspólnego kąta padania (AVA – Amplitude versus Angle), ponieważ modele teoretyczne łączą własności ośrodka ze współczynnikami odbicia (realnie rejestrowanymi amplitudami) w funkcji kąta padania fali na granicę (rys. 2.4.2).. 23.

(26) Rys. 2.4.2 Konwersja z domeny offsetu (po lewej) do domeny kąta padania (po prawej).. W praktyce dokonuje się konwersji offsetowej kolekcji WPG do kolekcji WPG (angle gather) i/lub do sumy wspólnego kąta (common angle stack). Konwersja wymaga znajomości rozkładu prędkości fali w ośrodku (w praktyce wykorzystuje się prędkości interwałowe i średnie). Dokładna konwersja wymaga również trasowania promienia sejsmicznego. W praktyce w systemach interpretacyjnych stosuje się przybliżenia (rys. 2.4.3). Najprostszym z nich jest metoda promienia bezpośredniego (straight ray method), która jest prawdziwa tylko dla pojedynczej poziomej warstwy. W jej przypadku wartość kąta padania θ wyraża się wzorem: X θ = arctan( ) (2.4.1) Vt gdzie: oznaczenia jak na rys. 2.4.3: W przypadku większej liczby warstw stosuje się metodę parametru promienia sejsmicznego (ray parameter method) wykorzystującą rozkład prędkości średnich kwadratowych (VRMS) oraz interwałowych (VINT). Wartość kąta padania oblicza się ze wzoru: XVINT θ = arcsin( 2 ) (2.4.2) tV RMS Powyższa metoda stanowi kompromis pomiędzy metodą promienia bezpośredniego i pełnym trasowaniem promienia sejsmicznego.. 24.

(27) Rys. 2.4.3 Idea konwersji offsetu (X) do kąta (θ) dla pojedynczej warstwy o prędkości V i miąższości h.. Wielkości amplitud są liczone różnie w zależności od tego czy wyjściem konwersji jest kątowa kolekcja WPG czy suma wspólnego kąta. W pierwszym przypadku dla każdej próbki na czasie t na danej trasie sejsmicznej znajdowane są dwie blisko położone trasy sejsmiczne i dla korespondujących próbek na czasach t obliczane są kąty θ1 i θ2 tak, aby została spełniona zależność: θ1 < kąt θ, dla którego chcemy obliczyć trasę w kolekcji kątowej < θ2 (2.4.3). Wartość amplitudy jest interpolowana z amplitud na dwóch wybranych trasach. W drugim przypadku procedura konwersji wygląda podobnie z tą różnicą, że brane są pod uwagę wszystkie trasy spełniające zależność: (kąt θ, dla którego chcemy obliczyć trasę w kolekcji kątowej) na trasie poprzedniej < kąt θ obliczony dla rozważanej próbki < (kąt θ, dla którego chcemy obliczyć trasę w kolekcji kątowej) (2.4.4). Sumowane są następnie wszystkie próbki spełniające powyższą zależność. Dzięki częściowemu sumowaniu suma wspólnego kąta jest gładsza i nie zawiera przerw w przeciwieństwie do kątowej kolekcji WPG. Superkolekcje WPG (supergathers) Technika tworzenia superkolekcji WPG ma za zadanie podniesienie stosunku sygnału do szumu na danych przed składaniem. Idea działania jest prosta i polega na sumowaniu tras odpowiadających temu samemu offsetowi, bądź wąskiemu przedziałowi offsetów w ramach kolekcji WPG, jak również często miesza się trasy o tym samym offsecie z sąsiednich WPG. Rezydualne poprawki kinematyczne Wprowadzenie rezydualnej poprawki kinematycznej (residual normal moveout –RNMO) ma na celu poprawienie wyrównania tras w obrębie kolekcji WPG po wprowadzeniu poprawki kinematyczej hiperbolicznej (NMO) . Zakłada się, że poprawka RNMO ma charakter paraboliczny, co daje dobre rezultaty dla struktur 25.

(28) o upadach mniejszych niż 45 stopni. Poprawkę zadaje się określając tabelę: czas rejestracji – t0i, resztkowa krzywizna – tRNMOi odpowiadająca temu czasowi, zakładając równocześnie offset odniesienia – xo, dla którego wprowadzamy poprawkę. Nowy czas – ti danej próbki dla offsetu x po wprowadzeniu poprawki oblicza się ze wzoru: x (2.4.5) t i = t 0i + t RNMOi ( ) 2 x0 Błędy poprawki kinematycznej mogą powodować duże błędy w analizie anomalii AVO Spratt (1987). Eliminacja refleksów wielokrotnych oraz tłumienie szumu przypadkowego. W systemie HRS powyższe procedury dostępne są w postaci filtru INVEST. Przy usuwaniu refleksów wielokrotnych, zakładany jest model reprezentujący refleksy wielokrotne w postaci wiązki parabol o predefiniowanym zakresie oraz o określonych nachyleniach ramion dla każdej próbki czasowej. Z uwagi na fakt, że filtracja jest aplikowana do kolekcji WPG po poprawce kinematycznej, odbicia jednokrotne powinny być zliniowane, natomiast wielokrotne wykazywać krzywiznę. Usuwanie refleksów wielokrotnych polega na minimalizacji funkcji celu J zdefiniowanej jako: J = Lm − d + λ m 2. p. (2.4.6). gdzie: m – model „parabol” w domenie Radona, L – operator transformacji Radona, d – analizowana kolekcja WPG, p – parametr określający „szorstkość” rozwiązania, λ – parametr wagowania. Dodatkowo jako wejście do procedury określamy zakres częstotliwości danych sejsmicznych z powodu przeprowadzania operacji filtracji Radona w dziedzinie spektralnej. Przy spełnieniu warunku minimalizacji funkcji celu J model odpowiadający refleksom wielokrotnym m jest odejmowany od danych d. Rezultatem jest kolekcja WPG zawierająca odbicia jednokrotne. Usuwanie szumów przypadkowych polega na odejmowaniu od zarejestrowanych WPG obliczonego szumu przypadkowego, obliczonego jako różnica zarejestrowanych kolekcji i kolekcji zawierających jedynie odbicia jednokrotne (otrzymywanych dzięki filtracji w domenie Radona). W celu uniknięcia zbyt wygładzonego wyjścia, szum przed odjęciem skaluje się, poprzez techniki doboru odpowiednich okien czasowych oraz założenia odnośnie spodziewanego stosunku sygnału do szumu. Skalowanie offsetów Procedura skalowania offsetów (offset scalling/balancing) (Ross i Beale, 1994) ma na celu korekcję zmian amplitud z offsetem, jeżeli taka jest niezbędna, wywołaną wpływem poprzedzających procedur przetwarzania, najczęściej takich, które zakładają stały poziom amplitud RMS na krótkich i długich offsetach. Procedura polega na. 26.

(29) określeniu trendu zmian amplitud z offsetem w oknie odpowiadającym warstwom nasyconym wodą złożową i następnym określeniu parametrów skalowania, definiujących trend zmian. Następnie zdefiniowane wartości aplikowane są do przeskalowania wybranego fragmentu wolumenu danych sejsmicznych w szerokim oknie czasowym. Statyka resztkowa (trim statics) Procedura ma na celu sprowadzenie do tego samego czasu refleksów dla różnych offsetów w obrębie kolekcji WPG (tzw. „zliniowanie”). Opiera się ona na określeniu wzajemnych przesunięć czasowych poszczególnych tras kolekcji WPG i trasy odniesienia (najczęściej trasy sumarycznej) poprzez procedurę korelacji krzyżowej. Parametryzacja procedury polega na zdefiniowaniu ilości, długości oraz kroku zstępowania okien korelacji krzyżowej. Obliczone wzajemne przesunięcia czasowe aplikowane są następnie do każdej trasy w kolekcji WPG. Odwrotna filtracja Q Odwrotna filtracja Q ma na celu kompensację zmian amplitud i fazy zapisu sejsmicznego wywołanych nieidealną sprężystością ośrodka (tłumienie wewnętrzne), która zakłóca odpowiedź AVO (Martinez, 1993). Luh (1993) zaproponował przybliżoną poprawkę do gradientu AVO – δG: fτ δG ≈ (2.4.7) Qe gdzie: f – częstotliwość dominująca sygnału sejsmicznego, τ – pionowy czas podwójny do reflektora, Qe – współczynnik dobroci dla całego nadkładu reflektora. Tłumienie wewnętrzne zakłóca w największym stopniu amplitudy przy padaniu promienia bliskim krytycznemu. Filtracja Q powinna być stosowana przed wprowadzeniem poprawek kinematycznych (Kasina, 1998). Filtrację Q dokonujemy poprzez zadanie krzywej zmian Q z głębokością (lub zakładamy stałe dla całego interwału), częstotliwości odcięcia (maksymalnej i minimalnej) oraz częstotliwości odniesienia. Filtracje można wykonać kompensując wpływ tłumienia wewnętrznego osobno na amplitudę lub fazę lub uwzględniając oba parametry jednocześnie.. 27.

(30) 3. Metody analiz zmienności amplitud z offsetem w kontekście identyfikacji akumulacji gazu Subtelny charakter zmian amplitud z offsetem (AVO) oraz duża objętość danych przed składaniem powoduje, że otrzymanie z nich informacji geologicznej wymaga wprowadzenia specjalnych metod interpretacyjnych. Polegają one na ekstrakcji z zapisu sejsmicznego atrybutów wyrażających charakter zmian amplitud z offsetem. Istnieją dwie generalne tendencje w sposobie dobierania atrybutów. Pierwsza z nich polega na obliczeniu atrybutu, który odzwierciedla interpretowalny parametr petrofizyczny. Druga zakłada obliczenie atrybutu, nie koniecznie fizycznego, którego wartości są np. maksymalne dla piaskowców nasyconych gazem, bliskie zeru lub ujemne dla tzw. tła czyli piaskowców nasyconych woda i skał ilastych. Przetestowano na prostych modelach seismogeologicznych metody interpretacji danych sejsmicznych wykorzystujące zmienność amplitud z offsetem. Celem testów i symulacji było wykazanie ich przydatności do identyfikacji akumulacji gazu ziemnego w silikoklastycznych utworach miocenu zapadliska przedkarpackiego. Dlatego też, w przeważającej ilości przypadków modele parametryzowano dobierając prędkość fal P i S oraz gęstość objętościową z dostępnych danych geofizyki otworowej. Przy analizach wzięto pod uwagę lokalne uwarunkowania geologiczne i petrofizyczne, przyjmując powszechnie stosowany w tego typu analizach model, zakładający występowanie gazu ziemnego w piaskowcu, który jest izolowany od góry nieprzepuszczalną warstwą iłowca. W rozdziale dużo miejsca poświęcono problematyce wpływu sposobu obliczeń atrybutów sejsmicznych charakteryzujących ośrodek geologiczny na ich wiarygodność. W tym celu posłużono się testami ekstrakcji uwzględniając uwarunkowania geologiczne oraz charakter dostępnych danych sejsmicznych. Zwrócono również uwagę na zależność analiz od dostępności danych geofizyki otworowej. Obliczenia zostały przeprowadzone w systemie MATLAB.. 28.

(31) 3.1 Atrybuty intercept, gradient, curvature i ich kombinacje 3.1.1 Ekstrakcja atrybutów A, B i C z kolekcji kątowej WPG. Drogę do matematycznego opisu zmian amplitud z offsetem otworzyły aproksymacje równań Knotta-Zoeppritza wiążące wielkość współczynnika odbicia – RPP z kątem padania fali na granicę – θ i wielkościami oddającymi własności sprężyste ośrodka geologicznego. Pierwszą przedstawił Bortfeld (1961) zakładając niewielkie zmiany własności sprężystych na granicy rozdziału. Podobnie rozumowanie przyjęli Richards i Frasier (1976) oraz Aki i Richards (1980) przedstawiając następującą aproksymację (zwaną najczęściej Akiego-Richardsa): RPP (θ ) ≈. ∆VP ∆VS  ∆ρ  1 1 2 2  + − 4 p p VS2 (1 − 4 p p VS2 ) 2 VS 2  ρ  2 cos θ VP. (3.1.1.1). gdzie: pp – parametr promienia sejsmicznego, ΔVP, ΔVS, Δρ – różnice parametrów na granicach, VP, VS, ρ – średnie wartości parametrów na granicy rozdziału. Po przekształceniach Wigginsa (za Hampson–Russell, 2008) aproksymacja (3.1.1) powiązała RPP(θ) z atrybutami A,B i C charakteryzującymi ośrodek: RPP (θ ) ≈ A + B sin 2 θ + C sin 2 θ tan 2 θ. (3.1.1.2). gdzie:. A= B=. 1  ∆VP ∆ρ  +  2  VP ρ . (3.1.1.3). V V ∆VS 1 ∆VP ∆ρ − 4( S ) 2 − 2( S ) 2 2 VP V P VS VP ρ. (3.1.1.4). 1 ∆VP 2 VP. (3.1.1.5). C=. Przyjęte nazwy atrybutów A – intercept, B – gradient, C – curvature wywodzą się z ich geometrycznej interpretacji. Amplitudy z kolekcji WPG odpowiadające refleksowi od danej granicy sejsmicznej naniesione na płaszczyznę sin2θ–RPP(θ) ułożą się zgodnie z określonym wzorem, wynikającym z własności sprężystych po obu stronach granicy sejsmicznej. Aproksymując punkty prostą, punkt jej przecięcia z osią RPP(θ) wskazuje wielkość intercept – A, jej nachylenie – gradient – B. Wziąwszy pod uwagę trzeci człon aproksymacji (3.1.1.2), chmurę punktów aproksymuje krzywa o krzywiźnie opisywanej przez atrybut C (rys. 3.1.1.1).. 29.

(32) Rys. 3.1.1.1 Aproksymacja Akiego-Richardsa dwu- (czarna) i trój- członowa (niebieska) wartości współczynników odbicia RPP wg równań Knotta-Zoeppritza (czerwona) na płaszczyźnie sin2θ-RPP(θ). Punkty obrazują wartości teoretyczne RPP z gaussowskim, nieskorelowanym z trasy na trasę, szumem przypadkowym.. Wartości atrybutów A, B i C można obliczyć rozwiązując równanie: (3.1.1.6). d = Gm. gdzie: d – wektor danych zawierający odczytane wielkości amplitud z kolekcji WPG dla danego czasu rejestracji t, m – wektor szukanych parametrów A, B i C, G –macierz opisującą model:. A(θ1 ) A(θ 2 ) d= M A(θ N ). 1 sin 2 θ1. A (3.1.1.7), m = B C. (3.1.1.8), G =. 1 sin 2 θ 2 M. M. 1 sin 2 θ N. sin 4 θ1 1 - sin 2 θ1 sin 4 θ 2 1 - sin 2 θ 2 M 4 sin θ N 1 - sin 2 θ N. (3.1.1.9). Macierzowy układ równań (3.1.1.6) opisuje zagadnienie z nadmiarem (overdetermined), bowiem ilość elementów wektora d jest większa od ilości estymowanych parametrów wektora m. Rozwiązanie równania wymaga minimalizacji funkcji celu fC zdefiniowanej w sensie najmniejszych kwadratów jako: f C = Gm − d. 2. (3.1.1.10). 30.

(33) Rozwiązanie powyższego równania stanowi wektor mest , którego elementami są szukane atrybuty A, B oraz C i wyraża się wzorem: mest = (G T G ) −1 G T d. (3.1.1.11). Rezultaty obliczeń atrybutów A, B i C z powyższego równania obarczone są wpływem szumu w wektorze danych d oraz długością kolekcji kątowej WPG, która wpływa na ilość elementów macierzy G. Przedstawiono zatem poniżej analizę wpływu wspomnianych czynników na wiarygodność estymat atrybutów. Metodykę analiz oparto na pracach Swan (1993) oraz teorii dotyczącej rozwiązywania zagadnień odwrotnych Menke (1984). Macierz kowariancji estymowanych parametrów modelu cov m definiujemy jako:. σ A2 cov m = nb na* nc na*. na nb*. na nc*. σ B2. nb nc*. * c b. nn. σ. (3.1.1.12). 2 C. gdzie: na – komponent szumu w estymacie atrybutu A, itd., σA2 – oznacza wariancję w estymacie atrybutu A, itd. Macierz cov m można wyznaczyć ze wzoru: cov m = (G T G ) −1 G T [cov n]G (G T G ) −1. (3.1.1.13). gdzie cov n jest macierzą kowariancji tras kątowej kolekcji WPG zdefiniowanej jako:. cov n =. σ 12 σ 12 L σ 1n σ 21 σ 22 L σ 2 n M. M. σ n1 σ n 2. (3.1.1.14). O M L σ n2. gdzie: σn2 – wariancja, σij – kowariancje poszczególnych tras. Założywszy że szum na poszczególnych trasach składowych kolekcji kątowej WPG jest nieskorelowany z trasy na trasę o średniej równej 0 i wariancji ν2 otrzymujemy: cov n = Iν 2. (3.1.1.15). Upraszcza to równanie (3.1.1.13) opisujące macierz kowariancji modelu: cov m = (G T G ) −1 G T Iν 2G (G T G ) −1 = (G T G ) −1ν 2. (3.1.1.16). Analiza elementów na diagonali składowych macierzy kowariancji modelu cov m umożliwia obliczenie względnego stosunku sygnału do szumu r(S/N) atrybutów A, B i C w zależności od rozpiętości kątowej kolekcji WPG. Względny stosunek definiujemy jako: r( S / N ) i =. e( S / N ) w( S / N ). =. A / cov mii A /ν. =. v cov mii. (3.1.1.17). 31.

(34) gdzie: |A| – wartość bezwzględna sygnału wejściowego, mii –element na diagonali macierzy kowariancji modelu cov m o indeksie ii, ν – odchylenie standardowe szumu w wektorze danych. Zgodnie z powyższą metodyką przeanalizowano względny stosunek sygnału do szumu r(S/N) w przypadku estymacji atrybutów A, B i C (rys. 3.1.1.2). Przyjmując amplitudę sygnału wejściowego jako Am, względny stosunek sygnału do szumu – r(S/N) dla atrybutu B wynosi 1% przy maksymalnym kącie zawartym w kolekcji równym ok. 9°. Wartość wzrasta do 10% przy długości kolekcji ok. 25° osiągając około 60% przy obecności kąta 50°. Dla przypadku atrybutu C stosunek r(S/N) osiąga wartość 1% dopiero przy ok. 21°, wzrastając do około 70% dla przy kącie 50°. Estymata atrybutu A posiada wysoki stosunek sygnału do szumu i wzrasta niewiele wraz ze wzrostem maksymalnego kąta obecnego w kolekcji. Estymata C osiąga większy stosunek sygnału do szumu niż B przy maksymalnym kącie około 46°. Zakładając estymację jedynie dwóch atrybutów A i B wartości stosunku r(S/N) rozkładają się odmiennie. Dlatego w przypadku długich kolekcji warto stosować estymację trzech parametrów, natomiast estymować parametr B stosując estymację dwóch parametrów (szerzej patrz rozdz. 3.3). Rys. 3.1.1.2 Teoretyczna wiarygodność estymacji atrybutów A, B i C (po lewej) oraz A, B (po prawej) w przypadku występowania nieskorelowanego z trasy na trasę szumu przypadkowego. Metodyka za Swan (1993).. Powyższe teoretyczne rozważania zakładały występowanie szumu gaussowskiego o wartości średniej równej zero, nieskorelowanego z trasy na trasę. W rzeczywistym przypadku kolekcje WPG mają skończoną długość oraz zakres kątowy, co powoduje, że założenia co do szumu nie są spełnione. W związku z powyższym przeanalizowano poniżej wpływ szumu na estymaty atrybutów A, B i C opierając się na symulacjach ekstrakcji atrybutów, odnosząc rozważania do konkretnego przypadku modelowego. Założono model składający się z dwóch kontaktujących się jednorodnych półprzestrzeni aproksymujących kontakt warstwy ilastej (górna) z piaskowcem nasyconym gazem (dolna) (zał. I – model I), parametryzowany na warunki III klasy AVO. Stosując rozwiązanie Knotta-Zoeppritza wygenerowano dla modelu współczynniki. 32.

(35) odbicia dla kątów padania od 0° do 50° co 1°. Otrzymany wektor współczynników odbicia dzoep wykorzystany został do dalszych analiz. Symulacja szumu przypadkowego, nieskorelowanego z trasy na trasę, polegała na losowaniu elementów wektora szumu – ni z rozkładu normalnego o wariancjach σ2 0,0001, 0,001 oraz 0,01 (rys. 3.1.1.3). Wygenerowano następnie 1000 wektorów odpowiadających szumowi n tworząc wektory danych d: d = d zoep + n (3.1.1.18) Dla każdego wektora d dokonano ekstrakcji atrybutów A, B i C. Otrzymane wyniki przedstawiono na wykresach krzyżowych A–B oraz A–C (rys. 3.1.1.4).. Rys. 3.1.1.3 Współczynniki odbicia dla modelu I (z gazem) z równania Zoeppritza – linia czarna oraz współczynniki zakłócone szumem przypadkowym o wariancji 0,0001(czerwone), 0,001(zielone), 0,01(niebieskie). Rys. 3.1.1.4 Wykresy korelacji krzyżowej A–B i A–C dla modelu I (zał. I). Ekstrakcja atrybutów A, B i C dla 1000 wektorów danych dla przedziału kątów [0° ; 50°] przy założeniu różnych poziomów szumu przypadkowego – rozwiązanie nietłumione.. 33.

(36) Analiza otrzymanych rezultatów pokazuje, że obecność szumu przypadkowego na kolekcjach CDP wpływa na wyniki obliczeń atrybutów A, B i C. Estymacja atrybutu A wykazuje najmniejszą czułość na obecność szumu przypadkowego w kolekcji WPG. Nawet przy średnim poziomie szumu odchylenie parametru od wartości teoretycznej jest rzędu 0,1. Estymaty atrybutu B i C zmieniają się w znacznie większym zakresie, rzędu 0,4 co jest nieakceptowalną wartością. W praktyce stosuje się stabilizację obliczeń atrybutów A, B i C poprzez dodanie do diagonali macierzy GTG jednostkowej przemnożonej przez element skalujący λ (tzw. miękkie więzy – soft constraints) będące rozwiązaniem tłumionym: mest = (G T G + λI ) −1 G T d. (3.1.1.19). Przeprowadzono identyczne doświadczenie obliczając atrybutu A, B i C, stosując powyższe podejście zakładając λ=1,0 (rys. 3.1.1.5) oraz λ=10,0 (rys. 3.1.1.6). Zwiększanie parametru λ powoduje zmniejszenie rozrzutu estymat atrybutów A, B i C, jednakże kosztem przesunięcia na płaszczyźnie ich wartości względem wartości teoretycznej, co szczególnie, dla powyższego przypadku, obserwuje się w przypadku atrybutu B.. Rys. 3.1.1.5 Wykresy korelacji krzyżowej A–B i A–C dla modelu I (zał. I). Ekstrakcja atrybutów A, B i C dla 1000 wektorów danych dla przedziału kątów [0°; 50°] przy założeniu różnych poziomów szumu przypadkowego – rozwiązanie tłumione λ=1,0.. 34.

(37) Rys. 3.1.1.6 Wykresy korelacji krzyżowej A–B i A–C dla modelu I (zał. I). Ekstrakcja atrybutów A, B i C dla 1000 wektorów danych dla przedziału kątów [0°; 50°] przy założeniu różnych poziomów szumu przypadkowego – rozwiązanie tłumione λ=10,0.. Przeanalizowano wpływ doboru λ na poprawienie wiarygodności estymat atrybutów i ich przesunięcie na płaszczyźnie względem wartości teoretycznych. Miary wielkości rozrzutu punktów zdefiniowano jako L oraz odchylenia estymaty od wartości teoretycznej zdefiniowano jako D (rys. 3.1.1.7):. L = ( Amaks − Amin ) 2 + ( Bmaks − Bmin ) 2. (3.1.1.20). D = ( Ateor − Aśr ) 2 + ( Bteor − Bśr ) 2. (3.1.1.21). Rys. 3.1.1.7 Elipsa szumu (żółta) o dłuższej półosi L związana z estymatami atrybutów A i B z zaznaczonym przesunięciem D związanym z wykorzystaniem tłumionej metody najmniejszych kwadratów.. Dokonano symulacji polegającej na 500-krotnej estymacji atrybutów A i B z wektorów współczynników odbicia d, wygenerowanych dla kontaktu iłowiecpiaskowiec nasycony wodą (model I – zał. I) i następnie zaszumionych. Dla każdej takiej symulacji założono λ z przedziału [10-4, ..., 102] zawierającego 100 parametrów λ 35.

(38) rozłożonych logarytmicznie. Dla każdego założonego λ obliczono wielkości D i L (rys. 3.1.8). Przy założonym λ z przedziału 10-4 do 3*10-2 długość elipsy szumu pozostaje stała, następnie zaczyna zmniejszać się osiągając wartość 0,1 przy λ=5*10-1 i 0,01 przy λ = 102. Odległość otrzymanych estymat atrybutów A i B względem ich wartości teoretycznych Ateor i Bteor – L nie zmienia się w przedziale 10-4 do 10-2 (anomalie wynikają ze skończonej (500) ilości powtórzeń), następnie wzrastając do 0,1 dla λ=1. Dobierając optymalne λ należy minimalizować długość elipsy szumu L nie pozwalając na znaczny wzrost odchylenia D (co w przypadku dysponowania danymi otworowymi można oszacować stosując modelowania). Dla analizowanego przypadku za optymalne λ wybrano wartość 0,5.. Rys. 3.1.1.8 Zależność długości półosi elipsy polaryzacji – L związanej z niejednoznacznością otrzymanej estymaty pary atrybutów A i B, a wielkością jej przesunięcia względem wartości teoretycznych – D.. 36.

(39) Stabilizacja gradientu. Powyższe rozważania pokazują, że estymacja gradientu – B obarczona jest wpływem szumu. Whitcombe et al. (2004) zaproponowali technikę stabilizacji obliczonego już atrybutu B, wykorzystującą obserwowalny fakt, że wielkości atrybutów A i B na ich wykresie korelacji krzyżowej, dla próbki niezakłóconego szumem pojedynczego odbicia sejsmicznego, układają się wzdłuż prostej (Keho et al., 2001). Szum wszelkiej natury powoduje odejście od tej reguły. Przedstawiono za Withcombem doświadczenie polegające na ekstrakcji atrybutów A i B z kolekcji kątowej od 0° do 30° z krokiem 1° o długości 300ms. Wartości dla każdego czasu t (co 1ms) kolekcji stanowiły wektor d, którego elementy były losowane z rozkładu normalnego o wariancji σ2=0,005 i μ=0. Dla kolekcji obliczono atrybuty A i B (rys. 3.1.1.9a). Wielkości atrybutów są ze sobą skorelowane i układają się wzdłuż linii trendu. Linia ortogonalna do linii trendu definiuje tzw. kierunek składania i jest nachylona do osi 0A pod kątem χ≈5°, który wynika z zależności R=A+0,087B (0,087≈sin2θśr[0-30º] =tan χ). Rotacja układu o kąt χ powoduje, że trend szumu układa się wzdłuż osi 0A (rys. 3.1.1.9b) na płaszczyźnie AR–BR. Następnie w ruchomym oknie zawierającym kilka próbek dokonano regresji liniowej wyznaczając współczynniki ar i br. Eliminacja szumu polega na wykorzystaniu wyliczonych współczynników do predykcji wartości atrybutu BT= ar+ ARbr, natomiast atrybut A pozostaje bez zmian AR= AT (rys. 3.1.9c). Ostatecznie wyliczone atrybuty AT i BT przedstawiono na płaszczyźnie AN–BN po odwrotnej rotacji układu współrzędnych o kąt – χ (rys. 3.1.9d). Elipsa szumu po stabilizacji ma mniejsze rozmiary.. 37.

(40) Rys. 3.1.1.9 Metoda stabilizacji gradientu dla kolekcji kątowej WPG zawierającej jedynie szum zaproponowana przez Withcomba; wykresy korelacji krzyżowej atrybutów A i B: a) dane wejściowe, b) dane po rotacji układu współrzędnych, c) dane po regresji w ruchomym oknie i predykcji atrybutu B, d) dane wyjściowe po cofnięciu rotacji układu współrzędnych.. Powyższą metodykę przetestowano na przykładzie teoretycznej kolekcji wyliczonej w oparciu o model IV (zał. I) z dodanym szumem przypadkowym (rys. 3.1.1.10). Dokonano ekstrakcji atrybutów A i B, poddając następnie atrybut B stabilizacji. Porównano atrybuty B ekstrahowane z kolekcji bez szumu, z szumem i z szumem po stabilizacji (rys. 3.1.1.11).. 38.

(41) Rys. 3.1.1.10 Kolekcja kątowa WPG wyliczona w oparciu o model IV (po lewej), ta sama kolekcja z dodanym szumem przypadkowym (po prawej).. Obliczenia dowiodły, że stabilizacja wyraźnie zmniejszyła rozrzut punktów odpowiadających refleksom od warstw nasyconych gazem oraz zminimalizowała trend odpowiadający szumowi przypadkowemu (rys. 3.1.1.11 a i d). Stabilizacji uległ także atrybut B obliczony z kolekcji WPG z szumem, co widać w interwałach czasowych odpowiadających obszarom występowanie jedynie szumu przypadkowego 20-50ms oraz 220-270ms (rys. 3.1.1.12). Efekt stabilizacji widać również dla interwału odpowiadającego refleksom z przedziału 180-205ms.. 39.

(42) Rys. 3.1.1.11 Metoda stabilizacji gradientu dla kolekcji kątowej WPG wyliczonej w oparciu o model IV z dodanym szumem; wykresy korelacji krzyżowej atrybutów A i B: a) dane wejściowe, b) dane po rotacji układu współrzędnych, c) dane po regresji w ruchomym oknie i predykcji atrybutu B, d) dane wyjściowe po cofnięciu rotacji układu współrzędnych.. 40.

(43) Rys. 3.1.1.12 Atrybuty A (na dole) i B (na górze) dla ekstrahowane z kolekcji dla modelu IV; krzywa niebieska – kolekcja WPG bez szumu, krzywa czarna – kolekcja WPG z szumem, krzywa czerwona – kolekcja WPG z szumem, atrybuty po stabilizacji.. 41.

(44) 3.1.2 Interpretacja atrybutów A, B i C Atrybut A, zakładając niewielkie zmiany w VP i ρ pomiędzy przyległymi warstwami, oddaje zmiany impedancji akustycznej. Przyjmuje on wartości dodatnie gdy na granicy zachodzi: ZN < ZN+1 i odzwierciedla refleksyjność fali podłużnej przy pionowym padaniu na granicę– RPP(0). Czuły jest na zmiany litologii i porowatości, a przede wszystkim rodzaj medium wypełniającego pory. Wraz ze wzrostem nasycenia gazem danej skały wartość atrybutu wzrasta, ale nie proporcjonalnie do wielkości nasycenia. Atrybut B jest trudniejszy w bezpośredniej interpretacji bowiem jego wartość jest funkcją złożoną. Ujemne wartości o dużym module atrybutu B w stropie piaskowca mogą wskazywać na nasycenie gazem w przypadku występowania I, II, IIp i III klasy AVO (rys.3.1.2.1).. Rys. 3.1.2.1 Kątowa kolekcja WPG wygenerowana w oparciu o model III (zał. I) i wyliczone atrybuty A, B oraz product – P (czerwone) z naniesionymi teoretycznymi wartościami A i B (czarne).. W przypadku dostępności wiarygodnej estymaty C odzwierciedlającego względne zmiany prędkości można zdefiniować atrybut:. A−C =. 1  ∆VP ∆ρ  1 ∆VP 1 ∆ρ + − = 2  VP ρ  2 VP 2 ρ. (3.1.2.1). 42.