DAS Earthquake Analysis

P. McGill
13Mar25

Import a function to read  Sintela DAS V3 format files, which was  developed in the notebook SintelaReadRawV3.nb according to the documentation at https://sintela.atlassian.net/wiki/spaces/DEMO/pages/336917/Phase+Data+Sensor+Data+Sintela+Binary+Format.

sintelaReadRawV3[fileName_String]:= Module[ (* v1.0.0 of this function *)
{fileSizeBytes,inStream,syncWord,hdrSize,hdrVers,sampCnt,sampTime,numChans,numSamps,sampRate,chanSpacing,gaugeLen,trigger,startChan,chanStep,gpsPPSTime,gpsPPSSampOffs,gpsStat,trigFlags,statFlags,packetSizeBytes,totalPackets,totalSamps,packet,header,rdBuf,data={}},

(* Decode header information *)
fileSizeBytes=UnitConvert[FileSize[filename],"Bytes"]//QuantityMagnitude;
inStream=OpenRead[fileName,BinaryFormat->True];
syncWord=BinaryRead[inStream, "UnsignedInteger32"] (* offset 0 *);
hdrSize=BinaryRead[inStream, "UnsignedInteger32"] (* offset 4 *);
hdrVers=BinaryRead[inStream, "UnsignedInteger32"] (* offset 8 *);
If[
syncWord==16^^11223344&&hdrSize==96,
Print["Sintela Data Type 3 - Raw 32 Bit."], 
(Print["Not a Sintela Type 3 File."];Close[inStream];Abort[])
] ;(* end if *)

sampCnt=BinaryRead[inStream, "UnsignedInteger64"] (* offset 12 *);
sampTime=BinaryRead[inStream,"UnsignedInteger64"] (* offset 20 *);
numChans=BinaryRead[inStream,"UnsignedInteger32"] (* offset 28 *);
numSamps=BinaryRead[inStream,"UnsignedInteger32"] (* offset 32 *);
sampRate=BinaryRead[inStream,"Real32"] (* offset 36 *);
chanSpacing=BinaryRead[inStream,"Real32"] (* offset 40 *);
gaugeLen=BinaryRead[inStream,"Real32"] (* offset 44 *);
trigger=BinaryRead[inStream,"UnsignedInteger32"] (* offset 48 *);
startChan=BinaryRead[inStream,"UnsignedInteger32"] (* offset 52 *);
chanStep=BinaryRead[inStream,"UnsignedInteger32"] (* offset 56 *);
gpsPPSTime=BinaryRead[inStream,"UnsignedInteger32"] (* offset 60 *);
gpsPPSSampOffs=BinaryRead[inStream,"Integer32"] (* offset 64 *);
gpsStat=BinaryRead[inStream,"UnsignedInteger32"] (* offset 68 *);
trigFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 72 *);
statFlags=BinaryRead[inStream,"UnsignedInteger8"] (* offset 73 *);

(* Display first header *)
Print["Packet Channels x Samples : ",numChans," x ",numSamps];Print["Start Channel: ",startChan];
Print["Fs: ",NumberForm[sampRate,{8,3}]," Hz"] ;Print["Channel Spacing: ",NumberForm[chanSpacing,{4,1}]," m"];
Print["Gauge Length: ",NumberForm[gaugeLen,{4,1}]," m"];

Print["Start Time: ",DateString[FromUnixTime[sampTime/10^9,TimeZone->"GMT"]],".",Mod[sampTime,10^9]," GMT"];

packetSizeBytes=(hdrSize+(numChans*numSamps*4));
totalPackets=Round[fileSizeBytes/packetSizeBytes];
totalSamps=totalPackets*numSamps;
Print["Total Samples Per Channel: ",totalSamps];
Print["Total Packets: ",totalPackets];
Print["Bytes per Packet: ",packetSizeBytes];

(* Read and extract data samples (skipping over header blocks) *)
SetStreamPosition[inStream,0];(* Rewind file *)
Monitor[ (* show a progress indicator as we read the packets *)
For[packet=1,packet<=totalPackets,packet++,
(* skip over header and read samples *)
SetStreamPosition[inStream,StreamPosition[inStream]+hdrSize];
rdBuf=BinaryReadList[inStream,"Real32",numChans*numSamps];
AppendTo[data,rdBuf];
] (* end For *),
ProgressIndicator[packet,{1,totalPackets}]
] (* end Monitor *);

(* data are now in sublists of packetSizeBytes length, and we must rearrange to sublists of numChans length *);
data=Partition[Flatten[data],numChans];

Close[inStream];
Return[data]
] (* end module *)

After reviewing the list of events that Aaron M identified, which can be found at

 ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Aarons\ Event\ List

I choose what looks like an earthquake event on 1Feb25 UTC. Display the screen capture of the event.

(*Import["/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Aarons\ Event\ List/Screenshots/31\ Jan\ -\ 13.png","PNG"]*)

Import["/Volumes/Geo-Sense/Geo-Sense\ observations_AM/Screenshots/31\ Jan\ -\ 13.png","PNG"]

The event appears to have occurred just after 2120 PST on 31Jan25, which is 0520 UTC on 1Feb25.

Import the DAS data file.

(*filename="/Volumes/ProjectLibrary/902401_GeoSense/Project.Data/Deployments/DAS_2025_01_28/Geo-Sense\ data/decimator_2025-02-01_05.15.23_UTC_000481.raw";*)

filename="/Volumes/Geo-Sense/2025\ Jan\ Data/decimator_2025-02-01_05.15.23_UTC_000481.raw";

dasDat=sintelaReadRawV3[filename];

Dimensions@dasDat

MinMax@dasDat

Make a data browser that shows the first 10 seconds of data for any channel.

Manipulate[ListPlot[dasDat[[1;;10000,chan]],PlotRange->{-5,5},Frame->True,ImageSize->Large],{chan,1,162,1}]

Make a data browser that shows a 10 second period containing the event. We know the data file starts at 051524 UTC, and we want to start the viewer at 052003 UTC.

delay=QuantityMagnitude[UnitConvert[DateDifference[{2025,2,1,5,15,24},{2025,2,1,5,20,3}],"Seconds"]]//Round

We need to skip 279 seconds into the file.

Manipulate[ListPlot[dasDat[[delay 1000;;(delay+10) 1000,chan]],PlotRange->{-4,6},Frame->True,GridLines->Automatic,ImageSize->Large],{chan,1,162,1}]

There is an event that starts about 4.5 seconds into the plot. Find the time.

eventDateWarp=DatePlus[DateObject[{2025,2,1,5,15,23.54},TimeZone->"UTC"],Quantity[delay+4.5,"Seconds"]]

Channel 18 appears to have a strong signal. Try to improve the SNR by applying a 2 Hz highpass filter.

ListPlot[HighpassFilter[dasDat[[(delay+2) 1000;;(delay+12) 1000,18]],Quantity[2,"Hertz"],SampleRate->1000],PlotRange->{-4,6},Frame->True,FrameLabel->{{"Phase (radians)",""},{"Time (ms)","Geo-Sense DAS channel 18\r2 Hz Highpass Filter"}},GridLines->Automatic,ImageSize->Large]

There appears to be an arrival of more energy about 5 s after the onset, typical of the S wave following the P wave of an earthquake.

Using information at

http://eqseis.geosc.psu.edu/cammon/HTML/Classes/IntroQuakes/Notes/waves_and_interior.html

we can find the distance to the earthquake.

dt=Quantity[5,"Seconds"]; (* time difference of P- and S-wave arrivals *)
Vp=Quantity[6,"Kilometers"/"Seconds"]; (* P-wave velocity in km/s *)
Vs=Quantity[4,"Kilometers"/"Seconds"]; (* S-wave velocity in km/s *)
dist=(dt/(((1/Vs)-(1/Vp))))

Apply the filter to all the channels.

Manipulate[ListPlot[HighpassFilter[dasDat[[(delay+2) 1000;;(delay+12) 1000,chan]],(Quantity[2, "Hz", hertz, "Hertz"]),SampleRate->1000],PlotRange->{-4,6},Frame->True,GridLines->Automatic,ImageSize->Large],{chan,1,162,1}]

Plot a periodogram for the time including 5 s before the event until 5 s after.

Periodogram[dasDat[[delay 1000;;(delay+10) 1000,18]],SampleRate->1000,PlotRange->{-80,20},Frame->True,FrameLabel->{"Frequency (Hz)","Power (dB)"},ImageSize->Large]

Overlay two periodograms for a 10 s span before the event and a 10 s span during the event, for the frequency range 0–10 Hz.

Periodogram[{dasDat[[(delay-10) 1000;;delay 1000,18]],dasDat[[(delay+4) 1000;;(delay+14) 1000,18]]},SampleRate->1000,Frame->True,FrameLabel->{"Frequency (Hz)","Power (dB)"},ImageSize->Large,PlotRange->{{0,10},{-40,30}},PlotLegends->{"Before","During"}]

ArrayPlot[Transpose[BlockMap[RootMeanSquare[#-Mean[#]]&,#,100]&/@Transpose[dasDat[[(delay+2) 1000;;(delay+22) 1000,All]]]],Frame->True,FrameLabel->{"Time (s)","DAS Channel"},FrameTicks->Automatic,PlotRange->{0,5},ColorFunction->"BlueGreenYellow",ImageSize->Full]

eventDateWarp

eventDateTrue=(DateObject["Fri 31 Jan 2025 17:44:34""UTC", DateObject[{2025,1,31,17,44,34.63891334854998`},"Instant","Gregorian","UTC"]])

timeDelta=Quantity[60,"Minutes"];

EarthquakeData[Here,1,{eventDateTrue-timeDelta,eventDateTrue+timeDelta}]

EarthquakeData[Here,1,{{2025,1,31,0},{2025,1,31,23,59,59}},TimeZone->"UTC"]

data=EarthquakeData[Here,1,{{2025,1,31,0},{2025,1,31,23,59,59}},"Position",TimeZone->"UTC"]["Values"];
GeoGraphics[{GeoStyling["ReliefMap"],Red,PointSize[Large],Point[data]}]