（地球探测与信息技术专业论文）enkf整合动态和地震数据的方法研究.pdf

摘要 随着我国油田进入中后期，提高采收率变得困难。目前，为了更进一步认识油藏 非均质性，在油田中后期进行了大规模的高精度3 d 地震测量；因此，如何结合重复 观测的地震数据和动态数据进一步认识油藏就显得更为重要。 本文在前人研究的基础上，研究了集合卡尔曼滤波整合地震数据和动态数据更新 储层参数。储层静态模型是储层表征的重要内容，通常利用已知信息，结合先验性认 识对未知储层空间进行估计。然而，这样得到的储层静态模型经常与实际生产观测数 据以及地震观测数据存在着较大差异。集合卡尔曼滤波能够通过观测数据反推系统模 型的状态向量，并使得观测数据和正演数据之间的差异达到最小。本文首先利用随机 导航点法的思路优化生成初始模型集合的算法。然后，介绍集合卡尔曼滤波在更新油 藏储层静态参数的应用。根据集合卡尔曼滤波产生的必要条件，介绍了集合卡尔曼滤 波整合动态数据和3 d 地震数据的基本原理。通过实例分析，验证了集合卡尔曼滤波 整合动态数据和3 d 地震数据在更新油藏模型参数方面的有效性，更新后的油藏模型 能够较好地反映储层非均质性，并且与3 d 地震数据有较好的一致性。还比较了利用 3 d 地震数据和4 d 地震数据进行更新油藏模型的差别。通过e n k f 方法将地震数据和 动态数据结合起来描述油藏特征，很好地拟合观测数据，并获得较好的模型。本文通 过试验模型的基础上验证了所采用方法的正确性及有效性。 关键词：集合卡尔曼滤波，静态模型，3 d 地震，整合 s t u d yo na s s i m i l a t i n gd y n a m i cd a t a a n ds e i s m i cd a t a u s i n gt h ee n s e m b l ek a i m a n f i l t e r c h e n y u d o n g ( g e o p h y s i c a lp r o s p e c t i n ga n di n f o r m a t i o nt e c h n o l o g y ) d i r e c t e db yp r o f e s s o rh u a n gx u r ia n dp r o f e s s o ry i nx i n g y a o a b s t r a c t a st h eo i lf i e l de n t e r si nt h em i da n dl a t e rp e r i o di no u rc o u n t r y ，i tb e c o m e sm o r ea n d m o r ed i f f i c u l tt oe n h a n c eo i lr e c o v e r ye f f i c i e n c y a tp r e s e n t ，w em a k eal a r g es c a l eo f3 d s e i s m i cm e a s u r e m e n t sw i t hh i g l lp r e c i s i o nf o rf u r t h e rs t u d y i n gt h er e s e r v o i r h e t e r o g e n e i t y t h e r e f o r e ，h o wt oc o m b i n et h ed u p l i c a t es e i s m i co b s e r v a t i o nd a t aa n dt h ed y n a m i cd a t at o f u r t h e rc h a r a c t e r i z et h er e s e r v o i ri sb e c o m i n gm o r ei m p o r t a n tn o w o nt h eb a s i so fp r e v i o u ss t u d i e s ，ir e s e a r c ht h a tt h ee n s e m b l ek a l m a nf i l t e ri ss t u d i e d t oi n t e g r a t es e i s m i cd a t aa n dd y n a m i cd a t af o ru p d a t i n gr e s e r v o i rp r o p e r t i e s r e s e r v o i r s t a t i cm o d e li sa l li m p o r t a n t p a r to fr e s e r v o i rc h a r a c t e r i z a t i o n i ng e n e r a l ，b a s e do n c o n d i t i o n e dd a t aa n dw i t hap r i o r ii n f o r m a t i o n ，t h ep r o p e r t i e sa tu n s a m p l e dl o c a t i o n sc a n b ee s t i m a t e d h o w e v e r ，t h i sp r o c e s si sn o ta b l et oe n s u r et h er e p r o d u c t i o no f o b s e r v e dd a t a i n c l u d i n gp r o d u c t i o na c t u a ls e i s m i cd a t a t h ee n s e m b l ek a l m a nf i l t e rc a ni n v e r tt h es t a t e v e c t o ro ft h es y s t e mm o d e lt h r o u g ht h eo b s e r v a t i o nd a t a ，a n dm i n i m i z et h ed i f f e r e n c e b e t w e e no b s e r v a t i o nd a t aa n dm o d e l e dd a t a i nt h i ss t u d y ，r a n d o mp i l o t p o i n tm e t h o di s u s e dt oo p t i m i z et h ei n i t i a lm o d e l t h e n ，w ew i l li n t r o d u c et h ea p p l i c a t i o no fe n s e m b l e k a l m a nf i l t e ro nu p d a t i n gr e s e r v o i rs t a t i cp a r a m e t e r s a c c o r d i n gt ot h ea p p l i c a t i o n c o n d i t i o n so ft h ee n s e m b l ek a l m a n ，w ei n t r o d u c et h eb a s i cp r i n c i p l e so fi n t e g r a t i o no f d y n a m i cd a t aa n d3 ds e i s m i cd a t au s i n gt h ee n s e m b l ek a l m a nf i l t e r t h r o u g he x a m p l e s ， w ec a ns e et h ev a l i d i t yo ft h ee n s e m b l ek a l m a nf i l t e ri nt h ei n t e g r a t i o no fd y n a m i cd a t a a n d3 ds e i s m i cd a t a ，w h i c hr e a l i z eab e t t e rc h a r a c t e r i z a t i o no ft h er e s e r v o i rh e t e r o g e n e i t y ， a n dab e t t e rc o n s i s t e n c yw i t h3 ds e i s m i cd a t a a l s ow ec a m p a r e dt h em o d e lu p d a t i n gu s i n g 3 ds e i s m i cd a t aa n d4 ds e i s m i cd a t a u s i n gt h ee n k fm e t h o dt oi n t e g r a t es e i s m i cd a t aa n d d y n a m i cd a t at od e s c r i b er e s e r v o i rc h a r a c t e r i s t i c s ，i tf i t so b s e r v e dd a t aw e l l ，a n do b t a i n sa b e t t e rm o d e l i nt h i sw o r k ，t h es y n t h e t i cs t u d ys h o w st h ea c c u r a c ya n dv a l i d i t yo ft h i s m e t h o d k e yw o r d s ：e n s e m b l ek a l m a nf i l t e r , s t a t i cm o d e l ，3 ds e i s m i c ，a s s i m i l a t e u 关于学位论文的独创性声明 本人郑重声明：所呈交的论文是本人在指导教师指导下独立进行研究工作所取得 的成果，论文中有关资料和数据是实事求是的。尽我所知，除文中已经加以标注和致 谢外，本论文不包含其他人已经发表或撰写的研究成果，也不包含本人或他人为获得 中国石油大学( 华东) 或其它教育机构的学位或学历证书而使用过的材料。与我一同 工作的同志对研究所做的任何贡献均已在论文中作出了明确的说明。 若有不实之处，本人愿意承担相关法律责任。 学位论文作者签名：酞玉拉! 日期：年月日 学位论文使用授权书 本人完全同意中国石油大学( 华东) 有权使用本学位论文( 包括但不限于其 印刷版和电子版) ，使用方式包括但不限于：保留学位论文，按规定向国家有关 部门( 机构) 送交学位论文，以学术交流为目的赠送和交换学位论文，允许学位 论文被查阅、借阅和复印，将学位论文的全部或部分内容编入有关数据库进行检 索，采用影印、缩印或其他复制手段保存学位论文。 保密学位论文在解密后的使用授权同上。 学位论文作者签 指导教师签名： 日期：年月日 日期： 年月 日 中国石油大学( 华东) 硕士学位论文 第一章前言 在我国许多油田的开采已经进入了中后期，在现有的技术条件下采收率只不过在 百分之三十左右；提高采收率是现在迫切需要解决的难题。目前，为了更进一步认识 油藏非均质性，在各大油田开发的中后期进行了大规模的高精度3 d 地震测量；因此， 如何结合观测的地震数据和动态数据进一步认识油藏就显得更为重要。在本文中，我 们研究集合卡尔曼滤波( e n s e m b l ek a l m a nf i l t e r ，简称“e n k f ) 整合地震数据和 动态数据的问题；并且希望利用此方法能够解决这些问题。 在地球流体中，各种运动的具体形式是随着时间和空间变化的。比如大气系统， 其真实的背景误差应该是很大程度上依赖于具体的形势，也就是说理想的背景误差应 该是时变的。最优插值和三维变分所用的背景误差统计量是一个静态的形式，一般被 认为是大体上均匀、各向同性、时间上变化少的。而k a l m a n 滤波( i ( f ) 是用一种系 统方法根据模式动力计算随时间演变的预测误差统计量，以便更准确地描述流体运动 的不确定性；集合卡尔曼滤波是传统的k a l m a n 滤波的一个m o n t ec a r l o 近似，其背景 误差协方差是由一组有限的预测集合估计得到的，也具有随着空间和时间变化的复杂 结构。 然而，k a l m a n 滤波是在假设误差是高斯和线性增长的前提下才可以达到最优结 果，并且其庞大的计算量，尤其是预测误差协方差矩阵的预测，因此在石油、大气、 海洋的多维数系统中的实际应用方面很难实现。e n k f 结合了k a l m a n 滤波和集合预报 的优点；即用一组短期集合预测来估计具有时变的预测误差统计量以整合新的观测信 息。虽然集合卡尔曼滤波计算量也很大；但是由于它的概念简单，不需要作线性假设 和求解模式的切线性。它还适合于并行计算等优点；因此集合卡尔曼滤波是一个目前 比较受关注的方法。 继1 9 9 4 年e v e n s e n 提出集合卡尔曼滤波之后，很多人在这方面做了进一步的研究。 并证明集合卡尔曼滤波能够有效处理强非线性动力系统和多维数的状态空间【1 翻。已 在石油、大气、海洋方面等得到实际应用；而且其理论上的优越性在数值试验中能够 显现出来。 1 1 研究目的与意义 早期的地震勘探主要是解决地下的构造问题。针对我国油田开发的中后时期，提 高采收率越来越难；许多油田也做了大量的重复地震测量；怎样更好地利用这些地震 第一章前言 数据是一个紧迫的问题。本文的主要目的是如何结合这些地震数据和动态数据挖掘出 更准确的地下信息。 油藏模型( 包括静态参数，如孔隙度和动态参数，如流体饱和度) 的建立是油气田 开发过程中的关键问题，将实时的开发动态数据和地震数据特别是延时地震数据整 合，更新油藏模型对开发方案的更新具有十分重要的指导意义。油藏模型模拟涉及到 多学科综合，包括地质、地球物理与油藏工程，从这方面看，本文研究的内容具有重 要的理论意义与实际应用价值。 用油藏数值模拟技术预测油田生产的未来动态，可以为制定最佳的油田开发方案 提供有利依据；然而预测结果的精度取决于历史拟合的结果。历史拟合是整个油藏数 值模拟工作的关键；而传统的历史拟合存在很多的不足。因此，需要一种新型的历史 拟合方法取代传统的历史拟合方法。 传统油藏历史拟合是通过调整油藏模型的参数( 例如渗透率、孔隙度和相渗关系 等) ，使得调整后的模型数值模拟的结果与实际观测的生产数据相匹配。它需要计算 敏感性系数和反复的数值模拟运算。在传统历史拟合过程中模拟器没有设计成为不断 与新得到的生产数据进行匹配的形式。当有新的生产观测数据时，历史拟合过程必须 用所有的观测数据( 包括之前的数据) 重新进行重复处理，这就限制了传统历史拟合技 术的效率和实用性。 集合卡尔曼滤波进行历史拟合是考虑了油藏模型数据概率分布的多实现整体( 集 合) ，通过吸收油藏的生产观测数据对油藏模型整体进行不断地更新，从而得到能够 反映油藏非均质性特征的模型。它是一种蒙特卡罗方法，用多个实现的集合的平均值 代表真实值。利用集合卡尔曼滤波更新得到的油藏模型与实际生产历史数据有较好的 一致性，才能更好地进行生产预测。 集合卡尔曼滤波可以利用岩石物理模型作为桥梁，将油藏的动态数据和地震数据 联系起来。在有地震数据约束下的油藏中，地震数据可以作为观测向量参与到油藏模 型的更新过程中。这样就可以将其所含有的油藏非均质特征的信息揭示出来。地震数 据可以单独地作为观测向量更新油藏模型，也可以与生产数据联合作为观测向量更新 模型，从而得到更能反映油藏的非均质性特征结果。 e n k f 在油藏模型参数更新方面被认为是很有应用前景的一种方法。本课题将通 过集合卡尔曼滤波整合动态和地震数据应用到油藏更新模型参数方面；在动态数据的 基础上加入地震数据，使集合卡尔曼滤波方法可以更好地应用在油藏藏模型参数更新 2 中国石油大学( 华东) 硕士学位论文 中，更加准确地预测油田生产未来动态。 1 2 国内外发展研究现状 卡尔曼滤波最初应用于线性系统的状态更新。它是基于最小均方差推导出来的， 也是维纳滤波器的进一步发展。卡尔曼滤波是一个不断预测修正和递推的过程。在地 震的资料处理中，用卡尔曼滤波进行反褶积，压缩地震子波来提高地震资料的分辨率。 在应用卡尔曼滤波时，首先要建立系统的状态模型；然后通过吸收随着时间变化的观 测向量，不断地更新系统模型中的状态矢量。 鉴于卡尔曼滤波方法在较强非线性的和模型较大的条件下应用的局限性，1 9 9 4 年e v e n s e n 首次引入集合卡尔曼滤波的方法用于更新非线性的海洋模型。此后集合卡 尔曼滤波方法在天气预报、海洋汇图、水文地理等领域得到了广泛的应用。后来，也 逐渐地应用到石油方面。2 0 0 4 年e v e n s e n 通过取样策略和平方根分析的方法，说明了 相同集合大小的情况下，通过初始模型集合的选择、观测系统的扰动方面可以提高集 合卡尔曼滤波的效果。2 0 0 5 年e v e n s e n 讨论了参数结合与状态估计的问题。 1 2 1 集合卡尔曼滤波在其他领域的发展 利用e n k f 进行数据整合，最初的工作是由e v e n s e n ( 1 9 9 4 ) 完成的；并且e v e n s e n 使用了一个试验个例证明了e n k f 能够有效地解决e k f 方法在对非线性系统进行数 据整合方面的问题。由于使用线性近似和闭合方案而带来协方差演变的计算不准确问 题【3 】；随后他们又进行了一些应用试验，e v e n s e n 和v a nl e e u w e n ( 1 9 9 6 ) 讨论了利用 e i l i ；口整合阿加勒斯流域测高数据的应用【4 1 ，后来v a nl e e u w e n 和e v e n s e n ( 1 9 9 6 ) 又 使用集合平滑( e n s e m b l es m o o t h 或e s ) 重新作了该试验并进行了对比分析，发现利用 e s 方法时，在线性系统中e s 的估计精度更高；但对于模式的有色噪声滤除效果比集 合卡尔曼滤波差【5 】。 e v e n s e n ( 1 9 9 7 ) 使用l o r e n z 系统作了应用举例；试验证明：即使是对于l o r e n z 系 统那样一个具有混沌和高度非线性的系统；e n k f 能够有效地追踪系统的状态转变， 计算的误差与实际的误差估计相一致【6 】。 h o u t e k a m e r 和m i t c h e l l ( 1 9 9 8 ) 介绍了一种新的e n k f 应用方法。该方法中将两个 模式状态的集合同时向前积分，用从一个集合预测中得到的误差统计来更新另一个集 合；之所以使用两个集合，是因为这样可以减少分析过程中所遇到的“不繁殖性 问 题。针对这一观点，v a nl e e u w e n ( 1 9 9 9 ) 提出了一些有争议的评论，并且h o u t e k a m e r 3 第一章前言 和m i t c h e l l ( 1 9 9 9 ) 做出了回答【7 1 。 m a d s e n 和c a n i z a r e r s ( 1 9 9 9 ) 在一个2 d 飓风湍流模式中，对比研究了e n k f 和e k f 的降秩方根( r r s q r t ) 整合方法。这是一个弱非线性问题，在此项试验中发现e n k f 和e k f 的结果十分一致【引。 e c h e v i n 等( 2 0 0 0 ) 使用普林斯顿海洋模式对e n k f 进行了研究：并且重点研究了多 元协方差函数的水平结构和海表面以上的垂直结构。试验结果表明，e n k f 能够跟踪 由于海岸线和海岸动力学引起的非各向同性的协方差函数；并且在这个领域中算法简 单是一个特殊的优点【9 1 。 h o u t e k a m e r 和m i t c h e l l ( 2 0 0 1 ) 继续研究了使用两个集合的整合方法；同时引入了 一个方案来进行含有大量观测数据的e n k f 的全球分析。同时他们也引入了一种方法 来整合由于有限的集合尺寸所导致的最终的虚假长期相关项【1 0 】。 r e i c h l e 等( 2 0 0 2 ) 成功的讨论了e n k f 的“描述解 。他们发现e n k f 能够成功的 向“描述解”逼近；但由于e n k f 中误差高斯分布的假设，在分析时刻会引起一定的 偏差。在所谓的描述方法中克服了这一缺点，并且描述方法能够获得最大可能性平滑 解【i l 】。 n a t v i k 和e v e n s e n ( 2 0 0 3 ) 第一次展示了e n k f 在海洋的生态模式中理想的三维应 用。他们的文章证明了如何通过整合w i f s 海洋的有色数据来控制海洋生态模式的演 变。除此之外，还引入几种诊断的方法；这些方法能够有效的用来检验集合的统计属 性t l z , t 3 】。 集合卡尔曼滤波方法可以充分考虑描述误差演变的各阶分量，可以极大的减小存 储代价，提高运算效率。另一方面，集合卡尔曼滤波方法的理论清晰易懂，算法简单 可行，无需构造切线性和伴随算子，可以同时考虑初值误差、观测误差以及模式误差 的存在。在实际应用当中使用真正的观测资料进行数据整合有很好的应用前景。 1 2 2 集合卡尔曼滤波在石油领域的发展 近几年，越来越多的参数敏感性计算方法应用到油田的开发生产中；用于监测油 藏的压力、温度或流动速度等参数。大量监测数据的积累为集合卡尔曼滤波在油藏开 发领域的应用创造了条件。 n a v d a l ( 2 0 0 2 ) 等利用集合卡尔曼滤波方法吸收生产观测数据更新油井附近的油 藏模型，然后进行油藏的生产预测。更新后的模型能够较好地与观测数据保持一致。 g u 和o l i v e r ( 2 0 0 4 ) 用对p u n q 3 油藏模型利用集合卡尔曼滤波方法进行了历史拟 4 中国石油大学( 华东) 硕士学位论文 合。n a v d a l ( 2 0 0 5 ) 对一个二维油藏模型用集合卡尔曼滤波方法吸收生产观测数据， 并连续进行油藏模型更新。更新的油藏模型数据包括动态数据( 压力和含水饱和度) 、 静态数据( 渗透率) 。r j l o r e n t z e n ( 2 0 0 5 ) 用集合卡尔曼滤波方法更新油藏模型的孔 隙度和渗透率；研究了集合卡尔曼滤波效果与初始集合的关系，并且用更新得到的油 藏模型进行了生产的预测。研究表明，随着生产数据不断的被吸收预测的效果逐渐变 好。n l i u ( 2 0 0 5 ) 利用集合卡尔曼滤波进行地质相约束的自动历史拟合。n l i u ( 2 0 0 5 ) 比较了利用基于梯度最小化方法和集合卡尔曼滤波方法在历史拟合中相边界的估计 问题。m z a f a r i ( 2 0 0 5 ) 对集合卡尔曼滤波方法在油藏描述和生产预测中的不确定性 问题进行了分析。g g a o 和r e y n o l d s ( 2 0 0 5 ) 在贝叶斯框架内定量分析比较了p u n q 3 油藏模型利用r m l 方法和集合卡尔曼滤波方法的不确定性问题。h a u g e n ( 2 0 0 6 ) 等 人利用集合卡尔曼滤波在北海油田进行历史拟合。e v e n s e n ( 2 0 0 7 ) 等人利用集合卡 尔曼滤波对北海油藏模型做辅助性历史拟合。 通过叠前地震波形反演估计岩石物理的性质是一个研究热点。由于叠前地震波形 反演的非线性和非唯一性，故它是一项富有挑战性的课题。近来几年，对于线性反演 和全局化反演方法的梯度的局限性，基于连续反演的e n k f 是很有效的工具。在2 0 0 8 年s e g 年会上，l o n gj i n 等人研究了基于e n k f 的一维叠前波形反演【l4 1 ，即与偏移 距或入射角度有关的信息按顺序反演，类似于在石油或地下水水文学方面连续地利用 与时间有关的信息。利用小波变换参数化也可减少模型参数的数目。基于e n k f 的一 维叠前地震波形反演主要方法是逐步拟合与角度相关的振幅。集合卡尔曼滤波状态向 量包括两部分：模型参数( 声阻抗和横波阻抗) 和数据；对于每一步需要计算一个新 的集合卡尔曼滤波因子，并利用它更新模型参数。 随着地震勘探技术的发展，时移地震技术也不断趋向成熟。如何将地震数据和生 产数据联合起来对油藏的特征进行描述是一个重要的研究课题。1 9 9 7 年x u r ih u a n g 首先提出了一个定量的方法，应用地质统计学方法使得4 d 真实振幅和从油藏模拟计 算得来的合成振幅之间的差异匹配值达到最小。该方法被应用到墨西哥湾干气田研究 中，提高了预测的可靠性。2 0 0 5 年，j k h a z a n e h d a r i 利用4 d 地震作为约束结合生产 数据对北海油田中的油藏进行了历史拟合和生产预测研究。 在国内，2 0 0 6 年郑占【i5 】通过利用e n k f 将时移地震的数据应用到油藏的历史拟 合中去，改善仅仅应用生产数据的e n k f 方法对油藏模型参数的估计，形成一种定量 的计算方法；并对方法的效果进行分析。2 0 0 7 年，季秀峰 1 6 】在研究地质统计学的基 5 第一章前言 础上，建立一个实现生成多个实现的方法，使e n k f 能更广泛地应用在油藏历史拟合 中。2 0 0 7 年，凌云、黄旭日、孙德胜等提出了3 5 d 地震勘探方法【1 7 】。3 5 d 地震勘探 方法应用油田开发中、晚期的高精度3 d 地震数据，结合油田开发动态信息解决油田 开发中存在的问题，寻找剩余油气分布区域。通过对中国西部某油田的开发实例研究， 认为：3 5 d 地震可以较有效地解决油田开发中的问题，发现剩余油气的分布，不但 解决了一些油田早期没有3 d 地震或早期3 d 地震数据质量存在严重非重复性噪声影 响的问题，同时还可降低时间推移地震勘探的风险，提高油田开发的经济效益。 在此基础上，本课题的主要目的：借助于油藏历史拟合，利用集合卡尔曼滤波结 合地震数据和油田开发动态信息解决油田开发中存在的问题。将地震数据作为观测向 量的模型更新方法，并设法把地震向量设为单一3 d 地震；探索油田开发中后期采集 的3 d 地震应用方法。 我国的许多油藏为砂、泥薄互层油藏，且地震的分辨率不足以满足我们生产的需 要。因此，如何综合地震和动态数据建立高分辨率的储层参数( 如渗透率) 模型是目 前迫切需要解决的难题。集合卡尔曼滤波方法是此难题的一种解决方法。从在国外的 实际应用效果来看，该方法能够很好地解决整合地震和动态数据中存在的问题，有着 较好的优越性。这种方法在国内尚未见到应用，一方面非常有必要研究此方法对我国 特定油藏的适应性；另一方面，有必要研究此方法的普遍性。 6 中国石油人学( 华东) 硕士学位论文 第二章卡尔曼滤波原理及应用 2 1 卡尔曼滤波的发展 最佳线性滤波理论起源于四十年代美国科学家w i e n e r 和前苏联科学家等人的研 究工作，后人统称为维纳滤波理论。从理论上说，维纳滤波的最大缺点是必须用到无 限过去的数据，不适用于实时处理。为了克服这一缺点，六十年代k a l m a n 把状态空 间模型引入滤波理论，并导出了一套递推估计算法，后人称之为卡尔曼滤波理论。卡 尔曼滤波是以最小均方误差作为估计的最佳准则，来寻求一套递推估计的算法。其基 本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现时刻的观测 值来更新对状态变量的估计，求出当前时刻的估计值。因此，它适合于实时处理和计 算机运算。 卡尔曼滤波是匈牙利数学家卡尔曼在1 9 6 0 年提出的这一种线性滤波和预测的新 算法。它通过递归算法高效地处理离散数据线性滤波问题；卡尔曼滤波器的本质是利 用递归算法获得贝叶斯规则最优解。其中，递归指的是时域内可以高效计算：贝叶斯 规则指的是在已知观测信号的分布的前提下由真实变量的先验估计推算其后验概率 分布：最优解指的是后验估计的误差协方差最小：即最小均方误差准则。 我们知道维纳滤波是根据当前x ( n ) 和过去全部的观测值x ( n 一1 ) ，x ( n 一2 ) ，来 估计信号的当前值s ( 甩) ，它的解形式是以均方误差最小为原则下的系统的传递函数 日( z ) 或单位脉冲响应办( 刀) 。而卡尔曼滤波不需要过去全部的观测值，它是根据前一 个估计值( h 一1 ) 和最近一个观测值x ( 咒) 来估计信号的当前值( 以) ，它是用状态方程 和递推方法进行估计的，因而卡尔曼滤波对信号的平稳性和时不变性不做要求。我们 利用维纳滤波的模型引入到卡尔曼滤波的信号模型。 2 2 卡尔曼滤波的信号模型 2 2 1 状态方程和观测方程 要给出卡尔曼滤波的信号模型，先来讨论状态方程和观测方程。在维纳滤波的模 型，信号可以认为是由白噪声激励一个线性系统的响应，假设响应和激励的时域关系 可以表示： s ( 咒) = a s ( n - 1 ) + w ( n - 1 ) ( 2 1 ) 上式( 2 1 ) 也就是一阶a r 模型。在卡尔曼滤波中信号s ( n ) 被称为是状态变量，用 7 第二章卡尔曼滤波原理及应用 矢量的形式表示为s ( k ) ，在k 时刻的状态用s ( k ) 表示，在k l 时刻的状态用s ( k 1 ) 表 示。激励信号w ( n ) 也用矢量表示为w ( k ) ，激励和响应之间的关系用传递矩阵a ( k ) 来 表示，它是由系统的结构确定的，与彳( z ) 有一定关系。有了这些假设后我们给出状 态方程： s ( k ) = a ( k ) s ( k 1 ) + w ( k - 1 ) ( 2 2 ) 上式表示的含义就是在k 时刻的状态s ( k ) 可以由它的前一个时刻的状态s ( k 一1 ) 来求 得，即认为k 一1 时刻以前的各状态都已记忆在状态s ( k 一1 ) 中。 卡尔曼滤波是根据系统的观测数据( 即量测数据) 对系统的运动进行估计的，所 以除了状态方程之外，还需要观测方程。还是从维纳滤波的观测信号模型入手，观测 数据和信号的关系为：x ( n ) = s 研) 十v ( 托) ，v ( n ) 一般是均值为零的高斯白噪声。在卡 尔曼滤波中，用x ( k ) 表示观测到的信号矢量序列，v ( k ) 表示观测时引入的误差矢量， 则观测矢量x ( k ) 与状态矢量s ( 七) 之间的关系可以写成： x ( k ) = s ( 尼) + ，( 七) ( 2 3 ) 上式( 2 3 ) 和维纳滤波的工0 ) = s o ) + y 0 ) 概念上是一致的，也就是说卡尔曼滤 波的一维信号模型和维纳滤波的信号模型是一致的。 可以把上式( 2 - 3 ) 推广就得到更普遍的多维观测方程 彳( 七) = c ( 七) s ( 尼) + v ( k ) ( 2 4 ) 式( 2 - 4 ) 中的c ( 露) 称为观测矩阵，它的引入原因是，观测矢量的维数不一定与状态 矢量的维数相同，因为我们不一定能观测到所有需要的状态参数。假如x ( k ) 是m l 的矢量，s ( k ) 是甩x l 的矢量，c ( k ) 就是肌，z 的矩阵，v ( k ) 是mx l 的矢量。 2 2 2 信号模型 有了状态方程s ( k ) = a ( k ) s ( k 一1 ) + w ( k - t ) 和观测方程x ( k ) = c ( k ) s ( k ) + 1 ，( 七) 后， 我们就给出卡尔曼滤波的信号模型，如图2 1 所示。 图2 1卡尔曼滤波的信号模型 f i 9 2 - 1s i g n a lm o d e lo fk a l m a nf i l t e r 8 x ( 膏) 孛蘑蟊演大学( 华东) 磺学位论文 2 3 卡尔曼滤波方法 建立好了卡尔曼滤波的信号模型、状态方程以及观测方程后，要解决的问题就 是要寻找在最小均方误差下信号s ( 惫) 的估计值蓉( 露) 。 2 3 1 卡尔曼滤波的一步递推法模型 把状态方程和观测方程重毅给出： s ( 后) 篇a ( k ) s ( k - 0 + 以奄一1 ) ( 2 5 ) r ( 后) 。c ( 七) s ( 豇) + y ( 七) ( 2 6 ) 式( 2 。5 ) 和( 2 - 6 ) 中童是) 和e ( 素) 是已知的，x ( 对是观测到的数据，也是已知 的，假设信号的上个估计值蜃( 彪一1 ) 已知，现在的问题就是如何来求当前时刻的估 计值s ( 意) 。 如果上两式中没有v ( k ) 与w ( k 1 ) ，可以立鼯求得s ( 素) ，估计问题的出现就是因 为信号与噪声的叠加。假设暂不考虑v ( 尼) 与w ( k 一1 ) ，用上两式得到的s ( 詹) 和又( 七) 分 别用季7 ) 和j ( 女) 表示，得： s ( 奄) = a ( k ) s ( k - o ( 2 7 ) x ( 尼) = c ( k ) s ( 七) = c ( k ) a ( k ) s ( k - 1 ) ( 2 8 ) 必然，鼹测值x ( 露) 稿估计值童( 露之间有一定的误差，它们之闻的差岩露) 称为 新息( i n n o v a t i o n ) ： x ( k ) = x ( 露) 一x ( 露) ( 2 9 ) 最然，新息的产生是由于我们前面忽略了v ( k ) 与w ( k 1 ) 所引起静，也就是说新 息里面包含了v ( k ) 与w ( k 一1 ) 的信息成分。因而我们用新息膏( 七) 乘以个修正矩阵 h ( k ) ，用它来代替式( 2 5 ) 的畎毒一1 ) 来对s ( k ) 进行估计： s ( 七) = 彳( 意) s ( 惫一1 ) + 日( 露) x ( 露) = a ( k ) s ( k 一1 ) + h ( 七) x ( 七) 一c ( k ) a ( k ) s ( k - 1 ) 1 ( 2 1 0 ) 囊式( 2 - 5 ) ( 2 - 1 0 ) 可以溺出卡尔曼滤波对s ( k ) 进行估计的递接模型，如图 2 2 所示，输入为观测值x ( k ) ，输出为信号估计值雪( 七) 。 9 第二章卡尔曼滤波原理及应用 。， - 1 、 ， ! r l 、 lll- 1l r ， i 1l r - 贾7 似r 丽1 冈一固 s ( k ) 图2 - 2 卡尔曼滤波的一步递推法模型 f i 9 2 - 2o n es t e pr e c u r s i v em o d e lo fk a i m a nf i l t e r 2 3 2 卡尔曼滤波的递推公式 从图2 2 中容易看出，要估计出( 七) 就必须要先找到最小均方误差下的修正矩阵 h ( k ) ，先后把式( 2 6 ) 和( 2 5 ) 分别代入式( 2 1 0 ) 得： s ( 后) = a ( k ) s ( k 一1 ) + 日( 尼) c ( 后) s ( 后) + v ( 后) 一c ( k ) a ( k ) s ( k - 1 ) 】 = a ( k ) s ( k - 1 ) + 日( 七) c ( 七) 【彳( 后) s ( 七一1 ) + w ( k 一1 ) + ，( 七) 一c ( 七) 彳( 七) s ( 七一1 ) 】 = 1 一月7 ( 后) c ( 后) 彳( 七) s ( 后一1 ) + 日( 尼) c ( 后) 彳( 七) s ( 尼- 1 ) + w ( k 一1 ) 】+ h ( k ) v ( k ) ( 2 11 ) 根据上式来求最小均方误差下的h ( k ) ，然后把求到的h ( k ) 代入( 2 - 1 0 ) 则可以 得到估计值s ( 七) 。 设真值和估计值之间的误差为：雪( 七) = s ( 尼) 一雪( 尼) 误差是个矢量，因而均方误差 是一个矩阵，用s ( 七) 表示。把式( 2 1 1 ) 代入得 s ( k ) = s ( d - s ( d = 4 ( 忌) 【s ( 七一1 ) 一s ( k 一1 ) 】+ 【，一日( 克) c - ( 庇) 】“七一1 ) - h ( k ) c ( 露) 彳( 足) s ( 七一1 ) 一s ( 七一1 ) 卜日( | ) 尼) = ，一日( 七) c ( 七) 彳( 尼) s ( 七- 1 ) 一s ( k 一1 ) 】+ w ( k - 1 ) 卜h ( k ) v ( k ) ( 2 1 2 ) 均方误差矩阵： s ( 尼) = e j ( 七) 雪( 七) r ( 2 1 3 ) 其中，丁表示对向量取共轭转置。为了计算方便，令 s ( 七) = e 【( s ( 七) 一s ( 七) ) ( s ( 七) 一雪( 尼) ) r 】 ( 2 14 ) 利用条件：v ( d 与w ( k ) 都是零均值的高斯白噪声，且它们互不相关，协方差矩 阵分别为研以尼) 以) r 】= q ( k ) 6 ( k - j ) 和研，( 七) 1 ，( ) r 】= r ( k ) 8 ( k 一歹) ； s ( k 一1 ) 与 w ( k 1 ) 不相关；( 足一1 ) 与w ( 庀一1 ) 及v ( 尼) 不相关。它们的相关矩阵分别为 哪删) = o 毁掌哪m ) = 愍岩 1 0 中国石油大学( 华东) 硕士学位论文 还假设状态的初始值与动态噪声、观测噪声，在k 0 时均不相关；并且动态噪 声向量与观测噪声向量也不相关，即有e ( w ( 七) v ( 甩) r ) = o ， f i n ，k 。 找到和均方误差矩阵的关系： s ( 七) = e 【( 彳( 七) s ( 七- 1 ) + w ( k - 1 ) - a ( k ) s ( k - 1 ) ) ( a ( d s ( k - 1 ) + w ( k - 1 ) - a ( k ) s ( k 一1 ) ) 1 】 = 彳( 七) e ( s ( 七一1 ) 一s ( k 一1 ) ) ( s ( 七- 1 ) 一s ( k 一1 ) ) 7 彳( 尼) 7 + e w ( k 一1 ) w ( k 一1 ) 1 】 = a ( k ) e ( k - 1 ) a ( k ) 7 + q ( k - 1 ) ( 2 1 5 ) 把式( 2 1 2 ) 代入式( 2 1 3 ) ，最后化简得： s ( 后) = e 雪( 七) 雪( 七) r 】 = ，一日( 七) c ( 七) 】 彳( 尼) s ( 七一1 ) 彳( 尼) r + q ( 七一1 ) ，一日( 尼) c ( 尼) 】7 + ( 七) 尺( 七) ( 后) r ( 2 1 6 ) 把式( 2 1 5 ) 代入式( 2 1 6 ) 得 c ( k ) = ，一h ( 尼) c ( 尼) s ( 七) ，一日( 尼) c ( 尼) 】r + 日( 七) r ( 七) 日( 尼) r = s ( 七) 一h ( k ) c ( k ) 6 ( 七) 一s ( 足) c ( 七) rh ( 七) r + h ( 七) 【c ( 足) s ( 七) c ( 忌) r + r ( 尼) 】h ( 七) 7 设存在矩阵m ，u ；使得m m r = c ( 尼) 占( 尼) c ( 七) r + r ( k ) ，u = s ( 七) c ( 尼) r ；然后 代入上式化简： 占( 七) = 占( 七) 一h ( k ) u r u h ( k ) r + h ( k ) m m r 日( 尼) r = f ( 七) + 日( 尼) m u ( m r ) 1 】 日( 七) m u ( m7 ) 卅】r u ( m m r ) 1 u r ( 2 1 7 ) 式( 2 1 7 ) 第一项和第三项与修正矩阵h ( k ) 无关，第二项是半正定矩阵，要使 得均方误差最小，则必须h ( k ) m u ( m r ) 一= 0 ，于是可以求得最小均方误差下的修 正矩阵日( 尼) 为： 日( 尼) = 占( 七) c ( 七) r 【c ( 七) 占7 ( 尼) c ( 豇) r + r ( 七) 】一1 ( 2 1 8 ) 把式( 2 - 1 8 ) 代入式( 2 - 1 0 ) 后，我们就可以得到均方误差最小条件下的雪( 七) 递推 公式。 相应的式( 2 1 7 ) 的第二项为零，得最小均方误差为： c ( k ) = 6 ( k ) - u ( m m7 ) 一1 u r = i x h ( 忌) c ( 足) 】s ( 足) ( 2 1 9 ) 综上所述，得到卡尔曼滤波的一步递推公式： 占( 七) = a ( k ) 6 ( k - 1 ) a ( k ) r + q ( k - 1 ) ( 2 2 0 ) h ( k ) = 占( 七) c ( 七) r 【c ( 七) s ( 尼) c ( 尼) r + 尺( 七) 一1 ( 2 - 2 1 ) 6 ( k ) = i - h ( k ) c ( k ) 8 ( 尼) ( 2 2 2 ) s ( 七) = 爿( 七) s ( 尼一1 ) + 日( 尼) x ( 尼) 一c ( k ) a ( k ) s ( k - 1 ) 】 ( 2 2 3 ) 第二章卡尔曼滤波原理及应用 有了上面四个递推公式后我们就可以得到雪( 尼) 和占( 尼) 。如果初始状态s ( o ) 的统计 特性已知，并且令雪( o ) = e s ( o ) 】，占( o ) = e ( s ( o ) 一j ( o ) ( s ( o ) 一j ( o ) r 】= v a r ( s ( o ) ) ，且 矩阵a ( k ) ，q ( k ) ，c ( k ) ，r ( k ) 都是已知的，以及观测量x ( k ) 也是已知的，就能用 递推计算j 。z 。，1 4 所有的( 庀) 和s ( 忌) ：将初始条件s ( 0 ) 代入式( 2 2 0 )