지구물리학 분야의 Wolfram Mathematica

블로그 작성자에게 감사드립니다. 안톤 에키멘코 그의 보고를 위해

소개

이 메모는 회의가 끝난 후 작성되었습니다. Wolfram 러시아 기술 컨퍼런스 그리고 내가 제공한 보고서의 요약이 포함되어 있습니다. 이 행사는 2016월 상트페테르부르크에서 열렸습니다. 제가 컨퍼런스 현장에서 한 블록 떨어진 곳에 근무한다는 점을 고려하면 이번 행사에는 참석하지 않을 수 없었습니다. 2017년과 XNUMX년에는 컨퍼런스 보고를 듣고, 올해는 발표를 했습니다. 첫째, 우리가 개발 중인 흥미로운 주제가 나타났습니다. 키릴 벨로프, 둘째, 제재 정책에 관한 러시아 연방 법률에 대한 오랜 연구 끝에 제가 일하는 기업에서 최대 XNUMX개의 라이센스가 나타났습니다. Wolfram Mathematica.

연설 주제로 넘어가기 전에, 행사의 좋은 조직에 대해 말씀드리고 싶습니다. 회의의 방문 페이지는 카잔 대성당의 이미지를 사용합니다. 대성당은 상트페테르부르크의 주요 명소 중 하나이며 회의가 열렸던 홀에서 매우 명확하게 보입니다.

상트 페테르부르크 주립 경제 대학 입구에서 참가자들은 학생들 중 조수들과 만났으며 길을 잃지 않았습니다. 등록하는 동안 작은 기념품이 제공되었습니다(장난감-번쩍이는 스파이크, 펜, Wolfram 기호가 있는 스티커). 점심 식사와 휴식 시간도 회의 일정에 포함되었습니다. 나는 이미 그룹 벽에 맛있는 커피와 파이에 대해 언급했습니다. 요리사는 훌륭합니다. 이번 소개 부분에서는 행사 자체와 형식, 장소가 이미 긍정적인 감정을 불러일으키고 있다는 점을 강조하고 싶다.

나와 Kirill Belov가 작성한 보고서는 "Wolfram Mathematica를 사용하여 응용 지구물리학 문제 해결"이라는 제목입니다. 지진 데이터 또는 "고대 강이 흐르던 곳"에 대한 스펙트럼 분석. 보고서의 내용은 두 부분으로 구성됩니다. 첫째, 다음에서 사용할 수 있는 알고리즘의 사용입니다. Wolfram Mathematica 지구물리학적 데이터를 분석하는 방법, 두 번째로 Wolfram Mathematica에 지구물리학적 데이터를 입력하는 방법입니다.

지진 탐사

먼저 지구물리학에 대한 짧은 여행을 해야 합니다. 지구물리학은 암석의 물리적 특성을 연구하는 과학입니다. 글쎄요, 암석은 전기적, 자기적, 탄성적 특성을 가지기 때문에 지구물리학에는 그에 상응하는 전기 탐사, 자기 탐사, 지진 탐사 등의 방법이 있습니다. 이 기사의 맥락에서는 지진 탐사에 대해서만 더 자세히 논의하겠습니다. 지진 탐사는 석유와 가스를 찾는 주요 방법입니다. 이 방법은 탄성 진동의 여기와 연구 지역을 구성하는 암석의 반응에 대한 후속 기록을 기반으로 합니다. 진동은 육지(다이너마이트 또는 비폭발성 탄성 진동원 사용) 또는 바다(공기총 사용)에서 발생합니다. 탄성 진동은 암반을 통해 전파되며, 서로 다른 특성을 갖는 층의 경계에서 굴절 및 반사됩니다. 반사된 파동은 표면으로 돌아와 육지의 지오폰(보통 코일에 매달린 자석의 움직임을 기반으로 하는 전기 역학 장치) 또는 바다의 수중음기(압전 효과 기반)에 의해 기록됩니다. 파도가 도달할 때쯤이면 지질층의 깊이를 판단할 수 있다.

지진 선박 견인 장비


에어건은 탄성 진동을 자극합니다.


파도는 암석 덩어리를 통과하며 수중음파로 녹음됩니다.


상트페테르부르크 블라고베셴스키 다리 근처 부두에 있는 지구물리학 연구 선박 "Ivan Gubkin"

지진 신호 모델

암석은 다양한 물리적 특성을 가지고 있습니다. 지진 탐사에서는 탄성 특성, 즉 탄성 진동의 전파 속도와 밀도가 주로 중요합니다. 두 레이어가 동일하거나 유사한 속성을 갖는 경우 파동은 두 레이어 사이의 경계를 "알지 못합니다". 층의 파동 속도가 다르면 층의 경계에서 반사가 발생합니다. 속성의 차이가 클수록 반사가 더 강해집니다. 강도는 반사 계수(rc)에 의해 결정됩니다.

여기서 ρ는 암석 밀도, ν는 파동 속도, 1과 2는 상층과 하층을 나타냅니다.

가장 간단하고 가장 자주 사용되는 지진 신호 모델 중 하나는 기록된 지진 추적이 프로빙 펄스와 일련의 반사 계수의 컨볼루션 결과로 표시되는 컨볼루션 모델입니다.

어디 s(t) — 지진 흔적, 즉 정해진 녹음시간 동안 수중음파나 지음기로 녹음한 모든 것, 승(티) - 공기총에 의해 생성된 신호, n(티) - 무작위 소음.

예를 들어 합성 지진 흔적을 계산해 보겠습니다. 지진 탐사에 널리 사용되는 리커 펄스(Ricker pulse)를 초기 신호로 사용하겠습니다.

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]

초기 지진 충격량


300ms와 600ms의 깊이에 두 개의 경계를 설정하고 반사 계수는 난수입니다.

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]]

반사 계수의 시퀀스


지진 흔적을 계산하여 표시해 봅시다. 반사 계수의 부호가 다르기 때문에 지진 추적에서 두 개의 교번 반사를 얻습니다.

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]]

시뮬레이션된 트랙


이 예에서는 예약이 필요합니다. 실제로 레이어의 깊이는 물론 미터 단위로 결정되며 지진 추적 계산은 시간 영역에 대해 발생합니다. 깊이를 미터 단위로 설정하고 레이어의 속도를 알고 도착 시간을 계산하는 것이 더 정확할 것입니다. 이 경우에는 즉시 시간축에 레이어를 설정합니다.

현장 조사에 관해 이야기하면 그러한 관찰의 결과로 유사한 시계열(지진 흔적)이 엄청나게 많이 기록됩니다. 예를 들어 길이 25km, 너비 15km의 사이트를 연구할 때 작업 결과 각 트레이스가 25x25미터 크기의 셀(이러한 셀을 빈이라고 함)의 특성을 나타내면 최종 데이터 배열에는 600000개의 트레이스가 포함됩니다. 샘플링 시간이 1ms이고 녹음 시간이 5초인 경우 최종 데이터 파일은 11GB 이상이 되며 원본 "원시" 자료의 용량은 수백 기가바이트에 이를 수 있습니다.

그들과 협력하는 방법 Wolfram Mathematica?

패키지 지질학IO

패키지 개발이 시작되었습니다 문제 러시아어를 사용하는 지원 그룹의 VK 벽에 있습니다. 커뮤니티의 반응 덕분에 해결책이 매우 빨리 발견되었습니다. 그리고 그 결과 심각한 발전으로 발전했습니다. 동 Wolfram 커뮤니티 벽 게시물 심지어 중재자에 의해 표시되었습니다. 현재 이 패키지는 지질학 산업에서 활발히 사용되는 다음 데이터 유형 작업을 지원합니다.

1. ZMAP 및 IRAP 형식의 지도 데이터 가져오기
2. LAS 형식 웰의 측정 가져오기
3. 지진 파일 형식의 입력 및 출력 SEGY

패키지를 설치하려면 조립된 패키지의 다운로드 페이지에 있는 지침을 따라야 합니다. 다음 코드를 실행하십시오. Mathematica 노트북:

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

그런 다음 패키지가 기본 폴더에 설치되며 경로는 다음과 같습니다.

FileNameJoin[{$UserBasePacletsDirectory, "Repository"}]

예를 들어 패키지의 주요 기능을 보여드리겠습니다. 호출은 전통적으로 Wolfram 언어의 패키지에 대해 수행됩니다.

Get["GeologyIO`"]

패키지는 다음을 사용하여 개발되었습니다. 볼프람 워크벤치. 이를 통해 패키지의 주요 기능과 함께 문서 형식을 제공할 수 있으며, 이 문서는 Wolfram Mathematica 자체 문서와 다르지 않으며, 처음 접하는 사람을 위해 패키지에 테스트 파일을 제공할 수 있습니다.

특히 이러한 파일은 "Marmousi.segy" 파일입니다. 이는 프랑스 석유 연구소에서 개발한 지질 단면의 합성 모델입니다. 이 모델을 사용하여 개발자는 파동장 모델링, 데이터 처리, 지진 추적 역전 등에 대한 자체 알고리즘을 테스트합니다. Marmousi 모델 자체는 패키지 자체가 다운로드된 저장소에 저장됩니다. 파일을 얻으려면 다음 코드를 실행하십시오.

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

가져오기 결과 - SEGYData 객체


SEGY 형식에는 관찰에 대한 다양한 정보를 저장하는 작업이 포함됩니다. 첫째, 이것은 텍스트 주석입니다. 여기에는 작업 위치, 측정을 수행한 회사 이름 등에 대한 정보가 포함됩니다. 우리의 경우 이 헤더는 TextHeader 키가 있는 요청에 의해 호출됩니다. 다음은 단축된 텍스트 헤더입니다.

Short[marmousi["TextHeader"]]

"Marmousi 데이터 세트는 연구소에서 생성되었습니다. 최대 속도는 1500m/s이고 최대 속도는 5500m/s입니다."

"trace" 키를 사용하여 지진 추적에 액세스하여 실제 지질 모델을 표시할 수 있습니다(패키지의 기능 중 하나는 키가 대소문자를 구분하지 않는다는 것입니다).

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

모델 마르무시


현재 이 패키지를 사용하면 대용량 파일의 데이터를 부분적으로 로드할 수 있으므로 크기가 수십 기가바이트에 달하는 파일을 처리할 수 있습니다. 패키지의 기능에는 데이터를 .segy로 내보내고 부분적으로 파일 끝에 추가하는 기능도 포함되어 있습니다.

이와 별도로 .segy 파일의 복잡한 구조로 작업할 때 패키지의 기능에 주목할 가치가 있습니다. 키와 인덱스를 사용하여 개별 트레이스와 헤더에 액세스할 수 있을 뿐만 아니라 이를 변경한 다음 파일에 쓸 수 있기 때문입니다. GeologyIO 구현에 대한 많은 기술적 세부 사항은 이 기사의 범위를 벗어나므로 별도로 설명할 가치가 있을 것입니다.

지진 탐사에서 스펙트럼 분석의 관련성

지진 데이터를 Wolfram Mathematica로 가져오는 기능을 통해 내장된 신호 처리 기능을 실험 데이터에 사용할 수 있습니다. 각 지진 추적은 시계열을 나타내기 때문에 이를 연구하는 주요 도구 중 하나는 스펙트럼 분석입니다. 지진 데이터의 주파수 구성을 분석하기 위한 전제 조건 중 예를 들어 다음을 들 수 있습니다.

1. 다양한 유형의 파동은 서로 다른 주파수 구성을 특징으로 합니다. 이를 통해 유용한 파동을 강조하고 간섭파를 억제할 수 있습니다.
2. 다공성 및 포화도와 같은 암석 특성은 주파수 구성에 영향을 미칠 수 있습니다. 이를 통해 최상의 특성을 지닌 암석을 식별할 수 있습니다.
3. 두께가 다른 레이어는 다양한 주파수 범위에서 이상 현상을 발생시킵니다.

세 번째 요점은 이 기사의 맥락에서 주요한 내용입니다. 다음은 두께가 다양한 레이어(쐐기 모델)의 경우 지진 흔적을 계산하기 위한 코드 조각입니다. 이 모델은 전통적으로 여러 층에서 반사된 파동이 서로 중첩될 때 간섭 효과를 분석하기 위해 지진 탐사에서 연구되었습니다.

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]]

핀치아웃 형성 모델


쐐기 내부의 파동 속도는 4500m/s이고 쐐기 외부의 파동 속도는 4000m/s이며 밀도는 2200g/cmXNUMX로 일정하다고 가정됩니다. 이러한 모델의 경우 반사 계수와 지진 흔적을 계산합니다.

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]

쐐기형 모델의 지진 흔적


이 그림에 표시된 일련의 지진 흔적을 지진 단면이라고 합니다. 보시다시피, 반사파의 기하학적 구조가 이전에 지정한 모델과 명확하게 일치하므로 해석은 직관적인 수준에서 수행될 수도 있습니다. 추적을 더 자세히 분석하면 1에서 약 30까지의 추적이 다르지 않다는 것을 알 수 있습니다. 구조물의 지붕과 바닥의 반사가 서로 겹치지 않습니다. 31번째 트레이스부터 반사가 간섭하기 시작합니다. 그리고 모델에서는 반사 계수가 수평으로 변경되지 않습니다. 지진 흔적은 지층의 두께가 변경됨에 따라 강도가 변경됩니다.

지층의 상부 경계로부터의 반사 진폭을 고려해 봅시다. 60번째 경로부터 반사 강도가 증가하기 시작하고 70번째 경로에서 최대가 됩니다. 이는 지층의 지붕과 바닥에서 발생하는 파도의 간섭이 나타나는 방식으로, 일부 경우에는 지진 기록에 심각한 이상 현상이 발생합니다.

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]

쐐기의 상단 가장자리에서 반사파의 진폭 그래프


신호의 주파수가 낮을수록 큰 형성 두께에서 간섭이 나타나기 시작하고, 고주파 신호의 경우 더 작은 두께에서 간섭이 발생하는 것이 논리적입니다. 다음 코드 조각은 35Hz, 55Hz, 85Hz 주파수의 신호를 생성합니다.

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]

35Hz, 55Hz, 85Hz의 주파수를 갖는 소스 신호 세트


지진 흔적을 계산하고 반사파 진폭의 그래프를 작성함으로써 서로 다른 주파수에 대해 서로 다른 형성 두께에서 이상이 관찰된다는 것을 알 수 있습니다.

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]

다양한 주파수에 대한 쐐기의 상단 가장자리에서 반사파의 진폭 그래프


지진 관측 결과로부터 지층의 두께에 대한 결론을 도출하는 능력은 매우 유용합니다. 왜냐하면 석유 탐사의 주요 작업 중 하나는 유정을 놓기 위한 가장 유망한 지점(즉, 지층이 있는 지역)을 평가하는 것이기 때문입니다. 더 두꺼운). 또한 지질 단면에는 생성으로 인해 형성 두께가 급격히 변하는 물체가 있을 수 있습니다. 이는 스펙트럼 분석을 연구하는 데 효과적인 도구로 만듭니다. 기사의 다음 부분에서는 그러한 지질학적 물체를 더 자세히 고려할 것입니다.

실험 데이터. 어디서 얻었으며 그 안에서 무엇을 찾아야 합니까?

기사에서 분석된 자료는 서부 시베리아에서 얻은 것입니다. 이 지역은 누구나 예외 없이 알고 있듯이 우리나라의 주요 산유지역이다. 지난 세기 60년대에 이 지역에서 활발한 퇴적물 개발이 시작되었습니다. 석유 매장지를 찾는 주요 방법은 지진 탐사입니다. 이 지역의 위성 이미지를 보는 것은 흥미 롭습니다. 작은 규모에서는 엄청난 수의 늪과 호수를 볼 수 있으며, 지도를 확대하면 클러스터 우물 시추 장소를 볼 수 있고, 지도를 한계까지 확대하면 지진이 발생한 프로필의 개간도 구분할 수 있습니다. 관찰이 수행되었습니다.

Yandex 지도의 위성 이미지 - Noyabrsk 도시 지역


들판 중 하나에 있는 우물 패드 네트워크


서부 시베리아의 석유 함유 암석은 1km에서 5km까지 다양한 깊이에서 발견됩니다. 석유를 함유한 암석의 대부분은 쥐라기와 백악기에 형성되었습니다. 쥬라기 시대는 아마도 같은 이름의 영화를 통해 많은 사람들에게 알려졌을 것입니다. 쥐라기 기후 현대와는 확연히 달랐습니다. 브리태니커 백과사전에는 각 생물 시대를 특징짓는 일련의 고지 지도가 있습니다.

요즘

쥐라기


쥐라기 시대에 서부 시베리아의 영토는 해안(강과 얕은 바다가 교차하는 땅)이었다는 점에 유의하세요. 기후가 편안했기 때문에 당시의 전형적인 풍경은 다음과 같았을 것으로 추측할 수 있습니다.

쥐라기 시베리아


이 사진에서 우리에게 중요한 것은 동물이나 새보다는 배경에 있는 강의 이미지입니다. 강은 우리가 앞서 멈췄던 것과 동일한 지질학적 대상입니다. 사실 강의 활동으로 인해 잘 분류된 사암이 축적되어 석유 저장소가 됩니다. 이러한 저수지는 기이하고 복잡한 모양(강바닥과 같은)을 가질 수 있으며 두께가 다양합니다. 제방 근처에서는 두께가 작지만 수로 중앙에 가까우거나 구불구불한 지역에서는 두께가 증가합니다. 따라서 쥐라기에 형성된 강은 현재 약 XNUMXkm 깊이에 있으며 석유 저장소를 찾는 대상입니다.

실험 데이터. 처리 및 시각화

기사에 표시된 지진 자료에 대해 즉시 예약하겠습니다. 분석에 사용된 데이터의 양이 상당하기 때문에 원본 지진 흔적 세트의 일부만 기사의 텍스트에 포함되어 있습니다. 이를 통해 누구나 위의 계산을 재현할 수 있습니다.

지진 데이터로 작업할 때 지구물리학자는 일반적으로 다양한 유형의 데이터를 분석할 수 있고 편리한 그래픽 인터페이스를 제공하는 특수 소프트웨어(Petrel 또는 Paradigm과 같이 개발이 적극적으로 사용되는 여러 업계 리더가 있음)를 사용합니다. 모든 편리함에도 불구하고 이러한 유형의 소프트웨어에는 단점도 있습니다. 예를 들어 안정적인 버전에서 최신 알고리즘을 구현하는 데는 많은 시간이 걸리고 계산 자동화 가능성은 일반적으로 제한됩니다. 이러한 상황에서는 광범위한 알고리즘 기반을 사용하는 동시에 많은 루틴을 수행할 수 있는 컴퓨터 수학 시스템과 고급 프로그래밍 언어를 사용하는 것이 매우 편리해집니다. 이것이 Wolfram Mathematica에서 지진 데이터를 작업하는 데 사용되는 원리입니다. 데이터를 사용한 대화형 작업을 위해 풍부한 기능을 작성하는 것은 부적절합니다. 일반적으로 허용되는 형식에서 로드하고, 원하는 알고리즘을 적용하고, 외부 형식으로 다시 업로드하는 것이 더 중요합니다.

제안된 방식에 따라 원본 지진 데이터를 로드하고 이를 표시합니다. 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]

이렇게 다운로드하고 가져온 데이터는 10 x 5km 크기의 영역에 기록된 경로입니다. XNUMX차원 지진조사 방식(파동이 개별 지구물리학적 단면이 아닌 전체 지역에 걸쳐 동시에 기록됨)을 사용하여 데이터를 얻으면 지진 데이터 큐브를 얻는 것이 가능해집니다. 이는 지질 환경에 대한 상세한 연구를 가능하게 하는 수직 및 수평 단면의 XNUMX차원 물체입니다. 고려된 예에서 우리는 XNUMX차원 데이터를 다루고 있습니다. 텍스트 헤더에서 다음과 같은 정보를 얻을 수 있습니다.

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

C 1 이것은 GEOLOGYIO 패키지 테스트용 데모 파일입니다.
C 2
C 3
C 4
C 5 날짜 사용자 이름: WOLFRAM 사용자
C 6 설문 조사 이름: 시베리아 어딘가
C 7 파일 유형 3차원 지진량
C 8
C 9
C10 Z 범위: 처음 2200M 마지막 2400M

이 데이터 세트는 데이터 분석의 주요 단계를 보여주기에 충분합니다. 파일의 트레이스는 순차적으로 기록되며 각 트레이스는 다음 그림과 유사합니다. 이는 수직 축(깊이 축)을 따른 반사파의 진폭 분포입니다.

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]

지진단면 흔적 중 하나


연구 영역의 각 방향에 몇 개의 트레이스가 있는지 알면 3차원 데이터 배열을 생성하고 ImageXNUMXD[] 함수를 사용하여 표시할 수 있습니다.

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]]

지진 데이터 큐브의 XNUMXD 이미지(세로축 - 깊이)


관심 있는 지질학적 특징이 강렬한 지진 이상 현상을 생성하는 경우 투명성을 갖춘 시각화 도구를 사용할 수 있습니다. 녹화의 "중요하지 않은" 영역은 보이지 않게 만들어 예외적인 부분만 보이게 할 수 있습니다. Wolfram Mathematica에서는 다음을 사용하여 이 작업을 수행할 수 있습니다. 불투명[] и 래스터3D[].

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]

Opacity[] 및 Raster3D[] 함수를 사용한 지진 데이터 큐브 이미지

합성 예에서와 같이 원래 입방체의 단면에서 가변 기복이 있는 일부 지질학적 경계(층)를 식별할 수 있습니다.

스펙트럼 분석의 주요 도구는 푸리에 변환입니다. 이를 통해 각 트레이스 또는 트레이스 그룹의 진폭-주파수 스펙트럼을 평가할 수 있습니다. 그러나 데이터를 주파수 영역으로 전송한 후에는 주파수가 언제(어떤 깊이에서 읽혀지는지) 변경되는지에 대한 정보가 손실됩니다. 시간(깊이) 축에서 신호 변화를 국지화할 수 있도록 윈도우 푸리에 변환 및 웨이블릿 분해가 사용됩니다. 이 기사에서는 웨이블릿 분해를 사용합니다. 웨이블릿 분석 기술은 90년대부터 지진탐사에 활발히 활용되기 시작했다. 윈도우 푸리에 변환에 비해 장점은 더 나은 시간 분해능으로 간주됩니다.

다음 코드 조각을 사용하면 지진 추적 중 하나를 개별 구성 요소로 분해할 수 있습니다.

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]

추적을 구성 요소로 분해


서로 다른 파동 도달 시간에 반사 에너지가 어떻게 분포되는지 평가하기 위해 스케일로그램(스펙트로그램과 유사)이 사용됩니다. 일반적으로 실제로 모든 구성 요소를 분석할 필요는 없습니다. 일반적으로 저, 중, 고주파 성분이 선택됩니다.

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]

스케일로그램. 함수 결과 웨이블릿스케일로그램[]


Wolfram 언어는 웨이블릿 변환 기능을 사용합니다. 연속웨이브렛변형[]. 그리고 전체 추적 세트에 이 기능을 적용하는 것은 다음 기능을 사용하여 수행됩니다. 테이블[]. 여기에서 Wolfram Mathematica의 강점 중 하나인 병렬화 기능을 주목할 가치가 있습니다. 병렬 테이블[]. 위의 예에서는 병렬화가 필요하지 않습니다. 데이터의 양은 크지 않지만 수십만 개의 추적이 포함된 실험 데이터 세트로 작업할 때 병렬화가 필요합니다.

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

기능을 적용한 후 연속웨이브렛변형[] 선택한 빈도에 해당하는 새 데이터 세트가 나타납니다. 위의 예에서 이러한 주파수는 38Hz, 33Hz, 27Hz입니다. 주파수 선택은 테스트를 기반으로 가장 자주 수행됩니다. 다양한 주파수 조합에 대한 효과적인 맵을 얻고 지질 학자의 관점에서 가장 유익한 맵을 선택합니다.

결과를 동료와 공유하거나 고객에게 제공해야 하는 경우 GeologyIO 패키지의 SEGYExport[] 기능을 사용할 수 있습니다.

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];

이러한 큐브 중 세 개(저주파수, 중간 주파수 및 고주파수 구성 요소)를 사용하면 일반적으로 RGB 블렌딩을 사용하여 데이터를 함께 시각화합니다. 각 구성 요소에는 빨간색, 녹색, 파란색 등 고유한 색상이 지정됩니다. Wolfram Mathematica에서는 다음 함수를 사용하여 이 작업을 수행할 수 있습니다. 색상결합[].

그 결과는 지질학적 해석이 가능한 이미지입니다. 해당 단면에 기록된 구불구불한 지형을 통해 저수지일 가능성이 높고 석유 매장량이 포함된 고대 수로의 윤곽을 그릴 수 있습니다. 그러한 하천 시스템의 현대 유사체를 검색하고 분석하면 구불구불한 구불구불한 부분의 가장 유망한 부분을 결정할 수 있습니다. 수로 자체는 잘 분류된 사암층이 두껍게 쌓인 것이 특징이며 석유를 저장하는 좋은 장소입니다. "레이스" 이상 현상 외부의 지역은 현대의 범람원 퇴적물과 유사합니다. 범람원 퇴적물은 주로 점토질 암석으로 이루어져 있으며 이러한 지역에 대한 시추 작업은 효과적이지 않습니다.

데이터 큐브의 RGB 슬라이스. 중앙(중앙에서 약간 왼쪽)에서는 구불구불한 강을 따라갈 수 있습니다.

데이터 큐브의 RGB 슬라이스. 왼쪽에는 구불구불한 강을 따라가 볼 수 있습니다.


어떤 경우에는 지진 데이터의 품질을 통해 훨씬 더 선명한 이미지를 얻을 수 있습니다. 이는 현장 작업 방법론, 소음 감소 알고리즘에 사용되는 장비에 따라 다릅니다. 이러한 경우 하천 시스템의 일부만 볼 수 있는 것이 아니라 확장된 고강 전체도 볼 수 있습니다.

지진 데이터 큐브(수평 슬라이스)의 세 가지 구성 요소를 RGB로 혼합합니다. 깊이는 약 2km.

사라토프 인근 볼가 강의 위성 이미지

결론

Wolfram Mathematica를 사용하면 지진 데이터를 분석하고 광물 탐사와 관련된 응용 문제를 해결할 수 있으며, GeologyIO 패키지를 사용하면 이 프로세스가 더욱 편리해집니다. 지진 데이터의 구조는 내장된 방법을 사용하여 계산 속도를 높이는 구조입니다(병렬 테이블[], 병렬 실행[],…)은 매우 효율적이며 많은 양의 데이터를 처리할 수 있습니다. 이는 대부분 GeologyIO 패키지의 데이터 저장 기능을 통해 촉진됩니다. 그런데 이 패키지는 응용 지진 탐사 분야에만 사용할 수 있는 것이 아닙니다. 지표 투과 레이더 및 지진학에서도 거의 동일한 유형의 데이터가 사용됩니다. 결과를 개선하는 방법에 대한 제안이 있거나 Wolfram Mathematica 무기고의 어떤 신호 분석 알고리즘이 이러한 데이터에 적용 가능하거나 비판적인 의견이 있는 경우 코멘트를 남겨주세요.

출처 : habr.com