面波频散曲线反演是一种地震波形反演方法，用于估计地下结构的物理参数。其原理基于面波频散现象，即地震波在地下传播时会由于地下结构的变化而导致波速的变化，从而在地震记录中形成不同频率的相位延迟。具体而言，面波频散曲线反演过程如下：
1. 收集面波频散数据：利用地震台站网络采集地震记录，在不同频率范围内获取到面波的到时信息。比如地震记录如图：
2. 计算频散曲线：根据收集到的地震记录，计算每个频率处的相位延迟。通常使用频率－时间分析方法，如多积分法或小波变换等。根据下述代码可以绘制频散曲线如图：
在这幅图中黑色的点连成的线就是频率和相速度的关系即频散曲线。
3. 构建初始模型：根据已知的地质信息构建初始模型，即地下结构的物理参数初始分布。初始模型代码设置如下：
% Initial values for model parameters
n = 10;
h_initial = [1 1 2 2 4 5 4 6 5 7]; % m (array(n))
beta_initial = [75 90 150 180 240 290 290 300 320 330 340];%array(n+1) % m/s
n_unsat = n+1;
nu_unsat = 0.3;
alpha_temp = sqrt((2*(1-nu_unsat))/(1-2*nu_unsat))*beta_initial; % m/s
alpha_initial = [alpha_temp(1) 1440 1440 1440 1440 1440 1440 ...1440 1440 1440 1440]; % m/s array(n+1)
rho = [1850 1850 1850 1850 1900 1900 1900 1900 1950 1950 1950]; % kg/m^34. 正演模拟：以初始模型为基础，使用地震波传播的数值模拟方法，如有限差分法或有限元法等，计算面波的频散曲线。这里使用的是快速delta法，可以快速计算初始模型的频散曲线，这里使用函数MASWaves_theoretical_dispersion_curve_FDMA实现。代码如下：
%%
% [c_t,lambda_t] = MASWaves_theoretical_dispersion_curve_FDMA...
% (c_min,c_max,c_step,lambda,n,alpha,beta,rho,h,delta_c)
%%
% The function MASWaves_theoretical_dispersion_curve computes the
% theoretical fundamental mode dispersion curve for the stratified
% soil model defined by n, alpha, beta, rho and h at wavelengths lambda.
%
%% Input
% c_min Minimum testing Rayleigh wave phase velocity [m/s]
% c_max Maximum testing Rayleigh wave phase velocity [m/s]
% c_step Testing Rayleigh wave phase velocity increment [m/s]
% lambda Wavelength vector [m]
% n Number of finite thickness layers
% alpha Compressional wave velocity vector [m/s] (array of length n+1)
% beta Shear wave velocity vector [m/s] (array of length n+1)
% rho Mass density vector [kg/m^3] (array of length n+1)
% h Layer thickness vector [m] (array of length n)
% delta_c Zero search initiation parameter [m/s]
% At wave number k_i the zero search is initiated at
% a phase velocity of max{c_(i-1)-delta_c , c_min}, where
% c_(i-1) is the theoretical Rayleigh wave phase velocity
% value at wave number k_(i-1)
%
%% Output
% c_t Rayleigh wave phase velocity vector (theoretical
% dispersion curve) [m/s]
% lambda_t Rayleigh wave wavelength vector (theoretical dispersion
% curve) [m]
%
%% Subfunctions
% MASWaves_FDMA
%%%
function [c_t,lambda_t] = MASWaves_theoretical_dispersion_curve_FDMA...(c_min,c_max,c_step,lambda,n,alpha,beta,rho,h,delta_c)%
% c_min= c_min
% c_max=c_max
% c_step=c_step
% lambda=lambda_OBS
% n=n
% alpha=alpha_initial
% beta=beta_initial
% rho=rho
% h=h_initial
% delta_c=delta_c% Determine testing phase velocity values
c_test = c_min:c_step:c_max;% Wave numbers that correspond to wavelengths lambda
k = (2*pi)./lambda;% Number of modes (modes = 1, fundamental mode)
modes = 1;% Initialization
D = zeros(length(c_test),length(k));
c_t = zeros(length(k),modes);
lambda_t = zeros(length(k),modes);% For each wave number k, compute the dispersion function using
% increasing values of c_test until its value has a sign change.
m_loc = 1;
delta_m = round(delta_c/c_step);
for j = 1:length(k)for m = m_loc:length(c_test)D(j,m) = MASWaves_FDMA(c_test(m),k(j),n,alpha,beta,rho,h);if m==m_locsign_old = sign(MASWaves_FDMA(c_test(1),k(j),n,alpha,beta,rho,h));elsesign_old = signD;endsignD = sign(D(j,m));if sign_old*signD == -1c_t(j) = c_test(m);lambda_t(j) = 2*pi/k(j);m_loc = m - delta_m;if m_loc <= 0m_loc = 1;endbreakendend
end5. 与观测数据匹配：将模拟的频散曲线与实际观测得到的频散曲线进行比较作差，并通过某种优化算法蒙特卡洛法调整模型参数，使得模拟的频散曲线与观测数据相匹配。
beta_test = beta_opt + unifrnd(-(b_S/100).*beta_opt,(b_S/100).*beta_opt);
h_test = h_opt + unifrnd(-(b_h/100).*h_opt,(b_h/100).*h_opt);
6. 更新模型参数：根据优化算法得到的最佳参数值，更新初始模型的物理参数分布。
7. 重复以上步骤：通过迭代的方式，一直进行正演模拟和模型参数更新，直到达到收敛条件。
这样就可以得到横波速度随深度的变化的剖面。
所有的代码包放在这个链接里，或者看我的分享资源。