Fault-Inversion-Auto/Fun_SetInvStations.m at main · Yuxiao-git2025/Fault-Inversion-Auto · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
function [ob, g, loca, mm] = Fun_SetInvStations( ...
    ob, g, loca, mm, epicenter, srate, grid, gridsize, source, muSubfault, ...
    Range1,Range2,fBands, mag, opts)
% =========================================================================
% Fun_PrepareInversionStations
%
% Stage 1:
%   Preliminary station screening based on a low-iteration inversion.
%
% Stage 2:
%   Remove redundant stations for very dense station networks.
% =========================================================================
%% ------------------------------------------------------------------------
% Stage 1: Preliminary station screening
% -------------------------------------------------------------------------
Nsta = size(loca, 1);
if Nsta >= opts.MinStationScreen
    fprintf('[Do] Preliminary Station Screening\n');
    [obCandidate, gCandidate, locaCandidate, mmCandidate, ~] = divi_inver( ...
        epicenter, loca, ob, g, mm, srate, grid, gridsize, ...
        source, muSubfault, fBands, opts.PreScreenIteration, mag);
    % Retain unique stations and preserve their original candidate order.
    nCandidateStation = size(locaCandidate, 1);
    [locaUnique, keepStation] = unique(locaCandidate, 'rows', 'stable');

    mmUnique = mmCandidate(keepStation, :);
    % Waveform/Green channel layout:
    % [E1...EN, N1...NN, U1...UN]
    keepChannel = [ ...
        keepStation; ...
        nCandidateStation + keepStation; ...
        2 * nCandidateStation + keepStation];
    keepChannel = keepChannel(:);
    ob = obCandidate(:, keepChannel);
    g  = gCandidate(:, keepChannel, :);

    loca = locaUnique;
    mm   = mmUnique;
    fprintf('[Check] Retained %d stations after preliminary screening\n', ...
        size(loca, 1));
end

%% ------------------------------------------------------------------------
% Stage 2: Remove azimuthally redundant stations in dense networks
% -------------------------------------------------------------------------
Nsta = size(loca, 1);
if Nsta >= opts.DenseStationLimit
    fprintf('[Check] Fault-Plane Pre-Processing\n');
    % Use the first station only as a reference for distance screening.
    da = da_zh(loca, loca(1, :), 1);
    da1 = da(:, 1);
    % First try a relatively strict azimuth-distance threshold.
    %Range1 = 12;
    %Range2 = 5;
    Ndel = azim_del(da1, Range1);
    % If too many stations would be removed, relax the threshold.
    if Nsta - numel(Ndel) < opts.MinRemainStation
        Ndel = azim_del(da1, Range2);
    end
    % Convert station indices to E/N/U waveform-channel indices.
    Cdel = [ ...
        Ndel; ...
        Nsta + Ndel; ...
        2 * Nsta + Ndel];

    Cdel = Cdel(:);
    loca(Ndel, :) = [];
    mm(Ndel, :)   = [];

    ob(:, Cdel)  = [];
    g(:, Cdel, :) = [];
    fprintf('   Removed %d redundant stations\n', numel(Ndel));
end

%% ------------------------------------------------------------------------
% Dimension check
% -------------------------------------------------------------------------
Nsta = size(loca, 1);
nChannel = size(ob, 2);
assert(nChannel == 3 * Nsta, ...
    'The number of observation channels must equal 3*Nsta.');
assert(size(g, 2) == nChannel, ...
    'The Green-function channel number must match observation channels.');
assert(size(g, 1) == size(ob, 1), ...
    'Green functions and observations must have the same time length.');
fprintf('   Remained %d stations\n\n', Nsta);

end