Wolfram Mathematica dalam Geofisika

Terima kasih kepada penulis blog Anton Ekimenko untuk ceramahnya

Pendahuluan

Catatan ini ditulis setelah Konferensi Teknologi Rusia Wolfram dan berisi sinopsis dari laporan yang saya berikan. Acara tersebut berlangsung pada bulan Juni di kota St. Petersburg. Mengingat fakta bahwa saya bekerja satu blok dari tempat konferensi, saya tidak bisa tidak menghadiri acara ini. Pada 2016 dan 2017, saya mendengarkan laporan konferensi, dan tahun ini saya membuat presentasi. Pertama, topik yang menarik (menurut saya) telah muncul yang sedang kami kembangkan bersama Kirill Belov , dan kedua, setelah penelitian panjang mengenai undang-undang Federasi Rusia mengenai kebijakan sanksi, dua lisensi Wolfram Mathematica telah muncul di perusahaan tempat saya bekerja.

Sebelum beralih ke topik presentasi saya, saya ingin mencatat pengaturan acara yang baik. Halaman konferensi menggunakan gambar Katedral Kazan. Katedral adalah salah satu atraksi utama St. Petersburg dan sangat jelas terlihat dari aula tempat konferensi diadakan.



Di pintu masuk ke Universitas Ekonomi Negeri St. Petersburg, para peserta disambut oleh asisten mahasiswa - mereka tidak membiarkan mereka tersesat. Saat registrasi, suvenir kecil diberikan (adhesi mainan, pena, stiker dengan simbol Wolfram). Makan siang dan rehat kopi juga termasuk dalam jadwal konferensi. Tentang kopi dan pai lezat, saya sudah mencatat di dinding grup - koki yang baik. Dengan bagian pengantar ini, saya ingin menekankan bahwa acara itu sendiri, formatnya dan tempat sudah membawa emosi positif.

Laporan, yang disiapkan oleh saya dan Kirill Belov, disebut “Menggunakan Wolfram Mathematica untuk menyelesaikan masalah geofisika terapan. Analisis spektral data seismik atau "di mana sungai purba mengalir." Isi laporan mencakup dua bagian: pertama, penggunaan algoritma yang tersedia di Wolfram Mathematica untuk analisis data geofisika, dan kedua, ini adalah cara menempatkan data geofisika di Wolfram Mathematica.

Eksplorasi seismik

Pertama, Anda perlu melakukan sedikit perjalanan ke geofisika. Geofisika adalah ilmu yang mempelajari sifat fisik batuan. Nah, karena batuan memiliki sifat yang berbeda: listrik, magnetik, elastis, ada metode geofisika yang sesuai: eksplorasi listrik, eksplorasi magnetik, eksplorasi seismik ... Dalam konteks artikel ini, kami hanya akan menyentuh eksplorasi seismik secara lebih rinci. Eksplorasi seismik adalah metode utama untuk menemukan minyak dan gas. Metode ini didasarkan pada eksitasi getaran elastis dan rekaman selanjutnya respon dari batuan yang menyusun area penelitian. Eksitasi getaran dilakukan di darat (dinamit atau sumber getaran elastis yang tidak meledak) atau di laut (senjata udara). Getaran elastis merambat melalui ketebalan batuan, dibiaskan dan tercermin pada batas-batas lapisan dengan sifat yang berbeda. Gelombang yang dipantulkan kembali ke permukaan dan direkam oleh geofon berbasis darat (biasanya perangkat elektrodinamik berdasarkan pergerakan magnet yang tergantung pada koil) atau hidrofon di laut (berdasarkan efek piezoelektrik). Pada saat kedatangan gelombang, seseorang dapat menilai kedalaman lapisan geologis.

Peralatan derek kapal seismik:



Pistol udara menggairahkan getaran elastis:



Gelombang melewati massa batuan dan direkam oleh hidrofon:



Kapal penelitian eksplorasi geofisika "Ivan Gubkin" di dermaga di jembatan Blagoveshchensky di St. Petersburg:

Model sinyal seismik

Batuan memiliki sifat fisik yang berbeda. Untuk eksplorasi seismik, sifat elastis paling penting - kecepatan rambat getaran dan kepadatan elastis. Jika dua lapisan memiliki sifat yang sama atau serupa, gelombang tidak akan "memperhatikan" batas di antara mereka. Jika kecepatan gelombang di lapisan berbeda, maka refleksi akan terjadi pada batas lapisan. Semakin besar perbedaan sifat, semakin kuat pantulannya. Intensitasnya akan ditentukan oleh koefisien refleksi (rc):



di mana ρ adalah kepadatan batuan, ν adalah kecepatan gelombang, 1 dan 2 menunjukkan lapisan atas dan bawah.

Salah satu model sinyal seismik yang paling sederhana dan paling umum digunakan adalah model konvolusi, ketika jejak seismik terdaftar disajikan sebagai hasil konvolusi dari sekuens koefisien refleksi dengan pulsa probe:



di mana s (t) adalah trek seismik, mis. semua yang direkam oleh mikrofon atau geophone selama waktu perekaman tetap, w (t) adalah sinyal yang dihasilkan oleh airgun, n (t) adalah derau acak.

Kami menghitung misalnya survei seismik sintetis. Kami akan menggunakan pulsa Ricker yang banyak digunakan dalam eksplorasi seismik sebagai sinyal awal.

length=0.050; (*Signal lenght*) dt=0.001;(*Sample rate of signal*) t=Range[-length/2,(length)/2,dt];(*Signal time*) f=35;(*Central frequency*) wavelet=(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)]; ListLinePlot[wavelet, Frame->True,PlotRange->Full,Filling->Axis,PlotStyle->Black, PlotLabel->Style["Initial wavelet",Black,20], LabelStyle->Directive[Black,Italic], FillingStyle->{White,Black},ImageSize->Large,InterpolationOrder->2] 

Sumber pulsa seismik



Kami menetapkan dua batas pada kedalaman 300 ms dan 600 ms, dan koefisien refleksi akan menjadi angka acak.

 rcExample=ConstantArray[0,1000]; rcExample[[300]]=RandomReal[{-1,0}]; rcExample[[600]]=RandomReal[{0,1}]; ListPlot[rcExample,Filling->0,Frame->True,Axes->False,PlotStyle->Black, PlotLabel->Style["Reflection Coefficients",Black,20], LabelStyle->Directive[Black,Italic]] 

Urutan koefisien refleksi:



Hitung dan tampilkan trek seismik. Karena koefisien refleksi memiliki tanda-tanda yang berbeda, kami memperoleh dua refleksi bergantian di jalur seismik.

 traceExamle=ListConvolve[wavelet[[1;;;;1]],rcExample]; ListPlot[traceExamle, PlotStyle->Black,Filling->0,Frame->True,Axes->False, PlotLabel->Style["Seismic trace",Black,20], LabelStyle->Directive[Black,Italic]] 

Lagu yang disimulasikan:



Untuk contoh ini, reservasi harus dilakukan - pada kenyataannya, kedalaman lapisan ditentukan tentu saja dalam meter, dan perhitungan jalur seismik terjadi untuk domain waktu. Akan lebih tepat untuk mengatur kedalaman dalam meter dan menghitung waktu kedatangan mengetahui kecepatan di strata. Dalam hal ini, saya segera mengatur layer pada sumbu waktu.

Jika kita berbicara tentang studi lapangan, sebagai hasil dari pengamatan tersebut, sejumlah besar rangkaian waktu tersebut (trek seismik) dicatat. Misalnya, ketika menjelajahi situs sepanjang 25 km dan lebar 15 km, di mana sebagai hasil kerja setiap trek mencirikan sel berukuran 25x25 meter (sel tersebut disebut bin), susunan data akhir akan berisi 600.000 trek. Dengan langkah pengambilan sampel waktu 1 ms, waktu perekaman 5 detik, file data akhir akan lebih dari 11 GB, dan jumlah bahan baku "asli" bisa ratusan gigabyte.

Bagaimana cara bekerja dengan mereka di Wolfram Mathematica ?

Paket GeologiIO

Pengembangan paket dimulai dengan pertanyaan di dinding VK dari kelompok pendukung berbahasa Rusia. Berkat tanggapan masyarakat, solusi ditemukan dengan sangat cepat. Dan sebagai hasilnya, itu tumbuh menjadi perkembangan yang serius. Posting yang sesuai di dinding Wolfram Comunity bahkan dicatat oleh moderator. Saat ini, paket mendukung bekerja dengan tipe data berikut yang secara aktif digunakan dalam industri geologi:

1. mengimpor data peta format ZMAP dan IRAP
2. impor pengukuran dalam sumur format LAS
3. input dan output file seismik dalam format SEGY

Untuk menginstal paket, Anda harus mengikuti instruksi pada halaman pengunduhan paket yang sudah dirakit, mis. jalankan kode berikut pada buku catatan Mathematica apa pun:

 If[PacletInformation["GeologyIO"] === {}, PacletInstall[URLDownload[ "https://wolfr.am/FiQ5oFih", FileNameJoin[{CreateDirectory[], "GeologyIO-0.2.2.paclet"}] ]]] 

Setelah itu, paket akan diinstal di folder default, jalur yang dapat diperoleh sebagai berikut:

 FileNameJoin[{$UserBasePacletsDirectory, "Repository"}] 

Sebagai contoh, kami menunjukkan fitur utama dari paket. Panggilan secara tradisional untuk paket dalam Bahasa Wolfram:

 Get["GeologyIO`"] 

Paket ini sedang dikembangkan menggunakan Wolfram Workbench . Ini memungkinkan Anda untuk menemani fungsionalitas utama paket dengan dokumentasi, yang tidak berbeda dalam format presentasi dari dokumentasi Wolfram Mathematica itu sendiri dan menyediakan paket dengan file uji untuk kenalan pertama.





File seperti itu, khususnya, adalah file Marmousi.segy, model sintetis dari bagian geologi yang dikembangkan oleh French Petroleum Institute. Dengan menggunakan model ini, pengembang menguji algoritme mereka sendiri untuk memodelkan bidang gelombang, pemrosesan data, inversi jejak seismik, dll. Model Marmousi sendiri disimpan di repositori, dari mana paket itu sendiri diunduh. Untuk mendapatkan file, jalankan kode berikut:

 If[Not[FileExistsQ["Marmousi.segy"]], URLDownload["https://wolfr.am/FiQGh7rk", "Marmousi.segy"];] marmousi = SEGYImport["Marmousi.segy"] 

Hasil impor adalah objek SEGYData:



Format SEGY melibatkan penyimpanan berbagai informasi tentang pengamatan. Pertama, ini adalah komentar teks. Informasi tentang tempat kerja, nama-nama perusahaan yang melakukan pengukuran, dll dimasukkan di sini. Dalam kasus kami, tajuk ini dipanggil oleh permintaan dengan kunci TextHeader. Berikut judul teks singkat:

 Short[marmousi["TextHeader"]] 

"Kumpulan data Marmousi dihasilkan di Institute ... kecepatan maksimum 1500 m / s dan maksimum 5500 m / s)"

Model geologi yang sebenarnya dapat ditampilkan dengan menghubungi jejak seismik menggunakan kunci "jejak" (salah satu fitur dari paket ini adalah tidak sensitifnya huruf pada tombol):

 ArrayPlot[Transpose[marmousi["traces"]], PlotTheme -> "Detailed"] 

Model Marmousi:



Saat ini, paket ini juga memungkinkan Anda untuk mengunduh data dalam bagian-bagian dari file besar, sehingga memungkinkan untuk memproses file yang ukurannya dapat mencapai puluhan gigabyte. Paket ini juga mencakup fungsi untuk mengekspor data ke .segy dan penulisan ulang sebagian hingga akhir file.

Secara terpisah, perlu diperhatikan fungsionalitas paket ketika bekerja dengan struktur kompleks file .segy. Karena memungkinkan tidak hanya mengakses kunci dan indeks ke jejak individu, header, tetapi juga mengubahnya dengan tulisan berikutnya ke file. Banyak detail teknis penerapan GeologyIO berada di luar cakupan artikel ini dan mungkin layak mendapatkan deskripsi terpisah.

Relevansi analisis spektral dalam eksplorasi seismik

Kemampuan untuk mengimpor bahan seismik ke Wolfram Mathematica memungkinkan Anda untuk menggunakan fungsionalitas pemrosesan sinyal bawaan untuk data eksperimen. Karena setiap lintasan seismik adalah deret waktu, salah satu alat utama untuk mempelajarinya adalah analisis spektral. Di antara prasyarat untuk analisis komposisi frekuensi data seismik, misalnya, adalah sebagai berikut:

1. Berbagai jenis gelombang ditandai oleh komposisi frekuensi yang berbeda. Ini memungkinkan kita untuk memilih gelombang yang berguna dan menekan gelombang interferensi.
2. Sifat batuan seperti porositas dan saturasi dapat mempengaruhi komposisi frekuensi. Ini memungkinkan Anda memilih batu dengan properti yang lebih baik.
3. Lapisan dengan ketebalan berbeda menyebabkan anomali dalam rentang frekuensi yang berbeda.

Paragraf ketiga sangat penting dalam konteks artikel ini. Di bawah ini adalah cuplikan kode untuk menghitung jejak seismik dalam kasus lapisan dengan ketebalan bervariasi - model irisan. Model ini secara tradisional dipelajari dalam eksplorasi seismik untuk analisis efek interferensi, ketika gelombang yang dipantulkan dari banyak lapisan saling tumpang tindih.

 nx=200;(* Number of grid points in X direction*) ny=200;(* Number of grid points in Y direction*) T=2;(*Total propagation time*) (*Velocity and density*) modellv=Table[4000,{i,1,ny},{j,1,nx}];(* P-wave velocity in m/s*) rho=Table[2200,{i,1,ny},{j,1,nx}];(* Density in g/cm^3, used constant density*) Table[modellv[[150-Round[i*0.5];;,i]]=4500;,{i,1,200}]; Table[modellv[[;;70,i]]=4500;,{i,1,200}]; (*Plotting model*) MatrixPlot[modellv,PlotLabel->Style["Model of layer",Black,20], LabelStyle->Directive[Black,Italic]] 

Model pembentukan semburan:



Kecepatan gelombang di dalam irisan adalah 4500 m / s, di luar irisan adalah 4000 m / s, dan kepadatan diasumsikan konstan 2200 g / cm³. Untuk model seperti itu, kami menghitung koefisien refleksi dan jejak seismik.

 rc=Table[N[(modellv[[All,i]]-PadLeft[modellv[[All,i]],201,4000][[1;;200]])/(modellv[[All,i]]+PadLeft[modellv[[All,i]],201,4500][[1;;200]])],{i,1,200}]; traces=Table[ListConvolve[wavelet[[1;;;;1]],rc[[i]]],{i,1,200}]; starttrace=10; endtrace=200; steptrace=10; trasenum=Range[starttrace,endtrace,steptrace]; traserenum=Range[Length@trasenum]; tracedist=0.5; Rotate[Show[ Reverse[Table[ ListLinePlot[traces[[trasenum[[i]]]]*50+trasenum[[i]]*tracedist,Filling->{1->{trasenum[[i]]*tracedist,{RGBColor[0.97,0.93,0.68],Black}}},PlotStyle->Directive[Gray,Thin],PlotRange->Full,InterpolationOrder->2,Axes->False,Background->RGBColor[0.97,0.93,0.68]], {i,1,Length@trasenum}]],ListLinePlot[Transpose[{ConstantArray[45,80],Range[80]}],PlotStyle->Red],PlotRange->All,Frame->True],270Degree] 

Trek seismik untuk model wedge:



Urutan jejak seismik yang ditunjukkan pada gambar ini disebut bagian seismik. Seperti yang Anda lihat, interpretasinya dapat dilakukan pada tingkat intuitif, karena geometri gelombang yang dipantulkan secara unik sesuai dengan model yang ditetapkan sebelumnya. Jika Anda menganalisis trek secara lebih terperinci, maka Anda dapat melihat bahwa trek dari tanggal 1 hingga sekitar tanggal 30 tidak berbeda - pantulan dari bagian atas formasi dan dari solnya tidak tumpang tindih. Mulai dari rute 31, refleksi mulai mengganggu. Dan, meskipun, dalam model, koefisien pantulan tidak berubah secara horizontal - jalur seismik mengubah intensitasnya dengan perubahan ketebalan reservoir.

Pertimbangkan amplitudo refleksi dari batas atas reservoir. Mulai dari rute ke-60, intensitas pantulan mulai meningkat dan pada rute ke-70 menjadi maksimum. Ini adalah bagaimana gangguan gelombang dari atap dan sol untuk formasi dimanifestasikan, yang mengarah pada beberapa kasus ke anomali yang signifikan dari rekaman seismik.

 ListLinePlot[GaussianFilter[Abs[traces[[All,46]]],3][[;;;;2]], InterpolationOrder->2,Frame->True,PlotStyle->Black, PlotLabel->Style["Amplitude of reflection",Black,20], LabelStyle->Directive[Black,Italic], PlotRange->All] 

Grafik amplitudo gelombang yang dipantulkan dari tepi atas irisan



Adalah logis bahwa ketika sinyal berfrekuensi lebih rendah, maka interferensi mulai muncul pada ketebalan besar formasi, dan dalam kasus sinyal frekuensi tinggi, interferensi terjadi pada ketebalan yang lebih kecil. Fragmen kode berikut menciptakan sinyal dengan frekuensi 35 Hz, 55 Hz dan 85 Hz.

 waveletSet=Table[(1.0-2.0*(Pi^2)*(f^2)*(t^2))*Exp[-(Pi^2)*(f^2)*(t^2)], {f,{35,55,85}}]; ListLinePlot[waveletSet,PlotRange->Full,PlotStyle->Black,Frame->True, PlotLabel->Style["Set of wavelets",Black,20], LabelStyle->Directive[Black,Italic], ImageSize->Large,InterpolationOrder->2] 

Satu set sinyal sumber dengan frekuensi 35 Hz, 55 Hz, 85 Hz



Setelah menghitung jejak seismik dan memplot amplitudo dari gelombang yang dipantulkan, kita dapat melihat bahwa untuk frekuensi yang berbeda suatu anomali diamati pada ketebalan formasi yang berbeda.

 tracesSet=Table[ListConvolve[waveletSet[[j]][[1;;;;1]],rc[[i]]],{j,1,3},{i,1,200}]; lowFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[1]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; medFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[2]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; highFreq=ListLinePlot[GaussianFilter[Abs[tracesSet[[3]][[All,46]]],3][[;;;;2]],InterpolationOrder->2,PlotStyle->Black,PlotRange->All]; Show[lowFreq,medFreq,highFreq,PlotRange->{{0,100},All}, PlotLabel->Style["Amplitudes of reflection",Black,20], LabelStyle->Directive[Black,Italic], Frame->True] 

Grafik amplitudo dari gelombang yang dipantulkan dari tepi atas irisan untuk frekuensi yang berbeda



Kemampuan untuk menarik kesimpulan tentang ketebalan reservoir berdasarkan hasil pengamatan seismik sangat berguna, karena salah satu tugas utama dalam pencarian ladang minyak adalah untuk mengevaluasi titik-titik yang paling menjanjikan untuk meletakkan sumur (yaitu, bagian-bagian di mana reservoir memiliki ketebalan yang besar). Selain itu, di bagian geologis, mungkin ada benda-benda yang, karena asal usulnya, menyebabkan perubahan tajam dalam ketebalan formasi. Ini menjadikan analisis spektral alat yang efektif untuk mempelajarinya. Pada bagian selanjutnya dari artikel ini, kami akan mempertimbangkan objek gelogical seperti itu secara lebih rinci.

Data eksperimen. Di mana mereka diterima dan apa yang harus dicari di dalamnya?

Bahan-bahan yang dianalisis dalam artikel itu diperoleh di wilayah Siberia Barat. Wilayah ini, seperti yang diketahui oleh semua orang tanpa kecuali, adalah wilayah penghasil minyak utama di negara kita. Pengembangan aktif kelahiran metro dimulai di wilayah ini pada tahun 60an abad terakhir. Metode utama untuk menemukan ladang minyak adalah eksplorasi seismik. Sangat menarik untuk mempertimbangkan gambar satelit dari wilayah ini. Dalam skala kecil, sejumlah besar rawa dan danau dapat dicatat, dengan meningkatkan peta, Anda dapat melihat sumur cluster untuk pengeboran sumur, dan dengan meningkatkan peta hingga batasnya, Anda juga dapat membedakan glade profil di mana pengamatan seismik dilakukan.

Gambar satelit peta Yandex - area kota Noyabrsk:



Jaringan situs cluster di salah satu bidang:



Batuan penghasil minyak Siberia Barat terjadi di berbagai kedalaman - dari 1 km hingga 5 km. Volume utama batuan yang mengandung minyak terbentuk pada masa Jurassic dan Cretaceous. Periode Jurassic mungkin dikenal oleh banyak orang dengan film dengan nama yang sama. Iklim periode Jurassic secara signifikan berbeda dari yang modern. Dalam Encyclopedia Britannica ada serangkaian paleocard yang menjadi ciri setiap era gelogik.

Waktu sekarang:



Periode Jurassic:



Harap dicatat bahwa di Jurassic, wilayah Siberia Barat adalah pantai laut (daratan yang dilintasi oleh sungai dan laut dangkal). Karena iklimnya nyaman, dapat diasumsikan bahwa lansekap khas waktu itu tampak sebagai berikut:

Siberia dari Jurassic:



Dalam gambar ini, tidak begitu banyak hewan dan burung yang penting bagi kita, seperti gambar sungai di latar belakang. Sebuah sungai adalah objek yang sangat geologis tempat kami berhenti sebelumnya. Faktanya adalah bahwa aktivitas sungai memungkinkan akumulasi batupasir yang diurutkan dengan baik, yang kemudian menjadi reservoir untuk minyak. Reservoir ini dapat memiliki bentuk yang aneh dan kompleks (seperti dasar sungai) dan memiliki ketebalan yang bervariasi - sepanjang pantai ketebalannya kecil, dan bertambah lebih dekat ke pusat saluran atau di bagian berkelok-kelok. Jadi, sungai-sungai yang terbentuk pada zaman Jurassic sekarang berada di kedalaman sekitar tiga kilometer dan menjadi objek pencarian reservoir minyak.

Data eksperimen. Pemrosesan dan visualisasi

Kami segera membuat reservasi mengenai bahan seismik yang ditampilkan dalam artikel - karena fakta bahwa jumlah data yang digunakan untuk analisis adalah signifikan - hanya sebagian kecil dari set asli jejak seismik yang ditempatkan dalam teks artikel. Ini akan memungkinkan semua orang mereproduksi perhitungan di atas.

Ketika bekerja dengan data seismik, ahli geofisika biasanya menggunakan perangkat lunak khusus (ada beberapa pemimpin industri yang secara aktif menggunakan pengembangan mereka, misalnya, Petrel atau Paradigm), yang memungkinkan Anda untuk menganalisis berbagai jenis data dan memiliki antarmuka grafis yang nyaman. Terlepas dari semua kenyamanan, jenis perangkat lunak ini juga memiliki kelemahan - misalnya, penerapan algoritma modern dalam versi stabil membutuhkan banyak waktu, dan kemungkinan otomatisasi perhitungan biasanya terbatas. Dalam situasi seperti itu, menjadi sangat nyaman untuk menggunakan sistem matematika komputer dan bahasa pemrograman tingkat tinggi yang memungkinkan Anda untuk menggunakan basis algoritmik yang luas dan, pada saat yang sama, mengambil banyak rutinitas. Menggunakan prinsip ini, bekerja dengan data seismik di Wolfram Mathematica dibangun. Tidak praktis untuk menulis fungsional yang kaya dari pekerjaan interaktif dengan data - lebih penting untuk memastikan pemuatan dari format yang diterima secara umum, menerapkan algoritma yang diinginkan untuk mereka dan mengunggahnya kembali ke format eksternal.

Mengikuti skema yang diusulkan, muat data seismik asli dan tampilkan dalam Wolfram Mathematica :

 Get["GeologyIO`"] seismic3DZipPath = "seismic3D.zip"; seismic3DSEGYPath = "seismic3D.sgy"; If[FileExistsQ[seismic3DZipPath], DeleteFile[seismic3DZipPath]]; If[FileExistsQ[seismic3DSEGYPath], DeleteFile[seismic3DSEGYPath]]; URLDownload["https://wolfr.am/FiQIuZuH", seismic3DZipPath]; ExtractArchive[seismic3DZipPath]; seismic3DSEGY = SEGYImport[seismic3DSEGYPath] 

Data yang diunduh dan diimpor dengan cara ini adalah trek yang terdaftar di situs berukuran 10 kali 5 kilometer. Jika data diperoleh dengan metode eksplorasi seismik tiga dimensi (gelombang tidak direkam sepanjang profil geofisika yang terpisah, tetapi di seluruh area pada saat yang sama), dimungkinkan untuk mendapatkan kubus data seismik. Ini adalah objek tiga dimensi, bagian vertikal dan horizontal yang memungkinkan studi rinci tentang lingkungan geologi. Dalam contoh yang dipertimbangkan, kita hanya berurusan dengan data tiga dimensi. Kami dapat memperoleh beberapa informasi dari tajuk teks, misalnya seperti ini

 StringPartition[seismic3DSEGY["textheader"], 80] // TableForm 

C 1 INI ADALAH FILE DEMO UNTUK UJI PAKET GEOLOGYIO
C 2
C 3
C 4
C 5 NAMA PENGGUNA TANGGAL: PENGGUNA WOLFRAM
C 6 NAMA SURVEI: DI MANA SAJA DI SIBERIA
C 7 FILE TYPE 3D SEISMIC VOLUME
C 8
C 9
C10 Z RANGE: PERTAMA 2200M TERAKHIR 2400M

Kumpulan data ini akan cukup bagi kami untuk menunjukkan tahap utama analisis data. Trek dalam file direkam secara berurutan dan masing-masing terlihat seperti gambar berikut - ini adalah distribusi amplitudo dari gelombang yang dipantulkan di sepanjang sumbu vertikal (sumbu kedalaman).

 ListLinePlot[seismic3DSEGY["traces"][[100]], InterpolationOrder -> 2, PlotStyle -> Black, PlotLabel -> Style["Seismic trace", Black, 20], LabelStyle -> Directive[Black, Italic], PlotRange -> All, Frame -> True, ImageSize -> 1200, AspectRatio -> 1/5] 

Salah satu bagian seismik:



Mengetahui berapa banyak jejak yang terletak di setiap arah area yang dipelajari, Anda dapat membentuk array data tiga dimensi dan menampilkannya menggunakan fungsi Image3D []

 traces=seismic3DSEGY["traces"]; startIL=1050;EndIL=2000;stepIL=2; (*        *) startXL=1165;EndXL=1615;stepXL=2; (* Y       *) numIL=(EndIL-startIL)/stepIL+1; (*    *) numXL=(EndXL-startXL)/stepIL+1; (*    Y*) Image3D[ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}],ViewPoint->{-1, 0, 0},Background->RGBColor[0,0,0]] 

Gambar tiga dimensi dari kubus data seismik. (Sumbu vertikal - kedalaman):



Dalam hal objek-objek geologis yang menarik menciptakan anomali seismik yang kuat, alat visualisasi dengan transparansi dapat digunakan. Bagian "tidak penting" dari rekaman dapat dibuat tidak terlihat, hanya menyisakan kelainan. Di Wolfram Mathematica, ini bisa dilakukan menggunakan Opacity [] dan Raster3D [] .

 data = ArrayReshape[Abs[traces/Max[Abs[traces[[All,1;;;;4]]]]],{numIL,numXL,101}]; Graphics3D[{Opacity[0.1], Raster3D[data, ColorFunction->"RainbowOpacity"]}, Boxed->False, SphericalRegion->True, ImageSize->840, Background->None] 

Gambar kubus data seismik menggunakan fungsi Opacity [] dan Raster3D []



Seperti dalam contoh sintetis, pada irisan kubus asli, beberapa batas geologis (lapisan) dengan relief variabel dapat dibedakan.

Alat analisis spektral utama adalah transformasi Fourier. Dengan menggunakannya, Anda dapat mengevaluasi spektrum frekuensi amplitudo dari setiap jalur atau grup jalur. Namun, setelah mentransfer data ke domain frekuensi, informasi hilang tentang pada waktu apa (baca pada kedalaman apa) frekuensi berubah. Agar dapat melokalisasi perubahan sinyal pada sumbu waktu (dalam), digunakan jendela Transformasi Fourier dan dekomposisi wavelet. Artikel ini menggunakan dekomposisi wavelet. Teknologi analisis wavelet mulai digunakan secara aktif dalam eksplorasi seismik di tahun 90-an. Keuntungan atas transformasi Fourier jendela dianggap sebagai resolusi waktu terbaik.

Menggunakan fragmen kode berikut, Anda dapat menguraikan menjadi komponen-komponen individual dari salah satu jejak seismik:

 cwd=ContinuousWaveletTransform[seismicSection["traces"][[100]]] Show[ ListLinePlot[Re[cwd[[1]]],PlotRange->All], ListLinePlot[seismicSection["traces"][[100]], PlotStyle->Black,PlotRange->All],ImageSize->{1500,500},AspectRatio->Full, PlotLabel->Style["Wavelet decomposition",Black,32], LabelStyle->Directive[Black,Italic], PlotRange->All, Frame->True] 

Dekomposisi komponen



Untuk menilai bagaimana energi pantulan didistribusikan pada waktu kedatangan gelombang yang berbeda, skalogr (analog dari spektogram) digunakan. Sebagai aturan, dalam praktiknya, tidak perlu menganalisis semua komponen. Biasanya memilih komponen frekuensi rendah, sedang dan tinggi.

 freq=(500/(#*contWD["Wavelet"]["FourierFactor"]))&/@(Thread[{Range[contWD["Octaves"]],1}]/.contWD["Scales"])//Round; ticks=Transpose[{Range[Length[freq]],freq}]; WaveletScalogram[contWD,Frame->True,FrameTicks->{{ticks,Automatic},Automatic},FrameTicksStyle->Directive[Orange,12], FrameLabel->{"Time","Frequency(Hz)"},LabelStyle->Directive[Black,Bold,14], ColorFunction->"RustTones",ImageSize->Large] 

Scalogram. Hasil dari WaveletScalogram [] Fungsi



Bahasa Wolfram menggunakan fungsi ContinuousWaveletTransform [] untuk transformasi wavelet . Dan penerapan fungsi ini ke seluruh rangkaian jejak dilakukan dengan menggunakan fungsi Table [] . Di sini yang perlu diperhatikan salah satu kekuatan Wolfram Mathematica adalah kemampuan untuk menggunakan paralelisasi ParallelTable [] . Dalam contoh yang diberikan, paralelisasi tidak diperlukan - volume data tidak besar, tetapi ketika bekerja dengan set data eksperimental yang berisi ratusan ribu jejak, ini merupakan keharusan.

 tracesCWD=Table[Map[Hilbert[#,0]&,Re[ContinuousWaveletTransform[traces[[i]]][[1]]][[{13,15,18}]]],{i,1,Length@traces}]; 

Setelah menerapkan fungsi ContinuousWaveletTransform [] , muncul array data baru yang sesuai dengan frekuensi yang dipilih. Dalam contoh di atas, ini adalah frekuensi: 38Hz, 33Hz, 27Hz. Pilihan frekuensi paling sering dilakukan berdasarkan pengujian - mereka memperoleh peta yang efektif untuk kombinasi frekuensi yang berbeda dan memilih yang paling informatif dari sudut pandang seorang ahli geologi.

Jika Anda perlu membagikan hasilnya dengan kolega atau memberikannya kepada pelanggan, Anda dapat menggunakan fungsi SEGYExport [] dari GeologyIO

 outputdata=seismic3DSEGY; outputdata["traces",1;;-1]=tracesCWD[[All,3]]; outputdata["textheader"]="Wavelet Decomposition Result"; outputdata["binaryheader","NumberDataTraces"]=Length[tracesCWD[[All,3]]]; SEGYExport["D:\\result.segy",outputdata]; 

Memiliki tiga kubus seperti itu tersedia (komponen frekuensi rendah, frekuensi menengah dan frekuensi tinggi), sebagai aturan, pencampuran RGB digunakan untuk visualisasi data bersama. Setiap komponen diberi warna tersendiri - merah, hijau, biru. Di Wolfram Mathematica, ini dapat dilakukan menggunakan fungsi ColorCombine [] .

Hasilnya, gambar yang diperoleh dapat digunakan untuk melakukan interpretasi geologis. Berliku-liku yang menggantung pada bagian memungkinkan Anda untuk menguraikan paleorus, yang lebih cenderung menjadi reservoir dan mengandung cadangan minyak. Pencarian dan analisis analog modern dari sistem sungai memungkinkan kita untuk menentukan bagian yang paling menjanjikan dari berkelok-kelok. Lapisan itu sendiri ditandai oleh lapisan tebal batu pasir yang tersortir dengan baik dan merupakan reservoir yang baik untuk minyak. Area di luar anomali "renda" mirip dengan endapan dataran banjir modern. Sedimen dataran banjir terutama diwakili oleh batuan lempung dan pengeboran ke zona ini tidak akan efektif.

Potongan RGB dari kubus data. Di tengah (agak ke kiri tengah) Anda dapat menelusuri sungai yang berkelok-kelok.



Potongan RGB dari kubus data. Di sisi kiri Anda dapat menelusuri sungai yang berkelok-kelok.



Dalam beberapa kasus, kualitas data seismik memungkinkan gambar yang jauh lebih tajam. Itu tergantung pada metodologi pekerjaan lapangan, peralatan yang diadopsi oleh algoritma pengurangan kebisingan. Dalam kasus seperti itu, tidak hanya fragmen dari sistem sungai yang terlihat, tetapi juga seluruh paleorec yang diperluas.

Pencampuran RGB dari tiga komponen kubus data seismik (irisan horizontal). Kedalamannya sekitar 2 km.



Citra satelit dari sungai Volga di wilayah Saratov

Kesimpulan

Di Wolfram Mathematica, Anda dapat menganalisis data seismik dan menyelesaikan masalah aplikasi yang terkait dengan pencarian mineral, dan paket GeologyIO memungkinkan Anda membuat proses ini lebih nyaman. Struktur data seismik sedemikian rupa sehingga penggunaan metode bawaan untuk mempercepat perhitungan ( ParallelTable [] , ParallelDo [], ...) sangat efisien dan memungkinkan Anda memproses data dalam jumlah besar. Sebagian besar, ini difasilitasi oleh fitur penyimpanan data paket GeologyIO. Omong-omong, paket itu dapat digunakan tidak hanya di bidang eksplorasi seismik terapan. Tipe data yang hampir sama digunakan dalam georadar dan seismologi. Jika Anda memiliki saran tentang cara meningkatkan hasilnya, algoritma analisis sinyal Wolfram Mathematica mana yang berlaku untuk data tersebut, atau jika Anda memiliki komentar kritis, tinggalkan komentar.