Fault-Inversion-Auto/Fun_SetInvFault.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
% =========================================================================
% Fun_IterativeFaultPlaneUpdate
%
% Perform repeated finite-fault inversion and fault-plane boundary updates.
%
% The function stops when:
%   1. The slip distribution has stable low-slip margins;
%   2. Maximum outer iteration number is reached;
%   3. The inversion gives non-positive slip everywhere.
% =========================================================================
function [ob, g, loca, ...
          locasub, locadep, muSubfault, ...
          grid, source, index, ...
          substf, substfa, syn, obr, res, ...
          slip] = Fun_SetInvFault( ...
          ob, g, loca, ...
          locasub, locadep, muSubfault, ...
          PATHg, dataVelFolder, ...
          SDR, grid, gridsize, source, ...
          epicenter, srate, IndexEQ, index, depth, ...
          fBands, mag, opts)
substf  = [];
substfa = [];
syn     = [];
obr     = [];
res     = [];
slip    = [];
for iter = 1:opts.MaxOuterIteration
    fprintf('[Do] Fault-Plane Processing [%d/%d]\n', iter, opts.MaxOuterIteration);
    % --------------------------------------------------------------------
    % Finite-fault inversion
    % ---------------------------------------------------------------------
    nChannel = size(ob, 2);
    [~, substf, substfa, syn, obr, res, ~, loca] = ids_data( ...
        ob, g, loca, srate, grid, gridsize, source, ...
        nChannel, muSubfault, fBands, opts.InnerIteration, mag);

    %% --------------------------------------------------------------------
    % Convert source-time functions to final subfault slip
    % ---------------------------------------------------------------------
    slipVector = sum(substf, 1) ./ 3e16 ./ prod(gridsize);
    slip = reshape(slipVector, grid);
    Areas = prod(gridsize) * 1e6;
    SeisMomentNow = sum(slip(:) .* muSubfault(:) .* Areas);
    MwNow = m2m(SeisMomentNow);

    fprintf('[Result] Current Moment Magnitude: Mw=%.2f\n', MwNow);

    %% --------------------------------------------------------------------
    % Test whether the fault boundary is stable
    % ---------------------------------------------------------------------
    [isStable, NormSlip] = Fun_SetInvFaultBoundary(slip, locadep, gridsize, opts.SlipCutoff);
    if isStable
        fprintf('[Result] The fault boundary is stable---Stop updating\n\n');
        break;
    end

    %% --------------------------------------------------------------------
    % Stop if the outer-loop limit is reached
    % ---------------------------------------------------------------------
    if iter == opts.MaxOuterIteration
        fprintf('[Result] Maximum outer iteration is reached\n\n');
        break;
    end

    %% --------------------------------------------------------------------
    % Update fault geometry from significant-slip region
    % ---------------------------------------------------------------------
    [grid, source, index] = Fun_SetInvFaultUpdate( ...
        NormSlip, grid, source, gridsize, SDR, loca, epicenter, depth, IndexEQ);
    fprintf(['   New Grid: [%d, %d]; New Source: [%d, %d]; ', ...
             'New Index: %d\n\n'], grid(1), grid(2), source(1), source(2), index);

    %% --------------------------------------------------------------------
    % Rebuild Green functions for the new fault geometry
    % ---------------------------------------------------------------------
    nTime = size(ob, 1);
    [g, locasub, locadep, muSubfault] = Fun_RebuildGreenModel( ...
        PATHg, SDR, grid, gridsize, source, loca, epicenter, ...
        srate, IndexEQ, index, depth, dataVelFolder, nTime);


end

end