CN101109821A - 基于系统辨识提高地震资料分辨率的方法 - Google Patents

基于系统辨识提高地震资料分辨率的方法 Download PDF

Info

Publication number
CN101109821A
CN101109821A CNA2007100170290A CN200710017029A CN101109821A CN 101109821 A CN101109821 A CN 101109821A CN A2007100170290 A CNA2007100170290 A CN A2007100170290A CN 200710017029 A CN200710017029 A CN 200710017029A CN 101109821 A CN101109821 A CN 101109821A
Authority
CN
China
Prior art keywords
model
seismic
parameter
resolution
section
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.)
Pending
Application number
CNA2007100170290A
Other languages
English (en)
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.)
SHENGLI PETROLEUM ADMINISTRATION OF SINOPEC GROUP
Original Assignee
SHENGLI PETROLEUM ADMINISTRATION OF SINOPEC GROUP
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 SHENGLI PETROLEUM ADMINISTRATION OF SINOPEC GROUP filed Critical SHENGLI PETROLEUM ADMINISTRATION OF SINOPEC GROUP
Priority to CNA2007100170290A priority Critical patent/CN101109821A/zh
Publication of CN101109821A publication Critical patent/CN101109821A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于系统辨识提高地震资料分辨率的方法。主要包括下述步骤:(一)高分辨剖面预处理;(二)地层吸收系统模型构建;(三)地层吸收模型的系统辨识实现,该部分包括①系统辨识基本结构,②系统辨识的非参数系统模型实现,③系统辨识的参数系统模型实现;(四)地层吸收模型的系统结构验证;(五)地层吸收模型的系统响应求取;(六)地层吸收模型响应的三维空间外推;(七)地面地震资料的高频拓展实现等。本发明综合应用测井、井间地震、地面地震等多种尺度的地球物理资料,利用系统辨识方法,对地面地震资料进行保真性高频恢复,这对于油田开发后期储层的精细描述和微构造的准确认识具有重要的意义。

Description

基于系统辨识提高地震资料分辨率的方法
技术领域
本发明涉及石油地球物理勘探领域多尺度地球物理资料的联合应用,并利用现代信号处理技术提高地面地震资料的分辨率。具体是一种基于系统辨识提高地震资料分辨率的方法。
背景技术
提高地面地震资料分辨率技术一直是国内外研究的热点和难点。它在油田的勘探开发中有着重大的理论和实用价值。由于地层介质的滞弹性和非均质性,地震波在地层中的传播要经历吸收,即衰减和频散作用,理解、估计、补偿地震波的这种吸收作用对于提高地震资料的分辨率至关重要,而仅仅依靠地面地震资料是很难得到准确的地层吸收信息。常用的提高分辨率技术(反褶积、谱白化等),由于缺乏相对可靠的高分辨率目标数据体,仅仅利用低分辨率的地面地震叠加剖面求取到有效反褶积算子,因此在方法的有效性及实用性方面有一定的局限性。
现阶段,针对与地面地震记录对应的地质目标,可以得到以下形式的高分辨率剖面,即高分辨率目标数据体:①井间地震反射波成像剖面。井间地震是一种在井中激发,井中接收地震波的地球物理技术。由于观测和采集方式与常规地震不同,井间地震技术避开了地表低降速带对地震高频信号的吸收,可以获得较高分辨率的成像资料。利用井间地震资料可以得到两井之间精细的层析速度剖面和纵向分辨能力较高的反射波成像剖面;②井间地震合成地震记录剖面。首先由井间地震技术得到的两井之间精细层析速度剖面,根据Gardner公式得到对应的密度剖面(或者直接给定密度剖面),然后利用Zoppritz方程得到两井之间的反射系数剖面,最后基于地震合成记录的基本褶积模型,将高频地震子波与两井之间的反射系数剖面进行褶积,得到纵向分辨能力较高的合成地震记录剖面;③测井合成地震记录。与井间地震合成地震记录剖面产生类似,利用时深标定后的测井声波曲线,根据Gardner公式得到密度曲线(或者直接测量的密度曲线),再利用Zoppritz方程得到井点的反射系数曲线,最后基于基本褶积模型,将高频地震子波与反射系数曲线进行褶积,得到井点高分辨能力的合成地震记录。测井合成地震记录是井间地震合成地震记录在两井间距离为零时的特殊状态;④其它方法得到的纵向分辨能力较高的目标数据体。尽管这些高分辨率剖面产生的形式或方法不同,但它们都是对同一目标地质体在不同尺度上的性质进行刻画,因此,通过对常规地面地震剖面和高频信息衰减较少的高分辨剖面的联合应用,研究地震波在地层传播中的衰减规律,并对常规地面地震衰减的高频信息进行补偿,从而提高其分辨率有重要的实际意义。
系统辨识是系统和控制学科的一项重要技术,它是根据系统的输入、输出信号利用现代信号处理方法来估计它的数学模型,其三个基本要素是:数据、模型结构和准则。系统辨识形成了较完善的理论和技术体系,在控制、预测、规划、仿真研究、优化、过程监视和故障诊断等方面有着广泛应用。把地层或者岩石的物理性质作为系统进行研究,同样有重要的实践意义。当前,系统辨识在油气勘探开发中的应用还很少且不是很成熟,其重要原因之一就是在很多情况下系统辨识的基本要素不完备,即缺乏输入或输出数据,只能把理论分析或者模板数据补充为系统的输入输出信息,从而限制了系统辨识的应用范围和效果。对于地层系统而言,也只能对地面地震信号本身进行相关高阶谱分析和盲反褶积研究。目前国内外尚没有利用系统辨识方法,联合应用高分辨率剖面以提高地面地震资料分辨率的报道。
发明内容
本发明的目的是针对地面地震资料在地层分辨能力方面的不足,联合应用同一目标地质体的高分辨率剖面和地面地震等多尺度地球物理资料,建立高分辨剖面、地面地震和地层系统的响应模型,然后利用系统辨识方法得到地层吸收系统的系统函数,即反褶积算子,并对地面地震资料进行高频拓展,以提高地面地震资料的分辨率,实现对地下目标构造的更精确描述,而提出的基于系统辨识提高地震资料分辨率的方法。
地面地震资料所面临的一个重要问题就是其频带窄,分辨率不高,不能够提取或理解隐藏其中的细微地层信息。仅仅依靠地震资料本身信息,即使是采用了面向目标地区的处理技术,对地震资料进行重新处理,在很多情况下,仍然难以解决实际问题。
由于地层对地震波高频部分的吸收,常规地面地震的叠加剖面对地下目标的垂向分辨率很低;高分辨剖面在垂向上有较高的分辨率。因此,本发明的思路就是通过对同一目标地质体高分辨剖面和常规地面地震剖面的联合分析和处理,构建地面地震、高分辨剖面和地层系统吸收特性的响应模型,利用现代信号处理技术体系的系统辨识技术,求取地层系统的吸收特性,并进行适当的空间外推,作用于地面地震剖面,补偿地面地震资料由于地层吸收所导致的高频损失,拓宽其频带,以提高地面地震资料的分辨率。
本发明的技术方案,即基于系统辨识提高地震资料分辨率的方法主要包括下述步骤:
(一)高分辨剖面预处理
测井、井间地震和地面地震等多尺度地球物理技术能够反映地下同一目标物体不同尺度的性质。由于测量环境、测量仪器、处理技术等方面的不同,由这些不同尺度地球物理技术得到的成果数据信息呈现各自不同的特点。这就需要对这些不同技术得到的高分辨剖面进行预处理,实现高分辨剖面和地面地震资料在时深空间、采样率空间的匹配等。高分辨剖面的预处理包括以下三个方面:①根据给定的时深关系把高分辨剖面从深度域转换到时间域;②把深时转换后的高分辨剖面分别在纵向和横向上重采样为地面地震模式,使两者在形式上相互匹配;③对高分辨剖面的各个道赋以地面坐标,实现高分辨剖面各个道与地面地震剖面道的一一对应。
(二)地层吸收系统模型构建
地层对地震波高频部分的吸收,导致地面地震的叠加剖面对地下目标的垂向分辨率很低;而高分辨剖面避开或者没有经过地表低降速带对地震波高频部分的吸收,其垂向上有较高的分辨率。
设定地面地震经高频衰减后的地震子波采样序列为wx(n),地层反射系数序列为r(n),则基于褶积模型的地面地震叠加记录采样序列x(n)可记为(n=1,2,……N):
Figure A20071001702900121
针对地下同一个地质目标,即同一个地层反射系数序列r(n),设定高分辨剖面没有经过高频衰减的地震子波采样序列为wy(n),则高分辨剖面叠加记录y(n)为:
地层对地震波的高频吸收相当于特定的地层系统对地震波的吸收滤波。设定地层吸收系统的单位脉冲响应为则地面地震的地震子波与高分辨剖面地震子波的关系为:
Figure A20071001702900124
因此,可得
Figure A20071001702900125
Figure A20071001702900126
由此,得到地面地震叠加剖面与高分辨剖面的关系,即地面地震剖面可认为是高分辨剖面经过地层吸收系统后的输出。所以问题转化为求得地层吸收系统的单位脉冲响应或者频率响应。设h(n)为地层吸收系统的逆脉冲响应,I为单位冲激信号,则
Figure A20071001702900131
式(5)可以转换为
Figure A20071001702900132
因此,以地面地震叠加信号为输入参数,以高分辨剖面信号为输出参数,就可以利用下述步骤(三)、(四)、(五)所描述的系统辨识方法,求得地层吸收系统的逆脉冲响应,即反褶积算子,进而通过步骤(六)中系统响应的三维空间外推,实现步骤(七)所描述的地面地震资料的高频拓展。
(三)地层吸收模型的系统辨识实现
①系统辨识基本结构
如式(6)所述,在单位采样间隔下,输入地面地震信号为x(n),
Figure A20071001702900133
输出高分辨剖面地震信号为y(n),
Figure A20071001702900134
则输入输出关系,即地层吸收系统的时域参数基本模型描述为:
Figure A20071001702900135
Figure A20071001702900136
其中,q为移动算子,H(q)是系统的传递函数,v(n)为不可预测干扰,e(n)为方差为③的白噪音,v(n)可以进一步表述为噪音系统G(q)对白噪音滤波后的结果。H(q)x(n)是控制学科领域线性时不变系统响应的表示方式,相当于两者的褶积关系,可表述为:
Figure A20071001702900137
Figure A20071001702900138
序列为系统的冲激响应;h(k)表示输出在k时刻的响应;H(q)是令得到的频率函数(频率响应函数),即
Figure A200710017029001311
v(n)的其功率谱
Figure A20071001702900141
表示为:
Figure A20071001702900142
这样,就可以利用输入地面地震道信号x(n),所对应的输出高分辨剖面地震道信号y(n),基于一定的误差准则v(n),求得该道地层吸收系统的逆脉冲响应,即反褶积算子。
对式(8)的系统辨识实现方法包括:非参数模型实现和参数模型实现。
②系统辨识的非参数系统模型实现
在非参数模型系统辨识中,通过输入信号和输出信号,直接估计系统的单位脉冲响应或者频率响应。
直接估计系统单位脉冲响应的方法是相关法。对于时域参数基本系统模型,可以设计一个高阶有限脉冲响应模型
Figure A20071001702900143
利用最小均方误差方法直接估计单位脉冲响应的各个采样值:
Figure A20071001702900144
直接估计系统频谱响应的方法是谱分析。定义x(n)和y(n)间的协方差函数
Figure A20071001702900145
其傅立叶变换为
Figure A20071001702900146
设定x(n),y(n)相互独立,则由方程可得如下的频谱关系:
Figure A20071001702900147
Figure A20071001702900148
此时,可得到系统的频谱响应为:
Figure A20071001702900149
Figure A200710017029001410
③系统辨识的参数系统模型实现
在参数模型系统辨识中,不是简单的直接指定系统以
Figure A20071001702900151
为变量的特征函数H,G,而是通过设定合理的分子和分母参数,把H,G描述为以
Figure A20071001702900152
变量的函数。
典型的描述系统特性的参数模型为ARX模型。ARX模型设定:
Figure A20071001702900153
Figure A20071001702900154
Figure A20071001702900155
Figure A20071001702900156
其中,B,A分别为延迟算了
Figure A20071001702900157
的多项式,na,nb为多项式阶数,nk为滞后阶数。
则该模型为:
Figure A20071001702900158
Figure A20071001702900159
针对系统描述模型,给定观测到的输入-输出数据x,y,模型预测误差可由下式计算:
Figure A200710017029001511
因此预测误差e(n)是H,G的函数,而H,G由可由上述多项式表示。则得到H,G中多项式参数的估计值的常用方法最小均方误差准则估计方法。该方法通过使式(24)达到最小得到参数的估计值,即:
Figure A200710017029001513
(四)地层吸收模型的系统结构验证
利用地面地震剖面和高分辨剖面进行系统辨识首先要确定系统的结构。在系统辨识中,合理选择系统模型的结构参数是极为重要的问题。基于给定的输入输出数据,根据一定的误差准则和评判标准选择最优的系统模型结构和最佳结构参数,得到满足我们需要的适宜模型,称之为系统模型验证。系统模型验证是系统辨识的关键技术之一。
对同一模型不同阶数参数的比较分析,是模型验证最常用的方法。对于给定的数据系列,它通过对同一模型不同结构参数时的误差参数的比较,寻求最优的结构参数。常用的误差参数为Akaike定义的最终误差(FPE,Final PredicationError)、信息能量准则(AIC,Akaike Information Theoretic Criterion)以及Rissanen定义的最小描述长度(MDL,Minimum Description Length)。
针对ARX模型,其次要确定ARX模型的结构参数,即na,nb,nk的值。对于给定的地层模型,我们只有一个输入数据系列(即地面地震叠加信息),和对应的输出数据系列(即井间地震叠加信息),所以无法用其它的输入-输出数据系列来对模型的应用效果进行验证,而只能用估计模型的参数进行系统的验证。对于此种单输入输出系统的系统验证,误差参数FPE,AIC一般是随着模型结构参数的增大而减小,因此,比较适宜用MDL参数对模型的结构参数进行评判。
对于给定的高分辨剖面,其每个叠加记录、对应的地面地震叠加记录以及所对应位置地层的吸收系统函数,构成了系统辨识的基本要素,共构成若干个系统模型。对这些系统的ARX模型,分别计算阶数na,nb,nk在1-10范围内的误差参数FPE,AIC,MDL,从中选择各误差参数最小的值为各个系统的最佳结构参数。
(五)地层吸收模型的系统响应求取
针对给定的高分辨剖面,以地面地震的各个共深度点(CDP)叠加记录为输入参数序列,与各个CDP记录距离最近的高分辨剖面叠加记录为输出参数,根据步骤(四)中误差参数MDL所确定的模型结构参数为最佳参数,利用最小均方误差准则估计方法,对给定的系统模型进行模型参数的求取,即系统辨识,得到各个CDP记录所对应地层的吸收模型的响应函数。
对该剖面的所有地层吸收系统响应进行分析、比较,确定所构建的地层吸收系统结构的适应性和所得到系统响应的可靠性。
(六)地层吸收模型响应的三维空间外推
针对所处理区块的所有高分辨剖面,首先得到各个剖面所对应的平均系统响应,即平均反褶积算子,然后对其进行距离加权的三维空间插值,得到地面地震每个叠加记录CDP处的反褶积算子。
(七)地面地震资料的高频拓展实现
利用步骤(六)所得到的各个CDP记录处的反褶积算子,基于式(9)即可对三维空间地面地震资料的各个道进行处理,再经过相位校正和道均衡,得到高频拓展后的三维数据体。
此外,上述步骤(一)中的高分辨剖面可以是①井间地震经过VSP-CDP成像或者偏移成像后的反射波叠加剖面;②井间地震合成地震记录剖面;③测井合成地震记录剖面;④其它技术所得到的高分辨剖面。高分辨剖面在给定的区域可以有很多个。
上述步骤(一)中的时深关系可以是根据所研究地区的测井合成记录标定的时深关系,也可以是在测井合成记录标定时深关系基础之上,参考井间地震层析速度所得到的时深关系。
上述步骤(二)中的地层吸收系统是把地层对地震波的高频吸收作用抽象为一个线性时不变系统。
上述步骤(三)中的系统辨识非参数模型实现和系统辨识参数模型实现是系统辨识的两种基本方法,在实际应用中,可以根据需要选择相应的方法。
上述步骤(三)中的系统辨识参数模型可以是常用结构模型的以下几类等:
Figure A20071001702900181
Figure A20071001702900183
其中,
Figure A20071001702900185
Figure A20071001702900186
Figure A20071001702900187
Figure A20071001702900189
上述步骤(四)中的系统结构参数选择是通过一定参数范围内的穷举,寻求误差最小的结构参数为相应模型的最佳结构参数,再通过对一定范围内所有系统最佳结构参数的分析比较,确定一组最佳结构参数为该区域的最佳系统结构参数。
上述步骤(五)中是利用最小均方误差准则法对系统模型的参数进行估计,其也可以利用最优化领域的其它参数估计方法。
上述步骤(六)中的距离加权是指把三维空间有限点的反褶积算子,通过按距离加权的准则,对其进行三维空间的插值。加权的原则可以是线性插值,也可以建立指数衰减方式的插值等。
上述步骤(七)中的相位校正是指地震道经过反褶积算子的作用后,再通过反褶积算子的逆向作用,实现反褶积算子对地震道作用的零相位化。
上述步骤(七)中的道均衡是指把反褶积作用后的地震道的能量恢复到其作用前的水平。
本发明的效果:综合应用测井、井间地震、地面地震等多种尺度的地球物理资料,利用系统辨识方法得到地层吸收系统响应,对地面地震资料进行保真性高频恢复,突破地面地震资料频带限制所导致的分辨能力的不足,这对于油田开发后期储层的精细描述和微构造的准确认识具有重要理论和实践意义。
本发明实现了两个方面的创新:
①技术路线创新。
提高地面地震资料分辨率是石油地球物理勘探领域持续的难点和热点。前人一般都是在对地面地震资料本身进行统计分析的基础之上,利用反褶积或谱白化等信号处理技术实现地面地震信号高频信息的增强。本研究则拟通过研究不同尺度地球物理资料之间的相互联系,提高地面地震资料的分辨率。
②实现方法创新
测井、井间地震、和地面地震资料是地下地质体在不同测量条件下的储层特征的不同分辨率的反映,即同一个系统对不同输入信号的特征响应,因此在地震波传播和线性系统理论指导下,利用系统辨识等现代信号处理技术,研究地面地震资料的高频恢复技术,得到地下储层更为精细和准确的三维数据体。
附图说明
图1基于系统辨识提高地面地震资料分辨率的方法框图
图2 ARX模型的结构参数分析
图3垦71井间地震剖面地层吸收系统的频率响应和单位脉冲响应
图4井间地震连线剖面高频拓展前后比较示意
图5某横向剖面高频拓展前后比较示意
具体实施方式
下面结合附图和实例的进一步说明,该实施实例只是本发明的一部分,不能作为限制本发明的范围。
针对胜利油田垦71井区,以下进行了联合应用井间地震资料和地面地震资料提高地面地震资料分辨率的方法实施,其实施目的是为了证明,采用本发明的原理对地面地震资料进行高频拓展能够综合高、低频剖面的信息,在保持地层产状和原有层位信息基本不变的条件下,恢复低频剖面中缺失的高频信息,提高其垂向分辨率。
实施方案是:
(一)高分辨剖面预处理
测井、井间地震和地面地震等多尺度地球物理技术能够反映地下同一目标物体不同尺度的性质。由于测量环境、测量仪器、处理技术等方面的不同,由这些不同尺度地球物理技术得到的成果数据信息呈现各自不同的特点。这就需要对这些不同技术得到的高分辨剖面进行预处理,实现高分辨剖面和地面地震资料在时深空间、采样率空间的匹配等。高分辨剖面的预处理包括以下三个方面:①根据给定的时深关系把高分辨剖面从深度域转换到时间域;②把深时转换后的高分辨剖面分别在纵向和横向上重采样为地面地震模式,使两者在形式上相互匹配;③对高分辨剖面的各个道赋以地面坐标,实现高分辨剖面各个道与地面地震剖面道的一一对应。
(二)地层吸收系统模型构建
地层对地震波高频部分的吸收,导致地面地震的叠加剖面对地下目标的垂向分辨率很低;而高分辨剖面避开或者没有经过地表低降速带对地震波高频部分的吸收,其垂向上有较高的分辨率。
设定地面地震经高频衰减后的地震子波采样序列为wx(n),地层反射系数序列为r(n),则基于褶积模型的地面地震叠加记录采样序列x(n)可记为(n=1,2,……N):
Figure A20071001702900201
针对地下同一个地质目标,即同一个地层反射系数序列r(n),设定高分辨剖面没有经过高频衰减的地震子波采样序列为wy(n),则高分辨剖面叠加记录y(n)为:
Figure A20071001702900202
地层对地震波的高频吸收相当于特定的地层系统对地震波的吸收滤波。设定地层吸收系统的单位脉冲响应为则地面地震的地震子波与高分辨剖面地震子波的关系为:
Figure A20071001702900211
因此,可得
Figure A20071001702900212
Figure A20071001702900213
由此,得到地面地震叠加剖面与高分辨剖面的关系,即地面地震剖面可认为是高分辨剖面经过地层吸收系统后的输出。所以问题转化为求得地层吸收系统的单位脉冲响应或者频率响应。设h(n)为地层吸收系统的逆脉冲响应,I为单位冲激信号,则式(5)可以转换为
Figure A20071001702900215
因此,以地面地震叠加信号为输入参数,以高分辨剖面信号为输出参数,就可以利用下述步骤(三)、(四)、(五)所描述的系统辨识方法,求得地层吸收系统的逆脉冲响应,即反褶积算子,进而通过步骤(六)中系统响应的三维空间外推,实现步骤(七)所描述的地面地震资料的高频拓展。
(三)地层吸收模型的系统辨识实现
①系统辨识基本结构
如式(6)所述,在单位采样间隔下,输入地面地震信号为x(n),
Figure A20071001702900216
输出高分辨剖面地震信号为y(n),
Figure A20071001702900217
则输入输出关系,即地层吸收系统的时域参数基本模型描述为:
Figure A20071001702900218
Figure A20071001702900219
其中,q为移动算子,H(q)是系统的传递函数,v(n)为不可预测干扰,e(n)为方差为③的白噪音,v(n)可以进一步表述为噪音系统G(q)对白噪音滤波后的结果。H(q)x(n)是控制学科领域线性时不变系统响应的表示方式,相当于两者的褶积关系,可表述为:
Figure A20071001702900221
Figure A20071001702900222
序列
Figure A20071001702900223
为系统的冲激响应;h(k)表示输出在k时刻的响应;H(q)是令
Figure A20071001702900224
得到的频率函数(频率响应函数),即
Figure A20071001702900225
v(n)的其功率谱
Figure A20071001702900226
表示为:
Figure A20071001702900227
这样,就可以利用输入地面地震道信号x(n),所对应的输出高分辨剖面地震道信号y(n),基于一定的误差准则v(n),求得该道地层吸收系统的逆脉冲响应,即反褶积算子。
对式(8)的系统辨识实现方法包括:非参数模型实现和参数模型实现。
②系统辨识的非参数系统模型实现
在非参数模型系统辨识中,通过输入信号和输出信号,直接估计系统的单位脉冲响应或者频率响应。
直接估计系统单位脉冲响应的方法是相关法。对于时域参数基本系统模型,可以设计一个高阶有限脉冲响应模型
Figure A20071001702900228
利用最小均方误差方法直接估计单位脉冲响应的各个采样值:
直接估计系统频谱响应的方法是谱分析。定义x(n)和y(n)间的协方差函数
Figure A200710017029002210
其傅立叶变换为设定x(n),y(n)相互独立,则由方程可得如下的频谱关系:
Figure A200710017029002212
Figure A200710017029002213
此时,可得到系统的频谱响应为:
Figure A20071001702900231
Figure A20071001702900232
③系统辨识的参数系统模型实现
在参数模型系统辨识中,不是简单的直接指定系统以
Figure A20071001702900233
为变量的特征函数H,G,而是通过设定合理的分子和分母参数,把H,G描述为以
Figure A20071001702900234
变量的函数。
典型的描述系统特性的参数模型为ARX模型。ARX模型设定:
Figure A20071001702900235
Figure A20071001702900236
Figure A20071001702900238
其中,B,A分别为延迟算子的多项式,na,nb为多项式阶数,nk为滞后阶数。
则该模型为:
Figure A200710017029002310
Figure A200710017029002312
针对系统描述模型,给定观测到的输入-输出数据x,y,模型预测误差可由下式计算:
Figure A200710017029002313
因此预测误差e(n)是H,G的函数,而H,G由可由上述多项式表示。则得到H,G中多项式参数的估计值的常用方法最小均方误差准则估计方法。该方法通过使式(24)达到最小得到参数的估计值,即:
Figure A20071001702900241
Figure A20071001702900242
(四)地层吸收模型的系统结构验证
利用地面地震剖面和高分辨剖面进行系统辨识首先要确定系统的结构。在系统辨识中,合理选择系统模型的结构参数是极为重要的问题。基于给定的输入输出数据,根据一定的误差准则和评判标准选择最优的系统模型结构和最佳结构参数,得到满足我们需要的适宜模型,称之为系统模型验证。系统模型验证是系统辨识的关键技术之一。
对同一模型不同阶数参数的比较分析,是模型验证最常用的方法。对于给定的数据系列,它通过对同一模型不同结构参数时的误差参数的比较,寻求最优的结构参数。常用的误差参数为Akaike定义的最终误差(FPE,Final PredicationError)、信息能量准则(AIC,Akaike Information Theoretic Criterion)以及Rissanen定义的最小描述长度(MDL,Minimum Description Length)。
针对ARX模型,其次要确定ARX模型的结构参数,即na,nb,nk的值。对于给定的地层模型,我们只有一个输入数据系列(即地面地震叠加信息),和对应的输出数据系列(即井间地震叠加信息),所以无法用其它的输入-输出数据系列来对模型的应用效果进行验证,而只能用估计模型的参数进行系统的验证。对于此种单输入输出系统的系统验证,误差参数FPE,AIC一般是随着模型结构参数的增大而减小,因此,比较适宜用MDL参数对模型的结构参数进行评判。
对于给定的高分辨剖面,其每个CDP记录、对应的地面地震叠加记录以及所对应位置地层的吸收系统函数,构成了系统辨识的基本要素,共构成若干个系统模型。对这些系统的ARX模型,分别计算阶数na,nb,nk在1-10范围内的误差参数FPE,AIC,MDL,从中选择各误差参数最小的值为各个系统的最佳结构参数。
(五)地层吸收模型的系统响应求取
针对给定的高分辨剖面,以地面地震的各个CDP记录为输入参数序列,与各个CDP记录距离最近的高分辨剖面叠加记录为输出参数,根据步骤(四)中误差参数MDL所确定的模型结构参数为最佳参数,利用最小均方误差准则估计方法,对给定的系统模型进行模型参数的求取,即系统辨识,得到各个CDP记录所对应地层的吸收模型的响应函数。
对该剖面的所有地层吸收系统响应进行分析、比较,确定所构建的地层吸收系统结构的适应性和所得到系统响应的可靠性。
(六)地层吸收模型响应的三维空间外推
针对所处理区块的所有高分辨剖面,首先得到各个剖面所对应的平均系统响应,即平均反褶积算子,然后对其进行距离加权的三维空间插值,得到地面地震每个叠加记录CDP处的反褶积算子。
(七)地面地震资料的高频拓展实现
利用步骤(六)所得到的各个CDP记录处的反褶积算子,基于式(9)即可对三维空间地面地震资料的各个道进行处理,再经过相位校正和道均衡,得到高频拓展后的三维数据体。
根据所述步骤(一),首先对7对井间地震叠加剖面,在该地区测井合成记录标定和井间地震层析速度的基础上,把深度域的井间地震叠加剖面转换为时间域,并进行纵向空间和横向空间的重采样(纵向上重采样为与地面地震采样率相一致的2ms,横向上重采样为与地面地震各个道对应的记录道,其间隔为10m),使之能够与对应的地面地震剖面匹配。
根据本发明步骤(二)所描述的原理,以各个井间地震叠加剖面的各道为输出信号,以对应的地面地震记录道为输入信号,以该道记录所对应地层的吸收响应为未知信号,建立三者之间的联系。
根据本发明步骤(三),首先确定系统辨识的结构模型。由于系统辨识的ARX模型与地震勘探中的褶积模型很相似,在模型结构参数的选择上具有很大的灵活性,因此,本实例利用ARX模型(式18-式25),基于最小均方误差准则得到地层吸收系统的响应。
根据本发明步骤(四),确定ARX模型的结构参数,即na,nb,nk的值。对于给定的地层模型,我们只有一个输入数据系列(即地面地震叠加信息),和对应的输出数据系列(即井间地震叠加信息)。对于此种单输入输出系统的系统验证,误差参数FPE,AIC一般是随着模型结构参数的增大而减小,因此,比较适宜用MDL参数对模型的结构参数进行评判。针对垦71.检41-垦71.108井剖面,地面地震记录共有19道CDP记录,联合步骤(一)得到的对应井间地震叠加记录,作为一个输入-输出数据系列。这样,每个CDP记录、其对应的井间地震叠加记录以及所对应位置地层的吸收系统函数,构成了系统辨识的基本要素,共构成19个系统模型。对这19个系统的ARX模型,分别计算阶数na,nb,nk在1-10范围内的误差参数FPE,AIC,MDL,从中选择各误差参数最小的值为各个系统的最佳结构参数。计算结果如图2所示。根据误差参数MDL对模型结构参数的选择结果表明,利用地面地震信息和井间地震信息对地层ARX系统进行辨识的最佳阶数为:
Figure A20071001702900261
根据本发明步骤(五),以地面地震的各个CDP记录为输入参数序列,步骤(一)所确定的其对应的井间地震叠加记录为输出参数,根据步骤(四)所确定的模型结构参数为最佳参数,利用最小二乘法对ARX模型进行模型参数的求取,即系统辨识,得到各井间地震剖面各个道所对应地层的吸收响应及其平均响应。图3为垦71.检41-垦71.108井剖面各个道所对应地层吸收系统响应经归一化后的结果以及该剖面的平均响应,其中(a)和(b)分别为系统频率响应和单位脉冲响应。
根据本发明步骤(六),对步骤(五)所得到的7对井间地震剖面的反褶积算子,通过指数衰减距离加权,得到了三维空间各个地震道的反褶积算子。
根据本发明步骤(七),针对三维数据空间的每个记录道,首先作用以步骤(六)所得到各个道的反褶积算子,拓展其频带,其次进行相位校正,使其反褶积后的相位特征不变,最后经过道均衡,使各个道返褶积后的能量与反褶积前一致,从而得到高频拓展后的三维数据体。图4为垦71井区若干对井间地震剖面连线高频拓展前后的比较示意图。图5为垦71井区某横向剖面高频拓展前后的比较示意图。高频拓展后,地面地震资料主频提高20-30Hz,在地层产状和原有层位信息保持基本不变的条件下,剖面的分辨率得到了较大程度提高。

Claims (10)

1.一种基于系统辨识提高地震资料分辨率的方法,其特征是包括下述步骤:
(一)高分辨剖面预处理
高分辨剖面预处理包括以下三个方面:①根据给定的时深关系把高分辨剖面从深度域转换到时间域;②把深时转换后的高分辨剖面分别在纵向和横向上重采样为地面地震模式,使两者在形式上相互匹配;③对高分辨剖面的各个道赋以地面坐标,将高分辨剖面各个道与地面地震剖面道一一对应;
(二)地层吸收系统模型构建
设定地面地震经高频衰减后的地震子波采样序列为wx(n),地层反射系数序列为r(n),n=1,2,……N,则基于褶积模型的地面地震叠加记录采样序列x(n)可记为:
Figure A2007100170290002C1
针对地下同一个地质目标,即同一个地层反射系数序列r(n),设定高分辨剖面没有经过高频衰减的地震子波采样序列为wy(n),则高分辨剖面叠加记录y(n)为:
Figure A2007100170290002C2
设定地层吸收系统的单位脉冲响应为
Figure A2007100170290002C3
则地面地震的地震子波与高分辨剖面地震子波的关系为:
Figure A2007100170290002C4
因此,可得
Figure A2007100170290002C5
Figure A2007100170290002C6
由此,得到地面地震叠加剖面与高分辨剖面的关系,即地面地震剖面是高分辨剖面经过地层吸收系统后的输出,设h(n)为地层吸收系统的逆脉冲响应,I为单位冲激信号,则
Figure A2007100170290003C1
式(5)可以转换为
(三)地层吸收模型的系统辨识
①系统辨识基本结构
如式(6)所述,在单位采样间隔下,输入地面地震信号为
Figure A2007100170290003C3
输出高分辨剖面地震信号为
Figure A2007100170290003C4
则输入输出关系,即地层吸收系统的时域参数基本模型描述为:
Figure A2007100170290003C5
其中,q为移动算子,H(q)是系统的传递函数,v(n)为不可预测干扰,e(n)为方差为的白噪音,v(n)可以进一步表述为噪音系统G(q)对白噪音滤波后的结果,H(q)x(n)是线性时不变系统响应的表示方式,相当于两者的褶积关系,可表述为:
Figure A2007100170290003C7
序列
Figure A2007100170290003C8
为系统的冲激响应;h(k)表示输出在k时刻的响应;H(q)是令
Figure A2007100170290003C9
得到的频率函数,即
v(n)的功率谱
Figure A2007100170290003C11
表示为:
Figure A2007100170290003C12
由此利用输入地面地震道信号x(n),所对应的输出高分辨剖面地震道信号y(n),基于一定的误差准则v(n),求得该道地层吸收系统的逆脉冲响应,即反褶积算子,
②系统辨识的非参数系统模型
在非参数系统模型辨识中,通过输入信号和输出信号,直接估计系统的单位脉冲响应或者频率响应。
直接估计系统单位脉冲响应的方法是相关法,对于时域参数基本系统模型,可以设计一个高阶有限脉冲响应模型
Figure A2007100170290004C1
利用最小均方误差方法直接估计单位脉冲响应的各个采样值:
Figure A2007100170290004C2
直接估计系统频谱响应的方法是谱分析,定义x(n)和y(n)间的协方差函数
Figure A2007100170290004C3
其傅立叶变换为
Figure A2007100170290004C4
设定x(n),y(n)相互独立,则由方程可得如下的频谱关系:
Figure A2007100170290004C5
Figure A2007100170290004C6
此时,可得到系统的频谱响应为:
Figure A2007100170290004C7
Figure A2007100170290004C8
③系统辨识的参数系统模型
在参数系统模型的辨识中,不是简单的直接指定系统以
Figure A2007100170290004C9
为变量的特征函数H,G,而是通过设定合理的分子和分母参数,把H,G描述为以
Figure A2007100170290004C10
变量的函数,典型的描述系统特性的参数模型为ARX模型,ARX模型设定:
Figure A2007100170290005C1
Figure A2007100170290005C2
Figure A2007100170290005C3
Figure A2007100170290005C4
其中,B,A分别为延迟算子
Figure A2007100170290005C5
的多项式,na,nb为多项式阶数,nk为滞后阶数,则该模型为:
Figure A2007100170290005C6
Figure A2007100170290005C7
Figure A2007100170290005C8
针对系统描述模型,给定观测到的输入-输出数据x,y,模型预测误差可由下式计算:
Figure A2007100170290005C9
因此预测误差e(n)是H,G的函数,而H,G由可由上述多项式表示,则得到H,G中多项式参数的估计值的常用方法最小均方误差准则估计方法,该方法通过使式(24)达到最小得到参数的估计值,即:
Figure A2007100170290005C10
Figure A2007100170290005C11
(四)地层吸收模型的系统结构验证
利用地面地震剖面和高分辨剖面进行系统辨识确定系统的结构。在系统辨识中,通过对同一模型不同结构参数时的误差参数的比较,寻求最优的结构参数。常用的误差参数为Akaike定义的最终误差FPE、信息能量准则AIC或者以Rissanen定义的最小描述长度MDL。对于给定的高分辨剖面,其每个叠加记录、对应的地面地震叠加记录以及所对应位置地层的吸收系统函数,构成了系统辨识的基本要素,共构成若干个系统模型,对这些系统的ARX模型,分别计算阶数na,nb,nk在1-10范围内的误差参数FPE,AIC,MDL,从中选择各误差参数最小的值为各个系统的最佳结构参数;
(五)地层吸收模型的系统响应求取
针对给定的高分辨剖面,以地面地震的各个共深度点(CDP)叠加记录为输入参数序列,与各个CDP记录距离最近的高分辨剖面叠加记录为输出参数,根据步骤(四)中误差参数MDL所确定的模型结构参数为最佳参数,利用最小均方误差准则估计方法,对给定的系统模型进行模型参数的求取,即系统辨识,得到各个CDP记录所对应地层的吸收模型的响应函数,对该剖面的所有地层吸收系统响应进行分析、比较,确定所构建的地层吸收系统结构的适应性和所得到系统响应的可靠性;
(六)地层吸收模型响应的三维空间外推
针对所处理区块的所有高分辨剖面,首先得到各个剖面所对应的平均系统响应,即平均反褶积算子,然后对其进行距离加权的三维空间插值,得到地面地震每个叠加记录CDP处的反褶积算子;
(七)地面地震资料的高频拓展
利用步骤(六)所得到的各个CDP记录处的反褶积算子,基于式(9)即可对三维空间地面地震资料的各个道进行处理,再经过相位校正和道均衡,得到高频拓展后的三维数据体。
2.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是上述步骤(一)中的高分辨剖面可以是如下的一种:①井间地震经过VSP-CDP成像或者偏移成像后的反射波叠加剖面;②井间地震合成地震记录剖面;③测井合成地震记录剖面。高分辨剖面在给定的区域可以有很多个。
3.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(一)中的时深关系是根据所研究地区的测井合成记录标定的时深关系,或者是在测井合成记录标定时深关系基础之上,参考井间地震层析速度所得到的时深关系。
4.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(二)中的地层吸收系统是把地层对地震波的高频吸收作用抽象为一个线性时不变系统。
5.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是上述步骤(三)中的系统辨识的参数系统模型还可以是以下几类等:
Figure A2007100170290007C1
Figure A2007100170290007C2
Figure A2007100170290007C3
Figure A2007100170290007C4
其中,
Figure A2007100170290007C5
Figure A2007100170290007C6
Figure A2007100170290007C7
6.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(四)中的系统结构参数选择是通过一定参数范围内的穷举,寻求误差最小的结构参数为相应模型的最佳结构参数,再通过对一定范围内所有系统最佳结构参数的分析比较,确定一组最佳结构参数为该区域的最佳系统结构参数。
7.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(六)中的距离加权是指把三维空间有限点的反褶积算子,通过按距离加权的准则,对其进行三维空间的插值,加权的原则可以是线性插值,也可以是建立指数衰减方式的插值。
8.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(七)中的相位校正是指地震道经过反褶积算子的作用后,再通过反褶积算子的逆向作用,实现反褶积算子对地震道作用的零相位化。
9.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:上述步骤(七)中的道均衡是指把反褶积作用后的地震道的能量恢复到其作用前的水平。
10.根据权利要求1所述的基于系统辨识提高地震资料分辨率的方法,其特征是:对式(8)的系统辨识实现方法包括:非参数模型实现和参数模型实现。
CNA2007100170290A 2007-08-16 2007-08-16 基于系统辨识提高地震资料分辨率的方法 Pending CN101109821A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CNA2007100170290A CN101109821A (zh) 2007-08-16 2007-08-16 基于系统辨识提高地震资料分辨率的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CNA2007100170290A CN101109821A (zh) 2007-08-16 2007-08-16 基于系统辨识提高地震资料分辨率的方法

Publications (1)

Publication Number Publication Date
CN101109821A true CN101109821A (zh) 2008-01-23

Family

ID=39041972

Family Applications (1)

Application Number Title Priority Date Filing Date
CNA2007100170290A Pending CN101109821A (zh) 2007-08-16 2007-08-16 基于系统辨识提高地震资料分辨率的方法

Country Status (1)

Country Link
CN (1) CN101109821A (zh)

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147322A (zh) * 2011-01-13 2011-08-10 北京工业大学 一种用于考虑土结相互作用的多维多点抗震实验方法
CN102183787A (zh) * 2011-03-07 2011-09-14 中国海洋石油总公司 一种基于地震记录变子波模型提高地震资料分辨率的方法
CN101630016B (zh) * 2008-07-16 2011-10-05 中国石油天然气集团公司 一种提高垂直地震剖面成像质量的方法
CN102520438A (zh) * 2011-12-15 2012-06-27 中国石油集团川庆钻探工程有限公司地球物理勘探公司 分析干扰影响范围的方法
CN102540891A (zh) * 2012-01-17 2012-07-04 中冶南方工程技术有限公司 基于递推增广最小二乘法的结晶器armax模型辨识方法
CN102707317A (zh) * 2010-10-27 2012-10-03 中国石油化工股份有限公司 一种利用地震波吸收衰减特征进行储层分析的方法
CN102937720A (zh) * 2011-08-15 2013-02-20 中国石油化工股份有限公司 井控提高地震资料分辨率的方法
CN103364834A (zh) * 2013-07-29 2013-10-23 成都晶石石油科技有限公司 一种利用叠前地震频散分析预测储层渗透率的方法
CN103412335A (zh) * 2013-08-20 2013-11-27 成都晶石石油科技有限公司 一种利用地震物相体预测储层的方法
CN103714252A (zh) * 2013-12-29 2014-04-09 中国地震局工程力学研究所 基于地震反应记录尾部数据的结构振动台模型模态参数识别方法
CN104407385A (zh) * 2014-12-09 2015-03-11 魏继东 恢复动圈式检波器低频数据的方法
WO2016011587A1 (zh) * 2014-07-21 2016-01-28 王雅苹 一种大倾角地区的垂直地震数据桥式标定方法
CN106610508A (zh) * 2016-11-15 2017-05-03 中国石油天然气股份有限公司 一种裂缝的地震属性判识系统及方法
CN112305610A (zh) * 2020-11-09 2021-02-02 中国石油天然气股份有限公司 地震剖面的分辨率处理方法、装置、设备及可读存储介质
CN112782755A (zh) * 2019-11-07 2021-05-11 中国石油天然气集团有限公司 构建近地表结构模型的方法和装置

Cited By (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101630016B (zh) * 2008-07-16 2011-10-05 中国石油天然气集团公司 一种提高垂直地震剖面成像质量的方法
CN102707317B (zh) * 2010-10-27 2014-08-06 中国石油化工股份有限公司 一种利用地震波吸收衰减特征进行储层分析的方法
CN102707317A (zh) * 2010-10-27 2012-10-03 中国石油化工股份有限公司 一种利用地震波吸收衰减特征进行储层分析的方法
CN102147322A (zh) * 2011-01-13 2011-08-10 北京工业大学 一种用于考虑土结相互作用的多维多点抗震实验方法
CN102183787B (zh) * 2011-03-07 2013-05-29 中国海洋石油总公司 一种基于地震记录变子波模型提高地震资料分辨率的方法
CN102183787A (zh) * 2011-03-07 2011-09-14 中国海洋石油总公司 一种基于地震记录变子波模型提高地震资料分辨率的方法
CN102937720A (zh) * 2011-08-15 2013-02-20 中国石油化工股份有限公司 井控提高地震资料分辨率的方法
CN102520438A (zh) * 2011-12-15 2012-06-27 中国石油集团川庆钻探工程有限公司地球物理勘探公司 分析干扰影响范围的方法
CN102540891A (zh) * 2012-01-17 2012-07-04 中冶南方工程技术有限公司 基于递推增广最小二乘法的结晶器armax模型辨识方法
CN102540891B (zh) * 2012-01-17 2015-01-28 中冶南方工程技术有限公司 基于递推增广最小二乘法的结晶器armax模型辨识方法
CN103364834A (zh) * 2013-07-29 2013-10-23 成都晶石石油科技有限公司 一种利用叠前地震频散分析预测储层渗透率的方法
CN103364834B (zh) * 2013-07-29 2016-08-10 成都晶石石油科技有限公司 一种利用叠前地震频散分析预测储层渗透率的方法
CN103412335A (zh) * 2013-08-20 2013-11-27 成都晶石石油科技有限公司 一种利用地震物相体预测储层的方法
CN103412335B (zh) * 2013-08-20 2015-12-09 成都晶石石油科技有限公司 一种利用地震物相体预测储层的方法
CN103714252A (zh) * 2013-12-29 2014-04-09 中国地震局工程力学研究所 基于地震反应记录尾部数据的结构振动台模型模态参数识别方法
WO2016011587A1 (zh) * 2014-07-21 2016-01-28 王雅苹 一种大倾角地区的垂直地震数据桥式标定方法
CN104407385A (zh) * 2014-12-09 2015-03-11 魏继东 恢复动圈式检波器低频数据的方法
CN106610508A (zh) * 2016-11-15 2017-05-03 中国石油天然气股份有限公司 一种裂缝的地震属性判识系统及方法
CN106610508B (zh) * 2016-11-15 2019-05-07 中国石油天然气股份有限公司 一种裂缝的地震属性判识系统及方法
CN112782755A (zh) * 2019-11-07 2021-05-11 中国石油天然气集团有限公司 构建近地表结构模型的方法和装置
CN112305610A (zh) * 2020-11-09 2021-02-02 中国石油天然气股份有限公司 地震剖面的分辨率处理方法、装置、设备及可读存储介质

Similar Documents

Publication Publication Date Title
CN101109821A (zh) 基于系统辨识提高地震资料分辨率的方法
US7492664B2 (en) Method for processing acoustic reflections in array data to image near-borehole geological structure
CN102508293B (zh) 一种叠前反演的薄层含油气性识别方法
CN102937720B (zh) 井控提高地震资料分辨率的方法
Wang et al. Current developments on micro-seismic data processing
CN106405651B (zh) 一种基于测井匹配的全波形反演初始速度模型构建方法
CN102393532B (zh) 地震信号反演方法
CN105319589B (zh) 一种利用局部同相轴斜率的全自动立体层析反演方法
CN104516018A (zh) 一种地球物理勘探中岩性约束下的孔隙度反演方法
CN105652316A (zh) 一种基于裂缝模型的智能优化地震多属性融合方法
WO2012139082A1 (en) Event selection in the image domain
CN104502997A (zh) 一种利用裂缝密度曲线预测裂缝密度体的方法
CN109738952A (zh) 基于全波形反演驱动的被动源直接偏移成像方法
Li et al. Skeletonized inversion of surface wave: Active source versus controlled noise comparison
CN102565852B (zh) 针对储层含油气性检测的角度域叠前偏移数据处理方法
CN112946752B (zh) 基于裂缝密度模型预测裂缝概率体的方法
EP3215872B1 (en) Compensating for space and slowness/angle blurring of reflectivity
CN110703329B (zh) 基于弱振幅地震反射形成机制的岩性油藏边界确定方法
CA2815211A1 (en) System and method for characterization with non-unique solutions of anisotropic velocities
CN108680957B (zh) 基于加权的局部互相关时频域相位反演方法
CN102353991B (zh) 基于匹配地震子波的物理小波的地震瞬时频率分析方法
CN109613615B (zh) 基于叠前地震响应分析的地质体尺度定量估算方法
CN112305595B (zh) 基于折射波分析地质体结构的方法及存储介质
CN116068663A (zh) 基于磁震联合低频建模的火成岩波阻抗反演方法
Nunziata et al. Active and passive experiments for S-wave velocity measurements in urban areas

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C02 Deemed withdrawal of patent application after publication (patent law 2001)
WD01 Invention patent application deemed withdrawn after publication

Open date: 20080123