地球物理學中的 Wolfram 數學

感謝博客作者 安東·埃基門科 為他的報告

介紹

這篇筆記是在會議結束後寫的 Wolfram 俄羅斯技術會議 並包含我提供的報告的摘要。 活動於六月在聖彼得堡舉行。 考慮到我的工作地點距離會場只有一個街區，所以我忍不住參加了這次活動。 2016年、2017年我聽了會議報告，今年我做了報告。 首先，出現了一個有趣的（在我看來）主題，我們正在開發這個主題 基里爾·別洛夫其次，經過長時間研究俄羅斯聯邦有關制裁政策的立法，在我工作的企業，出現了多達兩張許可證 Wolfram Mathematica.

在進入演講主題之前，我想先介紹一下這次活動的良好組織。 會議的訪問頁面使用了喀山大教堂的圖像。 大教堂是聖彼得堡的主要景點之一，從會議舉行的大廳可以清楚地看到。

在聖彼得堡國立經濟大學入口處，學生們的助手迎接了參與者，他們不讓他們迷路。 註冊期間，我們會贈送一些小紀念品（一個玩具 - 一根閃光的釘子、一支筆、帶有 Wolfram 符號的貼紙）。 午餐和茶歇也包括在會議日程中。 我已經在小組的牆上註意到了美味的咖啡和餡餅 - 廚師很棒。 透過這個介紹部分，我想強調的是，活動本身、它的形式和地點已經帶來了正面的情緒。

我和 Kirill Belov 準備的報告名為「使用 Wolfram Mathematica 解決應用地球物理學中的問題」。 地震資料或「古代河流流經的地方」的頻譜分析。 報告內容分為兩部分：一是現有演算法的使用。 Wolfram Mathematica 用於分析地球物理數據，其次，這是如何將地球物理數據放入 Wolfram Mathematica 中。

地震勘探

首先，您需要對地球物理學進行簡短的了解。 地球物理學是研究岩石物理性質的科學。 那麼，由於岩石具有不同的性質：電、磁、彈性，因此有相應的地球物理學方法：電法勘探、磁法勘探、地震勘探……在本文中，我們將只更詳細地討論地震勘探。 地震勘探是尋找石油和天然氣的主要方法。 該方法基於彈性振動的激發以及隨後記錄構成研究區域的岩石的響應。 在陸地上（使用炸藥或彈性振動的非爆炸性振動源）或在海上（使用氣槍）激發振動。 彈性振動透過岩體傳播，在不同性質的層的邊界處折射和反射。 反射波返回表面，並由陸地上的地震檢波器（通常基於懸掛在線圈中的磁鐵的運動的電動裝置）或海洋中的水聽器（基於壓電效應）記錄。 透過波浪到達的時間，可以判斷地質層的深度。

地震船拖曳設備


氣槍激發彈性振動


波穿過岩體並被水聽器記錄


地球物理調查研究船「伊凡古布金號」停泊在聖彼得堡布拉戈維申斯基大橋附近的碼頭

地震訊號模型

岩石具有不同的物理性質。 對於地震勘探，彈性特性至關重要—彈性振動的傳播速度和密度。 如果兩層具有相同或相似的屬性，那麼波「將不會注意到」它們之間的邊界。 如果層中的波速不同，則在層的邊界處會發生反射。 性質差異越大，反射越強烈。 其強度由反射係數 (rc) 決定：

式中，ρ為岩石密度，ν為波速，1和2表示上層和下層。

最簡單且最常用的地震訊號模型之一是卷積模型，當記錄的地震道表示為反射係數序列與偵測脈衝的捲積結果時：

其中 s(t) — 地震道，即水聽器或地震檢波器在固定記錄時間內記錄的所有內容， 重量（噸） - 氣槍產生的訊號， n(t) - 隨機噪音。

讓我們以計算合成地震道為例。 我們將使用地震勘探中廣泛使用的瑞克脈衝作為初始訊號。

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]

初始地震脈衝


我們將在300 ms和600 ms的深度設定兩個邊界，反射係數將是隨機數

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

模擬賽道


對於這個例子，有必要進行保留 - 實際上，層的深度當然是以米為單位確定的，並且地震道的計算是在時域中進行的。 以公尺為單位設定深度並在了解各層速度的情況下計算到達時間會更正確。 在這種情況下，我立即在時間軸上設定圖層。

如果我們談論實地研究，那麼由於這種觀察，大量相似的時間序列（地震痕跡）就會被記錄下來。 例如，當研究一個長 25 公里、寬 15 公里的地點時，工作結果是每條跡線都表徵一個 25x25 公尺的單元（這樣的單元稱為 bin），最終的數據陣列將包含 600000 個跡線。 採樣時間為1毫秒，記錄時間為5秒，最終的資料檔案將超過11GB，原始「原始」材料的體積可達數百GB。

如何與他們合作 Wolfram Mathematica?

包 地質學IO

開始開發包 問題 在俄語支援小組的 VK 牆上。 感謝社區的回應，很快就找到了解決方案。 結果，它發展成為一項重大的發展。 相應的 Wolfram 社區牆柱 它甚至被版主標記了。 目前，該軟體套件支援使用地質產業中積極使用的以下資料類型：

1. 導入 ZMAP 和 IRAP 格式的地圖數據
2. 以 LAS 格式孔匯入測量值
3. 地震檔案格式的輸入與輸出 賽格

要安裝軟體包，您必須按照已組裝軟體包的下載頁面上的說明進行操作，即在任意位置執行下列程式碼 數學筆記本:

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 資料集是在研究所產生的…最小速度為 1500 m/s，最大速度為 5500 m/s）”

您可以透過使用「traces」鍵存取地震道來顯示實際的地質模型（該套件的功能之一是鍵不區分大小寫）：

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

模型 Marmousi


目前，該軟體包還允許您從大檔案中分批載入數據，從而可以處理大小可達數十GB的檔案。 該套件的功能還包括將資料匯出到 .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]]

尖滅地層模型


楔內的波速為 4500 m/s，楔外的波速為 4000 m/s，假設密度恆定為 2200 g/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]

楔塊上邊緣反射波振幅圖


邏輯上，當訊號頻率較低時，幹擾開始出現在大的地層厚度處，而在高頻訊號的情況下，幹擾出現在較小的厚度處。 以下程式碼片段建立頻率為 35 Hz、55 Hz 和 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]

一組頻率為 35 Hz、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 地圖衛星圖 - 諾亞布爾斯克市區


其中一個油田的井場網絡


西西伯利亞的含油岩石的深度範圍很廣——從 1 公里到 5 公里。 含油岩石的主要體積形成於侏羅紀和白堊紀。 侏羅紀時期可能是許多人透過同名電影而了解的。 侏羅紀氣候 與現代有很大不同。 大英百科全書有一系列古地圖，描繪了每個神學時代的特徵。

現在時

侏羅紀時期


請注意，在侏羅紀時期，西西伯利亞的領土是海岸（河流和淺海穿過的土地）。 由於氣候舒適，我們可以假設當時的典型景觀是這樣的：

侏羅紀西伯利亞


在這張照片中，對我們來說重要的不是動物和鳥類，而是背景中河流的圖像。 這條河與我們之前停過的地質物體是同一個。 事實上，河流的活動使得分選良好的砂岩累積起來，然後成為石油儲層。 這些水庫可能具有奇異、複雜的形狀（如河床），並且厚度可變 - 靠近河岸的厚度較小，但靠近河道中心或在蜿蜒區域的厚度會增加。 因此，侏羅紀形成的河流現在深度約為三公里，是尋找油藏的目標。

實驗數據。 處理和可視化

讓我們立即對文章中顯示的地震材料進行保留 - 由於用於分析的數據量很大 - 文章文本中僅包含原始地震道集的片段。 這將允許任何人重現上述計算。

在處理地震資料時，地球物理學家通常使用專門的軟體（有幾個行業領導者的開發成果被積極使用，例如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×5公里的區域內記錄的路線。 如果使用三維地震勘測方法獲得資料（不是沿著單獨的地球物理剖面記錄波，而是在整個區域同時記錄波），則可以獲得地震資料立方體。 這些是三維物體，其垂直和水平截面允許對地質環境進行詳細研究。 在所考慮的範例中，我們正在處理三維資料。 我們可以從文字標題中獲取一些信息，就像這樣

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

C 1 這是地質包測試的示範文件
Ç2的
Ç3的
Ç4的
C 5 日期用戶名：WOLFRAM 用戶
C 6 調查名稱：西伯利亞某處
C 7 檔案類型 3D 地震體積
Ç8的
Ç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]

地震剖面軌跡之一


知道研究區域的每個方向上有多少條跡線，您可以產生三維資料數組並使用 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]]

地震資料立方體的 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 混合。 深度約2公里。

薩拉托夫附近伏爾加河的衛星圖

結論

Wolfram Mathematica 可讓您分析地震資料並解決與礦產勘探相關的應用問題，而 GeologyIO 軟體套件讓此過程更加方便。 地震資料的結構是這樣的，使用內建方法來加速計算（並行表[], 並行[],...) 非常高效，允許您處理大量資料。 在很大程度上，GeologyIO 套件的資料儲存功能促進了這一點。 順便說一句，該軟體包不僅可以用於應用地震勘探領域。 探地雷達和地震學中使用的數據類型幾乎相同。如果您對如何改進結果有建議，Wolfram Mathematica 武器庫中的哪些信號分析演算法適用於此類數據，或者您有任何批評意見，請發表評論。

來源： www.habr.com