CN105467451B - 基于全变差最小化约束的地震反射系数反演方法 - Google Patents

基于全变差最小化约束的地震反射系数反演方法 Download PDF

Info

Publication number
CN105467451B
CN105467451B CN201610020556.6A CN201610020556A CN105467451B CN 105467451 B CN105467451 B CN 105467451B CN 201610020556 A CN201610020556 A CN 201610020556A CN 105467451 B CN105467451 B CN 105467451B
Authority
CN
China
Prior art keywords
mrow
msub
mtd
mtr
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610020556.6A
Other languages
English (en)
Other versions
CN105467451A (zh
Inventor
吴秋波
张洞君
邹文
黄东山
王颀
刘开元
周晶晶
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Petroleum Corp
BGP Inc
Original Assignee
BGP Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGP Inc filed Critical BGP Inc
Priority to CN201610020556.6A priority Critical patent/CN105467451B/zh
Publication of CN105467451A publication Critical patent/CN105467451A/zh
Application granted granted Critical
Publication of CN105467451B publication Critical patent/CN105467451B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明提供了一种基于全变差最小化约束的地震反射系数反演方法。对叠后地震数据采用逐道逐时窗反演。对于单个时窗，首先对时窗内的叠后地震数据和事先提取的地震子波进行傅立叶变换，然后得到反射系数的频域表达式，再对时域反射系数进行傅里叶变换，在反射系数奇偶分解的基础上提取其实部和虚部，构建相应的求解方程，在传统共轭梯度算法基础上采取最小全变差约束进行方程求解，得到奇、偶反射系数并重构得到时窗内的原始反射系数；所有时窗依次进行反演，得到单道地震记录的反射系数。对所有地震道依次进行反演，得到每道地震记录的反射系数。本发明的方法可以有效地反演地震反射系数，利于提高地震数据分辨率和提高储层预测精度。

Description

基于全变差最小化约束的地震反射系数反演方法
技术领域
本发明属于地震油气勘探技术领域，具体为一种基于全变差最小化约束共轭梯度算法的地震反射系数反演方法。
背景技术
作为地震油气勘探中的重要步骤和环节，地震储层预测一直在油气勘探中扮演着不可或缺的角色。随着岩性油气藏勘探开发的不断深入，对小于地震波长四分之一的薄储层(薄层)预测问题日益突出。模型反演，稀疏脉冲反演等常规地震反演方法对薄储层识别能力有限，直接影响到后续的油气藏评价和井位设计。
地震反射系数反演是在频率域实施的，可获得高分辨率的时间域反射系数信息，实现对小于调谐厚度的薄储层进行有效识别，以达到提高地震数据分辨率和提高储层预测精度的目的。
地震反射系数反演是用有限的频域信息得到全域频域信息，是一个欠定反演问题，传统共轭梯度方法求解欠定问题效率低且易受异常值影响，会导致反演结果背离实际地层特性，与地震剖面匹配较差等不利后果，不能正确的反应出真实的地层信息。
目前应用于地震反射系数估计的反演算法包括：Charles I.Puryear和JohnP.Castagna采用的最小二乘共轭梯度法。其基本思想是把共轭性与最小二乘方法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜素，求出目标函数的极小点。
Rui Zhang和John P.Castagna采用的基追踪算法。它采用反射系数的范数作为稀疏性的度量，通过最小化L1范数将反射系数稀疏表示问题定义为一类有约束的极值问题，进而转化为线性规划问题进行求解。
袁三一和王尚旭等采用的粒子群算法和列文伯格-马夸尔特法的联合算法。其首先采用粒子群算法进行全局寻优，来反演出地层反射系数位置，然后采用列文伯格-马夸尔特法算法来精准地求取反射系数的数值大小。粒子群算法和列文伯格-马夸尔特法算法相联合，反演结果更精确，收敛性更快。
秦德文等采用的基于随机搜索的蒙特卡罗算法。它是以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。蒙特卡罗方法有很强的适应性，该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度。
Kelyn Paola和Germán Ojeda采用的模拟退火算法和遗传算法。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解。而遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。应用于反射系数估计问题时，遗传算法效果比模拟退火算法稍好。
柴新涛，李振春等采用的最小二乘QR分解算法。它是利用Lanezos方法求解最小二乘问题的一种投影法。由于在求解过程中用到QR因子分解，在对数据误差传递的压制和求解收敛效率上具有明显的优越性。总体来讲，这些算法都可以不同程度的解决反射系数的反演问题。
发明内容
针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的一个或多个问题。例如，本发明的目的之一在于解决现有技术中存在的反射系数反演效率低和易受异常值影响的问题，提供一种基于全变差最小化约束的共轭梯度地震反射系数反演方法。该方法利用地层连续性约束，在频率域提升共轭梯度反演的准确性和效率。通过全变差最小化约束增加对反射系数的识别能力，抑制共轭梯度方法易受到数据变化扰动以及反演结果地层跳变的问题，获得高分辨率的时间域反射系数，提升地震数据的分辨率，提高对小于调谐厚度的薄储层反射系数识别能力。
为了实现上述目的，本发明提供了一种基于全变差最小化约束的地震反射系数反演方法。所述反演方法包括以下步骤：
A、从叠后地震数据S中提取地震子波w并进行傅立叶变换得到地震子波w的频域表示W(f)，其中：
W(f)＝FFT(w) (1)
B、从所述叠后地震数据S中取一个地震道的一个时窗内的地震数据s进行傅立叶变换，得到地震数据s的频域表示S(f)，其中：
S(f)＝FFT(s) (2)
C、根据地震数据形成原理由所述地震数据s的频域表示S(f)和所述地震子波w的频域表示W(f)得到反射系数的频域表示R(f)，其中：
R(f)＝S(f)/W(f) (3)
D、以所述时窗中心为分析点，对所述时窗内的地震数据s的时域反射系数进行傅立叶变换，得到反射系数的频域表达R(f)＇，其中：
在等式(4)中，N为由点数表示的分析数据长度，Ti代表反射系数对之间的时间间隔，j为虚数单位，f为频率。
E、在所述步骤C得到的反射系数的频域表示R(f)和所述步骤D得到的反射系数的频域表达R(f)＇相等的基础上，根据反射系数奇偶分解原理构建所述时窗内反射系数反演的目标函数方程。
F、利用共轭梯度算法求解所述目标函数方程以得到所述时窗内反射系数的奇分量和偶分量，并重构出所述时窗内的反射系数，其中，在每次迭代过程中对共轭梯度算法所求得的解进行全变差最小化约束，同时将经全变差最小约束之后的解作为下一次迭代计算的初始解。
G、在所述地震道上按步长step滑动所述时窗，得到下一个时窗的地震数据，重复步骤B至F，直到所述时窗遍历所述地震道的地震数据，完成所述地震道的反演，得到所述地震道的反射系数。
H、取下一地震道的地震数据并按照步骤B到G进行处理，直到得到每一个地震道的反射系数。
根据本发明的基于全变差最小化约束的地震反射系数反演方法的一个实施例，所述步骤E包括：用欧拉公式将所述步骤D得到反射系数的频域表达R(f)＇中的指数项展开，得到等式(5)：
根据奇偶分解原理，等式(5)中实部是反射系数对ri与rN-i+1的偶分量，用re(i,N-i+1)来表示，虚部是反射系数对ri与rN-i+1的奇分量，用ro(i,N-i+1)来表示，得到等式(6)：
将所述等式(3)分析点移至时窗中点，并结合等式(6)建立目标函数方程，得到等式(7)：
在等式(7)中，flow为低频截止频率，fhigh为高频截止频率，αe为偶分量权重，αo为奇分量权重，Re代表着实部，Im代表着虚部，re为反射系数的偶分量，ro为反射系数的奇分量，Δt为相对于分析点的位移量，tw为时窗win的半窗长，T代表反射系数对之间的时间间隔，t为分析点的时间。
将等式(7)在频域和时域离散化并写为奇分量和偶分量的形式，得到下面的等式(8)和等式(9)：
在等式(8)和等式(9)中，f1是起始频率，取值为所述flow；fm是截止频率，取值为所述fhigh。
将所述等式(8)和等式(9)分别写为矩阵形式，得到下面的等式(10)和等式(11)：
be＝Ae×re (10)
bo＝Ao×ro (11)
建立目标函数方程：
在等式(10)、(11)和(12)中，αe为偶分量权重，αo为奇分量权重；be是频域反射系数的实部；bo是频域反射系数的虚部；Ae为偶分量变换矩阵，即：
Ao为奇分量变换矩阵，即：
将所述目标函数方程(12)写为更一般的优化问题，得到目标函数方程(15)：
b＝Ax(15)
在等式(15)中，x为待求的反射系数。
根据本发明的基于全变差最小化约束的地震反射系数反演方法的一个实施例，所述步骤F包括：
将所述目标函数方程(15)的求解问题转化为||Ax＝b||2优化问题，进一步写成二次函数形式：
在等式(16)中，Q＝ATA，y＝ATb，T表示转置运算。
初始化：x0＝0，r0＝b-Ax0，c＝||Q||2，k＝1，k为迭代次数。
循环步骤：
(1)更新梯度方向：
gk＝Qxk-1-y (17)
(2)更新共轭方向：
若k＝1，则pk＝r0；若k≠1，则
pk＝-gk+αkpk-1 (19)
(3)更新第k次迭代的解：
(4)对xk进行全变差最小化，
其中，D表示整体剖面集，采用对称边界条件，xir,jc表示第ir行第jc个元素。
(5)再次更新残差：
rk＝b-Axk (22)
在等式(22)中，rk为第k次迭代的残差。
停止条件：当残差小于预定阈值时收敛，停止迭代，输出xk；否则进入第k+1次迭代。
与现有技术相比，本发明的有益效果包括：本发明采用的基于全变差最小化约束的共轭梯度地震反射系数反演方法利用了地层的横向连续性特征，克服经典的共轭梯度算法求解时效性差，易受异常值影响的缺点，提升反演的精度与效率。
附图说明
通过下面结合附图进行的描述，本发明的上述和其他目的和特点将会变得更加清楚，其中：
图1是地震信号反射系数奇偶分解示意图。
图2是模型地震剖面。
图3是传统共轭梯度方法反演得到的反射系数。
图4是利用本发明示例性实施例的基于全变差最小化约束的地震反射系数反演方法对图2进行反演处理的反射系数剖面。
图5是原始实际地震剖面。
图6是利用本发明示例性实施例的基于全变差最小化约束的地震反射系数反演方法对图5进行反演处理后的结果剖面。
其中，图2至图6的横坐标为道，纵坐标为时间，单位是毫秒。
具体实施方式
在下文中，将结合附图和示例性实施例详细地描述根据本发明的基于全变差最小化约束的地震反射系数反演方法。需要说明的是，在本申请中，地震道简称为道。
本发明研究的是求解欠定方程的反演算法。这里采用了全变差最小化约束的正则化方法，对方程解的性质做一些先验约束，提高解的可靠性同时提高运行效率。
本发明最重要的创新之处在于采用了全变差最小化约束以利用地层横向连续性，并将其作为下一次迭代过程的初始解，如此达到求解稳定高效的目的。
本发明的反射系数反演方法的技术思路是：对叠后地震数据采用逐道逐时窗反演。对于单个时窗，首先对时窗内的叠后地震数据和事先提取的地震子波进行傅立叶变换，然后由地震数据形成原理得到反射系数的频域表达式，再对时域反射系数进行傅里叶变换，在反射系数奇偶分解的基础上提取其实部和虚部，构建相应的求解方程，在传统共轭梯度算法基础上采取全变差最小化约束进行方程求解，得到奇、偶反射系数，由奇、偶反射系数重构得到时窗内的原始反射系数，其中，奇反射系数又称为反射系数的奇部或奇分量，偶反射系数又称为反射系数的偶部或偶分量。所有时窗依次进行反演，便得到单道地震记录的反射系数。对所有地震道依次进行反演，最后便得到每道地震记录的反射系数。
根据本发明示例性实施例的基于全变差最小化约束的地震反射系数反演方法包括以下步骤：
步骤一，数据输入：读入叠后地震剖面数据S(又称为叠后地震数据，简称地震数据，地震记录)。
步骤二，初始参数设置：根据地震数据的时频分析结果设置优势频段[fhigh，flow]，其中，flow为低频截止频率，fhigh为高频截止频率；设置奇偶权重αe，αo，设置反演时窗长度win和步长step。
步骤三，反演步骤：
A，从叠后地震数据S中提取地震子波w，并对地震子波w进行傅立叶变换得到地震子波w的频域表示W(f)，其中，
W(f)＝FFT(w) (1)
B，从叠后地震数据S中取一个地震道的一个时窗内的地震数据s进行傅立叶变换，得到地震数据s的频域表示S(f)，其中，
S(f)＝FFT(s) (2)
C，根据地震数据形成原理，由地震数据s的频域表示S(f)和地震子波w的频域表示W(f)得到反射系数的频域表示R(f)：
R(f)＝S(f)/W(f) (3)
D，将时窗的中点设为分析点，则分析点上下每一对采样点的反射系数模型进行奇偶分解，图1是地震信号反射系数奇偶分解示意图，具体由时域反射系数傅立叶变换来实现，即下面的等式(4)：
在等式(4)中，N为由点数表示的分析数据长度，Ti为反射系数对之间的时间间隔，j为虚数单位。用欧拉公式将上式中的指数项展开，得等式(5)：
E、在所述步骤C得到的反射系数的频域表示R(f)和所述步骤D得到的反射系数的频域表达R(f)＇相等的基础上，根据反射系数奇偶分解原理构建所述时窗内反射系数反演的目标函数方程。
具体地，根据奇偶分解原理，等式(5)中实部即是反射系数对ri与rN-i+1的偶分量，用re(i,N-i+1)来表示，虚部即是反射系数对ri与rN-i+1的奇分量，用ro(i,N-i+1)来表示，得到等式(6)：
将等式(3)分析点移至时窗中点，由上所述，实部代表着偶分量，虚部代表着奇分量，并结合等式(6)建立目标函数，即等式(7)：
在等式(7)中，flow为低频截止频率，fhigh为高频截止频率，αe为偶分量权重，αo为奇分量权重，Re代表着实部，Im代表着虚部，re为反射系数的偶分量，ro为反射系数的奇分量，Δt为相对于分析点的位移量，tw为时窗win的半窗长，T代表反射系数对之间的时间间隔，t为分析点的时间。
将等式(7)在频域和时域离散化，并写为奇、偶分量的形式，如下面的等式(8)和等式(9)：
等式(8)和等式(9)中的f1是起始频率，取值为所述flow；fm是截止频率，取值为所述fhigh。将所述等式(8)和等式(9)简写为矩阵形式(10)和(11)，并建立目标函数方程，即等式(12)：
be＝Ae×re (10)
bo＝Ao×ro (11)
在等式(10)、(11)和(12)中，αe为偶分量权重，αo为奇分量权重；be是频域反射系数的实部；bo是反射系数的虚部。值得指出的是，这里的be和bo是由测量的地震数据和子波的傅立叶变换相除得到的。
Ae为偶部变换矩阵，即：
Ao为奇部变换矩阵，即：
将等式(12)写为更一般的优化问题，即等式(15)：
b＝Ax(15)
在等式(15)中，x为待求的反射系数。反射系数反演中的核心问题最终归结为求解线性方程组的形式。
F，由于求解之前并不知道反射系数的实际位置，通常先假定所有的时间点上都有反射系数。在实际的地层中，同一地层面的变化是连续的，所以邻近道的反射系数连续性应该较好，因此等式(15)中所求的解在不同道之间的偏差是较小的。同时在频域进行了截频处理，使得反射系数反演问题成为一个欠定问题。直接使用传统的共轭梯度方法求解等式(15)时，就不能够考虑到这方面的要求，因此本文提出在传统共轭梯度算法基础上采取全变差最小化约束，以达到解稳定性的目的。算法的基本思想是，在每次迭代过程中对共轭梯度法所求得的解进行全变差最小化约束，同时将其作为下一次迭代过程的初始解。算法流程如下：
将目标函数方程求解问题转化为||Ax＝b||2优化问题，进一步写成二次函数形式：
其中，Q＝ATA，y＝ATb，这里的T表示转置运算。
初始化：x0＝0，r0＝b-Ax0，c＝||Q||2，k＝1，k为迭代计数器，即k为迭代次数，开始迭代，直至收敛。
循环步骤：
(1)更新梯度方向：
gk＝Qxk-1-y (17)
(2)更新共轭方向：
若k＝1，则pk＝r0；若k≠1，则
pk＝-gk+αkpk-1 (19)
(3)更新第k次迭代的解：
(4)对xk进行全变差最小化，
其中，D表示整体剖面集，此处采用对称边界条件，xir,jc表示第ir行第jc个元素。
(5)再次更新残差：
rk＝b-Axk (22)
在等式(22)中，rk为第k次迭代的残差。
停止条件：当残差||r||2小于给定阈值时收敛，停止迭代，输出xk，否则进入第k+1次迭代。
经过上述方法得到的解x只是反射系数的奇分量ro(t)和偶分量re(t)，按照式(23)重构出时窗内的反射系数：
r(t)＝ro(t)+re(t) (23)
G，按照步长step滑动时窗，得到下一个时窗的地震数据，重复步骤B到F，直到时窗滑到该道地震数据结束(即时窗遍历该地震道的地震数据)，完成一道地震数据的反演，得到该道地震数据的反射系数。
H，取下一道地震数据按照步骤B到G进行处理，得到每一道的反射系数。
进一步地，为验证上述技术的效果，发明人采用了一个模型地震剖面进行了测试。图2是模型地震剖面，图3是由传统共轭梯度算法反演处理得到的反射系数剖面，图4是利用本发明示例性实施例的基于全变差最小化约束的反射系数反演方法对图2进行反演处理后的反射系数剖面，对比分析图3和图4可知，本发明反演处理的结果较好地解决了层位跳变的问题。更进一步，发明人选取了四川盆地某地区采集到的地震数据进行了应用测试，图5是原始地震剖面，图6是利用本发明示例性实施例的基于全变差最小化约束的反射系数反演方法对图5进行反演处理后的反射系数剖面。从图5与图6比较分析知道，利用本发明处理后的剖面分辨率明显提高。
综上所述，本发明提供了一种基于全变差最小化约束的共轭梯度反射系数反演方法，通过引入全变差最小化约束，利用了地层连续性约束，抑制共轭梯度方法易受数据变化扰动以及反演结果地层跳变问题，在频率域提升共轭梯度反演的准确性和效率，增加了对反射系数的识别能力，能够实现对小于调谐厚度的薄储层(例如，小于地震波长四分之一的薄储层)进行有效识别，利于提高地震数据分辨率和提高储层预测精度。
尽管上面已经通过结合示例性实施例描述了本发明，但是本领域技术人员应该清楚，在不脱离权利要求所限定的精神和范围的情况下，可对本发明的示例性实施例进行各种修改和改变。

Claims (3)

1.一种基于全变差最小化约束的地震反射系数反演方法，其特征在于，所述反演方法包括以下步骤：
A、从叠后地震数据S中提取地震子波w并进行傅立叶变换得到地震子波w的频域表示W(f)，其中，
W(f)＝FFT(w) (1)
B、从所述叠后地震数据S中取一个地震道的一个时窗内的地震数据s进行傅立叶变换，得到地震数据s的频域表示S(f)，其中，
S(f)＝FFT(s) (2)
C、根据地震数据形成原理由所述地震数据s的频域表示S(f)和所述地震子波w的频域表示W(f)得到反射系数的频域表示R(f)，其中，
R(f)＝S(f)/W(f) (3)
D、以所述时窗中心为分析点，对所述时窗内的地震数据s的时域反射系数进行傅立叶变换，得到反射系数的频域表达R(f)＇，其中，
<mrow> <mi>R</mi> <msup> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mfrac> <mi>N</mi> <mn>2</mn> </mfrac> </munderover> <mo>&amp;lsqb;</mo> <msub> <mi>r</mi> <mi>i</mi> </msub> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mfrac> <msub> <mi>T</mi> <mi>i</mi> </msub> <mn>2</mn> </mfrac> </mrow> </msup> <mo>+</mo> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mo>-</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mfrac> <msub> <mi>T</mi> <mi>i</mi> </msub> <mn>2</mn> </mfrac> </mrow> </msup> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
在等式(4)中，N为由点数表示的分析数据长度，Ti代表反射系数对之间的时间间隔，j为虚数单位，f为频率；
E、在所述步骤C得到的反射系数的频域表示R(f)和所述步骤D得到的反射系数的频域表达R(f)＇相等的基础上，根据反射系数奇偶分解原理构建所述时窗内反射系数反演的目标函数方程；
F、利用共轭梯度算法求解所述目标函数方程以得到所述时窗内反射系数的奇分量和偶分量，并重构出所述时窗内的反射系数，其中，在每次迭代过程中对共轭梯度算法所求得的解进行全变差最小化约束，同时将经全变差最小化约束之后的解作为下一次迭代计算的初始解；
G、在所述地震道上按步长step滑动所述时窗，得到下一个时窗的地震数据，重复步骤B至F，直到所述时窗遍历所述地震道的地震数据，完成所述地震道的反演，得到所述地震道的反射系数；
H、取下一地震道的地震数据并按照步骤B到G进行处理，直到得到每一个地震道的反射系数。
2.根据权利要求1所述的基于全变差最小化约束的地震反射系数反演方法，其特征在于，所述步骤E包括：用欧拉公式将所述步骤D得到反射系数的频域表达R(f)＇中的指数项展开，得到等式(5)：
<mrow> <mi>R</mi> <msup> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mfrac> <mi>N</mi> <mn>2</mn> </mfrac> </munderover> <mo>&amp;lsqb;</mo> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mo>-</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>cos&amp;pi;fT</mi> <mi>i</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>r</mi> <mrow> <mi>N</mi> <mo>-</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mi>j</mi> <mi> </mi> <msub> <mi>sin&amp;pi;fT</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
根据奇偶分解原理，等式(5)中实部是反射系数对ri与rN-i+1的偶分量，用re(i,N-i+1)来表示，虚部是反射系数对ri与rN-i+1的奇分量，用ro(i,N-i+1)来表示，得到等式(6)：
<mrow> <mi>R</mi> <msup> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mfrac> <mi>N</mi> <mn>2</mn> </mfrac> </munderover> <mo>&amp;lsqb;</mo> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <msub> <mi>cos&amp;pi;fT</mi> <mi>i</mi> </msub> <mo>+</mo> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>j</mi> <mi> </mi> <msub> <mi>sin&amp;pi;fT</mi> <mi>i</mi> </msub> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
将所述等式(3)分析点移至时窗中点，并结合等式(6)建立目标函数方程，得到等式(7)：
<mrow> <mi>O</mi> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>e</mi> </msub> <mo>,</mo> <msub> <mi>r</mi> <mi>o</mi> </msub> <mo>,</mo> <mi>T</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <msub> <mi>f</mi> <mrow> <mi>l</mi> <mi>o</mi> <mi>w</mi> </mrow> </msub> <msub> <mi>f</mi> <mrow> <mi>h</mi> <mi>i</mi> <mi>g</mi> <mi>h</mi> </mrow> </msub> </msubsup> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;alpha;</mi> <mi>e</mi> </msub> <mo>{</mo> <mi>Re</mi> <mo>&amp;lsqb;</mo> <mfrac> <mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>W</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>&amp;rsqb;</mo> <mo>-</mo> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <msub> <mi>t</mi> <mi>w</mi> </msub> </mrow> <msub> <mi>t</mi> <mi>w</mi> </msub> </msubsup> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>c</mi> <mi>o</mi> <mi>s</mi> <mo>&amp;lsqb;</mo> <mi>&amp;pi;</mi> <mi>f</mi> <mi>T</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mi>d</mi> <mi>t</mi> <mo>}</mo> <mo>+</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;alpha;</mi> <mi>o</mi> </msub> <mo>{</mo> <mi>Im</mi> <mo>&amp;lsqb;</mo> <mfrac> <mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>W</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mn>2</mn> <mi>&amp;pi;</mi> <mi>f</mi> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </msup> <mo>&amp;rsqb;</mo> <mo>-</mo> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <msub> <mi>t</mi> <mi>w</mi> </msub> </mrow> <msub> <mi>t</mi> <mi>w</mi> </msub> </msubsup> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>sin</mi> <mo>&amp;lsqb;</mo> <mi>&amp;pi;</mi> <mi>f</mi> <mi>T</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mi>d</mi> <mi>t</mi> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mi>d</mi> <mi>f</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
在等式(7)中，flow为低频截止频率，fhigh为高频截止频率，αe为偶分量权重，αo为奇分量权重，Re代表着实部，Im代表着虚部，re为反射系数的偶分量，ro为反射系数的奇分量，Δt为相对于分析点的位移量，tw为时窗win的半窗长，T代表反射系数对之间的时间间隔，t为分析点的时间；
将等式(7)在频域和时域离散化并写为奇分量和偶分量的形式，得到下面的等式(8)和等式(9)：
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>Re</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>Re</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>Re</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>e</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>Im</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>Im</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>Im</mi> <mo>&amp;lsqb;</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mi>r</mi> <mi>o</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
在等式(8)和等式(9)中，f1是起始频率，取值为所述flow；fm是截止频率，取值为所述fhigh；
将所述等式(8)和等式(9)分别写为矩阵形式，得到下面的等式(10)和等式(11)：
be＝Ae×re (10)
bo＝Ao×ro (11)
建立目标函数方程：
<mrow> <mi>O</mi> <mrow> <mo>(</mo> <msub> <mi>r</mi> <mi>e</mi> </msub> <mo>,</mo> <msub> <mi>r</mi> <mi>o</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>&amp;alpha;</mi> <mi>e</mi> </msub> <mo>(</mo> <msub> <mi>b</mi> <mi>e</mi> </msub> <mo>-</mo> <msub> <mi>A</mi> <mi>e</mi> </msub> <mo>&amp;times;</mo> <msub> <mi>r</mi> <mi>e</mi> </msub> <mo>)</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&amp;alpha;</mi> <mi>o</mi> </msub> <mo>(</mo> <msub> <mi>b</mi> <mi>o</mi> </msub> <mo>-</mo> <msub> <mi>A</mi> <mi>o</mi> </msub> <mo>&amp;times;</mo> <msub> <mi>r</mi> <mi>o</mi> </msub> <mo>)</mo> </mtd> </mtr> </mtable> </mfenced> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
在等式(10)、(11)和(12)中，αe为偶分量权重，αo为奇分量权重；be是频域反射系数的实部；bo是频域反射系数的虚部；Ae为偶分量变换矩阵，即：
<mrow> <msub> <mi>A</mi> <mi>e</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
Ao为奇分量变换矩阵，即：
<mrow> <msub> <mi>A</mi> <mi>o</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>1</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mn>2</mn> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mn>3</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;pi;f</mi> <mi>m</mi> </msub> <msub> <mi>T</mi> <mrow> <mi>N</mi> <mo>/</mo> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
将所述目标函数方程(12)写为更一般的优化问题，得到目标函数方程(15)：
b＝Ax (15)
在等式(15)中，x为待求的反射系数。
3.根据权利要求2所述的基于全变差最小化约束的地震反射系数反演方法，其特征在于，所述步骤F包括：
将所述目标函数方程(15)的求解问题转化为||Ax＝b||2优化问题，进一步写成二次函数形式：
<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msup> <mi>x</mi> <mi>T</mi> </msup> <mi>Q</mi> <mi>x</mi> <mo>-</mo> <msup> <mi>y</mi> <mi>T</mi> </msup> <mi>x</mi> <mo>+</mo> <mfrac> <mrow> <msup> <mi>b</mi> <mi>T</mi> </msup> <mi>b</mi> </mrow> <mn>2</mn> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>
在等式(16)中，Q＝ATA，y＝ATb，T表示转置运算；
初始化：x0＝0，r0＝b-Ax0，c＝||Q||2，k＝1，k为迭代次数；
循环步骤：
(1)更新梯度方向：
gk＝Qxk-1-y (17)
(2)更新共轭方向：
若k＝1，则pk＝r0；若k≠1，则
<mrow> <msub> <mi>&amp;alpha;</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>g</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>Qp</mi> <mi>k</mi> </msub> </mrow> <mrow> <msubsup> <mi>p</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>Qp</mi> <mi>k</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow>
pk＝-gk+αkpk-1 (19)
(3)更新第k次迭代的解：
<mrow> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>=</mo> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>+</mo> <mfrac> <mn>1</mn> <mi>c</mi> </mfrac> <msub> <mi>p</mi> <mi>k</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>
(4)对xk进行全变差最小化，
<mrow> <mi>X</mi> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mi>r</mi> <mo>,</mo> <mi>j</mi> <mi>c</mi> <mo>&amp;Element;</mo> <mi>D</mi> </mrow> </munder> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>r</mi> <mo>,</mo> <mi>j</mi> <mi>c</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>r</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mi>c</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>r</mi> <mo>,</mo> <mi>j</mi> <mi>c</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>x</mi> <mrow> <mi>i</mi> <mi>r</mi> <mo>,</mo> <mi>j</mi> <mi>c</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow>
其中，D表示整体剖面集，采用对称边界条件，xir,jc表示第ir行第jc个元素；
(5)再次更新残差：
rk＝b-Axk (22)
在等式(22)中，rk为第k次迭代的残差；
停止条件：当残差小于预定阈值时收敛，停止迭代，输出xk；否则进入第k+1次迭代。
CN201610020556.6A 2016-01-13 2016-01-13 基于全变差最小化约束的地震反射系数反演方法 Active CN105467451B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610020556.6A CN105467451B (zh) 2016-01-13 2016-01-13 基于全变差最小化约束的地震反射系数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610020556.6A CN105467451B (zh) 2016-01-13 2016-01-13 基于全变差最小化约束的地震反射系数反演方法

Publications (2)

Publication Number Publication Date
CN105467451A CN105467451A (zh) 2016-04-06
CN105467451B true CN105467451B (zh) 2018-05-15

Family

ID=55605348

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610020556.6A Active CN105467451B (zh) 2016-01-13 2016-01-13 基于全变差最小化约束的地震反射系数反演方法

Country Status (1)

Country Link
CN (1) CN105467451B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105842596B (zh) * 2016-05-24 2018-06-22 四川大学 一种高灵敏度电力电缆局部缺陷诊断方法
CN106054247B (zh) * 2016-05-25 2020-09-29 中国石油集团东方地球物理勘探有限责任公司 基于转换波地震数据求取高精度反射系数的方法
CN106707338B (zh) * 2016-11-18 2018-12-07 中国石油集团东方地球物理勘探有限责任公司 一种强屏蔽下储层高精度预测方法
CN110618450B (zh) * 2018-06-20 2021-07-27 中国石油化工股份有限公司 基于岩石物理建模的致密储层智能化含气性预测方法
CN109212602B (zh) * 2018-09-05 2019-11-08 湖南科技大学 一种改进地震数据分辨率的反射系数反演方法
CN110579800B (zh) * 2019-10-19 2021-06-01 西南石油大学 一种基于高精度同步挤压变换的地震数据数字处理方法
CN111708082B (zh) * 2020-05-29 2022-04-12 成都理工大学 一种随深度变化的深度域地震子波的提取方法
CN111812715B (zh) * 2020-07-22 2021-04-06 中国矿业大学(北京) 一种煤矿工作面围岩多方向全变分正则化层析成像方法
CN113970789B (zh) * 2020-07-24 2024-04-09 中国石油化工股份有限公司 全波形反演方法、装置、存储介质及电子设备
CN115327624B (zh) * 2022-08-02 2023-07-14 西安交通大学 一种地震子波和反射系数的反演方法及反演系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6263284B1 (en) * 1999-04-22 2001-07-17 Bp Corporation North America Inc. Selection of seismic modes through amplitude characteristics
CN102393532A (zh) * 2011-09-06 2012-03-28 电子科技大学 地震信号反演方法
CN104122588A (zh) * 2014-07-30 2014-10-29 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于谱分解提高叠后地震资料分辨率的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7561491B2 (en) * 2005-03-04 2009-07-14 Robinson John M Radon transformations for removal of noise from seismic data

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6263284B1 (en) * 1999-04-22 2001-07-17 Bp Corporation North America Inc. Selection of seismic modes through amplitude characteristics
CN102393532A (zh) * 2011-09-06 2012-03-28 电子科技大学 地震信号反演方法
CN104122588A (zh) * 2014-07-30 2014-10-29 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于谱分解提高叠后地震资料分辨率的方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Geman范数约束的频率域反射系数反演;王本锋等;《石油地球物理勘探》;20140831;第49卷(第4期);第667-671页 *
Reflectivity computation from separated wavefields: application on the Sigsbee2B example;Alba Ordo&ntilde;ez et al.;《2015 SEG New Orleans Annual Meeting》;20151231;第4563-4568页 *
Sparse reflectivity inversion for nonstationary seismic data;Xintao Chai et al.;《GEOPHYSICS》;20140630;第79卷(第3期);第V93-V105页 *
Thin-bed reflectivity inversion;Satinder Chopra et al.;《SEG/New Orleans 2006 Annual Meeting》;20061231;第2057-2061页 *
井间地震资料全变差正则化波形反演;王守东;《石油物探》;20110531;第50卷(第3期);第234-240页 *
预条件共轭梯度反褶积的改进及其应用;&#131428;晓宇;《地球物理学进展》;20061231;第21卷(第4期);第1192-1197页 *

Also Published As

Publication number Publication date
CN105467451A (zh) 2016-04-06

Similar Documents

Publication Publication Date Title
CN105467451B (zh) 基于全变差最小化约束的地震反射系数反演方法
Zhang et al. Parameter prediction of hydraulic fracture for tight reservoir based on micro-seismic and history matching
US11333783B2 (en) Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain
Reynolds et al. Memoir 71, Chapter 10: Reducing Uncertainty in Geostatistical Description with Well-Testing Pressure Data
US8670960B2 (en) Proxy methods for expensive function optimization with expensive nonlinear constraints
CN101206264A (zh) 高分辨率非线性地震波阻抗反演方法
CN107589448A (zh) 一种多道地震记录反射系数序列同时反演方法
Iturrarán-Viveros Smooth regression to estimate effective porosity using seismic attributes
CN103454677B (zh) 基于粒子群与线性加法器结合的地震数据反演方法
CN105388518A (zh) 一种质心频率与频谱比联合的井中地震品质因子反演方法
CN104122581B (zh) 一种叠后声波阻抗反演方法
Yue et al. A multi‐grid method of high accuracy surface modeling and its validation
CN111897006B (zh) 一种基于方位弹性阻抗差异奇异值分解的裂缝密度及方向预测方法及系统与应用
CN106443770A (zh) 一种页岩气地质甜点的预测方法
NO20101734A1 (no) Hastighetsmodeller for en enkeltbronner og for et sett av bronner
CN105005073A (zh) 基于局部相似度和评价反馈的时变子波提取方法
CN104516015B (zh) 一种确定煤层气纵波和横波速度的方法
US20230168405A1 (en) Deep learning architecture for seismic post-stack inversion
CN111897004B (zh) 一种基于大数据分析技术的测井预测方法
CN102353991B (zh) 基于匹配地震子波的物理小波的地震瞬时频率分析方法
Wang et al. Some theoretical aspects of elastic wave modeling with a recently developed spectral element method
US20150127262A1 (en) Method and apparatus for formation tester data interpretation with diverse flow models
CN104062680B (zh) 一种计算波阻抗反演目标函数梯度的方法
Jiang et al. Transdimensional simultaneous inversion of velocity structure and event locations in downhole microseismic monitoring
Bebbington et al. A robust heuristic estimator for the period of a Poisson intensity function

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
CB03 Change of inventor or designer information

Inventor after: Wu Qiubo

Inventor after: Zhang Dongjun

Inventor after: Zou Wen

Inventor after: Huang Dongshan

Inventor after: Wang Ken

Inventor after: Liu Kaiyuan

Inventor after: Zhou Jingjing

Inventor before: Wu Qiubo

Inventor before: Zhang Dongjun

Inventor before: Zou Wen

Inventor before: Huang Dongshan

Inventor before: Wang Xin

Inventor before: Liu Kaiyuan

Inventor before: Zhou Jingjing

COR Change of bibliographic data
CB02 Change of applicant information

Address after: No. 216, No. 216, Huayang Avenue, Tianfu New District, Sichuan, Sichuan

Applicant after: GEOPHYSICAL EXPLORATION COMPANY OF CNPC CHUANQING DRILLING ENGINEERING Co.,Ltd.

Address before: No. 216, No. 216, Huayang Avenue, Huayang Town, Shuangliu County, Shuangliu County, Sichuan

Applicant before: GEOPHYSICAL EXPLORATION COMPANY OF CNPC CHUANQING DRILLING ENGINEERING Co.,Ltd.

CB02 Change of applicant information
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20180330

Address after: No. 189, fan Yangxi Road, Zhuozhou City, Baoding, Hebei

Applicant after: BGP INC., CHINA NATIONAL PETROLEUM Corp.

Address before: No. 216, No. 216, Huayang Avenue, Tianfu New District, Sichuan, Sichuan

Applicant before: GEOPHYSICAL EXPLORATION COMPANY OF CNPC CHUANQING DRILLING ENGINEERING Co.,Ltd.

GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20201109

Address after: 100007 Beijing, Dongzhimen, North Street, No. 9, No.

Patentee after: CHINA NATIONAL PETROLEUM Corp.

Patentee after: BGP Inc., China National Petroleum Corp.

Address before: No. 189, fan Yangxi Road, Zhuozhou City, Baoding, Hebei

Patentee before: BGP Inc., China National Petroleum Corp.