利用共轭梯度法反演二度体密度界面的起伏

第28卷第4期2020年10月
Vol.28No.4
Oct.202059石油工业计算机应用
COMPUTER APPLICATIONS OF PETROLEUM
•基础科学•
利用共辄梯度法反演二度体密度界面的起伏
刘星
（陕西能源职业技术学院）
摘要:本文利用H—共觇的性质以及多个彼此共匏向量的构造方法，实现了二次极小化问题的共範梯度法算法。

对于二次函数的极小值问题，从初始点x（x为n维向量）出发,分别沿n个彼此共铤的n维向量，做一维搜索，最多n步可以达到极小点x*。

但实际应用中由于舍入误差影响，往往n步并不能达到极小点x*，此时应当取壬切为新的初始点，重复共觇化过程得到极小点。

利用共觇梯度法反演二度体密度界面的起伏，以观测异常与理论异常差的平方和为目标函数，通过给定不同初值以及反演的平均深度,研究初值与平均深度以及界面起伏对反演结果的影响。

反演结果表明，初值对反演结果影响较小，反演的平均深度与理论模型的平均深度不一致时,反演结果差异较大，界面起伏变大时,反演结果精度减小。

关键字：共辄梯度法；二次极小化；反演；二度体密度界面；收敛；精度
1引言
密度分界面与区域构造、储油构造有密切关系，因此计算密度分分界面的起伏或深度的变化在区域构造研究和石油勘探中具有重要意义。

求密度界面起伏的反演方法确定了密度分界面深度及起伏,实际上也是确定地质体的参数，二者的方法原理也是相同的叫目前主要的方法有线性回归法、压缩质面法、迭代法、频率域反演方法。

线性回归法前提是地下界面起伏平缓，则认为重力变化与界面起伏近似呈线性关系，但当界面起伏较大时,反演结果误差较大。

压缩质面法闵是将密度界面从最小深度h和最大深度H处向中间挤压，使之在界面的平均深度D=罟上压缩成沿X轴方向密度不均匀的物质面,然后将这个物质面以一定的宽度分成许多水平物质带，每一个物质带的面密度为卜=ff-AhnAhj,为第j带界面实际深度旦与平均深度D的差值。

同样压缩质面法也要求密度界面起伏较小。

其反演结果受平均深度取值影响。

压缩质面法实际上也是一种迭代法，迭代计算是提高计算精度的重要手段。

迭代法反演是通过地质与地球理资料建立地球物理模型,给定初值，计算模型初值的重力效应，根据此重力效应与观测值的差值对模型进行修正，继续计算新模型的重力效应，进行修正,直到达到精度要求为止。

频率域反演是Oldenburg™根据Parker^公式提出的一种频率域密度界面迭代反演方法，由于采用快速傅里叶变换，其运算速度较快。

实际上，只要平均深度取得适当,不管采用何种反演方法，都能得到达到一定精度而且相似的结果，关键问题是用于反演的异常必须是单纯由反演目标引起的。

本文根据压缩质面的原理，利用共辄梯度法反演地下二度体密度界面起伏,通过对不同参数试验,研究了反演的初值、平均深度以及界面起伏对反演结果的影响。

2H-共辘的概念及性质
2.1H-共辄的概念
H—共辄的定义:若n维向量R和匕满足
叫心：必（2-1）
则称n维向量R和R是H—共辄的。

若H矩阵取为单位矩阵E,则有
严旳=吓巧=旳=｛：*（2-2）
作者简介:刘星（1989-）,男，汉族，陕西西安,研究生，硕士学位，研究方向为助教。

60石油工业机应用
2020年10月
所以，关于单位矩阵共辄的向量事实上就是我 们熟知的正交向量。

2.2 H-共辄的性质
性质1：设有一组m 个n 维向量J\,P2……,P "彼
此共辄，艮卩
(/ * J ； =(i = j;
=
(2-3)
则Pl , P2……Pm —定是线性无关的O
性质2 ：若n 维向量R 和Pj 是H —共辄的，则新的
向量
Q, =a/>…a^O, Q,=pP ”B"
也是H —共辄的。

即
0*7)
(/=7)
(2-4)
3 H-共辄的构造方法
设山吕……,d "是n 维空间中线性无关的一组 向量，可以通过其线性组合得出一组m 个向量R,
P 2……，P ”彼此H —共辄。

⑴首先取PR 。

(2)为了产生P2,取P2为d2和R 的线性组合。

即有P2=ck + a P-由于P2与Pi H 共辄Pi T HP 2=
PjHPdrF aRTHR=0,得出將o 因此可以
得出
爺3)
+ cd\ 二 〃2 -
p;帆P ： H P 、
设已经求出了 P ” P2……,R,它们彼此是H —共 辄的，为了得出一个向量Pk"与P1,P2……，玖它们彼
此是H 一共辄的，记
r=l
为使P 魚与R , P2……,Pk 彼此H —共辄,则有
用网"+ Xa r ^T H^ =0 (i = 1,2,...,k)
r=l
可以求出 %= _ 7豊;1 (i = l,2”.,k)
由此可以得出
/J=d “+£a,C =d “ - J ：(等器■” (k =0,1,2，―m-l)(3_2)
至此，给出了由m 个线性无关的向量构造m 个
H-共辄的向量的方法。

由该构造过程可以看出，
这种结果是不唯一的。

4二次函数的共辄方向极小化问题
设二次函数/(x) = *M + Z/x + c 的极小点为X* ,
X®为任一给定的起始魚。

若P-P2……,P "彼此H —
共辄。

则它们构成了n 维空间的一个基，所以任何
一个n 维向量都可以唯一的表示成这一组向量的线 性组合，特别对向量有
x •一x<" = a£十+... + %/；
(4-1)
如果能求出Pl, P2……，P",则便能得到最小值
点£这样就把求f(x)的极小问题化为求n 个数色，
a 2……，a "的问题。

下面给出具体求解方法：
设二次函数/(© =討Hr+長+c 的极小点为
X* ,占为一个迭代点，现在从占出发沿给定的方向 R 来求二次函数f(x)极小点。

也就是用精确一维搜
索完成
niin /(**)+ AP k ) = niin F(2)
根据一维极值的必要条件，以及复合函数求导
规则，可以得出
FU )=[w (x “>+“)『人=[v (x ,t<,,)r p k =0(4-2)
对于一般函数均有此式成立。

显然对二次函数则有 W)=«v + Z>.代入上式得
[//(,*> + 才 A ) + A = [(Hr® +b) + XHP t ]p t
= [(/fr u, +Z>)f P t +XP…T H r P… = [v/'(^*))r h + 才"R = 0
求出才一卅
又X 堤f(x)的极小点,则有 W') = <> o 对二次
函数
V/'(x*)=Hr ，+h = 0
(4-3)
另外，由极小点的组成有
x*=x <1>+ a,/；-|-a 2/^ + ... + a…P… (4-4)
对上式两边左乘PJ 得出
+ a^+ a 2P 2 +... + a…P…)=用&“ + 凶用期P^Hx'+P^b=P^tix w +P l r b+a l P l r HP /]r (Hv* +Z>)=/；r (Hv (')+h)+a 1^r H/] = 0
由此得出
=_川(山"'+/，) = _旳(沁)
L 一 ~
7
这个ai 就是从x ⑴出发沿方向R 求二次函数f
(X )极小点的步长X 1*0对任一k,ak 就是从x®出发
沿方向巩求二次函数f(x)极小点的步长X 「。

表达
第28卷第4期刘星:利用共辄梯度法反演二度体密度界面的起伏61
[VcOM
式为勺一
对于二次函数极小值问题，若有了彼此H—共辄的n个方向R,P2……,P”,不管初始点x®如何取,则从初始点出发，分别沿R,P2……,玖做一维搜索,则最多做n次一维搜索便可得出极小点x*o n个单位向量e,(i=1,2”..,")是线性无关的，于是可以利用共辄化过程产生n个共辄方向R, P2……,P”依次便可建立二次函数的共辄方向极小化方法。

下面给出具体流程：
⑴给定初值x(1),并且执行^=e,;l=■k。

(2)严*> +"其中心满足
=_iW2)M
*P：HP*
(3)A+l^/c o若则转向⑸；若心”,则作⑷。

”
⑷求P^e k+ia r P r.%=_等箫.转向⑵。

(5)求出了极小点X=X(k)o
在二次函数的极小化过程中，原则上最多进行n步迭代就可以得到极小点坐标x=x(k+1)，但由于舍入误差的影响会造成号(*“"片0o此时就不能用上面的共共辄化过程做下去，只能取此为新的初始点x®,重复上面的共辄化过程来得到极小点。

5几种常见的共辄梯度方法
5.1Daniel共觇梯度法
Daniel共辄梯度法与前面介绍的二次函数共辄方向的主要区别就是共辄向量的求取上。

Daniel方法在构造共辄向量时用负梯度向量代替了原来的单位向量。

计算步骤如下：
⑴给定初值占，控制误差s0并且执行I斗
(2)如果k<n,则
(2.1)计算8*=盼")=诃+6，以及
r«l
(2.2)严>=尹+佔,其中儿满足
，一一gE
*P；H P*
(2.3)若||.严"-x纠|"•则最优解xW”;否则k+\=>k。

⑶若kMn,则k=l,x(1)=x(k+1),转向(2)。

在该方法中要计算共辄向量的组合系数a”(r=l,2,……,k-l)o另外，存贮共辄向量也要占用很大的内存,现在研究其有关性质。

性质：若&=%=12”.」-1)为共辄向量侧有
g"=W(?*>)T P；=0,(J=12....,*-1)(4-1)证明：因为新的迭代点严"=严+砧.&=W(F“)可以令
y k~Sk
对于二次函数有
”=g“-g*=W(xg")-W(x">)
=(Hr"-"+b)-(Hr">+6)=//(.严"-"*>)=H人P*
若P】,P2……，巩彼此共辄侧有
P：HP)=0,/#}\i,j=1,2”k
所以人/f眄=斤(的)=£工=0,i*；;i,j=1,2”k
又因为严是从严"出发沿方向P"做一维搜索得到的点，前面已经证明有：
=&为T=0
于是对k>j有
gjpj ngjpj+(-g【£+g:/)+(-g:/+g;#)+•••+(-“/+g;詞) =(g J弓-gid)+(g【再-g;/)+•••+©M-g;#)+g;M
将前面两个式子代入得到g"=0。

现用前面介绍的性质对该计算方法进行简化：
已知：g；P严0(丿=1,2...,—1)
y；Pi=P^yj-o(i h j\ij=1,2,…，斤一1)因为号=-g,+Ya,£或％=_弓+£久£(J=L2.....A-1)
所以g：&=g[(-C+^«Z)=-g*+£a,(g；/：)=«
又因为HP,==區产
卍up=g*(g~i-g」_g*g”i_g*g,
A人
这里r=l,2.........,k-l o当r=r=l,2...........,k-2时，g；HP「=0
由此霑出
Pk=一一g*+£a£+a k_,P k_x=一g*+a k_x P k_x
r«l
其中:
«*-!=一P^HPx
修改后的Danie哄辄梯度法计算步骤:
⑴给定初值占，控制误差£。

并且执行S
(2)如果k<n,则
(2.1)计算g*=W*')=H严+b,以及
0(Ar=
+=_g什昭」
T
62石油工业机应用2020年10月
(2.2)严"=卅+",其中\匕满足
—小P*
*P：HP*
卅卜£，，则最优解x*=x(®;
(2.3)若|严
否则k+n。

(3)若2”，则《=W>，转向⑵。

5.2Sorenson-Wolfe共辄梯度法
在前面Daniel法中已经将计算量作了简化。

即
在计算共辄向量R时，将原来的巩计算式变成了
3=-&+陷乩，«*.,=，篇1
从而减少了计算量。

现在讨论叭的计算问题。

a_gjHP*"_UHP-_g/Ha/i)
i Pj H P*“A^PjHP^人(兀/_)
伫严Z g/(g*-gi)"T芒-站g*5r
A t(g*—g*-i)A t^*-1
Sorenson-Wolfe共辄梯度法计算步骤:
(1)给定初值宀，控制误差£。

并且执行1=”o
(2)如；果k<n,则
(2.1)计算g*=W(x">)=H\">+6，以及
0(k=l)什=-&+%a=g/)，i=g/(g*-g—)化>])
(2.2)*“)=x w+\P t,其中X k满足
人一輕
P：HP*
(2.3)若卜z-門二，则最优解x'd";否则A-+1=>Zr o
⑶若kMn,则k=l,x(1,=x(k+1),转向(2)。

5.3Fletcher-Reeves共辄梯度法
在前面Sorenson-Wolfe法中共辄向量组合系数aK-i已经表示为
a_g/(g*_gz)_
^*-i(g*-g*-J A-i yt-i
该式还可简化。

因为由前面的共辄梯度与共辄向量的一个性质有
gJPj=0,(/=1,2.L*-1)
另外，由Daniel修改后的共辄梯度法方案中有
g【g丿=0,0=1,2丄,《-1)
这样HK-!的计算就可以表示为
a=g/(g*-g*-J=g/g*-g/g—=_gJSt 1一母/(&-&1厂岛加-纬/g-一纬Jg-
=_________g/g*________=g/g*=Pk
g*-l(ig*-l+a*_2P»-2)St-lSk-i Pk-\
可以看出Fletcher-Reeves的计算方案更加节省空间和计算量。

⑴给定初值占，控制误差£。

并且执行nk o
(2)如果《<",则
(2.1)计算©=W w)=曲“+b,以及
[0(k=l)
A=g：g”«=U(k>l,0"O)
(2.2)严=x(*>+AR.其中X i满足
_g"
*P：HP*
(2.3)若||严-門w&■则最优解x=x(ktI);否则k+\=>k°
⑶若kMn,则k=l,x(I)=x(k+1),转向⑵。

6共辄梯度法反演二度体密度界面起伏下面我们就共辄梯度法方法在二度体密度界面反演问题中的应用做一简单的研究。

如图6-1所示，已知实测剖面上的重力异常观测值细"心”0),求解地下密度分界面相对于深度D的高度Ah。

把地下界面起伏视作相对深度D的一系列水平长方柱体高度Ah的变化，当D比Ah大很多的时候,近似有
師x,0)=2问：([：
?乙2茗(6-1)
式中,f=6.672X10_11m3/(kg-s2)为引力常数，b 为上下界面密度差,D为某一深度；Ah(£)为基底起伏;x为测点横坐标；E为模型体的横坐标。

将式(6-1)离散化
本go)=2碍(廿胃“殂(>=1,2”.“M)(6-2
)
第28卷第4期刘星:利用共辄梯度法反演二度体密度界面的起伏63
^Ah=(Ah!,Ah2……，Ah m)T，显然Ah为决策变量。

因此，该问题以观测值与理论值之差的平方和建立目标函数：
0⑷)=-虫(兀,0))2
=》(电％,,0)-2./吃
(x.-^)2+D2一般情况下,结合工区的钻井资料，部分测点的△h；是已知的，即Ah=Ah™,那么对于反演结果△h“则这些已知点上有约束
£[JA,-<~]2->min(6-3)
/=!
故该问题最优化模型为：
.»/j//-D
Min:")=工他％”0)-2/吃一'g-4^)2(/=1.2.....M) w 気(*,-£)+D-
“送【型-她V-*min.其中m，WM
!若测点为间距a,模型宽度△&均匀剖分，也为a侧有
A/J/j.D,、
4g(x”0)=2/ir^(/=(6-4)问题的最优化模型变为：
•V y All-D
Min:倾価)=艺(缗"‘(兀,0)-2/aV7------鼻~—«)2(/=1,2，…,M)
A/r
st.艺[劭-Ah；M f t min.其中M匕M。

t=i
令
A=4f%”0);C=(c,)=(2巴“_鳥；+Q：a);h=(All,)则处弘)可以写成
<p(Ah)=(A-m(A-Ch)T=AA t-2ACh r+hCC r h T
做呦的梯度为
=-2AC+2CC r h r
利用共辄梯度法反演地下密度界面起伏中H= 2CC T,b=-2AC
显然该问题是二次函数极小化问题，图6-2为该问题的N-S图。

给定理论模型，正演深度Dl=100m;反演深度D2=100m,观测点距(模型剖分距离)a=lm;观测总点数n=101;剩余密度o'=103kg/m30
图6-3反演结果表明，在初值相等时,误差越小,反演次数越多,但随着误差越来越小，反演次数增加越来越少,实际上这是由于反演方法是共辄梯度法所致，共辄梯度法理论上可以在有限步精确达到二次函数的极小点，当反演次数接近未知量的维
X mi初tth£1挖创极小点收1^,1，大迭代次
&控铝檢度向■大小
迭代次救试初(ftL-l;兴耗向杯k n点两次迭代出杯的腔的敘敕赋初<ftcno«2c ululeierrO>mnil>
g-W(J》-g+6-2CC r x M-2AC',fl-g r g
k^n?
百-----*.1?些一"*"*"
or・0a'i
计算循的共魏方向P/»--«+«/»»
计算步长—"HP
P T IIP
计算新的楼小戍J=J+"
计«»■!：«的哉敢厶0■氏赋仇J
A.*+!；/./+!
极小点X“
图6-2共辘梯度法N-S图
图6-3不同误差限共辄梯度反演结果图
数时,迭代结果无限逼近于精确解,这也就导致了，误差变得非常小时，反演次数却增加不明显。

通过给不同初值(图6-4),显然共辄梯度法都会收敛到真解附近，反演次数比较少，收敛速度快。

而且迭代次数相差并不大，因此共辄梯度法的初值x_0对最终结果和计算时间影响较小。

理论模型与反演模型拟合程度也非常高。

其所产生的正演异常几乎一致。

这说明了共辄梯度法反演精度高。

误差：e=10-B m.(a)初值:h=1m;迭代次数l：37(b)初值：h= 100m;迭代次数l：37
图6-4
不同初值共辄梯度反演结果图
64石油工业机应用
2020年10月
显然，图6-5表明当反演的平均深度与真实的 平均深度差异较大时,反演结果与理论模型差异非 常大,实际上，当反演的平均深度与理论模型不一致
时,如反演的平均深度比理论模型大时，此时反演的 结果深部异常与浅部异常一致,那反演所得到的模
型必然与理论模型不一致。

图6-5不同反演深度共觇梯度反演结果图
图6-6是密度界面起伏较大时的反演结果，显 然，当界面起伏剧烈时，该方法的反演结果误差较大, 尤其是在极值处与边部误差较大，边部误差主要是由
于在模型两端延伸处观测数据的的缘故引起的。

图6-6剧烈密度界面起伏反演结果图7结论与建议
通过利用共觇梯度法反演地下密度界面起伏，
得出以下几点结论；
(1) 初值对共辄梯度法反演地下密度界面起伏
影响较小，尤其是单极值密度界面，初值对反演结果
无影响。

(2) 反演的平均深度对结果影响较大,当反演
的平均深度远偏离真实平均深度时，反演结果与真 实模型差异非常大。

(3) 地下密度界面起伏剧烈时,利用共辄梯度
法反演的模型，误差较大。

通过研究共辄梯度法反演地下密度界面起伏，
对以后的工作提出以下几点建议:
(1) 在反演过程中，建立目标函数时，尽量加入
约束条件，反演结果将更接近真实的地下结构。

(2) 二度体密度界面反映地下结构比较片面，
以后应着力于研究三维地下密度界面起伏反演。

参考文献：
[1] 曾华霖，重力场与重力勘探.北京:地质出版社.2005.[2] 刘元龙、王谦身，用压缩质面法反演重力资料以估算地壳
构造,地球物理学报,1977,20(1): 59-69.
[3] Oldenburg The Inversion and Interpretation of Gravity Anomalies.Geophysics. 1974, 39:526—536.
[4] P aker R L.The Rapid Calculation of Potential Anomalies.
Geophys.1973,31 ：447—455.
[5] B ott M H P.The use of Rapid Digital Computing Methods for Direct Gravity Interpretation of Sedimentary Basins.
Geophys. 1960,3:63—67.
⑹袁亚湘、孙文瑜,最优化理论与方法，北京：科学出版社,
1997.
[7]姚姚,地球物理反演基本理论与方法应用，北京：中国地质
大学出版社.2002.
•最新测井专利•(外国专利)
专利名称：METHODS TO SYNCHRONIZE SIGNALS AMONG ANTENNAS WITH DIFFERENT CLOCK SYSTEMS
专利申请号：US201715779461 公开粤：US2019376383(Al) 申请日：2017.06.27 公开日：2019.12.12
申请人:Halliburton Energy Services, Inc.
A method for synchronizing signals among transmitters and receivers of a logging tool positioned in a borehole is provided. Measurement signals generated from operating a transmitter in the borehole are acquired by a receiver. An operating frequency of the receiver is determined by a processing unit. The operating frequency of the receiver is different from an operating frequency of the transmitter. A sampling frequency of the receiver is determined based on the determined operating frequency. A phase delay of the receiver is determined by the processing unit. The acquired measurement signals are adjusted by the processing unit based on the determined sampling frequency and the phase delay of the
receiver.。