本发明涉及油气勘探开发领域，具体包含一种针对强各向异性vti介质的各向异性参数预测方法，以提高油气储层勘探精度。
背景技术：
：：1、随着油气勘探开发不断向深层、深水、非常规等领域推进，目标油气藏的各向异性强度也逐渐增强，特别是可等效为垂直横向各向同性（vti）介质的页岩等层状地层一般表现出较强的各向异性程度。叠前地震反演方法是各向异性储层参数预测的有效手段。2、目前针对各向异性vti介质参数预测的地震方法主要包括两类，一类是基于扰动理论的叠前地震反演方法，另一类是基于精确各向异性zoeppritz方程的叠前地震非线性反演方法。第一类主要利用介质扰动理论，将各向异性介质视为各向同性背景与各向异性扰动的叠加，通过舍掉与各向异性参数相关的高阶项，保留与各向异性参数相关的低阶项来建立线性或非线性近似反射/透射系数方程，基于所建立的方程，采用合理的优化算法实现vti介质各向异性参数预测，但该类方法舍掉了部分与各向异性参数相关的项，损失了各向异性参数的预测精度，尤其当反射系数方程仅保留各向异性参数一阶项，表现为线性形式时，该方法针对强各向异性介质参数预测将表现出较低的精度；第二类主要基于精确各向异性zoeppritz方程，采用合理的非线性求解算法，实现vti介质各向异性参数地震预测，该方法有效保留各向异性参数高阶信息，在一定程度上提高了参数预测精度，但参数间极强的耦合作用和反射系数方程极强的非线性程度，降低了求解算法的运行效率，增加了反演耗时，加剧了预测多解性。3、申请号为201710280917.5的发明专利公开了一种地震反演方法和装置，以达到提高vti参数的反演精度，指导油田高效开发的目的。虽然该专利在一定程度上提高了强各向异性vti介质的地震描述精度，但在推导纵波和横波垂直慢度近似表达式时，依旧采用弱各向异性假设，所以该专利在理论上不适用于具有任意各向异性强度的vti地层地震预测。4、申请号为202110533801.4的发明专利公开了基于vti介质的页岩储层avo反演表征方法及系统。该专利基于弱各向异性假设，分别推导了包含6个未知参数的饱和vti介质地震反射系数方程和包含7个未知参数的衰减vti介质地震反射系数方程，多参数同步反演会降低vti介质参数预测稳定性，所以该专利通过参数组合，将这两个反射系数方程均改写为包含4个参数的表达式，然后利用四参数方程开展了vti介质参数反演。但该专利没有考虑页岩等储层具有的强各向异性特征。5、申请号为202010053042.7的发明专利公开了一种vti介质的各向异性参数分步反演方法。申请号为202110218907.5的发明专利公开了一种基于反射系数精确式的页岩vti储层叠前混合反演方法。这两个专利所提技术相似，均采用vti介质精确反射系数方程开展储层地震预测，区别在于采用的求解算法不一样。所以适用于任意强各向异性vti介质，但vti介质精确反射系数方程的形式非常复杂，求解难度大、运行效率低，严重制约了该方法的工业化应用。6、总之，现有的各向异性储层参数叠前地震反演方法主要基于弱各向异性假设，降低了页岩等强各向异性vti介质参数的预测精度，所以亟需研究考虑强各向异性特征的vti介质参数预测方法，以提高页岩等层状地层的解释精度。技术实现思路1、为了解决上述技术问题，本发明提出一种更便于在实际工程中应用、理论精度更高的考虑vti介质强各向异性特征的储层参数地震预测方法，该方法在考虑vti介质强各向异性特征的基础上平衡各向异性参数求解速度和反演精度。本发明并通过案例分析验证了新方法的可行性。2、为解决上述技术问题，本发明采用的技术方案是：一种考虑强各向异性特征的vti介质参数预测方法，该方法包括以下步骤：3、步骤1，构建各向异性参数、纵波速度、横波速度和质量密度直接表达的强各向异性vti介质反射系数方程；4、步骤2，基于强各向异性vti介质反射系数方程和贝叶斯反演理论，构建待反演模型参数服从的后验概率密度函数；5、步骤3，构建建议分布函数、接受概率函数；6、步骤4，基于随机采样策略直接反演各向异性参数。7、进一步的，步骤1中，对于分离两个单斜各向异性介质的水平反射面，当界面两侧介质参数差异较小时，下行入射p波在该界面处的pp反射系数近似为：8、 （1）9、（1）式中，下角标d代表下行波，下角标p代表p波，上标t代表矩阵转置操作；代表强各向异性介质pp波反射系数；代表下行p波的标准化极化矢量；代表下行p波慢度矢量的第k分量；代表地层介质的四阶刚度张量，其元素可写为；符号代表差分算子，用于求取反射界面两侧介质属性的差异，且10、 （2a）11、 （2b）12、<mrow><mstyledisplaystyle="true"mathcolor="#000000"/><mstyledisplaystyle="true"mathcolor="#000000"><mi>=</mi></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mo>[</mo><mtable><mtr><mtd><mi>−</mi><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mi>−</mi><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd></mtr></mtable><mo>]</mo></mstyle></mrow> （2c）13、对于vti介质，p波垂直慢度为：14、 （3）15、（3）式中，16、17、18、上式中，代表水平慢度，代表介质刚度系数的voigt表达形式，代表质量密度；19、vti介质的p波极化矢量表达为：20、<mstyledisplaystyle="true"mathcolor="#000000"><msub><mi>v</mi><mi>dp</mi></msub><mi>=</mi><msup><mrow><mo>[</mo><mtable><mtr><mtd><msub><mi>ε</mi><mi>p</mi></msub><msubsup><mi>p</mi><mi>x</mi><mrow/></msubsup></mtd><mtd><msub><mi>ε</mi><mi>p</mi></msub><msubsup><mi>p</mi><mi>y</mi><mrow/></msubsup></mtd><mtd><msub><mi>ε</mi><mi>p</mi></msub><msub><mi>q</mi><mi>p</mi></msub></mtd></mtr></mtable><mo>]</mo></mrow><mi>t</mi></msup></mstyle> （4）21、（4）式中，和分别代表水平慢度的x和y分量，，代表极化矢量的标准化系数， ,22、vti介质下行p波相速度表达为：23、<mstyledisplaystyle="true"mathcolor="#000000"><msub><mi>v</mi><mi>p</mi></msub><mi>=</mi><msqrt><mfrac><mn>1</mn><mrow><mn>2</mn><mi>ρ</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msubsup><mi>c</mi><mn>11</mn><mrow/></msubsup><mo>+</mo><msubsup><mi>c</mi><mn>44</mn><mrow/></msubsup></mrow><mo>)</mo></mrow><msup><mi>sin</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><mrow><mo>(</mo><mrow><msubsup><mi>c</mi><mn>33</mn><mrow/></msubsup><mo>+</mo><msubsup><mi>c</mi><mn>44</mn><mrow/></msubsup></mrow><mo>)</mo></mrow><msup><mi>cos</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><msqrt><mi>d</mi></msqrt></mrow><mo>]</mo></mrow></msqrt></mstyle> （5）24、（5）式中，代表p波相速度；代表p波入射角，且25、<mstyledisplaystyle="true"mathcolor="#000000"><mi>d</mi><mi>=</mi><msup><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><msub><mi>c</mi><mn>11</mn></msub><mi>−</mi><msub><mi>c</mi><mn>44</mn></msub></mrow><mo>)</mo></mrow><msup><mi>sin</mi><mn>2</mn></msup><mi>θ</mi><mi>−</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mn>33</mn></msub><mi>−</mi><msub><mi>c</mi><mn>44</mn></msub></mrow><mo>)</mo></mrow><msup><mi>cos</mi><mn>2</mn></msup><mi>θ</mi></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>+</mo><mn>4</mn><msup><mrow><mo>(</mo><mrow><msubsup><mi>c</mi><mn>13</mn><mrow/></msubsup><mo>+</mo><msubsup><mi>c</mi><mn>44</mn><mrow/></msubsup></mrow><mo>)</mo></mrow><mn>2</mn></msup><msup><mi>sin</mi><mn>2</mn></msup><mi>θ</mi><msup><mi>cos</mi><mn>2</mn></msup><mi>θ</mi></mstyle>;26、将公式（2）~公式（5）带入公式（1），并引入thomsen各向异性参数、、及、、，得到界面两侧介质属性弱差异近似下的所述强各向异性vti介质反射系数方程：27、 （6）28、（6）式中，<mrow><mstyledisplaystyle="true"mathcolor="#000000"><mn>1</mn><mn>2</mn></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>g</mi><mo>+</mo><mi>fu</mi></mrow><mo>)</mo></mrow><msup><mi>tan</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><msup><mi>u</mi><mn>2</mn></msup></mrow><mo>]</mo></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><msup><mrow/><mrow><mi>−</mi><mn>1</mn></mrow></msup></mstyle></mrow>29、<mstyledisplaystyle="true"mathcolor="#000000"><msub><mi>h</mi><mi>vp</mi></msub><mi>=</mi><mi>g</mi><mi>(</mi><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi><msup><mi>)tan</mi><mn>4</mn></msup><mi>θ</mi><mo>+</mo><mrow><mo>[</mo><mrow><mn>2</mn><mrow><mo>(</mo><mrow><mi>ε</mi><mi>−</mi><mi>δ</mi></mrow><mo>)</mo></mrow><mi>−</mi><mi>s</mi><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><msup><mi>tan</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><mi>gs</mi><mrow><mo>(</mo><mrow><mn>1</mn><mi>−</mi><msubsup><mi>k</mi><mn>2</mn><mi>r</mi></msubsup><mi>s</mi></mrow><mo>)</mo></mrow></mstyle>30、31、<mstyledisplaystyle="true"mathcolor="#000000"><msub><mi>h</mi><mrow><mi>vs</mi><mn>0</mn></mrow></msub><mi>=</mi><mrow><mo>[</mo><mrow><mn>2</mn><mrow><mo>(</mo><mrow><mi>δ</mi><mi>−</mi><mi>ε</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>s</mi><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><msup><mi>tan</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><mi>s</mi><mrow><mo>(</mo><mrow><mn>1</mn><mi>−</mi><mi>s</mi></mrow><mo>)</mo></mrow></mstyle>32、33、34、<mstyledisplaystyle="true"mathcolor="#000000"><mi>=</mi><mfrac><mn>1</mn><msup><mi>f</mi><mn>2</mn></msup></mfrac><mrow><mo>{</mo><mtablecolumnalign="left"><mtr><mtd><mrow><mo>[</mo><mrow><mn>2</mn><mrow><mo>(</mo><mrow><mi>ε</mi><mi>−</mi><mi>δ</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>g</mi><mrow><mo>(</mo><mrow><mn>2</mn><mi>δ</mi><mo>+</mo><mn>2</mn><mi>f</mi><mo>+</mo><mn>2</mn><mi>−</mi><mi>g</mi></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi></mrow><mo>)</mo></mrow><msup><mi>tan</mi><mn>4</mn></msup><mi>θ</mi></mtd></mtr><mtr><mtd><mo>+</mo><mrow><mo>[</mo><mtablecolumnalign="left"><mtr><mtd><mn>4</mn><mi>ε</mi><mrow><mo>(</mo><mrow><mi>g</mi><mi>−</mi><mi>s</mi></mrow><mo>)</mo></mrow><mo>+</mo><mn>2</mn><mi>δ</mi><mrow><mo>(</mo><mrow><mn>2</mn><mi>g</mi><mi>−</mi><mi>s</mi></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mi>g</mi><mi>−</mi><mn>1</mn></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mn>2</mn><mi>gf</mi><mrow><mo>(</mo><mrow><mi>g</mi><mi>−</mi><mi>s</mi></mrow><mo>)</mo></mrow><mo>+</mo><mi>g</mi><mrow><mo>(</mo><mrow><mi>g</mi><mi>−</mi><mn>2</mn></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mi>s</mi><mi>−</mi><mn>2</mn><mi>g</mi></mrow><mo>)</mo></mrow><mi>−</mi><mi>s</mi></mtd></mtr></mtable><mo>]</mo></mrow><msup><mi>tan</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><msup><mrow><mo>(</mo><mrow><mi>g</mi><mi>−</mi><mi>s</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup></mtd></mtr></mtable><mo>}</mo></mrow></mstyle>35、，36、，37、，38、<mstyledisplaystyle="true"mathcolor="#000000"><mi>s</mi><mi>=</mi><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi><mo>+</mo><mi>g</mi></mrow><mo>)</mo></mrow><msup><mi>tan</mi><mn>2</mn></msup><mi>θ</mi><mo>+</mo><mi>(</mi><mn>1</mn><mo>+</mo><mi>g</mi><mi>)</mi><mo>+</mo><mi>ω</mi><msup><mi>sec</mi><mn>2</mn></msup><mi>θ</mi></mrow><mo>]</mo></mrow></mstyle>，39、<mstyledisplaystyle="true"mathcolor="#000000"><mi>ω</mi><mi>=</mi><msqrt><msup><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>2</mn><mi>ε</mi><mi>−</mi><mi>g</mi></mrow><mo>)</mo></mrow><msup><mi>sin</mi><mn>2</mn></msup><mi>θ</mi><mi>−</mi><mi>(</mi><mn>1</mn><mi>−</mi><mi>g</mi><msup><mi>)cos</mi><mn>2</mn></msup><mi>θ</mi></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>+</mo><mn>4</mn><msup><mi>f</mi><mn>2</mn></msup><msup><mi>sin</mi><mn>2</mn></msup><mi>θ</mi><msup><mi>cos</mi><mn>2</mn></msup><mi>θ</mi></msqrt></mstyle>，40、41、42、。43、进一步的，步骤2中，根据所述强各向异性vti介质反射系数方程，得到地震正演方程：44、 （7）45、（7）式中，d代表地面观测地震数据，m代表待反演模型参数，rpp(m)为各向异性vti介质反射系数，w指代子波矩阵，n为观测数据中包含的干扰噪音向量；；46、将公式（7）变换为频域正演方程：47、 （8）48、（8）式中，顶端带有弧线的字符代表频率域变量或算子，、、和分别代表频域地震数据、频域地震子波矩阵、反射系数的傅里叶变换矩阵、频域噪音，且。49、采用贝叶斯理论构建强各向异性vti介质待反演模型参数m服从的后验概率密度分布函数：50、 （9）51、（9）式中，为模型参数m服从的后验概率密度分布函数，为似然函数，p(m)代表待反演模型参数的先验概率密度分布，为地面观测地震数据的概率密度函数；52、选用高斯分布作为待反演参数的先验概率分布，待反演模型参数的先验概率分布为：53、<mrow><mstyledisplaystyle="true"mathcolor="#000000"/><mstyledisplaystyle="true"mathcolor="#000000"><mi>=</mi></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mfrac><mn>1</mn><msup><mrow><mo>(</mo><mrow><mn>2</mn><mi>π</mi></mrow><mo>)</mo></mrow><mrowbevelled="true"><mi>n</mi><mostretchy="true">/</mo><mn>2</mn></mrow></msup></mfrac></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mfrac><mn>1</mn><msup><mrow><mo>|</mo><msubsup><mi>c</mi><mi>m</mi><mrow/></msubsup><mo>|</mo></mrow><mrowbevelled="true"><mn>1</mn><mostretchy="true">/</mo><mn>2</mn></mrow></msup></mfrac></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mi>e</mi></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mi>x</mi></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mi>p</mi></mstyle><mstyledisplaystyle="true"mathcolor="#000000"><mo>[</mo><mrow><mi>−</mi><mfrac><mn>1</mn><mn>2</mn></mfrac><msup><mrow><mo>(</mo><mrow><mi>m</mi><mi>−</mi><msubsup><mi>m</mi><mn>0</mn><mrow/></msubsup></mrow><mo>)</mo></mrow><mi>t</mi></msup><msubsup><mi>c</mi><mi>m</mi><mrow><mi>−</mi><mn>1</mn></mrow></msubsup><mrow><mo>(</mo><mrow><mi>m</mi><mi>−</mi><msubsup><mi>m</mi><mn>0</mn><mrow/></msubsup></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mstyle></mrow> （10）54、（10）式中，n代表采样点个数；cm代表模型参数的协方差矩阵，利用测井数据进行估计；m0代表模型参数的低频约束；55、假设时间域和频率域的似然函数服从高斯分布，且干扰噪音具有随机性和相互独立性，则时频联合域反演的似然函数为：56、 （11）57、（11）式中，和分别代表随机干扰噪音和n的协方差矩阵；58、将公式（10）和公式（11）代入公式（9）中，得到待反演模型参数m服从的后验概率密度分布函数：59、 （12）60、（12）式中，；；。61、进一步的，步骤3中, 选用高斯分布函数作为建议分布函数：62、（13）63、（13）式中， 代表候选模型参数相对当前模型参数第k次迭代的协方差矩阵，通过测井数据等先验认识估计得到；64、接受概率函数表达为：65、（14）66、对公式（14）求对数得：67、（15）68、（15）式中：69、 。70、进一步的，步骤4包括以下步骤：71、步骤41利用测井数据获得纵波速度、横波速度、质量密度、各向异性参数ε、各向异性参数δ的测井解释结果；72、步骤42利用井点数据进行内插外推构建纵波速度、横波速度、质量密度、各向异性参数ε、各向异性参数δ的三维低频模型，将该低频模型作为模型参数的当前状态m；73、步骤43，设置循环次数最大为y，当前值为y=0；当前模拟的模型参数序号为x=0（x为整数）；实现次数最大为u，当前值为u=0；迭代次数最大为s，当前值为s=0；74、步骤44，判断当前循环次数y<y？如果否，结束反演，模型参数当前状态m为最终预测结果并输出；如果是，则y=y+1，将模型参数当前状态设置为本次循环的初始状态，m0=m，转至步骤45；75、步骤45，判断当前模拟的模型参数序号x<5？如果否，x=0，返回至步骤44；如果是，则x=x+1，将第x个模型参数视为待模拟参数，并固定其他模型参数，转至步骤46；76、步骤46，判断当前实现次数u<u？如果否，将u次实现结果的平均值作为第x个模型参数的当前值m(u)，u=0，返回至步骤45；如果是，则u=u+1，转至步骤47；77、步骤47，判断当前迭代次数s<s？如果否，统计分析本次实现的收敛序列，将收敛序列的平均值作为本次实现结果，s=0，返回至步骤46；如果是，转至步骤48；78、步骤48，根据建议分布产生第x个模型参数的候选状态，且建议分布为：79、将候选状态和第s次迭代状态带入公式（15）计算接受概率的对数，从均匀分布uniform(0, 1)中随机抽取变量j：如果ln(j)≤lnα，则将候选状态记为第x个模型参数的新状态，s=s+1，返回至步骤47；否则，,s=s+1，返回至步骤47。80、与现有技术相比，本发明具有以下有益效果：81、本发明构建的vti介质反射系数方程仅假设反射界面两侧属性差异较小、理论上适用于具有任意vti各向异性强度的情况，所以在强各向异性地层描述中，其精度高于传统弱各向异性理论；82、本发明构建的vti介质反射系数方程的复杂度远低于vti介质精确反射系数方程，所以相较于基于精确反射系数方程的vti介质参数反演方法，本发明求解难度小、运行效率高，为强各向异性介质地震预测技术的工业化应用提供了很好的支撑；83、本案所构建的反演流程相较于其他非线性反演流程，形式更简单、操作更方便，且可实现各向异性参数的直接反演。当前第1页12当前第1页12