根据克里斯托弗方程（Christoffel equation）计算相速度和偏振(极化)方向，进而计算群速度

一. 基础理论

对于各向异性介质，从基础理论上认识几大方程以及它们之间的关系有助于加深对各向异性性质的了解，首先要搞清楚以下几点：

1.弹性波动方程（Elastic wave equation）与克里斯托弗方程（Christoffel equation）之间的关系为何？

2.克里斯托弗方程（Christoffel equation）与程函方程（Eikonal equation）之间的关系为何？

3.几大方程的物理意义是什么？

要回答这些问题，还是从推导入手。弹性波方程是通过应力（stress）应变（strain）关系（广义胡克定律，或本构关系）+运动微分方程（动量守恒线性方程 Linearized equation of momentum conservation Navier 方程或牛顿第二定律）导出，描述了波的所有动力学和运动学传播特性。它们描述了弹性介质内部质点的位移、 应力 、 应变之间相互联系的普遍规律，是建立弹性波动方程的基础 。

而Christoffel 方程是各向异性弹性波动方程的平面波（单相波，对应相速度）解形式，可以用来独立描述三种类型的波的极化和相速度特征，即运动学（速度）和偏振（振幅以及方向）。程函方程也是将平面波解带入波动方程，在高频近似下导出，只描述运动学特征（走时）。

克里斯托弗方程（Christoffel equation）形式上与特征矩阵相似，相速度即为对应矩阵的特征值，极化方向即为对应矩阵的特征向量。因此，可以通过克里斯托弗方程（Christoffel equation）计算相速度和偏振方向。有了相速度就可以通过Berryman(1979)推导的相速度群速度，相角群角之间的关系求对应群角的群速度。

二. 代码：

%% Christoffel equation, Phase and Group Velocity %% 
%%     Written By Zhang Jianming，2021.10.10      %% 

clear ;clc
pi=4*atan(1);

rou=3000;
% c11=10.2e9;
% c66=2.3e9;
% c44=1.36e9;
% c33=7.0e9;
% c13=4.96e9;

c11=66.6e9;
c66=23.5e9;
c44=10.9e9;
c33=39.9e9;
c13=39.4e9;
c12=c11-2*c66;
fi=pi/4;
N=360;
pahse_angle=zeros(1,N);
vs1=zeros(1,N);
vs2=zeros(1,N);
vp=zeros(1,N);

g_vs1=zeros(1,N);
g_vs2=zeros(1,N);
g_vp=zeros(1,N);

g_angle_vs1=zeros(1,N);
g_angle_vs2=zeros(1,N);
g_angle_vp=zeros(1,N);

for i=1:1:N
    sita=pi/180*i;
    pahse_angle(i)=sita;
    %Christoffel矩阵
    tao11=c11*(sin(sita)^2)*(cos(fi)^2)+c66*(sin(sita)^2)*(sin(fi)^2)+c44*(cos(sita)^2);
    tao22=c66*(sin(sita)^2)*(cos(fi)^2)+c11*(sin(sita)^2)*(sin(fi)^2)+c44*(cos(sita)^2);
    tao33=c44*(sin(sita)^2)*(cos(fi)^2)+c44*(sin(sita)^2)*(sin(fi)^2)+c33*(cos(sita)^2);
    tao12=(c11-c66)*(sin(sita)^2)*sin(fi)*cos(fi);
    tao13=(c13+c44)*sin(sita)*cos(sita)*cos(fi);
    tao23=(c13+c44)*sin(sita)*cos(sita)*sin(fi);
    a=[tao11 tao12 tao13;tao12 tao22 tao23;tao13 tao23 tao33];

    %求特征值
    %[V,D] = eig(A) 返回特征值的对角矩阵 D 和矩阵 V，其列是对应的右特征向量，使得 A*V = V*D。
    [x,y]=eig(a);
    %特征值
    vs1(1,i)=sqrt(y(1,1)/rou);
    vs2(1,i)=sqrt(y(2,2)/rou);
    vp(1,i)=sqrt(y(3,3)/rou);
    %特征值对应的矩阵
    p_vs1(1,i)=x(1,1);p_vs1(2,i)=x(2,1);p_vs1(3,i)=x(3,1);
    p_vs2(1,i)=x(1,2);p_vs2(2,i)=x(2,2);p_vs2(3,i)=x(3,2);
    p_vp(1,i)=x(1,3);p_vp(2,i)=x(2,3);p_vp(3,i)=x(3,3);
end


for i=1:1:N-1
    g_vs1(1,i)=sqrt(vs1(1,i)^2+(vs1(1,i+1)-vs1(1,i))^2/(pi/180)^2);
    g_vs2(1,i)=sqrt(vs2(1,i)^2+(vs2(1,i+1)-vs2(1,i))^2/(pi/180)^2);
    g_vp(1,i)=sqrt(vp(1,i)^2+(vp(1,i+1)-vp(1,i))^2/(pi/180)^2);

    g_angle_vs1(1,i)=atan(1/vs1(1,i)*(vs1(1,i+1)-vs1(1,i))/(pi/180))+pi/180*i;
    g_angle_vs2(1,i)=atan(1/vs2(1,i)*(vs2(1,i+1)-vs2(1,i))/(pi/180))+pi/180*i;
    g_angle_vp(1,i)=atan(1/vp(1,i)*(vp(1,i+1)-vp(1,i))/(pi/180))+pi/180*i;
end
g_vs1(1,N)=vs1(1,N);
g_vs2(1,N)=vs2(1,N);
g_vp(1,N)=vp(1,N);
g_angle_vs1(1,N)=pi/180*N;
g_angle_vs2(1,N)=pi/180*N;
g_angle_vp(1,N)=pi/180*N;

figure();
%画相角相速度
polar(pahse_angle,vp,'--b');
hold on;
polar(pahse_angle,vs1,'--r');
hold on;
polar(pahse_angle,vs2,'--g');

%画群角群速度
polar(g_angle_vp,g_vp,'-b');
hold on;
polar(g_angle_vs1,g_vs1,'-r');
hold on;
polar(g_angle_vs2,g_vs2,'-g');
hold on;

%%特征向量计算
for i=1:5:N
    x_vs1(i)=vs1(1,i)*cos(pi/180*i);
    y_vs1(i)=vs1(1,i)*sin(pi/180*i);
    u_vs1(i)=sin(acos(p_vs1(3,i)))*cos(acos(p_vs1(1,i)/norm([p_vs1(1,i) p_vs1(2,i)])));
    v_vs1(i)=sin(acos(p_vs1(3,i)))*sin(acos(p_vs1(1,i)/norm([p_vs1(1,i) p_vs1(2,i)])));

    x_vs2(i)=vs2(1,i)*cos(pi/180*i);
    y_vs2(i)=vs2(1,i)*sin(pi/180*i);
    u_vs2(i)=norm([p_vs2(1,i) p_vs2(2,i) p_vs2(3,i)])*sin(fi)*cos(pi/180*i);
    v_vs2(i)=norm([p_vs2(1,i) p_vs2(2,i) p_vs2(3,i)])*sin(fi)*sin(pi/180*i);

    x_vp(i)=vp(1,i)*cos(pi/180*i);
    y_vp(i)=vp(1,i)*sin(pi/180*i);
    u_vp(i)=norm([p_vp(1,i) p_vp(2,i) p_vp(3,i)])*sin(fi)*cos(pi/180*i);
    v_vp(i)=norm([p_vp(1,i) p_vp(2,i) p_vp(3,i)])*sin(fi)*sin(pi/180*i);
end
%画偏振方向(特征向量)
%quiver(x,y,u,v,scale)在(x,y)位置处画出(u,v)的向量箭头,scale为矢量长度
quiver(x_vs1,y_vs1,u_vs1,v_vs1,'r')
quiver(x_vs2,y_vs2,u_vs2,v_vs2,'g')
quiver(x_vp,y_vp,u_vp,v_vp,'b')

legend('Phase-Vp','Phase-Vs1','Phase-Vs2','Group-Vp','Group-Vs1','Group-Vs2',Location='southeastoutside')

三. 计算实例

下面给出两种典型VTI介质的计算结果：

刚度系数1：； 刚度系数2：

四. 致谢

该博客内容是在同济大学地球物理系程玖兵教授的《各向异性与多分量地震学》课程上学习完成，部分内容截图自程老师课件。更多优质内容和参考学习资料欢迎讨论互相学习。