וולפרם מתמטיקה בגיאופיזיקה

תודה לכותבת הבלוג אנטון אקימנקו לדוח שלו

מבוא

פתק זה נכתב בעקבות הוועידה ועידת טכנולוגיה רוסית של וולפרם ומכיל תקציר של הדוח שנתתי. האירוע התקיים ביוני בסנט פטרבורג. בהתחשב בכך שאני עובד בלוק מאתר הכנס, לא יכולתי שלא להשתתף באירוע הזה. ב-2016 וב-2017 הקשבתי לדוחות הכנס, והשנה העברתי מצגת. ראשית, הופיע נושא מעניין (נראה לי) שאיתו אנו מפתחים קיריל בלוב, ושנית, לאחר מחקר ארוך של החקיקה של הפדרציה הרוסית בנוגע למדיניות הסנקציות, במפעל שבו אני עובד, הופיעו עד שני רישיונות וולפרם מתמטיקה.

לפני שנעבור לנושא דברי, ברצוני לציין את הארגון הטוב של האירוע. עמוד הביקור של הכנס משתמש בתמונה של קתדרלת קאזאן. הקתדרלה היא אחת האטרקציות המרכזיות של סנט פטרסבורג והיא נראית היטב מהאולם בו התקיים הכנס.

בכניסה לאוניברסיטה הכלכלית הממלכתית של סנט פטרבורג פגשו את המשתתפים עוזרים מקרב הסטודנטים - הם לא אפשרו להם ללכת לאיבוד. במהלך ההרשמה חולקו מזכרות קטנות (צעצוע - ספייק מהבהב, עט, מדבקות עם סמלי וולפרם). גם הפסקות צהריים וקפה נכללו בלוח הזמנים של הכנס. כבר ציינתי לגבי קפה טעים ופשטידות על קיר הקבוצה - השפים נהדרים. עם חלק ההיכרות הזה, אני רוצה להדגיש שהאירוע עצמו, הפורמט והמיקום שלו כבר מעלים רגשות חיוביים.

הדו"ח שהוכן על ידי ועל ידי קיריל בלוב נקרא "שימוש ב-Wolfram Mathematica לפתרון בעיות בגיאופיזיקה יישומית. ניתוח ספקטרלי של נתונים סיסמיים או "היכן זרמו נהרות עתיקים". תוכן הדוח מכסה שני חלקים: ראשית, השימוש באלגוריתמים הזמינים ב וולפרם מתמטיקה לניתוח נתונים גיאופיזיים, ושנית, כך מכניסים נתונים גיאופיזיים ל-Wolfram Mathematica.

חיפוש סיסמי

ראשית עליך לעשות טיול קצר לגיאופיזיקה. גיאופיזיקה היא המדע החוקר את התכונות הפיזיקליות של סלעים. ובכן, מכיוון שלסלעים יש תכונות שונות: חשמליות, מגנטיות, אלסטיות, קיימות שיטות מתאימות לגיאופיזיקה: חיפוש חשמלי, חיפוש מגנטי, חיפוש סיסמי... בהקשר של מאמר זה, נדון רק בחיפוש סיסמי ביתר פירוט. חיפוש סיסמי הוא השיטה העיקרית לחיפוש נפט וגז. השיטה מבוססת על עירור תנודות אלסטיות ורישום לאחר מכן של התגובה מהסלעים המרכיבים את אזור המחקר. רעידות מתרגשות ביבשה (עם דינמיט או מקורות רטט לא נפיצים של רעידות אלסטיות) או בים (עם רובי אוויר). רעידות אלסטיות מתפשטות דרך מסת הסלע, נשברות ומשתקפות בגבולות של שכבות בעלות תכונות שונות. גלים מוחזרים חוזרים לפני השטח ומתועדים על ידי גיאופונים ביבשה (בדרך כלל מכשירים אלקטרודינמיים המבוססים על תנועת מגנט התלוי בסליל) או הידרופונים בים (על בסיס האפקט הפיאזואלקטרי). עד למועד הגעת הגלים, ניתן לשפוט את עומקן של שכבות גיאולוגיות.

ציוד לגרירת כלי שיט סיסמיים


רובה האוויר מעורר רעידות אלסטיות


הגלים עוברים דרך מסת הסלע ומתועדים באמצעות הידרופונים


ספינת מחקר גיאופיזית "איבן גובקין" במזח ליד גשר בלגוובשצ'נסקי בסנט פטרסבורג

מודל אות סייסמי

לסלעים תכונות פיזיקליות שונות. עבור חקר סיסמי, תכונות אלסטיות חשובות בעיקר - מהירות ההתפשטות של רעידות וצפיפות אלסטיות. אם לשתי שכבות יש תכונות זהות או דומות, אז הגל "לא יבחין" בגבול ביניהן. אם מהירויות הגלים בשכבות שונות, אזי תתרחש השתקפות בגבול השכבות. ככל שההבדל בתכונות גדול יותר, ההשתקפות חזקה יותר. עוצמתו תיקבע על ידי מקדם ההחזר (rc):

כאשר ρ היא צפיפות הסלע, ν היא מהירות הגל, 1 ו-2 מציינים את השכבה העליונה והתחתונה.

אחד המודלים הפשוטים והנפוצים ביותר של אותות סיסמיים הוא מודל הקונבולציה, כאשר העקיבה הסייסמית המוקלטת מיוצגת כתוצאה מהקונבולולוציה של רצף של מקדמי השתקפות עם דופק גישוש:

איפה s(t) - עקבות סיסמיים, כלומר. כל מה שהוקלט על ידי הידרופון או גיאופון בזמן הקלטה קבוע, w(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 ו-600 אלפיות השנייה, ומקדמי ההשתקפות יהיו מספרים אקראיים

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 ms וזמן הקלטה של ​​5 שניות, קובץ הנתונים הסופי יהיה יותר מ-11 GB, ונפח חומר הגלם המקורי יכול להיות מאות גיגה-בייט.

איך לעבוד איתם וולפרם מתמטיקה?

חבילה גיאולוגיהIO

החל פיתוח החבילה שאלה על קיר VK של קבוצת התמיכה דוברת הרוסית. הודות לתגובות הקהילה, נמצא פתרון מהר מאוד. וכתוצאה מכך, זה גדל להתפתחות רצינית. תוֹאֵם פוסט קיר של קהילת וולפראם זה אפילו סומן על ידי מנחים. נכון לעכשיו, החבילה תומכת בעבודה עם סוגי הנתונים הבאים הנמצאים בשימוש פעיל בתעשייה הגיאולוגית:

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 מ"ש ומקסימום של 5500 מ"ש)"

אתה יכול להציג את המודל הגיאולוגי בפועל על ידי גישה לעקבות הסייסמיות באמצעות מקש "עקבות" (אחת התכונות של החבילה היא שהמפתחות אינם רגישים לאותיות גדולות):

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

דגם של תצורת צביטה


מהירות הגלים בתוך הטריז היא 4500 מ"ש, מחוץ לטריז 4000 מ"ש, וצפיפות ההנחה היא קבועה של 2200 גרם/ס"מ³. עבור מודל כזה, אנו מחשבים מקדמי השתקפות ועקבות סיסמיות.

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 הרץ, 55 הרץ ו-85 הרץ.

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 הרץ, 55 הרץ, 85 הרץ


על ידי חישוב עקבות סיסמיות ושרטוט גרפים של משרעות גל מוחזר, אנו יכולים לראות כי עבור תדרים שונים נצפית אנומליה בעובי תצורה שונים.

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. לא מתאים לכתוב פונקציונליות עשירה לעבודה אינטראקטיבית עם נתונים – חשוב יותר להקפיד על טעינה מפורמט מקובל, החלת האלגוריתמים הרצויים עליהם והעלאתם חזרה לפורמט חיצוני.

בעקבות הסכימה המוצעת, נטען את הנתונים הסיסמיים המקוריים ונציג אותם וולפרם מתמטיקה:

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 זהו קובץ הדגמה לבדיקת חבילת GEOLOGYIO
ג 2
ג 3
ג 4
C 5 תאריך שם משתמש: WOLFRAM USER
C 6 שם הסקר: אי שם בסיביר
C 7 נפח סיסמי תלת מימדי מסוג קובץ
ג 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]]

תמונה תלת מימדית של קוביית נתונים סיסמיים. (ציר אנכי - עומק)


אם המאפיינים הגיאולוגיים המעניינים יוצרים חריגות סיסמיות אינטנסיביות, ניתן להשתמש בכלי הדמיה עם שקיפות. ניתן להפוך אזורים "לא חשובים" בהקלטה לבלתי נראים, ולהשאיר רק חריגות גלויות. ב-Wolfram Mathematica ניתן לעשות זאת באמצעות אֲטִימוּת[] и 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]

תמונת קוביית נתונים סיסמיים באמצעות פונקציות 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]

סקלוגרמה. תוצאת פונקציה WaveletScalogram[]


שפת וולפרם משתמשת בפונקציה לטרנספורמציה של גלים ContinuousWveletTransform[]. והיישום של פונקציה זו על כל סט העקבות יתבצע באמצעות הפונקציה שולחן[]. כאן כדאי לשים לב לאחת מנקודות החוזק של Wolfram Mathematica - היכולת להשתמש בהקבלה טבלה מקבילה[]. בדוגמה שלעיל אין צורך בהקבלה – נפח הנתונים אינו גדול, אך כאשר עובדים עם מערכי נתונים ניסיוניים המכילים מאות אלפי עקבות, זהו הכרחי.

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

לאחר החלת הפונקציה ContinuousWveletTransform[] ערכות נתונים חדשות מופיעות התואמות לתדרים שנבחרו. בדוגמה למעלה, התדרים הללו הם: 38Hz, 33Hz, 27Hz. בחירת התדרים מתבצעת לרוב על בסיס בדיקות - הם משיגים מפות יעילות לשילובי תדרים שונים ובוחרים את האינפורמטיבית ביותר מנקודת מבטו של גיאולוג.

אם אתה צריך לשתף את התוצאות עם עמיתים או לספק אותן ללקוח, אתה יכול להשתמש בפונקציה SEGYExport[] של חבילת 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];

עם שלוש מהקוביות הללו (רכיבי תדר נמוך, תדר בינוני ותדר גבוה), מיזוג RGB משמש בדרך כלל כדי להמחיש את הנתונים יחד. לכל רכיב מוקצה צבע משלו - אדום, ירוק, כחול. ב-Wolfram Mathematica ניתן לעשות זאת באמצעות הפונקציה ColorCombine[].

התוצאה היא תמונות שמהן ניתן לעשות פרשנות גיאולוגית. הפיתולים המתועדים בקטע מאפשרים לשרטט ערוצים פליאו-ערוצים, שסביר יותר שהם מאגרים ומכילים מאגרי נפט. החיפוש והניתוח של אנלוגים מודרניים של מערכת נהר כזו מאפשרים לנו לקבוע את החלקים המבטיחים ביותר של הפיתולים. הערוצים עצמם מתאפיינים בשכבות עבות של אבן חול ממוינת היטב ומהווים מאגר טוב לנפט. אזורים מחוץ לחריגות ה"תחרה" דומים למרבצי מישור ההצפה המודרניים. מרבצי מישור שיטפונות מיוצגים בעיקר על ידי סלעים חרסיתיים וקידוח לתוך אזורים אלה לא יהיה יעיל.

פרוסת RGB של קוביית הנתונים. במרכז (קצת משמאל למרכז) אפשר להתחקות אחר הנהר המתפתל.

פרוסת RGB של קוביית הנתונים. בצד שמאל ניתן לעקוב אחר הנהר המתפתל.


במקרים מסוימים, איכות הנתונים הסיסמיים מאפשרת תמונות ברורות יותר באופן משמעותי. זה תלוי במתודולוגיית העבודה בשטח, בציוד המשמש את אלגוריתם הפחתת הרעש. במקרים כאלה, לא רק שברי מערכות נהרות נראים, אלא גם נהרות פליאו מורחבים שלמים.

ערבוב RGB של שלושה מרכיבים של קוביית נתונים סיסמיים (פרוסה אופקית). עומק כ-2 ק"מ.

תמונת לוויין של נהר הוולגה ליד סרטוב

מסקנה

Wolfram Mathematica מאפשר לך לנתח נתונים סיסמיים ולפתור בעיות יישומיות הקשורות לחקר מינרלים, וחבילת GeologyIO הופכת את התהליך הזה לנוח יותר. המבנה של נתונים סיסמיים הוא כזה ששימוש בשיטות מובנות כדי להאיץ את החישובים (טבלה מקבילה[], ParallelDo[],...) יעיל מאוד ומאפשר לך לעבד כמויות גדולות של נתונים. במידה רבה, זה מקל על ידי תכונות אחסון הנתונים של חבילת GeologyIO. אגב, החבילה יכולה לשמש לא רק בתחום של חקר סייסמי יישומי. כמעט אותם סוגי נתונים משמשים ברדאר חודר קרקע ובסיסמולוגיה. אם יש לך הצעות כיצד לשפר את התוצאה, אילו אלגוריתמים לניתוח אותות מארסנל וולפרם Mathematica ישימים לנתונים כאלה, או אם יש לך הערות קריטיות כלשהן, בבקשה השאר תגובה.

מקור: www.habr.com