ω循环型边界条件
第４６卷第６期
２００３年１１月
地球物理学报
ＣＨＩＮＥＳＥ
Ｊ０ｕＲＮＡＬ
０Ｆ
Ｖ０１．４６．Ｎｏ．６
Ｎｏｖ．．２００３
ＧＥＯＰＨＹＳＩＣＳ
＝———＝＝＝———＝＝＝；——＝＝＝——＝—＝—————＝＝————＝＝＝＝——＝＝＝一
国循环型边界条件
梅金顺
刘
洪
１０００２９
中国科学院地质与地球物理研究所，北京
摘要
通过利用预条件共轭梯度法对对称正定’ｒｏｅｐｌｉｔｚ矩阵系统进行分析，重点介绍了一种新的嵌入式预条件
矩阵构造方法，证明了以前的预条件矩阵构造方法大都是这种方法的特例．提出了ｍ循环型边界条件，并将其与普通循环型及螺旋型边界条件作了分析、比较后得到了一种新的边界条件即混合型边界条件．
关键词
Ｔｏｅｐｌｉｔｚ矩阵
预条件共轭梯度法（ＰｃＧ）。循环边界条件
收稿日期２００２—０４—０８，２００３—０６—０６收修定稿
文章编号０００ｌ一５７３３（２００３）０６—０８３５—０ｒ７中图分类号Ｐ６３１
∞－ＣＩＲＣＵＬＡＮＴＢｏＵＮＤＡＲＹＣｏＮＤＩＴＩｏＮ
ＭＥＩ
ＪＩＮｓＨｕＮ
ＬＩＵＨｏＮＧ
胁ｍ眦ｅ
Ｄ，＆Ｄ蛐Ⅱ，以＆ｏｐ伽如，∞ｉ心ｅＡｃ础啊旷＆妇，螂，曰ｅ舛ｎｇ１０００２９，伽胍
Ａｂｓｔｒａｃｔ
ｗｅ
ｕｓｅ
ｔｈｅｐｒｅｃｏｎｄｉｔｉｏｎｅｄ
ｗｅ
ｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔ（ＰＣＧ）ｍｅｔｈｏｄ
ａ
ｔｏ
ａＩｌａｌｙｚｅｔｈｅｓｙｍｍｅｔｒｉｃｐｏｓｉｔｉｖｅ－ｄｅｆ－
ｉｎｉｔｅＴｏｅｐｌｉｔｚｏｔｈｅｒｗａｙｓｏｆ
ｓｙｓｔｅｍｓ．Ｔｈｅｎｐｒｅｓｅｎｔ
ｎｅｗ
ｅｍｂｅｄｄｉｎｇｃｏｎｓｔｍｃｔｉｎｇｐｒｅｃｏｎｄｉｔｉｏｎｅｒｗａｙａｎｄｐｍｖｅｔｈａｔｍａｎｙ
ｇｅｎｅｒａｌｌｙｔｈｅ
ｃｏｎｓｔｍｃｔｉｎｇｐｒｅｃｏｎｄｉｔｉｏｎｅｒｓ
ａｒｅ
ｓｐｅｃｉａｌ
ｃａｓｅｓ
ｏｆｔｈｉｓ
ｍｅｔｈｏｄ．Ｗｅｐｍｐｏｓｅｔｈｅ∞一
ｃｉｒｃｕｌａｎｔｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎａｎｄｍａｋｅｉｘ
ｏｎｅ．
ｃｏｍｐａＪｉｓｏｎｗｉｔｈｔｙｐｅｏｆｂｏｕｎｄａｒｙ
ｂｏｔｈｔｈｅｏ耐ｉｎａｒｙｃｉｒｃｕｌａｎｔｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎａｒＩｄｈｅｌ．ｃｏｎｄｉｔｉｏｎ：ｔｈｅｒｎｉｘｅｄｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎ．
Ｆｕｎｈｅ瑚ｏｌ＿｝ｅ，
ｗｅ
ｏｂｔａｉｎ
ａｎｅｗ
Ｋｅｙｗｏｒ山；Ｔｏｅｐｌｉｔｚｍａｔｒｉｃｅｓ，Ｐｒｅｃｏｎｄｉｔｉｏｎｅｄｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔｍｅｔｈｏｄ，∞－ｃｉｒｃｕｌａｎｔ，Ｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎｓ．
泛。３刈．许多地球物理学家提出了很多好的想法来
１
引言
求解Ｔ０ｅｐｌｉｔｚ矩阵，如Ｃｌｅａｒｂｏｕｔ７Ｊ于１９９８年首先提出了螺旋坐标系方法；此后，国内外发表了很多有关螺旋坐标系方法的论文．这些论文大都没有对螺旋坐标系的稳定性及其与其他方法的比较进行分析．
本文从数学和信号分析中求解对称正定Ｔｏ，ｅｐｌｉｔｚ矩阵系统的预条件共轭梯度法这一新方法出发，系统地讨论了该方法在实际应用时与边界条件的结合问题，得到了一种扩展的ｃｕ循环边界条件（∞＝ｌ即为普通的循环边界条件），目的是能够快速、稳定地处理各类边界条件问题．
一个Ⅳ阶‰ｐｌｉｔｚ矩阵…是ＡⅣ＝（ｏ¨）．ｖ。。＝
（口。一ｉ）～。Ⅳ，记为ＡⅣ＝ｒⅣ（ｏ一＿Ⅳ＋１，…，ｏ一１，ｎｏ，口ｌ，
…，ｏⅣ一，）（其中ｎ表示Ⅳ阶‰ｐｌｉｔｚ矩阵）．７ｒｏｅｐｌｉｔｚ
矩阵是信号分析和数据处理中最为常见的一类矩阵．作为Ｔｏｅｐｌｉｔｚ矩阵Ａ。的延伸，若有８Ⅳ一★＝∞ｏ一。（０＜矗＜Ⅳ），则称其为∞循环矩阵ｕ。，记作Ａ。［∞］＝ｃｉｒｃ（ｎｏ，ｏ—ｌ，－．－，ｏ一Ⅳ＋１）。；特别地，∞＝±１时所对应的即为循环和斜循环矩阵．
‰ｐｌｉｔｚ矩阵在信号分析和数据处理中应用广
基金项目
国家自然科学基金委和大庆石油管理局联合资助（４９８９４１９０）．
作者简介梅金顺，男，１９６７年生，２００２年毕业于中国科学院地质与地球物理研究所固体地球物理专业，获博士学位，现在中国科学院地质与地
球物理所作博士后研究工作．主要从事地震波成像及相关的应用数学物理方法研究．Ｅ—ｍａｉｌ：ｍｙｆｌｏｗｅｒ６８＠ｓｏｈｕ．ｃｏｍ
８３６
地球物理学报（ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．）
４６卷
矩阵Ａ。中部分靠近主对角线元素组成的或按照一
２叫循环型边界条件
１９８６年，ｓｔｒａｎ∥１和ｏｌｋｉｎ‘９１提出了预条件共轭梯度法，它借助于共轭梯度迭代过程，利用预条件来求解对称正定Ｔ０ｅｐｌｉｔｚ系统．理论和实际计算结果表明旧一”，这种方法具有很快的收敛性．预条件共轭梯度法求解对称正定Ｔｏｅｐｌｉｔｚ系统时预条件矩阵Ｐ的构造方法可分为三种类型：
第一类方法是采用一定的准则直接从Ｔｏｅｐｌｉｔｚ
系数矩阵Ａ得到预条件矩阵Ｐ旧’１２““，不同的准则
定规则由Ａ。生成的对称正定’ｒｏｅｐｌｉｔｚ矩阵，即ｒ。＝ｎ（ｆ一＾，＋ｌ＇＿一，￡一ｌ，ｆｏ，￡ｌ，．一＇ｆ、一１），玑＝ｎ（￡ｌ’．一，￡Ⅳ一ｌ，０，…，０），Ｌ＿ｖ＝ｎ（０，…，Ｏ，￡一Ⅳ＋ｌ，…，￡一１），ＯⅣ为ｊｖ阶零矩阵，Ｅ。为Ⅳ阶单位矩阵Ｉ。，Ⅳ阶次对角线上元素为１，其余均为零的矩阵，。，则有
ｃ＿｜ｌｆ．Ⅳ［∞肘】ｘ肘，Ⅳ（（ｃＪＥ～）＝曰＾，，Ⅳ（∞露～）．
（２）选择层。＝，。，得到
（ｒⅣ＋∞ＬⅣ＋∞＿１【，＾）工Ⅳ：参Ⅳ，
我们可以选择ｎ＝ｎ＋ｃｕｋ＋ｃｕ“玑作为预条件
矩阵，并且有ｎ＝ｎ［叫］＝ｃｉｒｃ（￡。，￡一。＋（ｃＪ。１
ｆ。一。，
推导出不同的预条件矩阵Ｐ；第二类方法是函数拟合法ｎ２““１８。，预条件矩阵Ｐ及其逆Ｐ。１均可表示为一系列∞循环矩阵之和的形式；第三类方法为嵌入
式预条件矩阵构造方法．２．１∞预条件矩阵的构造
…，￡一。＋。＋∞。１￡。）。，也就是说，预条件矩阵ｎ是ｃＵ循环的．容易验证，玳［１］＝ｘ，及Ｐ。［一１］＝Ｋ：．
（３）选择层。＝＾且∞＝±１，由于ｎ为对称正
定矩阵，采用类似的方法可以得到如下两种预条件矩阵
ＱⅣ（±１）＝７～±．，Ⅳ（ＬⅣ＋ｃ，Ⅳ）
Ｋｕ等¨９。利用嵌入式方法构造了４种不同的预
条件矩阵Ｋ。（ｎ＝１，２，３，４），其中ｘ。为循环矩阵，
ｘ：为斜循环矩阵，Ｂ和Ｋ。既不是循环矩阵，也不是斜循环矩阵，但它们可以利用快速余弦和正弦变换来计算．ｘ。（ｎ＝１，２，３，４）具有一定的普遍性，很多第一类方法构造出的预条件矩阵都是它们的特
例．
容易验证，ＱⅣ（１）＝置，，Ｑ。（一１）＝Ｋ。．
从以上的分析我们不难看出，Ｋｕ等ｎ刊引入的４种不同的预条件矩阵眉。（凡＝ｌ，２，３，４）只是ｎ［ｃＵ］
和纵（土１）的４个特例．可以认为，利用这种嵌入
本文提出了一种更为普通的构造Ｔｏｅｐｌｉｔｚ矩阵
式方法构造出的预条件矩阵Ｐ。［∞］和Ｑ。（±１）更具有普遍性，其余各种方法构造出的预条件矩阵均为它们的特例．
叫
系统Ａ∥。＝６。预条件的方法：
（１）首先构造如下的分块矩阵
～～
Ⅳ
Ｕ
２．２∞循环型边界条件
ＴｏｅＤｌｉｔｚ矩阵系统属于褶积类运算ｈ０。．多数情况下需要分析的有效信号大都是能量有限的，在边
ｃ肝．Ⅳ［∞］＝
～
Ⅳ
ｎ“叽～“ｎ“～
吣
队“ｎ～虬
ｘＮ
１％叽～
．．．～
一
界附近信号的数值趋近于零，即吸收型边界条件．另一方面，如果采用环形网格进行观测（见图１）（方格网可以认为是环形网的一个特例，当环形网格的半径趋于无穷大时其周边部分网格（Ｆｉｇ．１ｃ）即
幽
Ｒ
，Ⅳ
∞Ｅｄ
Ｎ
构成方格网（Ｆｉｇ．１ｄ）），则会产生循环或部分循环型边界条件：单环环形网格（Ｆｉ异．１ａ）所对应的褶积运算生成的即为循环型边界条件；多环环形网格（ｎｇ．１ｂ）所对应的褶积运算生成的即为部分循环
边界条件．
Ｘ＾，．＾，（叫ＥⅣ）＝（∞Ｅ＿ｖ）２工Ⅳ
（∞ＥⅣ）肛１工＿Ⅳ
６Ⅳ
ｍＥ两Ｎ
曰＾ｆ．Ⅳ（∞ＥⅣ）＝
（∞层Ⅳ）２６Ⅳ
方格网所对应的一般是吸收类型边界条件；将
一个‰ｐｌｉｔｚ矩阵延拓成一个循环矩阵，其对应的吸
收型边界条件也就变成了循环型边界条件．其目的是引入Ｆｏｕｒｉｅｒ变换方法对原来的Ｔｏｅｐｌｉｔｚ矩阵系统作近似的、快速的分析，而原来的吸收型边界条件又能保证这种近似性具有较好的可靠性．
（ｃＵＥ。）肛１
６，
其中』ｌｆ≥３，叫＝ｅｘｐ（ｉ８）（８∈［０，２７ｃ）），ｎ是由系数
６期
梅金顺等：ｍ循环型边界条件８３７
题珍．篇鼢．黼
螺◇‘蹩渺’
Ｆｉｇ．１
雪亘
（ｃ）
（ｄ）
８ｙｓｔｅｍｓ
图ｌ环形网格观测系统．
（ａ）为单环环形网格（其中＋表示观测点位置）；（ｂ）为多环环形网格；（ｃ）为环形网格半径
趋于无穷大时其周边部分网格；（ｄ）为方格网．
Ｃｙｃｌｉｃｇｒｉｄｄｉｎｇ
采用Ｆ０ｕｒｉｅｒ变换方法对‰ｐｌｉｔｚ矩阵系统作分
析这一过程，包含了部分ＰｃＧ的思想，其对应的预条件矩阵即为扩展后的循环矩阵．但是，扩展后的循环系统有时并不理想：简单循环边界条件中通过延拓原Ｔｏｅｐｌｉｔｚ矩阵为循环矩阵来进行求解，然而得到的循环系统有时却变得非常奇异．
本文希望得到一种更具普遍性的新系统，并且容易对该系统进行简单而又快速地处理．我们采用的办法是引入ｃｕ循环系统，即利用∞循环矩阵替代简单循环矩阵，其对应的边界条件即为叫循环型边界条件：普通的简单循环边界条件只是叫循环型边界条件的一个特例（∞＝１）．３
阵，０。为四阶零矩阵）．
引入螺旋型边界条件得到的扩展后的矩阵为
曰１６＝ｒ１６（０，…，０，一１，０，０，一１，４，一１，０，Ｏ，一ｌ，０，
…，０）．在文献［７］中，利用ＩＪａＤｌａｃｅ算子的对称正定性，借助于自相关函数的谱分解方法，得到下三角矩
ｌ浑日１６＝Ｌ】６（１．８，一０．６，０，一０．２，一０．６，０，…，０），
且有日。。日‘ｒ６＝Ｂ，。．这样，求解二维ｈｐｌａｃｅ算子的
逆就变成了求解下三角形矩阵日。。的逆，原Ｔｏｅｐｌｉｔｚ矩阵系统的求解就变得很简单了．而正是由于原Ｔ０ｅｐｌｉｔｚ矩阵系统的对称正定性，保证了最终结果的有效性旧１’２２１（一种具有最小二乘意义的解）．
采用循环边界条件时情况又如何？采用分块循环的延拓方法，可以得到一个分块循环矩阵
Ｌ［口］
一，ａ
一，ｔ
ｃＵ循环型边界条件与螺旋边界条件的比较
引入螺旋坐标系方法ｎ１可以有效地进行数字信
仉
一，ａ
Ｃ。。［ａ，卢］＝
仉
一Ｊｔ
４
Ｌ［ａ］
Ｌ［口］
一Ｊ。
号处理．经螺旋坐标系处理后得到的数据体具有螺旋型边界条件，它是将上下或左右相邻两点进行错位对接．螺旋边界条件类似循环边界条件，只不过在循环边界条件中采用的是平行对接而非错位对接．为了便于进行两类边界条件的对比分析，本研究以４×４方格网上二维【ａｐｌａｃｅ算子
Ａ。一ｊ。
Ｖ２＝
一，。Ａ４一，。
一母ｌ
仉
∥仉“小
其中ｌ
ａ
Ｉ＝Ｉ
ｐ
４
Ｉ＝１，
—１４—１０
０一ｌ４—１
一ｄ一１
Ｌ［ａ］＝
—１０
一口
０—１４
Ｏ。
一，。
也即在该分块循环矩阵中，其内部是口循环，而其
作为分析对象（其
外部则是口循环．
从表面上看，两种边界条件有其相似性，矩阵曰。。和Ｃ。。［ａ，ｐ］都是从原ｈｐｌａｃｅ算子Ｖ２采用延拓的
ｏｄ
Ｏ。
４
Ａ。一ｌ
ａ
ｏｄ
一ｌ４—１０
０—１４—１
饥仉＿Ｉ钆ｏ
ｏ
方法变换而来，只不过具体延拓的方式有所差异．
，Ｊ。为四阶单位矩
采用螺旋边界条件时，其最终目的是得到一个形式
—１
中Ａ。＝
００
简单、阶数较大的‰ｐｌｉｔｚ矩阵曰。。；而采用循环边
界条件时得到的是分块∞循环矩阵，其内部每个分
一４
８３８
地球物理学报（ｃｈｉｎｅｓｅＪ．Ｇｅ叩ｈｙｓ．）
４６卷
块矩阵皆为循环类型的．
０，０，一ｌ，０，…，一（ｕ～，０，Ｏ，一ｃｕ‘１）。．很显然，可以
经过螺旋边界条件延拓后得到的对称正定Ｔｏ一印ｌｉｔｚ矩阵曰，。，采用谱分解的方法将其表示为一个下三角Ｔｏｅｐｌｉｔｚ矩阵日。。与它的转置的乘积．所有计算都是在空间（或时间）域内完成的．但是，由于这种方法是以函数的自相关为基础，因而只能应用于
方便地引入Ｆｏｕｒｉｅｒ变换理论对Ｄ。。［ｃｕ］进行处理和进一步的分析．５
实例分析
此处所有计算都是在ＩＢＭ卟ｉｎｋＰａｄ６００ｘ上完
对称正定型‰ｐｌｉｔｚ矩阵系统，这是它的局限所在；
并且由于扩展后的Ｔｂｅｐｌｉｔｚ矩阵阶数较大，该矩阵的条件数往往就会变得很坏（系统变得不稳定）．
经过循环边界条件延拓后得到的正定循环型矩阵Ｃ，。［乜，卢］，可以利用Ｆｏｕｒｉｅｒ变换理论对其进行处理，所有计算都是在频率域中完成的；并且由于ＦＦＴｒ的存在，使得能够对延拓后的任何数据体进行快速计算．最重要的一点是，在频率域中，可以很方便地对数据体进行各种有效的附加运算（如对数据体进行滤波等），这是空间（或时间）域中任何方法都无法
成的，所有图形都是用Ｍａｔｂａｌ软件包绘制的．实例ｌ：对于二维Ｌａｐｌａｃｅ算子Ｖ２＝ｎ（…，０，一ｌ，０，
０，０，０，０，０，一１，４，一１，Ｏ，０，０，０，Ｏ，０，一ｌ，０，…），我
们以Ｖ２ｘ．ｖ＝（１，ｌ，…，１，１）：作为分析对象并且选择
ｎ［∞］＝ｎ＋∞Ｌ，＋ｃｕ一玑作为预条件矩阵．本文
选择三种形式的ｎ，即ｎ＝４，。、Ｔｄｄｉａｇ（一ｌ，４，一１）和铲，构造出预条件矩阵Ｐ，［∞］．最终迭代次数是以满足｜｜ｎＩＩ：川，０Ｉｌ：＜１０’７为条件的，其中＾为第露次迭代后的剩余向量．
对于对称正定Ｔｏｅｐｌｉｔｚ矩阵系统Ｖ２ｘ。＝（１，１，…，１，１）‘：来说，采用ｎ＝４，。、Ｔ耐ｉａｇ（一１，４，一１）
比拟的．因此，采用循环边界条件对‰ｐｌｉｔｚ矩阵系
统进行处理时具有更大的灵活性和可操作性．
４混合型边界条件
进一步地，我们引入了混合型边界条件．本研究仍以４×４方格网上二维ＩＪａｐｌａｃｅ算子Ｖ２作为分析对象，此时混合型边界条件过程如下：首先对原对称正定Ｔｏｅｐｌｉｔｚ矩阵系统进行螺旋边界处理得到一个高阶Ｔｏｅｐｌｉｔｚ矩阵曰。。；其次再对口，。作循环边界处理得到一个新的ｃＵ循环矩阵Ｄ。。［叫］＝ｃｉｒｃ（４，一１，
表２
ＴａｂＩｅ
和铲时，其对应的预条件矩阵ｎ［∞］的条件数及满足｜｜＾Ｉｌ：川ｒ。｜｜：＜１０‘７的ＰｃＧ最终迭代次数等结果分别列于表１、２和３．ＰＣＧ的收敛速度很快，我们
表１
Ｔａｂｌｅ１
ｎ＝４“时ＰＣＧ的迭代次数
ＰＣＧ
Ｉｔ哪６ｖｅｎ岫ｂｅ璐ｗｊ恤ｎ＝４厶
注：此时预条件矩阵的条件数等于１
ｎ＝７蹦ｄｉａｇ（一ｌ，４。一１）时ｎ［∞】的条件数与ＰＣＧ的迭代次数
ＰｃＧ
６４
２
ｃｏｎｍｔｉ蚰ｎ岫ｂｅｒｓ粕ｄ
３２
Ｉｔｅｒｓ１２１５１５１５１５１５１２１５１５１５１５１５
ｉｔｅｎ廿ｖｅ
ｎ唧ｂｅ璐ｏｆ血ｅ
１２８
ｐｒ∞伽ｍ廿。眦ｒＰⅣ［∞］ｗｉｔｌｌｎ＝％ｄｉａｇ（一１．４。一１）
２５６
Ｉｔｅｒｓ２５２８２９２９２９２９２５２９２９２９２９２８
ＣｏｎｄＮ３．０００３．０００３．０００３．０００３．０００３．０００３．０００３．０００３．Ｏｏｏ３．０００３．０００３．０００
Ｉｔｅｒｓ３４４４４５４５４５４４３４４４４５４５４５
ＣｏｎｄＮ３．０００３．０００３．Ｏｏｏ３．ｏ（）０３．０００３．０００３．０００３．０００３．０００３．０００３．０００３．０００
５１２
Ｉｔｅ璐５０５０５０５０５０５０４９５０５０５０５０５０
日
ＣｏｎｄＮ
００３０。６００９００１２０。１５００１８００２１００２４０。２７００３０００３３００
３．０００２．９９９２．９９８２．９９５２．９９１２．９８７２．９８１２．９８７２．９９ｌ２．９９５２．９９８２．９９９
ＣｏｎｄＮ３．０００３．０００２．９９９２．９９９２．９９８２．９９７２．９９５２．９９７２．９９８２．９９９２．９９９３．Ｏｏｏ
Ｉｔｅｒｓ１８２０２０２０２０２０１７２０２０２０２０２０
ＣｏｎｄＮ３．０００３．０００３．ｏｏＯ３．０００２．９９９２．９９９２．９９９２．９９９２．９９９３．ｏ（）０３．０００３．０００
“
注：ｃｏｎｄＮ表示条件数，ＩｔｅＨ表示迭代次数，ｍ＝ｅｘｐ（ｉ日）
６期
梅金顺等：ｍ循环型边界条件
８３９
以Ⅳ＝３２时ｒ，：＝Ｖ２所对应的预条件矩阵Ｐ盟［∞］时该方法前１０次迭代后的剩余向量值列于表４，从表４很容易看出，迭代５、６次后剩余向量值就已经变得很小了．
很显然，当我们取ｎ＝铲且ｐ＝００时其对应的预条件矩阵ＰⅣ［１］即为普通Ｆｏｕｒｉｅｒ变换方法中所采用的系数矩阵；也即是说，口＝０。时对应的是普通循环边界条件．从表３、４中不难看出，此时ＰｃＧ根
件的必要性和重要性．
实际工作中，我们完全没有必要采用满足Ｉｌ＾ＩＩ：川ｒ。Ｉ｜：＜１０。为条件作为最终迭代次数．在做信号分析和数据处理时，由于实际数据体本身就包含有一定的误差，我们可以采用不同的标准终止迭代过程．表５显示的是Ⅳ＝３２且迭代过程满足
１１＾｜ｌ２川，ｏ｜｜２＜１０～（ｎ＝１，２，３，４，５，６，７）时ＰＣＧ
的最终迭代次数列表．从该表不难看出，只需有限
次迭代，ＰｃＧ即可得到很好的精度．本不收敛．这进一步地说明了引入叫循环型边界条
裹３ｎ＝铲时＾［国］的条件数与ＰＣＧ的迭代次数
ＴａＭｅ３
口
ＣｏｎｄＮ
００３０。６０。９００１２０。１５００１８０。２１０。
ｉｎｆ４５ｌ１１３５０２８１８１２１８２８５０１１３４５１
３２
Ｉｔｅ硌
＊
ｃ彻棚Ｈ帆ｎ咖ｂｅ璐ａｎｄ
“
ＣｏｎｄＮｉＩｌｆ１８０５４５ｌ２０ｌ１１３７２５０７２１１３２０ｌ４５１１８０５
ＰｃＧ溉ｎｄｖｅ咖ｂｅ璐ｏｆｕ地ｐｒｅｃ佃ｄｎｉｏ毗ｒＰⅣ［∞］们ｔｈ
１２８
Ｉｔｅ硌
＊
ｎ＝Ｖ２
５１２
２５６
Ｉｔｅｒｓ
＊
ＣｏｎｄＮｉＩｌｆ７２１８１８０５８０２４５ｌ２８９２０１２８９４５１８０２１８０５７２１８
ＣｏｎｄＮｉＩｌｆ２９０００７２１８３２０８１８０５１１５５８０２１１５５１８０５３２０８７２１８２９０００
Ｉｔｅ鸺
＊
ＣｏｎｄＮｉｒＩｆ１．２２９０００１３０００７２１８４６２０３２０８４６２０７２１８１３０（）ｏ２９０００１２００００
ｌｔｅ璐
＊
１３１３１２１２１２９１２１２１２１３１３
１５１４１４１４１２９１３１４１４１４１５
１５１４１４１４１４１０１４１４１４１５１６
＊＊
１５１４１４１４１０１４１５１５１５
＊
＊
１５１５１５１０１５１５１５
＊
擀
２７００３００。３３００
ｏ
注：ｃｏｎｄＮ表示条件数，Ｉｔｅｒｓ表示迭代次数，ｍ＝“ｐ（ｉ日），ｉＩｌｆ代表无穷大，＊表示不收敛
表４
ＴａＭｅ４
ｒ，：＝Ｖ３且预条件为ｎ［∞］时ＰＣＧ的前１０次迭代后的剩余向量模值列表
Ｎｏ珊ｍⅧｌｕ鹤ｏｆ也ｅ他湖眦Ｉｖｅｃｔｏ璐∞ｒｒ幅舯哪ｌｄｉ雌ｔｏ
ｗｉｔｈ
ｒ３２＝Ｖ２ａｎｄ
ｍｅｎｒｓｔ１０ｉｔｅ聃ｄ帆ｓｂｙＰＣＧ
ＰⅣ［∞】ａｓ也ｅｐｌ船ｃ咖ｍ廿伽ｅｒ
注：ｍ＝ｅｘｐ（ｉ∞，Ｎ∞表不结果过大、溢出．
实例２：我们采用一个１６×１６方格网上二维ＬａｐｌａｃｅＤ。。［∞］．本文利用Ｍ砒ａｌ软件包中的现成软件ｐｅａｋｓ
产生了三组数据体：用ｐｅａｌ【ｓ（１６）产生一个１６×１６三维规则网格数据体ｚ。，将ｚ，周边网格上的（６０个点）数值修改为零得到另一１６×１６三维规则网格数据体ｚ：；另外，用ｐｅａｌ【ｓ（１４）产生一个１４×１４三维规则网
算子Ｖ２作为分析对象，其对应的对称正定‰ｐｌｉｔｚ矩
阵系统为铲ｚ＝Ｒ，目的是对三种边界条件（循环型、螺旋型及混合型）所处理的结果进行比较，三种边界条件所对应的系数矩阵分别为ｃ。。［ａ，ｐ］、曰，。及
地球物理学报（ＣｈｉｎｅｓｅＪ．Ｇｅｏｐｈｙｓ．）
４６卷
表５采用三种边界条件对三组数据体Ｚ。、ｚ：和Ｚ，
进行反演计算的均方差误差结果
Ｔａｂｌｅ
５
Ｒ。、四。。ｚ。，：＝Ｒ。及Ｄ。６［∞］Ｚ¨＝Ｒ。得到ｚ¨、ｚ㈦和ｚ。．，；对数据体ｚ：和ｚ，的计算过程完全相同．此处与ＰＣＧ方法不同的是，所有计算结果都是一次性的（没有使用迭代方法）．对三组数据体的计算结果分别显示于图２～４及表５，其中表５的均方差计
ｅｒｒｏｒｓ
ｉｎｖｅ嗽Ｉｒｏｏｔ－ｍｅａｎ－ｓｑｕａｒｅ（ＲＭＳ）
ｏｆ妇ｄａｔａＺＩ、Ｚ２８ｎｄＺ３ｂｙｔｈｅ瞅ｏｆｔｌｌｒｅｅｍｆｆｅｎｎｔｂ０恤～ｄａｒｙｃ伽脚廿吣
Ｔｈｅ
厂１ｒ—————————～——
１６
算方法为
格数据体，然后将其周边进行扩展（用零进行充填），得到第三组１６×１６三维规则网格数据体ｚ，．
对于数据体ｚ。而言，我们先正演计算出Ｖ２ｚ。＝Ｒ，得到Ｒ。，然后分别反演计算ｃ。。［ａ，卢］ｚ¨＝
●５
√磊（ｚ…Ｑ∽。撕∥
１，２，３）．
从计算结果不难看出，采用ｃｕ循环型边界条件时反演结果误差最小；其次是混合型边界条件；螺旋型边界条件所对应的反演结果误差最大，最不稳定．
¨Ｏｘ５
，
８
Ｏ
蓦一
０
５
Ｏ
＠。多‘
５
，
徊、‘
＼
ｏ
／
１０１５０５
，
１０１５
图３对数据体ｚ：进行反演计算的结果比较
（ａ）为原始数据体ｚＩ；（ｂ）、（ｃ）、（ｄ）分别为循环型（ｃ。６［一ｌ，一１］）、
螺旋型（口。。）及混合型（Ｄ。。［一１］）边界条件的反演
计算结果ｚ¨、ｚ１．２和ｚ１．３．
Ｆｉｇ．２
（ａ）为原始数据体ｚ２；（ｂ）、（ｃ）、（ｄ）分别为循环型（ｃ。６［一ｌ，一１］）、
螺旋型（口。。）及混合型（Ｄ。６［一１］）边界条件的反演
计算结果ｚ２．１、ｚ２。２和ｚ２．３．
Ｉｎｖｅｒｓａｌｒｅｓｕｌｔｓｏｆｔｈｅｄａｔａ
ＺｌＦｉｇ．３
Ｉｎｖｅｒｓａｌｒｅｓｕｌｔｓｏ“ｈｅｄａｔａ
Ｚ２
６
总结
本文通过引入一种新的嵌入式预条件矩阵构造
方法对对称正定‰ｐｌｉｔｚ矩阵系统作分析，得到了两
类不同的预条件矩阵ｎ［∞］和Ｑ。（±１）；证明了其他几种预条件矩阵是这两类矩阵的特例．进一步导出并分析了∞循环型边界条件，讨论了其与螺旋边界条件的异同点，比较了它们在计算方法上各自的特点：叫循环型边界条件可以在频率域中实现，该方法在反演计算的同时亦可对数据体进行其他附加处
图４对数据体ｚ，进行反演计算的结果比较
（ａ）为原始数据体ｚ３；（ｂ）、（ｃ）、（ｄ）分别为循环型（ｃ１６［一ｌ，一１］）、
螺旋型（口，。）及混合型（Ｄ。。［一１］）边界条件的反演
计算结果ｚ３．１、ｚ３．２和ｚ３，３
理（如滤波等）；螺旋边界条件只能在空间域内进行．最后提出了一种混合型边界条件，它可在频率域中实现，是ｃｕ循环型和螺旋型两个边界条件的一种折衷．数据计算验证了它们的可行性．实际计算表明，采用∞循环型边界条件优于混合型边界条件，
Ｆｉｇ．４
ｌｎｖｅｒｓ８ｌｒｅｓｌ】１ｔｓｏｆｍｅｄａｔａＺ３
６期
梅金顺等：ｃｔ，循环型边界条件
８４ｌ
而混合型边界条件又优于螺旋边界条件参考文献
［１］
徐仲，张凯院，陆全．‰ｐ１汜矩阵类的快速算法．西安：
西北工业大学出版社，１９９９
ＸＵＺｈｏｎｇ，ＺＨＡＮＧ
Ｋａｉｙｕａｎ，叫Ｑｕ锄．ＱｕｉｃｋＡｌｇｏ打山ｍ０ｆｔｌｌｅＴ０－
ｅｐｌｉｔｚ
Ｍ栅ｃｅｓ．Ｘｉ’ａｎ：ＮｏｎｌｌｗｅｓｔＩｎｄｌｌｓ田ｕＩｌｉｖｅｒｓｉ‘ｙ
Ｐｒｅｓｓ，１９９９
［２］
蒋增荣，曾泳泓，余品能．快速算法．长沙：国防科技大学出
版社，１９９３ＪＩＡＮＧ
ｚｅ“ｇｒｏｒＩｇ，ＺＥＮＧ
Ｙｏｎ曲ｏｎｇ，ＹｕＰｉｎｎｅ“ｇ．Ｆａｓｔ
Ａｌ耐ｔｌｌｍ．
Ｃｈａｎ挚ｈａ：Ｎ撕ｏｎａｌ
Ｄｅｆｅｎｃｅｓｃｉｅｎｔｉ６ｃ
ＵｎｉＶｅ玲ｉｔｙＰｒｅｓｓ，１９９３
【３］
Ｂｕｎｃｈ
Ｊ．蛐ｉｌｉｔｙ
ｏｆｍ乩ｏｄｓｆｏｒ“ｖｉｎｇＴｏｅｐｌｉ乜ｓｙｓｔｅ鹏ｏｆ８ｑｕ撕。啮．
肋ＷＪ．尉．＆ｍ．Ｃｏ胛山．，１９８５，６：３４９～３“
【４］
Ｇ．ｍｉ啪Ｄ，Ｍ枷ｎ
Ｃ，ＬｕｎｄＪ．Ａｎａｌ徊ｃ蛐ｄｎｕｍ商ｃａｌ鲫ｐｅｃｔｓ
０ｆｔｌｌｅ
ｏｂｓｅｒｖ“０ｎ
ｏｆｔ｝Ｉｅ
ｈｅａｔ
ｅｑＩｌ撕ｏｎ．Ｐｒｏｃ．２６ｎ腰雎ｃｏ旷．０ｎＤ即ｂ如ｎ
ｏｎｄｃ０眦ｒ甜．１９８７。２：９７５～９７６
［５］
Ｈ８ｙｋｉｎＳ．Ａｄａ叫ｖｅ
ＦｉＩｔｅｒ
１１ｌｅｏｒｙ（２一ｅｄ）．ＮｅｗＪｅＩｓｅｙ：Ｐｒｅｒｎｉｃｅ
Ｈａｌｌ，Ｅｎｄｅｗ００ｄＣｌｉ雎，１９９１．１６７～１６８
［６］
硒ｎｇ
Ｒ，Ａｈ咖ｄｉＭ，Ｇｏｑ卿一ＮａｇＩＩｉｂ
Ｒ，ｅｔ
ａＩ．Ｄｉｇｉ“Ｆｄｔ谢“ｇ
ｉｎ０ｎｅ
舳ｄ个ｗｏＤｉｍｅｎｓｉｏｎｓ：Ｄ惜ｉ即ａＩｌｄＡｐｐｌｉｃ嘶。璐．Ｎｅｗ
Ｙｏｒｋ：ｐｌｅｎ啪
Ｐｒ∞ｓ．１９８９．３６８—３７９
［７］
ｃｌｅｍｂｏｕｔ
Ｊ．Ｍｕｌｄ・ｄｉ眦肥ｉｏｎａｌｒｅｃ岫ｉｖｅ
６ｌｔｅ璐ｖｉａ
ａ
ｌｌｅｌｉｘ．＆叩协－
ｊ［ｃｓ，１９９８，６３（５）：１５３２～１５４１
［８］
ｓｔｍｎｇ
Ｇ．Ａ
ｐ叩ａｌ
ｆ曲Ｔ０ｅｐｌ池ｍ８ｔｒｉｘｃａｌｃｕｌ血。璐．＆［１ｄ．ＡｐＰｆ．
Ｍｍｈ，．１９８６，７４：１７１—１７６
［９］Ｏｌｌ【ｉｎｊ．“ｎｅ盯锄ｄＮｏ“ｎｅ盯Ｄｅｃｏｎｖｏｌ“ｏｎＰｒｏｂｌｅ雌［ＰＨ．Ｄ．ｔｌｌｅ—
ｓｉｓｊ．Ｈｏｕｓｔｏｎ：彤ｃｅ
ＵＩｌｉｖ．，１９８６［１０］
Ａｘｅｌ蟠ｏｎ
Ｏ，Ｂ斌ｅｒＶ．Ｆｉ耐ｔｅ
Ｅｌｅ眦ｎｔ鼬ｌｎｉ叩０ｆ
Ｂ叫ｎｄａｒｙ
Ｖａｌｕｅ
ＰｒｏｂｌｅｌＩｌｓ，Ｔｈｅｏｒｙ肌ｄ
Ｃｏｍｐｕ喇ｏｎ．ｏｄ舯ｄｏ：Ａｃ８ｄｅＩｌｌｉｃ
Ｐｒｅｓｓ，１９８４
［１１］
ｓｔ啪ｇ
Ｇ，ＥｄｅｌＩｎａｎ
Ａ．Ｔｈｅ
Ｔ０ｅＰｌｉｔｚ－ｃｉｒｃＩｌｌ蚰ｔ
ＥｉｇｅｎｖａｌｕｅＰＩＤｂｌｅｍＡｘ
＝ＡＣｘ．Ｉｎ：ＬＢｍｇｇ，ＪＤｅｔｔｍａｎ，ｅｄｓ０瑚帅ｄ
Ｃｏｎｆ．ＯｎＰＤＥａｎ（１
Ａｐｐｌ．Ｍａ山．．Ｎｅｗ
Ｙｏｒｋ：ｈｎｇｍ肌ｓｃｉＴｅｃｈ．．１９８７
［１２］
ＣｈａｎＲ．ＣｉｒｃＩｌＩａｎｔｐｒｅｃｏｎｄｉⅡｏｎｅ璐ｆｏｒｈ－ｍｄｄａｎ
Ｔ０ｅｐｌｉｔｚｓｙｓｔｅＨｌｓ．ｓ，一
Ａ＿】ｌｆ，．＾‰血Ａ删．ＡｎＰｌ．，１９８９，１０：ｓ４２～５５０
［１３］
ＣｈａｎＴ．Ａｎｏｐｄｍａｌ
ｃｉｒｃＩｌＩ帅ｔ
ｐｒｅｃｏｎｄｉ“ｏｎｅｒｆ矗Ｔｏｅｐｌｉｔｚｓｙｓｔｅｎｌｓ．Ｓ，一
ＡＭ
Ｊ．＆ｉ．＆吡．（’０，印吡．，１９８８，９：７６６—７７ｌ
［１４］
啊ｙｓｈＩｌｉｋｏｖＥ．Ｏ砸ｍ８ｌ
ａｎｄｓｕｐｅｒｏｐｔｉｌＴｌａＩｃｉｒｃｕｌａｎｔｐｒｅｃｏｎｄｉ６０ｎｅｒｓ．
Ｓ玉ｔｊｌｆ
Ｊ．ＭⅡ打讧Ａ删．ＡＰＰＺ．，１９９２，１３：４５９～４７３
［１５］ｃｈａｎＲ，Ｎｇ
Ｍ．Ｔ０ｅｐｌｉｔｚ
ｐｒｅｃｏｎｄｉｄｏｎＰｒｓ
ｆｏｒＨｅ阳血ｔｉａｎ
７ｒｏｅｐｌ缸８ｙｓ・
ｔｅｒｎｓ．丘，ｌｅ甜Ａ咖６ｍＡ芦叫．，１９９３，ｌ如：１８１—２０８
［１６］
ｃｈａｎ
Ｒ，ｓ咖ｇＧ．１ｋｐ】池ｅｑｕ８“ｏｆｌｓ
ｂｙ
ｃｏｎｊ“９８ｔｅｇｒａｄｉｅｎ￡ｓｗｉｃｈｃ＊
ｃｕｌ蚰ｔ
ｐｒｅｃｏｎｄｉⅡｏｎｅｒｓ．Ｓ＾¨ｆＪ．蹦．ｓｃ眦．Ｃｏ唧眦．，１９８９，１０：１０４一
１１９
［１７］
ＣｈａｎＲ，ＹｅＩｌｆＩｇＭ．ＣｉｒｃＩｌｌａｎｔ
ｐｒｅｃ佣ｄｉｔｉｏｎｅ飓ｃｏｎｓ帅ｃｔｅｄ如ｍ
Ｋｅｒ．
ｎｅｌｓ．ｓ，ＡＭ
Ｊ．肌榭ｒ．Ａ删．。１９９２，２９（４）：１０９３一１１０３
［１８］
蹦ｅｄｌａｎｄｅｒＢ，ＭｏｒｆＭ，ＫａｉｌａｔｌｌＴ，ｅｔ
ａ１．Ｎｅｗｉｎｖｅ璐ｉｏｎ
ｆｏ咖Ｌｌｌ北ｆ０’
瑚ｍｃ鹧ｄ嘲ｉｎｅｄ
ｉｎ
ｔｅｍ毽０ｆ
ｔｈｅｉｒ
ｄｉｓＬ蛐ｃｅ
ｈｍ‰ｐｌ汜ｍａｔＩｉｃｅｓ．
厶ｗ口ｒ
Ａ妇６ｍ却一．，１９７９，２７：３ｌ一６０
［１９］
ＫｕＴ，Ｋｕｏｃ．Ｄ酷ｉｇｎ
ａｒＩｄ帅ａｌｙｓｉｓ
ｏｆ
Ｔ∞ｐｌｉｔｚ
ｐｒｅｃｏｎｄｉｎｏｎｅ璐．衄Ｅ
ｎ删玑０ｎＳ溏，越Ｐ眦醛，１９９２，４０：１２９～１４ｌ
［２０］程乾生．信号数字处理的数学原理．北京：石油工业出版社，
１９７９
ＣＨＥＮＧ
Ｑｉ咖ｈｅ“ｇ．Ｍ８ｔｌｌｅ咖ｄｃａｌ
Ｐｒｉｎｔ五ｐｌｅｓ
ｏｆＳｉｇｎａｌＤ８ｔａ
Ｐｔｏｃｅｓｓ－
ｉｎｇ．Ｂｅ妇ｉ“ｇ：Ｐｅｔｍｌｅｕｍｌｎｄｕｓｔｒｙ
Ｐｒｅｓｓ，１９７９．
［２１］
徐萃薇．计算方法引论．北京：高等教育出版社，１９９７
ｘｕ
ＣＩｌｉｗｅｉ．Ｉｎｔｒｏｄｕ舐ｏｎｔｏｔｌｌｅ
ｃｏ呷ｕｔａｄｏｎａｌ
Ｍｅｔｌｌｏｄ．Ｂｅ巧ｉｎｇ：Ｈｉ鼍；｝ｌ－
ｅｒ
Ｅｄｕｃ出０ｎａｌ
Ｐｒｅ∞．１９９７
［２２］
ｎｅ８ｓ
Ｗ，Ｔｅｕｋ０１８ｋｙ
Ｓ，Ｖｅｎｅｄｉ“ｇ
Ｗ，ｅｔ
８Ｉ．ＮｕｍｅｒｉｃａｌＲｅｃｉｐ船ｉｎＣ
（２叫Ｅｄｉ虹ｏｎ），Ｌｏｒｌｄｏｎ：ＣａｒｎｈｒｉｄｇｅＵｎｉｖｅ商ｔｙ
Ｐｋ踌，１９９２
ω循环型边界条件
作者：作者单位：刊名：英文刊名：年，卷(期)：被引用次数：
梅金顺， 刘洪
中国科学院地质与地球物理研究所,北京,100029地球物理学报
CHINESE JOURNAL OF GEOPHYSICS2003,46(6)2次
参考文献(22条)
1. 蒋增荣;曾泳泓;余品能 快速算法 1993
2. Ku T;Kuo C Design and analysis ofToeplitz preconditioners[外文期刊] 1992(1)
3. Friedlander B;Morf M;Kailath T New inversion formulae for matrices classified in terms oftheirdistance from Toeplitz matrices[外文期刊] 1979
4. Olkin J Linear and Nonlinear Deconvolution Problems 19865. Strang G A proposal for Toeplitz matrix calculations 1986
6. Clearbout J Multi-dimensional recursive filters via a helix 1998(05)
7. King R;Ahmadi M;Gorgui Naguib R Digital Filtering in One and Two Dimensions:Design andApplications 1989
8. HaykinS Adaptive Filter Theory (2nd ed) 1991
9. Gilliam D;Martin C;Lund J Analytic and numerical aspects of the observation of the heat equation1987
10. Bunch J Stability of methods for solving Toeplitz systems of equations[外文期刊] 198511. Press w;Teukolsky S;Vetterling W Numerical Recipes in C (2nd Edition) 199212. 徐萃薇 计算方法引论 1997
13. 程乾生 信号数字处理的数学原理 1979
14. Chan R;Yeung M Circulant preconditioners constructed from Kernels[外文期刊] 1992(04)
15. Chan R;Strang G Toeplitz equations by conjugate gradients with circulant preconditioners[外文期刊] 1989
16. Chan R;Ng M Toeplitz preconditioners for Hermitian Toeplitz systems[外文期刊] 199317. Tyrtyshnikov E Optimal and superoptimal circulant preconditioners[外文期刊] 199218. Chan T An optimal circulant preconditioner for Toeplitz systems 1988
19. Chan R Circulant preconditioners for hermitian Toeplitz systems[外文期刊] 198920. Strang G;Edelman A The Toeplitz-Circulant Eigenvalue Problem Ax =λ Cx 1987
21. Axelsson O;Barker V Finite Element Solution of Boundary Value Problems,Theory and Computation1984
22. 徐仲;张凯院;陆全 TOEPLITZ矩阵类的快速算法 1999
引证文献(2条)
1. 陈宝书. 梅金顺. 王润秋 褶积与反褶积的一种计算方法[期刊论文]-地球物理学进展 2009(2)2. 梅金顺. 刘洪 预条件方程组及其应用[期刊论文]-地球物理学报 2004(4)
本文链接：http://d.g.wanfangdata.com.cn/Periodical_dqwlxb200306017.aspx