CN103675916B - 一种三分量检波器埋置方向高精度校正的方法 - Google Patents

Abstract

本发明是一种利用地面多分量地震数据的偏振特性三分量检波器埋置方向高精度校正的方法。拾取初至时间，计算每一道三分量地震记录给定初至时窗的协方差矩阵，找出绝对值最大的非对角线元素，确定的旋转角度，构造一次旋转矩阵，对协方差矩阵的三个特征值按照从大到小的顺序进行排序，循环计算每一道的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角、方位角，计算检波器偏离角度，根据偏离角度对采集的三分量地震数据进行校正。本发明有很高的精度，有效地减小了随机误差，提高了计算结果的可靠性和校正的计算效率。

Description

一种三分量检波器埋置方向高精度校正的方法

技术领域

本发明属于石油物探地震资料处理过程中提高信噪比的技术，是一种利用地面多分量地震数据的偏振特性三分量检波器埋置方向高精度校正的方法。

背景技术

随着油气勘探开发的不断深入，常规地震勘探技术难以解决某些岩性油气藏、裂缝性油气藏等复杂油气藏的勘探问题，于是勘探工作者纷纷把目光投向多分量地震勘探技术。多分量地震勘探是指用纵/横波震源激发，用三分量检波器记录地震纵波、横波(包括快、慢横波)和转换波，从而使野外记录的地震数据信息更为丰富，为地质构造成像、裂隙和孔道确定、储层岩性解释等提供特定的信息。

根据弹性波理论，在各向同性介质中，在空间的任意方向上存在着3种不同类型的波，分别为P波、SV波和SH波，其中P波引起的质点振动在炮检连线所在的垂向平面内，质点振动方向与地震波的传播方向一致，SV波引起的质点振动也在炮检连线所在的垂向平面内，但质点振动方向垂直地震波的传播方向，SH波引起的质点振动垂直炮检连线所在的垂向平面内。由于不同类型的波有不同的质点振动方向，所以在多分量地震数据的采集过程中，需要严格控制三分量检波器的埋置质量，也就是利用三分量检波器上的方向标记和调平水泡将三分量检波器的垂直分量调整到垂直状态，将三分量检波器的X分量平行观测系统设计的测线方向，只有在这样的前提下，通过对两个水平分量数据做旋转，才能将多分量地震勘探中所要得到的SV波的能量集中在以炮检连线的径向分量上，将SH波能量集中在垂直炮检连线方向的切向分量上，垂直分量上则主要是P波能量。但是，在实际多分量数据采集过程中，由于野外地质施工条件的复杂性，很难保证每个三分量检波器都完全按照要求统一埋置，埋置的三分量检波器调平水泡未调平，方向标记未对准规定的方向或接收道极性反转的情况很难完全避免。这样就会使不同类型的地震波在三个分量上都有能量投影，也就是在接收的多分量数据中，Z分量的记录中含有横波成分，X分量的记录中含有纵波和SH波的能量，Y分量的记录中含有纵波的和SV波的成分。

针对这一问题，现有的文献大都是根据多波数据的偏振特征，利用多分量数据的初至能量对两个水平分量数据进行旋转，将多分量地震勘探中所要得到的SV波的能量集中在以炮检连线的径向分量上，将SH波能量集中在垂直炮检连线方向的切向分量上，但是这种做法只能解决三分量检波器的X分量与炮检连线方向不完全一致所带来的问题，由Z分量不垂直所造成的垂直分量和水平分量上的能量互相投影所带来的问题并没有解决。

常用的地震数据处理方法是利用多分量地震数据初至能量扫描的方法对三分量检波器进行重定向，其实现过程为：对于任意两个分量的地震记录，给定不同的旋转角，对这两个分量的地震记录进行旋转并计算给定初至时窗内的能量，通过角度扫描得到极小能量对应的角度，将该角度当作三分量检波器需要校正的角度。该方法原理简单，易于操作。但是，由于该方法是利用能量极小的规则扫描检波器三个分量的偏离角度，在0至360度的范围内存在两个相等的极值，这需要利用水听器记录分量的初至来判断偏离角度的极性，如果没有水听器分量，该方法的使用会受到限制。并且角度离散扫描的计算方法难以兼顾效率和精度。

发明内容

本发明目的在于提供一种使用不受到限制，兼顾效率和精度的三分量检波器埋置方向高精度校正的方法。

本发明通过以下步骤实现：

1)利用正交三分量检波器采集三分量地震记录；

2)拾取地震记录初至时间；

步骤2)所述的初至时间拾取是在Z分量或水平分量坐标旋转后的X分量地震记录上拾取每一道地震记录的初至时间fbt(i)，i是地震数据道顺序号，i＝1,2,…,Ntrace，Ntrace是每一个分量地震记录的总道数；

3)计算每一道三分量地震记录给定初至时窗的协方差矩阵

步骤3)所述的三分量地震记录给定初至时窗的协方差矩阵的计算公式为：

$c_{ij}^{\left(0\right)}=\frac{1}{n}\underset{i=i_1}{\overset{t_2}{\Σ}}{S}_i\left(t\right)S_j\left(t\right)$

是协方差矩阵的各个元素，i，j是协方差矩阵的行列序号，i，j＝1,2,3，S_i(t)和S_j(t)均是三分量地震记录，i和j是分量序号，i＝1,2,3，j＝1,2,3，t₁和t₂是计算时窗的起始和终止样点，t₁＝fbt(i)/dt，t₂＝(fbt(i)+twindow)dt，dt是地震记录采样间隔，twindow是给定的计算时窗长度，n＝t₂-t₁-1，根据计算的协方差矩阵元素得到3×3阶实对称矩阵，也就是三分量地震记录给定时窗的协方差矩阵

${\stackrel{\‾}{C}}^{\left(0\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(0\right)}& c_{12}^{\left(0\right)}& c_{13}^{\left(0\right)}\\ c_{21}^{\left(0\right)}& c_{22}^{\left(0\right)}& c_{23}^{\left(0\right)}\\ c_{31}^{\left(0\right)}& c_{32}^{\left(0\right)}& c_{33}^{\left(0\right)}\end{array}\right].$

4)在协方差矩阵中找出绝对值最大的非对角线元素p,q是协方差矩阵的行列序号，1≤p≤3，1≤q≤3，然后利用相应的矩阵元素计算出sinθ和cosθ，其中，θ是由协方差矩阵元素和确定的旋转角度，sinθ和cosθ相应的计算过程变量；

步骤4)所述的sinθ和cosθ的计算公式为：

$\left\{\begin{array}{c}tg2\mathrm{\θ}=\frac{2c_{pq}^{\left(0\right)}}{c_{pp}^{\left(0\right)}-c_{qq}^{\left(0\right)}}\\ \cos 2\mathrm{\θ}=\frac{1}{\sqrt{1+{\mathrm{tg}}^{2}2\mathrm{\θ}}}\\ \sin 2\mathrm{\θ}=\cos 2\mathrm{\θ}\·tg2\mathrm{\θ}\\ \cos \mathrm{\θ}=\sqrt{\frac{1+cos2\mathrm{\θ}}{2}}\\ \sin \mathrm{\θ}=\frac{\sin 2\mathrm{\θ}}{2\cos \mathrm{\θ}}\end{array}\right..$

5)利用步骤4)得到的sinθ和cosθ以及步骤3)得到的协方差矩阵计算一次旋转后的协方差矩阵的元素i,j＝1,2,3；

步骤5)所述的一次旋转后的协方差矩阵的各个元素的计算公式为：

$\left\{\begin{array}{c}{c}_{pj}^{\left(1\right)}=c_{pl}^{\left(0\right)}\cos \mathrm{\θ}+c_{ql}^{\left(0\right)}\sin \mathrm{\θ}=c_{jp}^{\left(1\right)},(l\≠i,j)\\ c_{qj}^{\left(1\right)}=-c_{pl}^{\left(01\right)}\sin \mathrm{\θ}+c_{ql}^{\left(0\right)}\cos \mathrm{\θ}=c_{jp}^{\left(1\right)}\\ c_{pp}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\cos }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\sin }^2\mathrm{\θ}+2c_{pl}^{\left(0\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ c_{qq}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\sin }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\cos }^2\mathrm{\θ}-2c_{pl}^{\left(01\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ c_{lm}^{\left(1\right)}=c_{ml}^{\left(0\right)}=c_{ml},(l,m\≠i,j)\\ c_{pq}^{\left(1\right)}=c_{qp}^{\left(1\right)}=\frac{1}{2}(c_{qq}^{\left(0\right)}-c_{pp}^{\left(0\right)})\sin 2\mathrm{\θ}+c_{pq}^{\left(0\right)}({\cos }^2\mathrm{\θ}-{\sin }^2\mathrm{\θ})\end{array}\right.$

其中，l,m是协方差矩阵的行列序号，l＝1,2,3,m＝1,2,3；

根据计算的各个矩阵元素得到一次旋转后的协方差矩阵为：

${\stackrel{\‾}{C}}^{\left(1\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(1\right)}& c_{12}^{\left(1\right)}& c_{13}^{\left(1\right)}\\ c_{21}^{\left(1\right)}& c_{22}^{\left(1\right)}& c_{23}^{\left(1\right)}\\ c_{31}^{\left(1\right)}& c_{32}^{\left(1\right)}& c_{33}^{\left(1\right)}\end{array}\right].$

6)构造一次旋转矩阵R⁽¹⁾(p,q,θ)；

步骤6)所述的构造取一次旋转矩阵R⁽¹⁾(p,q,θ)的过程为：利用步骤4)得到的sinθ和cosθ，在3×3阶单位阵I的p行p列，q行q列的交叉位置上置入 其它对角线元素为1，其它非对角线元素为0，得到的一次旋转矩阵为：

7)计算一次旋转后的协方差矩阵的非对角线元素的平方和；

步骤7)所述的一次旋转后的协方差矩阵的非对角线元素的平方和的计算公式为：

$E^{\left(1\right)}\left(A\right)=\underset{i\≠j}{\overset{3}{\underset{i,j=1}{\Σ}}}{c}_{ij}^{\left(1\right)}\·c_{ij}^{\left(1\right)}$

其中是步骤5)计算的一次旋转后的协方差矩阵的矩阵元素，i，j是该矩阵的行列序号。

8)旋转变换迭代计算；

步骤8)所述的旋转变换迭代计算过程为：判断步骤7)得到的非对角线元素的平方和E⁽¹⁾(A)是否小于给定的阀值δ，如果不小于δ，用步骤5)计算的一次旋转后的协方差矩阵更新步骤3)中的初始协方差矩阵重复步骤4)-7)，如此类推直到非对角线元素的平方和小于给定的阀值δ，迭代计算结束。

9)构造协方差矩阵的特征值和特征向量；

步骤9)所述的协方差矩阵的特征值和特征向量的构造过程为：设经过m次旋转变换计算后，迭代计算终止，则可得到m次旋转后的协方差矩阵以及每次迭代后的旋转矩阵R⁽¹⁾,R⁽²⁾,…,R^(m)，记变换矩阵P^(m)＝R⁽¹⁾R⁽²⁾…R^(m)，则可以得到协方差矩阵的特征值为是协方差矩阵的各个元素，特征值λ_j对应的特征向量P_j为变换矩阵P^(m)的第j列的各个元素，既，j＝1,2,3。

10)确定接收点处的质点振动的主方向；

步骤10)所述的质点振动主方向的确定过程是：对步骤9)得到的协方差矩阵的三个特征值按照从大到小的顺序进行排序，使λ₁≥λ₂≥λ₃，然后，得到排序后的特征值对应的特征向量V₁，V₂，V₃，其中，V₂＝(v₁₂ v₂₂ v₃₂)^T，V₃＝(v₁₃ v₂₃ v₃₃)^T，v_ij是特征向量的元素，最大特征值λ₁对应的特征向量V₁即为质点振动的主方向；

11)特征向量规则化；

步骤11)所述的特征向量规则化所用的公式为：

$W_j=\frac{V_j}{\[V_j,V_j\]}=\frac{V_j}{\underset{i=1}{\overset{3}{\Σ}}{v}_{ij}^2},j=1,2,3$

其中，V_j为由步骤10)得到的排序后的特征向量，v_ij为特征向量的元素。得到的规则化后的三个特征向量分别W₁，W₂，W₃，W₁＝(w₁₁ w₂₁ w₃₁)^T，W₂＝(w₁₂ w₂₂ w₃₂)^T，W₃＝(w₁₃ w₂₃w₃₃)^T，w_ij是规则化后的特征向量的元素；

12)计算特征向量W₁在XY平面内的投影与X轴的夹角以及特征向量W₁与Z轴的夹角，其中特征向量W₁为最大特征值λ₁对应的特征向量V₁规则化后的特征向量；

步骤12)所述的特征向量W₁在XY平面内的投影与X轴的夹角的计算公式为：

特征向量W₁与Z轴的夹角α的计算公式为：

α＝arccos(w₃₁)

其中，w₁₁，w₂₁，w₃₁为由步骤11)计算得到的规则化后的主特征向量W₁的元素。

13)循环计算每一道的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角；

步骤13)所述的计算过程为：设采集的地震记录有N炮，每炮有M个接收点，每个接收点三分量接收。

按照步骤3)至步骤12)的实现过程，计算出第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角α_ij，i＝1,…,N，j＝1,…,M。

14)循环计算每一道地震数据的方位角；

步骤14)所述的每一道地震数据的方位角的计算公式为：

其中，(x_i,y_i)为第i炮的坐标，(x_j,y_j)为第j个接收点处理的坐标，i＝1,…,N，j＝1,…,M。

15)计算每一个接收点处的三分量检波器的X分量与X坐标轴的偏离角度；

步骤15)所述的每一个接收点处三分量检波器的X分量与X坐标轴的偏离角度Ψ_j的计算公式为：

其中，是由步骤14)得到的第i炮，第j个接收点处的方位角，是由步骤13)得到的第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角，N是总炮数，M是总接收点数。

16)计算每一个接收点处的三分量检波器的Z分量与垂直方向的偏离角度；

步骤16)所述的每一个接收点处三分量检波器的Z分量与垂直方向的偏离角度的计算公式为：

其中，α_ij是由步骤13)得到的第i炮，第j个接收点处的主特征向量与Z轴的夹角，N是总炮数。

17)根据由步骤15)和步骤16)计算得到的三分量检波器的X分量与X坐标轴的偏离角度Ψ和Z分量与垂直方向的偏离角度Φ，对采集的三分量地震记录进行校正。

步骤17)所述的对三分量地震记录进行校正的公式为

$\left[\begin{array}{c}{C}_1\left(t\right)\\ C_2\left(t\right)\\ C_3\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}{L}_1& L_2& L_3\\ M_1& M_2& M_3\\ N_1& N_2& N_3\end{array}\right]\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ S_3\left(t\right)\end{array}\right]$

其中，

$L_2=\frac{sin\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)}{\mathrm{\Δ}}$

$L_3=\frac{-sin\mathrm{\Ψ}}{\mathrm{\Δ}}$

$M_1=\frac{cos\mathrm{\Ψ}cos\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)+sin\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_3=\frac{-cos\mathrm{\Ψ}}{\mathrm{\Δ}}$

$N_1=\frac{cos\mathrm{\Ψ}cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$N_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}\sin \mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

N₃＝0

$\mathrm{\Δ}=\left|\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\\ n_1& n_2& n_3\end{array}\right|$

l₁＝cosΨsinΦ m₁＝sinΨsinΦ n₁＝cosΦ

l₂＝cosΨcosΦ m₂＝sinΨcosΦ n₂＝-sinΦ

l₃＝m₁n₂-n₁m₂ m₃＝l₁n₂-n₁l₂ n₃＝l₁m₂-m₁l₂

S₁(t)、S₂(t)和S₃(t)是采集的三分量地震记录的X、Y和Z分量，C₁(t)、C₂(t)和C₃(t)是校正后的三分量地震记录，t＝1,2,…,N，N为每道地震记录的采样点数，Ψ是步骤17)计算得到的三分量检波器的X分量与X坐标轴的偏离角度，Φ是Z分量与垂直方向的偏离角度。通过对三分量数据进行校正，将水平分量中的纵波能量投影到垂直分量上，将垂直分量中的横波能量投影到水平分量上。

本发明具有如下特点：

(1)本发明利用地面地震多分量数据的偏振特性，通过计算三分量数据给定时窗的协方差矩阵的特征值和特征向量，来确定接收点处的质点振动的主方向，进而确定三分量检波器垂直分量偏离垂直方向的夹角以及X分量偏离X方向的夹角，在计算协方差矩阵的特征值和特征向量时，采用了收敛快、算法稳定的雅可比计算方法，得到的结果具有很高的精度。

(2)本发明利用多个炮点的地震数据来统计一个接收点上的三分量检波器垂直分量偏离垂直方向的夹角以及X分量偏离X方向的夹角，有效地减小了随机误差，提高了计算结果的可靠性。

(3)本发明推导了由三分量检波器垂直分量偏离垂直方向的夹角以及X分量偏离X方向的夹角表示的三分量地震记录校正的计算公式，提高了三分量地震记录校正的计算效率。

本发明利用地面地震多分量数据的偏振特性，通过计算三分量数据给定时窗的协方差矩阵的特征值和特征向量，来确定接收点处的质点振动的主方向，进而确定三分量检波器垂直分量偏离垂直方向的夹角以及X分量偏离X方向的夹角。通过利用推导的三分量地震记录标定公式对三分量数据进行校正，可以有效的将水平分量中的纵波能量投影到垂直分量上，以及将垂直分量中的横波能量投影到水平分量上。

附图说明

图1三分量检波器埋置方向校正角度示意图；

图2用本发明以前的共偏移距三分量合成记录；

图3用本发明以后的共偏移距三分量合成记录；

图4应用本发明以前的实际地面地震炮集记录(X分量)；

图5应用本发明以后的实际地面地震炮集记录(X分量)。

具体实施方式

以下结合附图详细说明本发明。

本发明的具体实施方式为：

1)利用正交三分量检波器采集三分量地震记录；

2)拾取地震记录初至时间；

拾取X和Y两个水平分量坐标旋转后的X分量地震记录的初至时间fbt(i)，i是地震数据道顺序号，i＝1,2,…,Ntrace，Ntrace是X分量地震记录的总道数；

3)计算三分量地震记录给定初至时窗的协方差矩阵协方差矩阵元素的计算公式为：

$c_{ij}^{\left(0\right)}=\frac{1}{n}\underset{t=t_1}{\overset{t_2}{\Σ}}{S}_i\left(t\right)S_j\left(t\right)$

是协方差矩阵的各个元素，i，j是协方差矩阵的行列序号，i，j＝1,2,3，S_i(t)和S_j(t)均是三分量地震记录，i和j是分量序号，i＝1,2,3，j＝1,2,3，t₁和t₂是计算时窗的起始和终止样点，t₁＝fbt(i)/dt，t₂＝(fbt(i)+twindow)/dt，dt是地震记录采样间隔，twindow是给定的计算时窗长度，n＝t₂-t₁-1。根据计算的协方差矩阵元素得到3×3阶实对称矩阵，也就是三分量地震记录给定时窗的协方差矩阵

${\stackrel{\‾}{C}}^{\left(0\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(0\right)}& c_{12}^{\left(0\right)}& c_{13}^{\left(0\right)}\\ c_{21}^{\left(0\right)}& c_{22}^{\left(0\right)}& c_{23}^{\left(0\right)}\\ c_{31}^{\left(0\right)}& c_{32}^{\left(0\right)}& c_{33}^{\left(0\right)}\end{array}\right].$

4)在协方差矩阵中找出绝对值最大的非对角线元素p,q是协方差矩阵的行列序号，1≤p≤3，1≤q≤3，然后利用相应的矩阵元素计算出sinθ和cosθ；

步骤4)所述的sinθ和cosθ的计算公式为：

$\left\{\begin{array}{c}tg2\mathrm{\θ}=\frac{2c_{pq}^{\left(0\right)}}{c_{pp}^{\left(0\right)}-c_{qq}^{\left(0\right)}}\\ \cos 2\mathrm{\θ}=\frac{1}{\sqrt{1+{\mathrm{tg}}^{2}2\mathrm{\θ}}}\\ \sin 2\mathrm{\θ}=\cos 2\mathrm{\θ}\·tg2\mathrm{\θ}\\ \cos \mathrm{\θ}=\sqrt{\frac{1+cos2\mathrm{\θ}}{2}}\\ \sin \mathrm{\θ}=\frac{\sin 2\mathrm{\θ}}{2\cos \mathrm{\θ}}\end{array}\right..$

5)利用步骤4)得到的sinθ和cosθ以及步骤3)得到的协方差矩阵计算一次旋转后的协方差矩阵的元素i,j＝1,2,3；

步骤5)所述的一次旋转后的协方差矩阵的各个元素的计算公式为：

$\left\{\begin{array}{c}{c}_{pj}^{\left(1\right)}=c_{pl}^{\left(0\right)}\cos \mathrm{\θ}+c_{ql}^{\left(0\right)}\sin \mathrm{\θ}=c_{jp}^{\left(1\right)},(l\≠i,j)\\ c_{qj}^{\left(1\right)}=-c_{pl}^{\left(01\right)}\sin \mathrm{\θ}+c_{ql}^{\left(0\right)}\cos \mathrm{\θ}=c_{jp}^{\left(1\right)}\\ c_{pp}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\cos }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\sin }^2\mathrm{\θ}+2c_{pl}^{\left(0\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ c_{qq}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\sin }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\cos }^2\mathrm{\θ}-2c_{pl}^{\left(01\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ c_{lm}^{\left(1\right)}=c_{ml}^{\left(0\right)}=c_{ml},(l,m\≠i,j)\\ c_{pq}^{\left(1\right)}=c_{qp}^{\left(1\right)}=\frac{1}{2}(c_{qq}^{\left(0\right)}-c_{pp}^{\left(0\right)})\sin 2\mathrm{\θ}+c_{pq}^{\left(0\right)}({\cos }^2\mathrm{\θ}-{\sin }^2\mathrm{\θ})\end{array}\right.$

其中，l,m是协方差矩阵的行列序号，l＝1,2,3,m＝1,2,3；

根据计算的各个矩阵元素得到一次旋转后的协方差矩阵为：

${\stackrel{\‾}{C}}^{\left(1\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(1\right)}& c_{12}^{\left(1\right)}& c_{13}^{\left(1\right)}\\ c_{21}^{\left(1\right)}& c_{22}^{\left(1\right)}& c_{23}^{\left(1\right)}\\ c_{31}^{\left(1\right)}& c_{32}^{\left(1\right)}& c_{33}^{\left(1\right)}\end{array}\right].$

6)构造一次旋转矩阵R⁽¹⁾(p,q,θ)；

步骤6)所述的构造取一次旋转矩阵R⁽¹⁾(p,q,θ)的过程为：利用步骤4)得到的sinθ和cosθ，在3×3阶单位阵I的p行p列，q行q列的交叉位置上置入 其它对角线元素为1，其它非对角线元素为0，得到的一次旋转矩阵为：

7)计算一次旋转后的协方差矩阵的非对角线元素的平方和；

步骤7)所述的一次旋转后的协方差矩阵的非对角线元素的平方和的计算公式为：

$E^{\left(1\right)}\left(A\right)=\underset{i\≠j}{\overset{3}{\underset{i,j=1}{\Σ}}}{c}_{ij}^{\left(1\right)}\·c_{ij}^{\left(1\right)}$

其中是步骤5)计算的一次旋转后的协方差矩阵的矩阵元素，i，j是该矩阵的行列序号。

8)旋转变换迭代计算；

步骤8)所述的旋转变换迭代计算过程为：判断步骤7)得到的非对角线元素的平方和E⁽¹⁾(A)是否小于给定的阀值δ，如果不小于δ，用步骤5)计算的一次旋转后的协方差矩阵更新步骤3)中的初始协方差矩阵重复步骤4)，5)，6)，7)，如此类推直到非对角线元素的平方和小于给定的阀值δ，迭代计算结束。

9)构造协方差矩阵的特征值和特征向量；

步骤9)所述的协方差矩阵的特征值和特征向量的构造过程为：设经过m次旋转变换计算后，迭代计算终止，则可得到m次旋转后的协方差矩阵以及每次迭代后的旋转矩阵R⁽¹⁾,R⁽²⁾,…,R^(m)，记变换矩阵P^(m)＝R⁽¹⁾R⁽²⁾…R^(m)，则可以得到协方差矩阵的特征值为是协方差矩阵的各个元素，特征值λ_j对应的特征向量P_j为变换矩阵P^(m)的第j列的各个元素，既j＝1,2,3。

10)确定接收点处的质点振动的主方向；

步骤10)所述的质点振动主方向的确定过程是：对步骤9)得到的协方差矩阵的三个特征值按照从大到小的顺序进行排序，使λ₁≥λ₂≥λ₃，然后，得到排序后的特征值对应的特征向量V₁，V₂，V₃，其中，V₂＝(v₁₂ v₂₂ v₃₂)^T，V₃＝(v₁₃ v₂₃ v₃₃)^T，v_ij是特征向量的元素，最大特征值λ₁对应的特征向量V₁即为质点振动的主方向；

11)特征向量规则化；

步骤11)所述的特征向量规则化所用的公式为：

$W_j=\frac{V_j}{\[V_j,V_j\]}=\frac{V_j}{\underset{i=1}{\overset{3}{\Σ}}{v}_{ij}^2},j=1,2,3$

其中，V_j为由步骤10)得到的排序后的特征向量，v_ij为特征向量的元素。得到的规则化后的三个特征向量分别W₁，W₂，W₃，W₁＝(w₁₁ w₂₁ w₃₁)^T，W₂＝(w₁₂ w₂₂ w₃₂)^T，W₃＝(w₁₃ w₂₃w₃₃)^T，w_ij是规则化后的特征向量的元素；

12)计算特征向量W₁在XY平面内的投影与X轴的夹角以及特征向量W₁与Z轴的夹角，其中特征向量W₁为最大特征值λ₁对应的特征向量V₁规则化后的特征向量；

步骤12)所述的笛卡尔坐标系下特征向量W₁在XY平面内的投影与X轴的夹角的计算公式为：

特征向量W₁与Z轴的夹角α的计算公式为：

α＝arccos(w₃₁)

其中，w₁₁，w₂₁，w₃₁为由步骤11)计算得到的规则化后的主特征向量W₁的元素。

13)循环计算每一道的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角；

步骤13)所述的计算过程为：设采集的地震记录有N炮，每炮有M个接收点，每个接收点三分量接收。按照步骤3)～步骤12)的实现过程，计算出第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角α_ij，i＝1,…,N，j＝1,…,M。

14)循环计算每一道地震数据的方位角；

步骤14)所述的第i炮，第j个接收点处的方位角的计算公式为：

其中，(x_i,y_i)为第i炮的坐标，(x_j,y_j)为第j个接收点处理的坐标，i＝1,…,N，j＝1,…,M。

15)计算每一个接收点处的三分量检波器的X分量与X坐标轴的偏离角度,如图1所示；

步骤15)所述的第j个接收点处三分量检波器的X分量与X坐标轴的偏离角度Ψ_j的计算公式为：

其中，是由步骤14)得到的第i炮，第j个接收点处的方位角，是由步骤13)得到的第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角，N是总炮数，M是总接收点数。

16)计算每一个接收点处的三分量检波器的Z分量与垂直方向的偏离角度，如图1所示；

步骤16)所述的第j个接收点处三分量检波器的Z分量与垂直方向的偏离角度的计算公式为：

${\mathrm{\Φ}}_j=\frac{\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\α}}_{ij}}{N},j=1,...,M$

其中，α_ij是由步骤13)得到的第i炮，第j个接收点处的主特征向量与Z轴的夹角，N是总炮数。

17)根据由步骤15)和步骤16)计算得到的三分量检波器的X分量与X坐标轴的偏离角度Ψ和Z分量与垂直方向的偏离角度Φ，对采集的三分量地震记录进行校正。

步骤17)所述的对三分量地震记录进行校正的公式为

$\left[\begin{array}{c}{C}_1\left(t\right)\\ C_2\left(t\right)\\ C_3\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}{L}_1& L_2& L_3\\ M_1& M_2& M_3\\ N_1& N_2& N_3\end{array}\right]\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ S_3\left(t\right)\end{array}\right]$

其中，

$L_2=\frac{sin\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)}{\mathrm{\Δ}}$

$L_3=\frac{-sin\mathrm{\Ψ}}{\mathrm{\Δ}}$

$M_1=\frac{\cos \mathrm{\Ψ}\cos \mathrm{\Φ}(l_1{m}_2-m_1{l}_2)+\sin \mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_3=\frac{-cos\mathrm{\Ψ}}{\mathrm{\Δ}}$

$N_1=\frac{cos\mathrm{\Ψ}cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$N_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}sin\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

N₃＝0

$\mathrm{\Δ}=\left|\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\\ n_1& n_2& n_3\end{array}\right|$

l₁＝cosΨsinΦ m₁＝sinΨsinΦ n₁＝cosΦ

l₂＝cosΨcosΦ m₂＝sinΨcosΦ n₂＝-sinΦ

l₃＝m₁n₂-n₁m₂ m₃＝l₁n₂-n₁l₂ n₃＝l₁m₂-m₁l₂

S₁(t)、S₂(t)和S₃(t)是采集的三分量地震记录的X、Y和Z分量，C₁(t)、C₂(t)和C₃(t)是校正后的三分量地震记录，t＝1,2,…,N，N为每道地震记录的采样点数，Ψ是步骤17)计算得到的三分量检波器的X分量与X坐标轴的偏离角度，Φ是Z分量与垂直方向的偏离角度。通过对三分量数据进行校正，将水平分量中的纵波能量投影到垂直分量上，将垂直分量中的横波能量投影到水平分量上。

图2为三分量合成地震记录的共偏移距道集，可以看出，由于受检波器放置方向的影响，本来只应当记录到X分量上的初至波也被投影到Z分量和Y分量上，图3为该三分量合成地震记录应用本发明以后的共偏移距道集，可以看出，通过检波器埋置方向校正，已成功地将Z分量和Y分量上的能量还原到X分量上。图4为实际地面地震炮集记录的X分量，图5为应用本发明对该地震记录进行检波器埋置方向校正后的X分量，可以看出，在校正后的X分量上，P波得到了消弱，转换波同相轴的连续性得到的改善。

Claims (13)

1.一种三分量检波器埋置方向高精度校正的方法，特点是采用以下步骤：

1)利用正交三分量检波器采集三分量地震记录；

2)拾取地震记录初至时间；

3)计算每一道三分量地震记录给定初至时窗的协方差矩阵

4)在协方差矩阵中找出绝对值最大的非对角线元素p，q是协方差矩阵的行列序号，1≤p≤3，1≤q≤3，然后利用相应的矩阵元素计算出sinθ和cosθ；其中，θ是由协方差矩阵元素和确定的旋转角度，sinθ和cosθ相应的计算过程变量；

5)利用步骤4)得到的sinθ和cosθ以及步骤3)得到的协方差矩阵计算一次旋转后的协方差矩阵的元素i,j＝1,2,3；

6)构造一次旋转矩阵R⁽¹⁾(p,q,θ)；

7)计算一次旋转后的协方差矩阵的非对角线元素的平方和E⁽¹⁾(A)；

8)旋转变换迭代计算；

9)构造协方差矩阵的特征值和特征向量；

10)对步骤9)得到的协方差矩阵的三个特征值按照从大到小的顺序进行排序，使λ₁≥λ₂≥λ₃，然后，得到排序后的特征值对应的特征向量V₁，V₂，V₃，其中，V₁＝(v₁₁v₂₁v₃₁)^T，V₂＝(v₁₂v₂₂v₃₂)^T，V₃＝(v₁₃v₂₃v₃₃)^T，v_ij是特征向量的元素，最大特征值λ₁对应的特征向量V₁即为质点振动的主方向；

11)特征向量规则化；

12)计算特征向量W₁在XY平面内的投影与X轴的夹角以及特征向量W₁与Z轴的夹角α，其中特征向量W₁为最大特征值λ₁对应的特征向量V₁规则化后的特征向量；

13)循环计算每一道的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角；

14)循环计算每一道地震数据的方位角；

15)计算每一个接收点处的三分量检波器的X分量与X坐标轴的偏离角度，该偏离角度Ψ_j的计算公式为：

其中，是由步骤14)得到的第i炮，第j个接收点处的方位角，是由步骤13)得到的第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角，N是总炮数，M是总接收点数；

16)计算每一个接收点处的三分量检波器的Z分量与垂直方向的偏离角度，该偏离角度Φ_j的计算公式为：

${\mathrm{\Φ}}_j=\frac{\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\α}}_{ij}}{N},j=1,...,M$

其中，α_ij是由步骤13)得到的第i炮，第j个接收点处的主特征向量与Z轴的夹角，N是总炮数；

17)根据由步骤15)和步骤16)计算得到的三分量检波器的X分量与X坐标轴的偏离角度Ψ和Z分量与垂直方向的偏离角度Φ，对采集的三分量地震记录进行校正，进行校正的公式为：

$\left[\begin{array}{c}{C}_1\left(t\right)\\ C_2\left(t\right)\\ C_3\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}{L}_1& L_2& L_3\\ M_1& M_2& M_3\\ N_1& N_2& N_3\end{array}\right]\left[\begin{array}{c}{S}_1\left(t\right)\\ S_2\left(t\right)\\ S_3\left(t\right)\end{array}\right]$

其中

$L_2=\frac{sin\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)}{\mathrm{\Δ}}$

$L_3=\frac{-sin\mathrm{\Ψ}}{\mathrm{\Δ}}$

$M_1=\frac{cos\mathrm{\Ψ}cos\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)+sin\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{m}_2-m_1{l}_2)-cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$M_3=\frac{-cos\mathrm{\Ψ}}{\mathrm{\Δ}}$

$N_1=\frac{cos\mathrm{\Ψ}cos\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}cos\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

$N_2=\frac{cos\mathrm{\Ψ}sin\mathrm{\Φ}(l_1{n}_2-n_1{l}_2)-sin\mathrm{\Ψ}sin\mathrm{\Φ}(m_1{n}_2-n_1{m}_2)}{\mathrm{\Δ}}$

N₃＝0

$\mathrm{\Δ}=\left|\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\\ n_1& n_2& n_3\end{array}\right|$

l₁＝cosΨsinΦ m₁＝sinΨsinΦ n₁＝cosΦ

l₂＝cosΨcosΦ m₂＝sinΨcosΦ n₂＝-sinΦ

l₃＝m₁n₂-n₁m₂ m₃＝l₁n₂-n₁l₂ n₃＝l₁m₂-m₁l₂

S₁(t)、S₂(t)和S₃(t)是采集的三分量地震记录的X，Y和Z分量，C₁(t)、C₂(t)和C₃(t)是校正后的三分量地震记录，t＝1,2,…,N，N为每道地震记录的采样点数，Ψ是步骤15)计算得到的三分量检波器的X分量与X坐标轴的偏离角度，Φ是Z分量与垂直方向的偏离角度，通过对三分量数据进行校正，将水平分量中的纵波能量投影到垂直分量上，将垂直分量中的横波能量投影到水平分量上。

2.根据权利要求1的方法，特点是步骤2)所述的初至时间拾取是在Z分量或水平分量坐标旋转后的X分量地震记录上拾取每一道地震记录的初至时间fbt(i)，i是地震数据道顺序号，i＝1,2,…,Ntrace，Ntrace是每一个分量地震记录的总道数。

3.根据权利要求1的方法，特点是步骤3)所述的三分量地震记录给定初至时窗的协方差矩阵的计算公式为：

$c_{ij}^{\left(0\right)}=\frac{1}{n}\underset{t=t_1}{\overset{t_2}{\Σ}}{S}_i\left(t\right)S_j\left(t\right)$

是协方差矩阵的各个元素，i，j是协方差矩阵的行列序号，i，j＝1,2,3，S_i(t)和S_j(t)均是三分量地震记录，i和j是分量序号，i＝1,2,3，j＝1,2,3，t₁和t₂是计算时窗的起始和终止样点，t₁＝fbt(i)/dt，t₂＝(fbt(i)+twindow)/dt，dt是地震记录采样间隔，twindow是给定的计算时窗长度，n＝t₂-t₁-1，根据计算的协方差矩阵元素得到3×3阶实对称矩阵，也就是三分量地震记录给定时窗的协方差矩阵

${\stackrel{\‾}{C}}^{\left(0\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(0\right)}& c_{12}^{\left(0\right)}& c_{13}^{\left(0\right)}\\ c_{21}^{\left(0\right)}& c_{22}^{\left(0\right)}& c_{23}^{\left(0\right)}\\ c_{31}^{\left(0\right)}& c_{32}^{\left(0\right)}& c_{33}^{\left(0\right)}\end{array}\right].$

4.根据权利要求1的方法，特点是步骤4)所述的sinθ和cosθ的计算公式为：

$\left\{\begin{array}{c}tg2\mathrm{\θ}=\frac{2c_{pq}^{\left(0\right)}}{c_{pp}^{\left(0\right)}-c_{qq}^{\left(0\right)}}\\ \cos 2\mathrm{\θ}=\frac{1}{\sqrt{1+{\mathrm{tg}}^{2}2\mathrm{\θ}}}\\ sin2\mathrm{\θ}=\cos 2\mathrm{\θ}\·tg2\mathrm{\θ}\\ \cos \mathrm{\θ}=\sqrt{\frac{1+\cos 2\mathrm{\θ}}{2}}\\ sin\mathrm{\θ}=\frac{sin2\mathrm{\θ}}{2\cos \mathrm{\θ}}\end{array}\right..$

5.根据权利要求1的方法，特点是步骤5)所述的一次旋转后的协方差矩阵的各个元素的计算公式为：

$\left\{\begin{array}{c}\begin{array}{cc}{c}_{pj}^{\left(1\right)}=c_{pl}^{\left(0\right)}\cos \mathrm{\θ}+c_{pl}^{\left(0\right)}\sin \mathrm{\θ}=c_{jp}^{\left(1\right)}& (l\≠i,j)\end{array}\\ c_{pj}^{\left(1\right)}=-c_{pl}^{\left(01\right)}\sin \mathrm{\θ}+c_{pl}^{\left(0\right)}\cos \mathrm{\θ}=c_{jp}^{\left(1\right)}\\ c_{pp}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\cos }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\sin }^2\mathrm{\θ}+2c_{pl}^{\left(0\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ c_{qq}^{\left(1\right)}=c_{pp}^{\left(0\right)}{\sin }^2\mathrm{\θ}+c_{qq}^{\left(0\right)}{\cos }^2\mathrm{\θ}-2c_{pl}^{\left(01\right)}\cos \mathrm{\θ}\sin \mathrm{\θ}\\ \begin{array}{cc}{c}_{lm}^{\left(1\right)}=c_{ml}^{\left(0\right)}=c_{ml}& (l,m\≠i,j)\end{array}\\ c_{pq}^{\left(1\right)}=c_{qp}^{\left(1\right)}=\frac{1}{2}(c_{qq}^{\left(0\right)}-c_{pp}^{\left(0\right)})\sin 2\mathrm{\θ}+c_{pq}^{\left(0\right)}({\cos }^2\mathrm{\θ}-{\sin }^2\mathrm{\θ})\end{array}\right.$

其中，l,m是协方差矩阵的行列序号，l＝1,2,3,m＝1,2,3；

根据计算的各个矩阵元素得到一次旋转后的协方差矩阵为：

${\stackrel{\‾}{C}}^{\left(1\right)}=\left[\begin{array}{ccc}{c}_{11}^{\left(1\right)}& c_{12}^{\left(1\right)}& c_{13}^{\left(1\right)}\\ c_{21}^{\left(1\right)}& c_{22}^{\left(1\right)}& c_{23}^{\left(1\right)}\\ c_{31}^{\left(1\right)}& c_{32}^{\left(1\right)}& c_{33}^{\left(1\right)}\end{array}\right].$

6.根据权利要求1的方法，特点是步骤6)所述的构造取一次旋转矩阵R⁽¹⁾(p,q,θ)的过程为：利用步骤4)得到的sinθ和cosθ，在3×3阶单位阵I的p行p列，q行q列的交叉位置上置入其它对角线元素为1，其它非对角线元素为0，得到的一次旋转矩阵为：

7.根据权利要求1的方法，特点是步骤7)所述的一次旋转后的协方差矩阵的非对角线元素的平方和的计算公式为：

$E^{\left(1\right)}\left(A\right)=\underset{i\≠j}{\overset{3}{\underset{i,j=1}{\Σ}}}{c}_{ij}^{\left(1\right)}\·c_{ij}^{\left(1\right)}$

其中是步骤5)计算的一次旋转后的协方差矩阵的矩阵元素，i，j是该矩阵的行列序号。

8.根据权利要求1的方法，特点是步骤8)所述的旋转变换迭代计算过程为：判断步骤7)得到的非对角线元素的平方和E⁽¹⁾(A)是否小于给定的阈值δ，如果不小于δ，用步骤5)计算的一次旋转后的协方差矩阵更新步骤3)中的初始协方差矩阵重复步骤4)-7)，如此类推直到非对角线元素的平方和小于给定的阈值δ，迭代计算结束。

9.根据权利要求1的方法，特点是步骤9)所述的协方差矩阵的特征值和特征向量的构造过程为：设经过m次旋转变换计算后，迭代计算终止，则可得到m次旋转后的协方差矩阵以及每次迭代后的旋转矩阵R⁽¹⁾,R⁽²⁾,…,R^(m)，记变换矩阵P^(m)＝R⁽¹⁾R⁽²⁾…R^(m)，则可以得到协方差矩阵的特征值为 是协方差矩阵的各个元素，特征值λ_j对应的特征向量P_j为变换矩阵P^(m)的第j列的各个元素，既

10.根据权利要求1的方法，特点是步骤11)所述的特征向量规则化所用的公式为：

$W_j=\frac{V_j}{\[V_j,V_j\]}=\frac{V_j}{\underset{i=1}{\overset{3}{\Σ}}{v}_{ij}^2},j=1,2,3$

其中，V_j为由步骤10)得到的排序后的特征向量，v_ij为特征向量的元素；得到的规则化后的三个特征向量分别W₁，W₂，W₃，W₁＝(w₁₁w₂₁w₃₁)^T，W₂＝(w₁₂w₂₂w₃₂)^T，W₃＝(w₁₃w₂₃w₃₃)^T，w_ij是规则化后的特征向量的元素。

11.根据权利要求1的方法，特点是步骤12)所述的夹角和夹角α的计算公式为：

α＝arccos(w₃₁)

其中，w₁₁，w₂₁，w₃₁为由步骤11)计算得到的规则化后的主特征向量W₁的元素。

12.根据权利要求1的方法，特点是步骤13)所述的循环计算过程为：设采集的地震记录有N炮，每炮有M个接收点，每个接收点三分量接收；

按照步骤3)-步骤12)的实现过程，计算出第i炮，第j个接收点处的主特征向量在XY平面内的投影与X轴的夹角以及主特征向量与Z轴的夹角α_ij，i＝1,…,N，j＝1,…,M。

13.根据权利要求1的方法，特点是步骤14)的第i炮，第j个接收点处的方位角的计算公式为：

其中，(x_i,y_i)为第i炮的坐标，(x_j,y_j)为第j个接收点处理的坐标，i＝1,…,N，j＝1,…,M。