（应用数学专业论文）多电极成象测井的数学模型和正则化.pdf

摘要 本文主要讨论多电极成象测井中的偏微分方程的反问题，将其化约为一个多 参数辩识问题，用拟牛顿方法和遗传算法进行求解。 多电极成象测井是一种新的电阻率测井技术，这一测井技术的电极系中含有 多个测量电极，可提供较多的测量信息，从而有助于以较高的分辨率确定地层电 阻率参数。多电极成象测井反问题可以归结为一个非线性多参数优化问题，其中 待确定的参数较多，计算量较大，而且一般有多个局部极值，是非线性规划中比 较复杂和困难的一类问题。本文中对这个问题建立了连续和离散化的数学模型， 运用非线性优化和演化计算等数学方法对多电极成象测井反演问题提出了数值 求解方法。并用一些数值结果证实了这些算法的有效性。 非线性多参数优化问题的数值方法通常是迭代法，但是一般的迭代法收敛 慢，因此我们采用计算h e s s i a r l 矩阵近似矩阵的拟牛顿方法，其中h e s s i a n 矩阵 的计算采用b f o s 迭代公式，从而避免了大量复杂的二次偏导数的计算。由于 每次迭代都需要计算目标函数的梯度，通常的差分法求梯度至少需要n + 1 次正问 题方程，计算量比较大，误差也比较大。本文采用伴随状态法求目标函数的梯度， 这一方法比一般的差分法更快更精确且不依赖参数的个数n ，只需解一次方程即 可求的梯度。 另外偏微分方程的反问题的一般是不适定的，对于测量资料的微小误差，导 致反演出的值和真值的误差比较大。在实际的测井工作中测量电极上的电位信息 ( 测量资料) 是测量出来的，测量误差不可避免，所以如果仍然用原来的目标函 数，去反演地层的电阻率，反演出来的电阻率的误差就比较大，在本文中我们在 模拟数据加上一定的误差，利用正则化方法再来反演，所得的结果也比较好，这 样就真正符合了实际的需求。此外我们还利用了比较定理改进了优化算法，节省 了计算的时间。 另外，我再利用一种无须计算导数的遗传算法来求解这一问题。遗传算法 模拟生物进化，采用简单的编码技术来表示复杂的结构，通过对一组编码进行遗 传操作和优胜劣汰的自然选择来指导学习和确定搜索方向。遗传算法是搜索和优 化问题中很有用的工具，特别是对一些复杂结构、多维问题尤其有效，因此我们 把遗传算法用于多电极成象测井的反问题的数值求解中。对于这种多参数优化问 题，一般方法可能会收敛到一个局部极值，而遗传算法可以收敛到一个全局最优 解。 关键词：多电极成象测井，反问题，优化算法，搜索，正则化，遗传算法 a b s t r a c t i nt h i st h e s i sw ed e a lw i t ha ni n v e r s e p r o b l e ma r i s i n g i no i l i n d u s t r y ： m a t h e m a t i c a lm o d e la n dm e t h o d sf o rm u l t i - e l e c t r o d e sr e s i s t i v i t yw e l l l o g g i n g t h i si s am u l t i p l ep a r a m e t e r si d e n t i f i c a t i o np r o b l e m w es o l v ei tu s i n gt w od i f f e r e n tw a y s ： q u a s i n e w t o nm e t h o da n dg e n e t i ca l g o r i t h m s m u l t i e l e c t r o d e r e s i t i v i t yw e l l l o g g i n g a n d i m a g i n g i san e w r e s i s t i v i t y w e l l - l o g g i n gt e c h n i q u e ，i n w h i c ht h ee l e c t r o d e a r r a y c o n t a i n s m a n ym e a s u r i n g e l e c t r o d e sa n dp r o v i d e sal o to fm e a s u r e m e n t s i ta l l o w su st oi d e n t i f i e rt h ef o m a t i o n r e s i s t i v i t y w i t h h i g h r e s o l u t i o n t h ei n v e r s e p r o b l e m o f m u l t i p l ep a r a m e t e r s o p t i m i z a t i o n i sad i f f i c u l ta n d c o m p l i c a t e dp r o b l e m i nn o n l i n e a r p r o g r a m m i n g b e c a u s eo f m a n yp a r a m e t e r st ob ed e t e r m i n e da n d a g r e a td e a lo f c a l c u l a t i o na n dw i t h m a n y 1 0 c a lo p t i m u m s i nt h i st h e s i sw eb u i l tt h em a t h e m a t i c a lm o d e lf o rt h ei n v e r s i o n p r o b l e m a n dg a v et h en u m e r i c a lm e t h o d s t os o l v et h ei n v e r s ep r o b l e m c o r r e s p o n d i n g t om u l t i - e l e c t r o d e r e s i s t i v i t yw e l l l o g g i n gp r o b l e m w i t hn o n l i n e a r p r o g r a m m i n g t e c h n i q u e sa n dg e n e t i ca l g o r i t h m s s o m en u m e r i c a lr e s u l t ss h o w nt h ee f f i c i e n c yo f t h e s em e t h o d s g e n e r a l l yw es o l v et h ep r o b l e m sf o rn o n l i n e a rm u l t i p l ep a r a m e t e r so p t i m i z a t i o n w i t hi t e r a t i o nm e t h o d b e c a u s eo ft h es l o wc o n v e r g e n e eo ft h e g e n e r a l i t e r a t i o n m e t h o d ，w ei n t r o d u c eq u a s i - n e w t o nm e t h o du s i n gh e s s i a nm a t r i x w ja r ef r e eo f c a l c u l a t i n g ag r e a td e a lo fc o m p l i c a t e d2 - o r d e rp a r t i a ld i f f e r e n t i a lc o e 街c i e n t b y a d o p t i n gb f g s i t e r a t i o nf o r m u l at oa p p r o x i m a t eh e s s i a nm a t r i x w js h o u l dc o m p u t e t h eg r a d i e n to ft h eo b j e c t i v ef u n c t i o ni ne a c hi t e r a t i o n s ow em u s ts o l v et h eh i g h o r d e re q u a t i o nn + 1t i m e sa tl e a s tw i t ht h eo r d i n a r i l yd i f f e r e n c em e t h o d i nt h et h e s i s w ei n t r o d u c ea d j o i n ts t a t em e t h o dt oc a l c u l a t et h eg r a d i e n t i ti sf a s t e rt 1 1 e nd i f f e r e n c e m e t h o da n dw i t hh i g h e rp r e c i s i o nb e c a u s ei ti si n d e p e n d e n to ft h ep a r a m e t e r sn u m b e r n i na d d i t i o n ，t h ei n v e r s ep r o b l e mi si l l - p o s e dp r o b l e m s t h e r ei sal a r g e re r r o r b e t w e e nt h er e c o n s t r u c t e dv a l u ea n de x a c ts o l u t i o no w i n gt ot h ee x c e p t i o n a l l ys m a l i e r r o ro fm e a s u r i n gd a t a i nf a c t ，t h ee l e c t r i c p o t e n t i a l o ft h ee l e c t r o d e ，w h i c hi s m e a s u r i n g ，i sm e a s u r e db yw e l l - l o g g i n gt 0 0 1 t h e r ei sn ow a yt oa v o i dt h em e a s u r i n g e r r o r i fw es t i l lu s et h eo r i g i n a l l yo b j e c tf u n c t i o nt or e c o n s t r u c t e dt h er e s i s t i v i t y , t h e r e m a yb eag r e a t e re r r o r i nt h i sp a p e r ，w eu s et h er e g u l a r i z a t i o nt or e c o n s t r u c t e dt h e r e s i s t i v i t y , a n dw eg o tab e t t e rr e s u l t ，w h i c ha c c o r d sw i t ht h en e e d si np r a c t i c e w e u s e dc o m p a r i s o nt h e o r e m sa n da m e l i o r a t e do p t i m i z a t i o na l g o r i t h m s i na d d i t i o n ，w ei n t r o d n c ean e wm e t h o d ，g e n e t i c a l g o r i t h m st os o l v et h i si n v e r s e p r o b l e mw h i c hn e e dn o tc o m p u t et h ed e r i v a t i v e g e n e t i ca l g o r i t h m sb a s e do nt h e f u n d a m e n t a lp r i n c i p l e so fd a r w i n i a ne v o l u t i o n t h e s ea l g o r i t h m se n c o d eap o t e m i a l s o l u t i o nt o s p e c i f i cp r o b l e mo ns i m p l ec h r o m o s o m e 1 i k ed a t as t r u c t u r ea n da p p l y r e c o m b i n a t i o no p e r a t o r st ot h e s es t r u c t u r es oa sp r e s e r v ec r i f i c a l i n f o r n q a t i o n g e n e t i c a l g o r i t h m sa r eap o w e r f u lt e c h n i q u ef o rs e a r c ha n do p t i m i z a t i o np r o b l e m ，a n da r e p a r t i c u l a r l yu s e f u li nt h eo p t i m i z a t i o no fc o m p o s i t es t r u c t u r e s s ow eu s et h e mi nt h e n u m e r i c a le v a l u a t i o ns o l v i n gt h ei n v e r s ep r o b l e m t h e yc a ns e a r c ha g l o b a lo p t i m u m s o l u t i o nd i f f e r e n tt h a nl o c a lo p t i m at h a to t h e rm e t h o d p e r h a p sg o t k e y w o r d s ：m u l t i e l e c t r o d e w e l l l o g g i n g ，i n v e r s e p r o b l e m ，o p t i m i z a t i o n a l g o r i t h m ，s e a r c h ，r e g u l a r i z a t i o n ，g e n e t i ca l g o r i t h m s 4 第一章引言 电阻率测井是石油勘探和开发中的重要技术之一，它主要被用来确定目的地 层的电阻率参数，从而为油田的油气储量和分布的评价和估计提供可靠的地质信 息。在众多的电阻率测井方法中，三侧向测井、七侧向测井、双侧向测井及微球 形聚焦测井等方法以被广泛地用来对油田地层进行解释分析，与其相适应的各种 数学方法也得到了相应的发展( 见 7 ) 。随着油田勘探开发工作的不断发展，对 油( 气，水) 层的综合解释要求越来越高。 3 ， 8 等文献讨论了用通常的测井技 术( 如双侧向测井等) 来反演一至二个参数的问题，并提出了很好的结果。但由于 地层及勘探钻井情况的复杂性，一方面，目的层的电阻率可以不是常数( 比如非 均匀复杂介质目的层) ，另一方面，目的层的电阻率也随泥浆侵入的深度和程度 不同面发生变化，这就要求我们能通过测井数据来同时反演多个地层参数，为此 需要收集更多的测量信息。根据这一实际问题，提出了一种新的直流测井技术 一一多电极成象测井技术。这种技术的主要特征是在电极系上配置了许多个测量 电极( 通常为1 0 至2 0 个) ，用这些测量电极可以测得沿井轴多个位置的电位信息。 从而通过利用众多的测量信息来反演多个地层参数，达到以较高的分辨率确定地 层的电阻率，即达到对地层电阻率成象的目的。 多电极成象测井研究工作目前己受到国际上的许多学者的关注，l i u q i n g - h u 等研究了使用八个测量电极的电极系来反演多个电阻率参数的问题( 见 1 0 ，1 7 ，1 9 ) 通过运用变形波昂迭代法( d b i m ) 和共轭梯度法等方法来进行数值 反演，得到了一些结果。另外s u g i oi m a m u r a 等也开展了这方面的研究，他们对 一个具有十个测量电极的多电极测井问题，通过利用有限元素法进行正演和利用 晟小二乘原理进行反演，亦取得了一定的进展( 参见 1 0 ) 。p mv a l l i n g a 等也 研究了具有八个以上测量电极的多电极成象问题，他们所用的反演方法主要是 r i d g e 回归或m a r q u a r t 方法等迭代算法。张建峰和李晶研究了具有十八和二十 个测量电极的多电极成象问题。但上面的研究中基本上没有考虑测量数据的误 差，反演所用的测量数据都是有正问题计算出来的，是没有误差的。在实际的测 井工作中测量数据的误差不可避免。我们知道反演问题一般都具有不适定性，测 量资料的误差对反演结果有比较大的影响。测量数据微小误差可能导致反演结果 和真值的误差比较大，也可能严重失真。 在本文中，我们介绍了多电极成象测井反问题的数学模型和数值方法：通 过运用具有等位面边值条件的偏微分方程，我们先对多电极成象测井正演问题建 立合适的数学模型，并运用有限元素法进行数值求解，然后我们应用反问题的理 论和技巧把反演问题转化为个带有约束的非线性优化问题。为了解决反问题的 不稳定性，我们利用正则化方法来处理这个问题，我们在原来的目标函数上增加 一项正则化项。这样即使在测量数据有误差的情况下我们也能得到比较好的结 果。 首先我们通过运用拟牛顿迭代法和求计算目标函数梯度的伴随状态法，对 非线性优化问题进行数值求解。在反问题计算中由于每次迭代都需要求解一个 高阶线性方程组( 正问题方程) ，计算量很大。一般的迭代( 如最速下降法) 收敛 缓慢，本文引入一种收敛较快的拟牛顿方法，通过运用b f g s 迭代公式计算 h e s s i a n 矩阵，从而避免了大量的二阶导数的计算。显然，每次迭代需要计算目 标函数的梯度，而一般的差分法至少需要求解n + 1 次的正问题方程，计算量太 大，且误差较大，因此我们采用快速计算梯度的伴随状态法，只需求解一次的 正问题方程即可得出梯度，不仅大大节约了计算时间，而且提高了精度。由于 伴随状态法不依赖于参数的个数n ，当n 比较大是，计算的效果也比较好。i f ：l k f 我们将比较定理运用到两个参数的反演问题中，这样比原来的算法节省了计算 时间。 另外，我们利用一种比较新的方法一遗传算法引入到求解这个反问题中。遗 传算法可以解决一般方法无法解决的复杂问题，尤其对于多维极值问题可以达 到更好的效果。一般方法容易收敛到一个局部极值，很难达到全局最优解：而 遗传算法从一组可能解出发，可以搜索到一个全局最优解。 本文的第二章中详细的讨论了多电极成象测井正演问题的数学模型和求解 的数值方法以及反问题的数学模型；第三章中我们介绍了一种对解非线性优化 问题比较有效的数值方法：拟牛顿算法，并利用它和计算目标函数梯度的伴随 状态法进行了数值求解。在这一章中我们对反问题的不稳定性做了一些数值实 验，并利用正则化的方法来处理测量数据有误差的情况。；第四章介绍了利用遗 传算法来求解多电极成象的反演问题。并给出了一些数值结果这些数值结果说 明了多电极成象测井工具和我们的正则化算法的可行性和有效性。在报告的最 后，我们列出了一些有待作进一步探讨的问题。 第二章正问题和反问题的数学模型 多电极成象测井是电法测井的一种，常见的电法测井问题的数学模型及其 理论和方法在 1 中已作了详尽的说明和讨论，这里我们针对多电极成象测井的 具体情形进行说明。 在直流电阻率测井的一些常见方法中，如三侧向测井，七侧向测井中，测井 仪器的电极系的构造是比较简单的，电极的个数通常比较少，而其中测量电极 的个数一般只有几个，因而在沿井轴上某一给定的位置，用这些仪器测量出来 的数据也就只有几个，要想获得较多的测量数据就必须在多个位置进行许多次 的测量。目前地质勘探及开发的发展要求我们对地层的电阻率进行高分辨率成 象，为此我们必须能够收集比较多的测量信息。实际应用中的这一需要引发了 油田测井工程师开发多电极成象测井工具的探索和对应的反演方法的研究。我 们研究的课题中，多电极成象测井方法中的电极系是有多个电极构成的复杂电 极系，它包括三个等距分布的发射电极，和二十个非等距分布的测量电极( 参 见图2 1 ) 。测井时将电极系沿井轴放置在井内某一固定位置，并由发射电极发 射低频电流。由于低频电流可以近似地看成为直流电流，因而在下面我们将运 用直流电阻率测井中的模型来对地层中的位场进行模拟。 图2 1 电极系参考图：a l ，a 2 ，a 3 为发射电极，l ，2 0 为测量电极 2 1 1 地层模型 2 。1 正问题的数学模型 在多电极成像测并中通常设底层是规则、层状的。图2 2 是多电极成像测 井中最常用的地层模型。我们假设地层模型是沿井轴对称的，它由井眼中的泥 浆带，上下围岩及目的地层构成。又假设目的地层又被细分为f 层， k 列( 其 中包含各种侵入带) 。电阻率参数在地层内的每一个区域中均可设为常数。 p m ， 量 p 。 p ， - 2 p 锄 图2 2 地层模型示意图 2 1 2多电极成像测井问题的微分方程和边界条件 在直流电阻率测井时，地层内因直流电的作用形成了一个稳定的电场( 与时间 的变化无关) 。对这样一个稳定电场，我们可以利用一些物理定律( 如欧姆定律等) 来建立合适的数学模型。由于测量电极尺寸较小，我们可以忽略其几何尺寸对地 层中电场的影响。多电极成象测井的正问题，即要求这些测量电极上的电位之值。 为此，我们必须求出地层中整个电场的电位分布状况。 ? ? j 【= 凡n “= u ( x ，y ，z ) 。如文献 1 ( p 卜3 ) 所述，在每一个电阻率等于常数的区域( 例如 泥浆，围岩，目的层内各分块) 内的任一点，电位函数“= u ( x ，y ，z ) 均满足偏微 分方程 瓦8 万1 面8 u ) + 茜( 吉雾) + 瓦8 万l 瓦a u ) = 。 ( 2 1 ) 其中电阻率p = ( 几，p 。n ：，n ，几) 显然为分块常数，p 。为泥浆电阻率， p 。p 。：分别为上下围岩电阻率，p ，m 为目的层内各分区域的电阻率。我们 可在柱面坐标系( r ，妒，：) 下重写这一方程。由于有轴对称的假设，方程应和妒无 关，利用变换r = 、可，妒= t a n 一，上，式( 2 1 ) 可以表示为 石o 万r 面o u ) + 瓦o 万r 瓦o u ) = 。 ( 2 2 ) 此外在地层各区域的交界面上，根据电位和电流的连续性，电位“亦必须 满足下面的交界面条件 “一i ，= “+ l ， ( 2 3 ) 古( 笔) 一1 ，= 去 豢) + 1 ， c z a ， 其中1 3 为交界面上的单位法线向量，它在交界两侧规定取同一方向，角标一， + 分别表示在交界砸两侧取值。 下面我们来考察多电极成象测井所形成电场所应具备的边界条件。一方面， 在距电极无穷远处，由于电极的影响已非常之小，所以其上的电位也应该非常 之小，因此求解方便起见，我们可用足够远处边界r 取代无穷远处边界( 图2 2 右端及上下边界) ，且将这些边界上的电位值取为0 ，就有 “卜= 0 ( 2 5 ) 另一方面，井轴及连接各个电极的绝缘棒形成了求解区域q 的内边界r 2 ，由 于已假设地层是沿井轴线轴对称的，故对称轴上无电流从其上流出或流入( 即 不是电流的源与汇) ，而绝缘棒上亦无电流进入，故在边界r 2 上成立如下的边界 条件 ! 丝：0 p 锄i ( 2 6 ) 其中n 为r ，上的单位外法线向量。我们再来看一看余下的边界，即三个发射电 极的表面r o ，ur o ：ur 0 ，由于这些电极均有金属导体组成，而金属导体的电阻 为o ，故在电场中构成短路，使每个电极上各个点间的电位都相等，即每一个电 极都是一个等位面。但在电极上的电位是未知的，由于三个电极间的彼此绝缘， 它们之间的电位也不一定相等，因此在每个发射电极上应有 u l r o ，= c ，( 未知常数) ，i = 1 , 2 ，3 ； ( 2 7 ) 而发射电极上的电流却是已知的，假设分别为，0 l ，0 2 ，0 3 ，则 土兰= k f 1 ，2 ，3 ； ( 2 8 ) ”1 o0 n 上面式( 2 3 ) 与( 2 4 ) 即所谓的等位面边值条件，它们构成了发射电极上的 边界条件。 综上所述，联合方程以上的边界条件我们可以对多电极成象测井所形成的电 场建立如下的数学模型 昙曙詈 + 鲁( j 罢) = 。 u l 2 0 蚓一o ， “i r ；= c ：( 未知常数) ，f = 1 ，2 ，3 ( 2 9 ) p 譬= ，o ( 已知常数) ，f - 1 ，2 ，3 。口咖 1 4 _ ，= “+ l ， 古( 景) 一l ，= 古( 象) + ， 利用具等位面边值条件的偏微分方程边值问题来对多电极成象测井的正 问题建模的一个主要特点是我们考虑到了发射电极的几何尺寸对电场的影 响，从而有别于其它的模型( 其中发射电极近似为线电极) 。 电法测井正演问题的数值解可以应用有限元素法( 见 1 等) 或数值模式 匹配法( 见 7 ) 等方法来计算。而对上面的多电极成象测井正演问题( 2 9 ) 我们可以证明下节所述的变分原理，从而可以运用有限元素法来求解。 2 2 1 变分模型 2 2 正问题的求解 为求解问题( 2 9 ) ，我们记矿为下面的集合 v = 0 0 ，z ) j 在q 上连续且适当光滑，“i l = o ，“h ，= 常数，i = 1 ，2 ，3 j ， 并记泛函： 喇吲1 k 艄r o u 2 俐卜磐叱。，汜 则变分原理可表述如下： 变分原理：若“= “+ ( o z ) 是定解问题( 2 9 ) 的连续且适当光滑的解，则在 函数集合v 中 = “( ，z ) 使泛函( 2 1 0 ) 达到最小值；反之，在函数集合v 中使 泛函( 2 1 0 ) 达到最小值的函数“= “+ ( r ，z ) 一定是定解问题的解。 此变分原理的证明类似于 1 中p 9 7 上的变分原理的证明，此处不再证明。 根据上述的变分原理，定解问题( 2 9 ) 就可等价地化为在集合v 中求泛函( 2 ，1 0 ) 的极值问题。在这种变分处理中，问题( 2 9 ) 中包含未知函数法向导数竺的边 d 胛 界条件均是自由边界条件，即对v 中使泛函( 2 1o ) 达到极值的函数来说是自动 满足的，从而避免了( 2 6 ) 和( 2 8 ) 式中的求导与积分计算。 2 2 2 有限元素法的计算格式 刚度矩阵的形成 根据上节的变分原理，我们对问题( 2 9 ) 进行离散处理，即可得到有限元素 法的计算格式。 由于定解问题( 2 9 ) 的求解区域q 比较简单，其边界线或交界线均是水平线 或垂直线，因此我们可以采用直交网格来处理，即用两组分别平行和垂直于水 平线的直线来形成直交网格，并在此基础上连接适当的对角线构成相应的直角 三角形有限元素分割( 洋见1 1 ) ，p 5 8 ) 。假设节点总数为n ，下面我们用一些几 何量来给出刚度矩阵的i 素表达式。假设节点为区域q 的内点( 见图2 3 ) 。 j 一1 j 十n ji 母卜l8 3 3 2 j + l 则刚度矩阵的第行的元素为( 详见 1 2 ) j n 0 ，。，o ，k j , ) - n , 0 ，k ，1 ，j ，j ，k j ，j + 1 ，0 ，。，o ，k j , j + n , 0 ，- ，0 ( 2 1 1 ) ，一= 一( 石r lc 。t d - + 2 f 见6 c 。t 届) ， k 川一( 去c o t 叫云c o t 黝， 女n 川2 一丽r 5c s + 丽r 4 c 。t 屈) ， ( 2 1 2 ) 尼“矿一( z 岛 4c o t l 2 , 4 - i - 丽r 3 c 。t p 4 ) ， k j ，= 一( f ，。，+ k j ，一l + k j 、+ t + 。) ， 式中l 是邻近节点的第f 个三角形的重心的，坐标，p 、p 。则为临近_ ， 的四个象限区域的电阻率值，而口，7 , ( i = 1 ，- ，6 ) 为j 点边上各单元的内角( 见图 2 3 ) 。 对q 的边界点亦可作类似处理，从而形成整个刚度矩阵 等位面处理 k = ( k 。) ( 2 1 3 ) 在对模型( 2 9 ) 形成刚度阵后，除了一般的约束处理( 如( 2 。6 ) 的约束) 外， 我们还须对等位面上的边界条件进行处理。不失一般性，我们假设只有个等 位面r 0 ，其上的节点序号为，l + 1 ，埘一l ，m ，为简便起见我们还假设l ：0 。 令要求的未知向量为 由等值面边界条件有 我们作如下面的变量替换 显然 而变换( 2 1 6 ) 的逆变换为 “，= “，+ l = 。= “一i = “ ( 2 1 5 ) v ，= “，一“1 v ，+ l = “f + 1 “f 十2 v m l = “m l 一“m v m = “ ( 2 16 ) v = v j + j 一一v = 0 ， ( 2 。1 7 ) “f = v ，+ v “l + + v m ， “i = v 十! + + v m 甜m l = v m 一】+ v m “m = v 应用上面的变换，泛函( 2 1 0 ) 的离散形式，目f j - - 次型 也可作相应的变换。令 ( 2 1 8 ) 巾= 三2 万7 k 占一垃4 a ：“h ( 2 1 9 ) 9 同时令 其中 又设 占a = 0 ，“，“。) 1 瓦= ( v 。，v ：，v 。) 1 隅 d = j 占。 ，曩 l j ：j 4 = ( “，。) t d 2 = 0 。，) 7 耻卜0 刘 为似一，+ 1 ) 一，+ 1 ) 阶上三角矩阵，由( 2 1 8 ) 我们有 因此 其中 占a = r 艿 8 = r 6 ( 2 2 0 ) ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) i 。，i ：分别是阶数为( f 一1 ) ，0 一m ) 的单位矩阵，故二次型( 2 1 9 ) 可化为 其中 m = 三2 j k 占一生4 1 r “。 月l ! 占，文占一i o l 。 24 r e ” k = 月1 k r ( 2 2 9 ) ( 2 3 0 ) 令 占= 0 l ，一，“。，“一，“。) 1 ， ( 2 3 1 ) 中= 妻占1 最占一占卢， ( 2 3 2 ) 户= ( ，o 加，o ) t 3 3 ) k 。对称】 ik ：k ， ( 2 3 4 ) 【k 3k ：k 2 i 其中 2 k = k ，对称 k ：巨l 。 k 3k ： k 2 k ：= ( 女k j = li = i i ，= k ， k ：= ( 呐，。， = i = ( 2 3 5 ) ( 2 3 6 ) ( 2 3 7 ) 这样，我们完成了等位面的处理，得到了下面的方程 k 占= p 其中j 6 由( 2 3 3 ) 给出。 综上所述，运用有限元素法于方程( 2 9 ) ( 2 3 9 ) 。为简化记号我们将方程( 2 3 9 ) 重写为 k ) “= p ， ( 2 3 9 ) 我们得到了一线性方程组 ( 2 4 0 ) 其中k ( p ) ：乏是刚度矩阵( 对称正定) ，依赖于电阻率p ，p ：芦。有关有限元方 法的细节问题可参阅 1 、 2 等文献。为简便起见，我们在下文中就直接以( 2 9 ) 的离散形式( 2 4 0 ) 作为正问题的支配方程。 2 3 反问题的数学模型 2 3 1 微分方程的反问题 在偏微分方程中，通常研究的是正问题，也就是给定了微分方程形式以及其 解满足的某些条件，如初始条件，边值条件或者混合初、边值条件，求满足给定 条件的解及研究解的一些性质。 然而，在实际的应用中常常会遇到与求解正问题相反的情况，作为微分方程 的解，我们不仅知道它应取的初、边值，而且还可以观测到某些其它附加信息， 但是某些物理参数或几何参数却作为未知量出现在微分方程的系数中，或出现在 方程的右端部分，或初、边值中，要求我们利用解的附加信息去反解未知参数， 这就是反问题。 在实际中更重要的令人更关心的是反问题，计算机强大的计算能力为反问题 开辟了广阔的前景。一般来说反问题是不适定的，即它的唯一性、稳定性都不能 解决，假使观测数据稍稍有些误差，则它的解可能会发生极大的偏差( 我们在第 三章中会给出这样的数值例子) 。如果这样，反演参数就没有什么意义了，长期以 来，曾以为不适定问题不能有实际意义。但是加上适当的条件，如正则化，可以 使它变成适定的问题。因此正则化被广泛采用，尤其是数值试验方面，正则化方 法显示出其强大的作用，很多反演问题都是借助正则化方法实现的，吉洪诺夫开 创了不适定问题的先河，他的不适定问题的解法参见文献【6 是一本被广泛阅 读的书。迄今为止正则化方法一致被人们认为是一种有效的求解不适定问题的方 法。本文主要讨论了数据几乎没有偏差的模型和数据有一定偏差的模型。对于数 据有一定偏差的模型是用正则化方法处理的，即在目标函数后面加上一项正则化 参数。 反问题是一门新兴的学科，发展只有2 0 年时间，理论还不成体系。文献 6 3 1 给出了数理方程反问题的数学结构分类、历史发展状况和主要研究方法。 文献 4 4 归纳了研究反问题的理论方法和数值方法，并指出研究反问题的发展 前景 对反问题存在性研究的主要方法有： 1 反s t u r m l i o u v i u e 法 2 。紧性原理方法 3 积分变换方法 4 非线性分析的一些方法，如不动点理论，压缩映象原理等 5 混合方法，如半群方法，g a l e r k i n 逼近方法。 对唯一性研究的主要方法有： 1 h o l d e r 连续方法 2 单调性方法 3 压缩映象方法 在数值研究方面涌现出许多不同的算法，如：量子反射反演法，迭代方法 单调同伦法，利用脉冲谱技术( p s t ) ，有限元反演法。 2 3 2 多电极成像测井的反问题 在石油勘探开发中，钻井要耗费许多的人力和物力，所以井的个数通常不 多，电法测井的反问题是要在较少的探井中利用特定的测井仪器和相应的分析 解释手段，来估计地层中的一些物理参数。具体来说，就是利用测井仪器在地 层中产生稳定电场，并在井轴的某些位置上测量电场的电位信息，通过测量的 这些电位信息来对地层中的电阻率分布进行解释。而多电极成象测井反问题则 是通过运用测量到的s 个电位信息( 比一般的测井方法所获的信息要多) 来确 定目的地层中的电阻率，即已知 “? ，f _ 1 ，s ) ，来求目的地层的电阻率。在实 际当中，我们所得到的电位信息和电场中真实的电位信息之间是有测量误差的， 由于反问题一般是不适定，初始数据的微小误差，就可能导致解的失真，为了 解决这一问题，许多人就针对不同的问题采取了各种不同的正则化方法。 首先我们针对利用测量所得到的s 个电位信息来识别地层中的电阻率参数 p 的反演问题来建立合适的数学模型。 我们使用加权最小二乘法技巧，可将已知恤? l ，求解p 的反问题化成_ f 述求 最小值问题： m i n o ( p 1 ， ( 2 4 ，1 ) p e t 其中c 互月“为p 的允许取值集合，通常可以取成闭集合兀墨。b f m ，p “j ，即每个 电阻率参数p ，必须约束在一个有界区间 p p ，“ 中；而且目标函数 g ( p ) ，和q ( 户) 则定义为 q ( p ) = - ，2 w ，1 “j ( p ) 一? 1 2 ， ( 2 4 2 ) l s ( p ) = 丢窆w ，岍p ) 一“，f 2 + d 州， ( 2 4 3 ) 其中 o ( i = 1 ，s ) 为加权常数，可根据测量误差等因素进行选取，例如 w i 可 取为，其中q 为测量“? 的估计误差。a 为正则化参数，一般取值比较小。i p l i ，r 为d 的半范数。 式( 2 4 2 ) 和( 2 4 3 ) 中的“j ( p ) 是当p 为已知时用支配方程( 2 4 0 ) 计算出 来的电场中对应与第i 测量电极的模拟电位值，即( p ，“? ( p ) ) 满足于支配方程 ( 2 4 0 ) 。易知q ( p ) ，亘( p ) 是p 的非线性函数，上面的问题( 2 4 1 ) 实质上是一个 非线性优化问题，解这样一个问题的常用方法是梯度迭代法。我们注意到，对 单个电极的具有等位面边界条件的边值问题， 8 证明了电位对电阻率的导数的 存在性。这个结果可以推广至我们这个问题，即对问题( 2 9 ) ，兰，f _ 1 ，n 存 o p ， 在，故目标函数9 的梯度v q 也存在，从而可以运用涉及梯度的迭代法来求解 ( 2 4 1 ) 。 第三章拟牛顿法求解反问题 在实际中更重要的令人更关心的问题是反问题。在石油勘探开发中，钻井 要耗费许多的人力和物力，所以井的个数通常不多，电法测井的反问题是要在 较少的探井中利用特定的测井仪器和相应的分析解释手段，来估计地层中的一 些物理参数。具体来说，就是利用测井仪器在地层中产生稳定电场，并在井轴 的某些位置上测量电场的电位信息，通过测量的这些电位信息来对地层中的电 阻率分布进行解释。而多电极成象测井反问题则是通过运用测量到的s 个电位 信息( 比般的测井方法所获的信息要多) 来确定目的地层中的电阻率，即已 知( “? ，f _ 1 ，s ) ，来求目的地层的电阻率。在实际当中，我们所得到的电位信 息和电场中真实的电位信息之间是有测量误差的，由于反问题一般是不适定， 初始数据的微小误差，就可能导致解的失真，为了解决这一问题，许多人就真 对不同的问题采取了各种不同的正则化方法。 为了方便起见，我们在下文中记待确定的电阻率参数为 p = ( p p 。) ( = k f ) 。换句话讲，反问题是要找出一组p 的值，使得方程 ( 2 4 0 ) 的对应于这组p 的解辩f 一，与测量电位 “? f 一，最大程度地接近。这一章 我们将讨论这一反问题的数学模型并用拟牛顿方法求解。 3 1 拟牛顿方法( q u a s i n e w t o n ) 几乎所有的解无约束非线性优化问题的数值方法都是迭代方法，目口给出目 标函数q ( 尸) 最优值p + 的一个初始估计p o ，迭代出一个数列 p 。 使得在某些有 关p o 及q 的条件的约束下， p 。 收敛于p 。而大部分迭代方法是下降算法，即 在几乎每一个迭代中都有q ( p “1 ) q ( p 。) 。为了说明下降算法，我们先定义下降 方向：d 为函数q 在点p 的下降方向( d o w n h i l l ) ，当且仅当v q ( p ) 7 d o 。 由t a y l o r 展开我们容易证明，当a 足够小时， q ( p 十口d ) q ( p ) 。一般 的下降算法可以写成下面的形式： ls 1 令ko ； 2 计算下降方向d 使得v q ( p ) 7d 0 ； 3 计算a 。使得 q ( p 。+ 口d k ) q ( p 。) 4 计算p “：p “= p + 口d ： 5 若p “1 满足给定的收敛准则，即停止： 6 令k = k + l ，转至2 通常算法的收敛准则( 或称停止法则) 是由下面的一些条件组成： 渺“一p k f | q ， ( 3 1 ) 粉p “) - q ( p ) j l h 。( p + ) 。在构造的许多种方法中，被证明行之有效的一种方法是 ：峨一型鱼塑型 其中 s = p “1 一p 。 y = v q ( p “) 一v q ( p ) 具体算法可写成 + f 1 + 三! ，鱼2 1 1 初始化：k = 0 ，p 。= p o ，计算q ( p 。) ，v g ( p o ) 2 确定下降方向d 。= 一日i 1 v q ( p ) ( 3 5 ) 3 利用a r m i j i o 规则确定步长吼，( 见文献 6 ) ，使得q ( p “1 ) g ( p ) ，其中 p “= p + 0 1 k d 女； 4 验证给定停止准则：满足即结束，否则转至5 ； 5 取k = k + 1 ，计算q ( p ) ，v q ( p ) ，并用b f g s 公式计算日。，转至2 3 2 用伴随状态法计算目标函数的梯度 由上面所介绍的算法可见，在每一次迭代中都需要计算目标函数q 及其梯度 v q ，而计算梯度最常用的办法是差分法，利用的差分公式为( 3 6 ) 但是这一方 法计算量大，每求一次梯度要解n + 1 次方程( 2 4 0 ) ，而且微扰量难确定，误差比较 大，我们介绍一种新的快速计算梯度的方法伴随状态法。 _ d q 。q ( p + a _ p j