（１）利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量；
（２）利用DEM模型将孔隙加进系统，并计算干燥岩石骨架的体积模量和剪切模量；
（３）利用Hudson理论和Schoenberg＆Sayers理论，在碳酸盐岩介质各向同性背景中加入裂缝系统，修正计算出的模量；
（４）利用Wood公式将孔隙流体进行混合，计算出混合流体的体积模量；
（ 5 ）计算裂缝型碳酸盐岩的弹性参数和各向异性参数．
1.利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量
利用岩石基质组分，计算混合矿物的弹性模量
井的矿物成分含量和所需的井曲线如下：
经过Hashin-Shtrikman计算后获得的等效矿物模型结果如下：
//////////////////////////////////////////////////////////////////////////////////////////////
void HashinShtrikman(float f[], int m, int n, float K[], float G[], float y[])//对应第一步 等效岩石基质的建立对应公式1和公式2
{
float *K_HS1 = new float[m];
float *G_HS1 = new float[m];
float *K_HS2 = new float[m];
float *G_HS2 = new float[m];
float Kmax = max(K, n);
float Kmin = min(K, n);
float Gmax = max(G, n);
float Gmin = min(G, n);
if (n == 2)
{
float K1 = K[0];
float K2 = K[1];
float G1 = G[0];
float G2 = G[1];
for (int i = 0; i < m; i++)
{
float f1 = f[2 * i];
float f2 = f[2 * i + 1];//组分体积比归一化
K_HS1[i] = K1 + f2 / (1.0f / (K2 - K1) + f1 / (K1 + 4.0f / 3.0f * Gmin));
G_HS1[i] = G1 + f2 / (1.0f / (G2 - G1) + f1 / (G1 + Gmin / 6 * (9 * Kmin + 8 * Gmin) / (Kmin + 2 * Gmin))); //公式1
K_HS2[i] = K1 + f2 / (1.0f / (K2 - K1) + f1 / (K1 + 4.0f / 3.0f * Gmax));
G_HS2[i] = G1 + f2 / (1.0f / (G2 - G1) + f1 / (G1 + Gmax / 6 * (9 * Kmax + 8 * Gmax) / (Kmax + 2 * Gmax))); //公式2
y[4 * i + 0] = K_HS1[i];
y[4 * i + 1] = K_HS2[i];
y[4 * i + 2] = G_HS1[i];
y[4 * i + 3] = G_HS2[i];
}
}
else
{
//float Kmax = max(K, n);
//float Kmin = min(K, n);
//float Gmax = max(G, n);
//float Gmin = min(G, n);
float Emax = Gmax / 6.0f * (9 * Kmax + 8 * Gmax) / (Kmax + 2 * Gmax);
float Emin = Gmin / 6.0f * (9 * Kmin + 8 * Gmin) / (Kmin + 2 * Gmin);
for (int i = 0; i < m; i++)
{
K_HS1[i] = 0;
K_HS2[i] = 0;
G_HS1[i] = 0;
G_HS2[i] = 0;
for (int j = 0; j < n; j++)//
{
K_HS1[i] = K_HS1[i] + f[i * n + j] / (K[j] + 4.0f / 3.0f * Gmax);//
K_HS2[i] = K_HS2[i] + f[i * n + j] / (K[j] + 4.0f / 3.0f * Gmin);
G_HS1[i] = G_HS1[i] + f[i * n + j] / (G[j] + Emax);
G_HS2[i] = G_HS2[i] + f[i * n + j] / (G[j] + Emin);
}
K_HS1[i] = 1.0f / K_HS1[i] - 4.0f / 3.0f * Gmax;
K_HS2[i] = 1.0f / K_HS2[i] - 4.0f / 3.0f * Gmin;
G_HS1[i] = 1.0f / G_HS1[i] - Emax;
G_HS2[i] = 1.0f / G_HS2[i] - Emin;
y[4 * i + 0] = K_HS1[i];
y[4 * i + 1] = K_HS2[i];
y[4 * i + 2] = G_HS1[i];
y[4 * i + 3] = G_HS2[i];
}
}
delete[]K_HS1;
delete[]K_HS2;
delete[]G_HS1;
delete[]G_HS2;
}
２.利用DEM模型将孔隙加进系统，并计算干燥岩石骨架的体积模量和剪切模量；
孔隙空间分成主成分孔隙和次成分孔隙
利用广义xuwhite模型计算岩石的弹性模量：
其中：
采用微分等效介质模型加入孔隙。
\[(1-y)d[K^*(y)]/dy=(K_2-K^*)P^{*2}(y)
\]
\]
\[(1-y)d[u^*(y)]/dy=(u_2-u^*)Q^{*2}(y)
\]
\]
初始条件是
K(0)=K1和μ（0）=μ1；
其中：
K1，μ1为初始主相材料的体积模量和剪切模量
K2，μ2为逐渐加入包含物材料的体积模量和剪切模量
y=相2的含量
对于流体包含物和空包含物，y等于孔隙度φ，项P和Q是给定的几何因素，P和Q的上标*2指的是几何因素具有等效的背景介质中的包含物材料（具体背景材料参数和系数参考岩石物理手册P110）。
DEM计算后的结果：
紫色为基质正演结果，红色为干骨架正演结果，绿色为实测数据；
void XuWhite_DEM2(float Por, float Sg, float So, float Vsh, int MatrixMixMethod, int PorGeometry, int NAspect, float *y)//DEM微分等效模型公式3
{
float Km = 71.07;
float Gm = 28.3833;
float Pm = 2.71;
float Ksh = 42.05;
float Gsh = 19.08;
float Psh = 2.62;
float Kw = 2.25;
float Gw = 0;
float Pw = 1;
float Kg = 0.71415;
float Gg = 0;
float Pg = 0.54;
float Ko = 0.91415;
float Po = 0.77;
float ratio[2] = { 0.12,0.05 };
float Vsand, Vsand1, Vsh1, Por_s, Por_c;
//matrix
Vsand = 1 - Vsh - Por;
Vsh1 = Vsh / (1 - Por);
Vsand1 = Vsand / (1 - Por);
Por_s = Vsand1 * Por;//体积比 Por_s 对应Tiijj0[0],Por_c对应Tiijj[1]即：组分1为砂岩，组分2为泥岩
Por_c = Vsh1 * Por;
//float f[2] = { 1 - Vsh, Vsh };//这个参数传递给HashinShtrikman()
float f[2] = { Vsand1, Vsh1 };
float K[2] = { Km, Ksh };
float G[2] = { Gm, Gsh };
float P[2] = { Pm,Psh };
float *y1 = new float[4];
////*********Dry matrix model (using DEM model)
// Xu-white coefficient
//float A,B,R,Tiijj,FF,alpha,theta,f,F1,F2,F3,F4,F5,F6,F7,F8,F9;
int NP = 5000;//NP是积分的个数，确定积分间隔
// int N=length(Fraction);
float Fraction[2] = { Por_s, Por_c };//第一孔隙是主岩性，第二孔隙是副岩性：例如砂泥组合第一孔隙是砂岩，第二孔隙是泥岩
float Aspect[2] = { ratio[0],ratio[1] };//第一孔隙是主岩性，第二孔隙是副岩性；
float p, dp;
float PP, QQ, Em, bm;
float K_in, Mu_in, Ein;
float K_m, Mu_m, K_tmp, Mu_tmp, Kd, Gd;
int i, j;
K_in = 0;//填充空气
Mu_in = 0;
K_m = Km;
Mu_m = Gm;
p = 0;
//for(i=0;i<NAspect;i++)//几种孔隙
for (i = NAspect - 1; i >= 0; i--)//几种孔隙
{
dp = Fraction[i] / NP;
K_tmp = K_m;
Mu_tmp = Mu_m;
for (j = 1; j <= NP; j++)
{
float A = -1;
float B = 0;
//float R = 3 * Gm / (3 * Km + 4 * Gm);
//float *Tiijj= new float[2];
//float *FF = new float[2];
float R = 3 * Mu_tmp / (3 * K_tmp + 4 * Mu_tmp);//K_tmp,Mu_tmp是每次的增加孔隙中间结果
float Tiijj, FF;//不需要记录每个孔隙
//for (int i = 0; i < 2; i++)
//{
float alpha = Aspect[i];
float theta = alpha / powf(1 - alpha * alpha, 1.5) * (acos(alpha) - alpha * sqrt(1 - alpha * alpha));
float f = alpha * alpha / (1 - alpha * alpha) * (3 * theta - 2);
float F1 = 1 + A * (1.5 * (f + theta) - R * (1.5 * f + 2.5 * theta - 4.0 / 3.0));
float F2 = 1 + A * (1 + 1.5 * (f + theta) - R / 2 * (3 * f + 5 * theta)) + B * (3 - 4 * R) + A / 2 * (A + 3 * B) * (3 - 4 * R) * (f + theta - R * (f - theta + 2 * theta * theta));
float F3 = 1 + A * (1 - (f + 1.5 * theta) + R * (f + theta));
//float F3 = 1 +A/2*( R*(2-theta)+(1+alpha*alpha)/alpha*f*(R-1) );
float F4 = 1 + A / 4 * (f + 3 * theta - R * (f - theta));
float F5 = A * (-f + R * (f + theta - 4.0 / 3.0)) + B * theta * (3 - 4 * R);
float F6 = 1 + A * (1 + f - R * (f + theta)) + B * (1 - theta) * (3 - 4 * R);
float F7 = 2 + A / 4 * (3 * f + 9 * theta - R * (3 * f + 5 * theta)) + B * theta * (3 - 4 * R);
float F8 = A * (1 - 2 * R + f / 2 * (R - 1) + theta / 2 * (5 * R - 3)) + B * (1 - theta) * (3 - 4 * R);
float F9 = A * ((R - 1) * f - R * theta) + B * theta * (3 - 4 * R);
Tiijj = 3 * F1 / F2;
FF = 2 / F3 + 1 / F4 + (F4 * F5 + F6 * F7 - F8 * F9) / (F2 * F4);
//printf("======%d=====%d=====%f=====%f===\n",i,j,alpha,theta,f);
//}
PP = 1.0 / 3.0*Tiijj;
QQ = 1.0 / 5.0*FF;
p = p + dp;
//p=dp;
K_tmp = (K_in - K_tmp)*PP*dp / (1 - p) + K_tmp;
Mu_tmp = (Mu_in - Mu_tmp)*QQ*dp / (1 - p) + Mu_tmp;
}
K_m = K_tmp;
Mu_m = Mu_tmp;
}
Kd = K_m;
Gd = Mu_m;
y[0] = Kd;
y[1] = Gd;
y[2] = Pm;
}
3.利用Ｈｕｄｓｏｎ理论和Ｓｃｈｏｅｎｂｅｒｇ ＆ Ｓａｙｅｒｓ理论，在碳酸盐岩介质各向同性背景中加入裂缝系统，修正计算出的模量。
从前面1,2步建立的等效岩石背景基质中，获得基质的体积模量和剪切模量。在其基础上加入裂缝的影响，修正计算的模量。
目前有两种应用广泛的裂缝模型，Ｈｕｄｓｏｎ 模型和Ｓｃｈｏｅｎｂｅｒｇ线性滑动模型。Ｈｕｄｓｏｎ理论是建立在整体平均的基础上。该模型假设弹性介质的内部分布着薄硬币形状的椭球缝隙。Ｈｕｄｓｏｎ模型利用裂缝密度和缝隙纵横比来描述裂缝系统的结构。裂缝密度e的定义如下：
\[e = \frac{{3\phi }}{{4\pi a}}
\]
\]
其中，φ代表孔隙度，α 是裂隙纵横比．则Ｈｕｄｓｏｎ模型的弹性矩阵犆C表示为：
\[C^{eff}_{ij}=C^{0}_{ij}+C^{1}_{ij}+C^{2}_{ij}
\]
\]
式中，C0为各向同性背景弹性矩阵，C1,C2为引入裂隙的各向异性一阶校正和二阶校正。
Ｓｃｈｏｅｎｂｅｒｇ模型的弹性矩阵由各向同性参数λ 和μ，以及垂直于裂缝面和平行于裂缝面的弹性差值ΔＮ 和ΔＴ 进行描述。表示裂缝对骨架岩石各向同性特征的影响。在裂缝诱导等效各向异性介质中，裂缝等效介质弹性参数矩阵C****表示为背景各向同性介质弹性参数Cb与裂缝介质弹性参数Cf之差；Ｓｃｈｏｅｎｂｅｒｇ线性滑动模型的弹性矩阵为：
△N和△T分别表示裂缝的法向弱度和切向弱度
下面列出Ｓｃｈｏｅｎｂｅｒｇ模型的计算代码：
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////计算裂缝密度，ΔＮ 和ΔＴ
for (int i = 0; i < nsample; i++)
{
Sg[i] = 1 - Sw[i];
e[i] = 3 * FraPOR[i] / (4 * PI*arf);
}
for (int i = 0; i < nsample; i++)
{
g[i] = (Vs_model[i] * Vs_model[i] / Vp_model[i] / Vp_model[i]);
deltaN[i] = 4 * e[i] / (3 * g[i] * (1 - g[i])*(1 + (k_f + 1.3333*miu_f) / (PI*(1 - g[i])*miu[i] * arf)));
deltaT[i] = 16 * e[i] / (3 * (3 - 2 * g[i])*(1 + 4 * miu_f / (PI*(3 - 2 * g[i])*miu[i] * arf)));
}
4.利用Ｗｏｏｄ公式将孔隙流体进行混合，计算出混合流体的体积模量；
wood公式进行的平均等效模型表达式为：
\[\begin{array}{l}
V = \sqrt {\frac{{{K_R}}}{\rho }} \\
\frac{1}{{{K_R}}} = \sum\limits_{i = 1}^N {\frac{{{f_i}}}{{{K_i}}}} \\
\rho = \sum\limits_{i = 1}^N {{f_i}} {\rho _i}
\end{array}\
\]
V = \sqrt {\frac{{{K_R}}}{\rho }} \\
\frac{1}{{{K_R}}} = \sum\limits_{i = 1}^N {\frac{{{f_i}}}{{{K_i}}}} \\
\rho = \sum\limits_{i = 1}^N {{f_i}} {\rho _i}
\end{array}\
\]
其中f，K，ρ分别为成分的体积含量，体积模量和密度；
// wood流体混合模式
{
if (nFluidCom == 0)
{
Kf = 1 / (So / Ko + (1 - So) / Kw);//水油
}
if (nFluidCom == 1)
{
Kf = 1 / (Sg / Kg + (1 - Sg) / Kw);//水气
}
else
{
Kf = 1 / (So / Ko + Sg / Kg + (1 - Sg - So) / Kw);//气油水
}
}
（ 5 ）计算裂缝型碳酸盐岩的弹性参数和各向异性参数．
由于研究区域只含有气，不用进行流体替换。
各向异性参数的表达式为：
\[\begin{array}{l}
\varepsilon = \left( {C[1][1] - C[3][3]} \right)/\left( {2*C[3][3]} \right);\\
\delta = {\rm{((C[1][3] + C[5][5])*(C[1][3] + C[5][5]) - (C[3][3] - C[5][5])*(C[3][3] - C[5][5])) / (2 * C[3][3] * (C[3][3] - C[5][5]));}}\\
\gamma {\rm{ = (C[6][6] - C[4][4]) / (2 * C[4][4]);}}
\end{array}
\]
\varepsilon = \left( {C[1][1] - C[3][3]} \right)/\left( {2*C[3][3]} \right);\\
\delta = {\rm{((C[1][3] + C[5][5])*(C[1][3] + C[5][5]) - (C[3][3] - C[5][5])*(C[3][3] - C[5][5])) / (2 * C[3][3] * (C[3][3] - C[5][5]));}}\\
\gamma {\rm{ = (C[6][6] - C[4][4]) / (2 * C[4][4]);}}
\end{array}
\]
对应的计算代码：
epsilon[i] = (CHTI[0][0] - CHTI[2][2]) / 2 / CHTI[2][2];
delta[i] = ((CHTI[0][2] + CHTI[4][4])*(CHTI[0][2] + CHTI[4][4]) - (CHTI[2][2] - CHTI[4][4])*(CHTI[2][2] - CHTI[4][4])) / (2 * CHTI[2][2] * (CHTI[2][2] - CHTI[4][4]));
gamma[i] = (CHTI[5][5] - CHTI[3][3]) / (2 * CHTI[3][3]);
计算所得的各向异性参数与实测的吻合度较高。
参考及其引用文献：
张广智，陈怀震，王琪等．基于碳酸盐岩裂缝岩石物理模型的横波速度和各向异性参数预测．地球物理学报，２０１３，５６（５）：１７０７
１７１５，ｄｏｉ：１０．６０３８／ｃｊｇ２０１３０５２８．
岩石物理手册.Gray Mavko and so on.
