% Processing of Low-Frequency DAS Displacement in a Bedrock Borehole
% Live Script in support of Article Submitted to Water Resources Research:
% Distributed Acoustic Sensing (DAS) as a Distributed Hydraulic Sensor in Fractured Bedrock
% M. W. Becker(1), T. I. Coleman(2), and C. C. Ciervo(1)
% 1 California State University, Long Beach, Geology Department, 1250 Bellflower Boulevard, Long Beach, California, 90840, USA.
% 2 Silixa LLC, 3102 W Broadway St, Suite A, Missoula MT 59808, USA.
% Corresponding author: Matthew W. Becker (matt.becker@csulb.edu).
%
% As described in the article, low-frequency DAS strain rate was measured in fiber optic cables coupled to boreholes completed in crystalline bedrock.  This example shows the results of displacement measured along borehole FSE9.
% The instrument used for this study was a Silixa iDAS.  The iDAS measures dynamic strain, i.e. strain rate through iphase demodulation of a laser signal to along the cable.  The settings used for the study was a gauge length of 10 m, a spatial sampling of 0.25 m, and a temporal sampling of 1 kHz.  Because these datasets are very large, the data were decimated to produce a sampling rate of 100 Hz.  Decimation was accomplished using the Matlab "decimate" command which applies a lowpass Chebyshev Type I infinite impulse response (IIR) filter of order 8.
% For this example, load the file FSE-9b_240SecP_Decimated_StepTest.mat which carries the decimated displacement rate (nm/sec). This is the data collected in response to a hydraulic pulse at FSE 6 with a period of 240 seconds (4 min).  A proprietary code from Silixa converts their raw data to displacement rate in nm/sample.   The format of the file a matrix with rows representing distance along the fiber (0.25 m sampling) and columns representing time (0.01 sec sampling) where each value in the matrix is displacement rate (nm/sample).

load('FSE-9b_240SecP_Decimated_StepTest.mat')

% The parameters of the test are specified here.

P=240;                           %<==  ENTER PERIOD OF OSCILLATION (s)
Fs = 100;                        %<==  ENTER SAMPLING FREQUENCY (samp/s=Hz)
Dz = 0.25;                        %<==  ENTER SPATIAL SAMPLING (m)
C1=50;                           %<==  ENTER FIRST CHANNEL TO EVALUATE (Each channel = 0.25m length)
C2=277;                          %<==  ENTER LAST CHANNEL TO EVALUATE(Each channel = 0.25m length)
TestName='FSE9b ';               %<==  ENTER NAME OF TEST

%Calculate initial parameters
Freq=1/P;                                   % Find frequency from period (Hz)
T = 1/Fs;                                   % Sample time (s)
Chan=C1:C2-1;                               % Create array of channels to process
C1_depth=0.2;                               % Set the reference depth for the first channel (m)
D=C1_depth+C1*Dz+Chan*Dz;                   % Convert channels to depths (m)
Dfrac=43;                                   % Expected depth of strained fracture (m)
Cfrac=round((Dfrac-C1_depth-C1*Dz)/Dz);     % Channel where the peak displacement is expected

% Process strings to report test parameters
PerStr=num2str(P/60);
formatSpec = '%s%s%s%s';
A2 = PerStr; A3=' min period';
name = sprintf(formatSpec,TestName,A2,A3);

% The native displacement rate (rawDAS) is scaled and converted to displacement through integration.

rawDAS=decmatrix(:,C1:C2-1);     % Subset raw data to channels of interest
L = length(rawDAS);             % Length of signal in number of samples
t = (0:L-1)*T;                  % Calculate time increments from samples (s)

%Integrate Raw data and convert to Displacement  
engDAS=116*rawDAS;              %convert to nm/sample (eng units)
disprate=engDAS*Fs;             %convert nm/sample*sample/s=nm/s
disp=cumtrapz(t,engDAS);        %integrate to displacement (nm)
disp=detrend(disp,4);        %detrend for envelope

% Plot the displacement matrix.  Displacement is in nm.  Notice the periodicity of the displacement along the borehole corresponding to the 240 min period of the test.

imagesc(disp'); c=colorbar; xlabel('Samples'); ylabel('Channels'); grid on;
c.Label.String='Displacement (nm)'; c.Label.FontSize = 12;

% The amplitude is extracted from the displacement signal by using the difference between a moving maximum and minimum with a window equal to the period of the hydraulic stress, in this case, 240 sec or 24000 samples.

dispFrac=disp(:,Cfrac);         % Isolate displacement at the fracture
NMov=P*Fs*1.25;                      % The window of the moving min/max
MMax=movmax(dispFrac,NMov);     % Calculate moving maximum envelope
MMin=movmin(dispFrac,NMov);     % Calculate moving minimum envelope
MovAmp=(MMax-MMin)/2;           % Calculate the displacement amplitude through time
%calculate the mean amplitude at the fracture ignoring early and late time data:
MeanMovAmp=mean(MovAmp(floor(L*.25):floor(L*.75)));
% Plot the displacement at the channel near fracture (Cpeak)
figure(1);  
        plot(t/60,dispFrac,'LineWidth',1.5); hold on;
        plot(t/60,MMax(:),t/60,MMin(:)); hold off;
         %axis([0,inf,0,inf]) ; grid on;
         xlabel('Minutes');
         ylabel('FFT Displacement (nm/m)');
         legend('Displacement','Upper Envelope','Lower Envelope')
         title(['Displacement: ' name]);

% We can also extract the amplitude by using a fast Fourier transform (FFT) algorthm and finding the amplitude at the known applied frequency.

NFFT = 2^nextpow2(L);               % Next power of 2 from length of y
f = Fs/2*linspace(0,1,NFFT/2+1);    % Make an array of frequencies for the spectrum
CF=find(f>Freq,1);                  % Find index at imposed frequency
YAmpD = fft(dispFrac,NFFT)/L;       % Calculate FFT
AmpD=2*abs(YAmpD(1:NFFT/2+1,:));    % Find amplitude from FFT
FFTAmpPeak=max(AmpD(CF-2:CF+2));    %Find peak amplitude from the FFT in vicinity of imposed frequency (CF)

% Plot the FFT spectrum

figure(2);  
    plot(f,AmpD,'LineWidth',1.5); hold on;
    plot(f(CF),FFTAmpPeak,'*'); hold off;
    axis([0,2*Freq,0,inf]) ; grid on;
    xlabel('Frequency (Hz)');
    ylabel('FFT disp (nm)');
    legend('FFT Amplitude','Amplitude Peak Value')
    title(['FFT Spectrum Displacement: ' name]);

% We note that the two methods do not produce the same amplitude.  The FFT may experience spectral leaking.

fprintf(' Mean FFT Amp from Moving Avg  = %.3f nm\n Mean FFT Amp from Integration = %.3f nm\n', MeanMovAmp, FFTAmpPeak);