时间：2023-03-05 14:29:20
Abstract： The application of poisson integral in the wave equation is introduced and proved detailedly by simple introduction of poisson integral and wave equation and by processing of two-dimensional sional and three-dimensional wave equation. In the introduction of application, we mainly used some of typical examples to discuss the application of poisson integral in the wave equation by linking theory with practice.
关键词： 泊松积分；波动方程；初值问题；调和函数
Key words： poisson points；Wave Equation；Initial Value Problems；Harmonic Function
中图分类号：G31 文献标识码：A文章编号：1006-4311（2011）11-0220-03
0引言
自17世纪，牛顿，莱布尼兹发明微积分后，科学家在利用微积分处理力学，物理学中各种问题的过程中导出大量的微分方程。在这些微分方程中，有些是常微分方程，比如力学中质点的运动方程m=f，但更多的是偏微分方程。欧拉，拉格朗日等科学家在研究流体力学，声音传播和膜振动等问题时，拉普拉斯在研究势函数和潮汐理论时，傅立叶在研究传导以及麦克斯韦在研究电磁理论时都导出一些偏微分方程，近代量子力学中出现的波动方程也是偏微分方程。我们把物理研究中出现的偏微分方程称为数学物理方程。
本文我们只涉及到两种最常见的偏微分方程―泊送方程和波动方程。19世纪最大的分析家之一，第一流的物理学家S.Poisson吸取Fourier的方法，他不仅从事热的理论研究，而且是弹位数学理论的奠基人之一，他最现将引力位势理论引入静电磁学。Poisson甚至相信所有的偏微分方程都可以用级数展开来求解，从1815年起，Poisson按照三角级数Logendre多项式，laplace曲面调和函数展开式求解了许多热传导问题所提出的泊松积分在三大典型的数学物理方程中就有很广泛的应用，在接下来的文章中我们主要探讨一下泊松积分在波动方程中的应用。
1基础知识
1.1 泊松方程与泊松积分的简介
1.1.1 泊松方程的定义方程-u=f（x，y，z）与方程u=0分别称为泊松方程和调和方程（或拉普拉斯方程），其中u=++，符号=++称为拉普拉斯算子。
特别的，当我们把这些偏微分方程与实际生活中的物理知识联系起来时，三维的泊松方程式非常常见的。如果经过相当长的时间后，区域G内各点的温度随时间的改变所发生的变化已不显著，在数学上可近似看做ut=0（即温度函数u与时间t无关，仅为x，y，z的函数）。这时，我们说温度分布趋于定常，方程可写为
u+u+u=0（1）
方程（1）通常称为三维拉普拉斯方程（laplace方程）。在有热源（与时间无关）而且温度分布定常的情况下，方程可写为
u+u+u+f（x，y，z）=0（2）
它也可以写成
u+u+u=f（x，y，z）
其中
f（x，y，z）=-（x，y，z）=-，F（x，y，z），
是热源强度，通常称方程（2）为三维泊松方程（Poisson方程）。三维的拉普拉斯方程经常写作u=0，这个方程我们刚刚也提到过。类似地我们也可以写出二维的拉普拉斯方程：u+u=0和二维的泊松方程：u+u=f（x，y）。至此，我们就结合物理知识导出了一类典型的数学物理方程。由于拉普拉斯方程和泊松方程的关系密切，在一定的程度上可以相互转化，所以在这里我们同时把这两个方程拿出来对比着介绍一下，一是便于记忆，二是为下面泊松积分在波动方程中的应用的叙述奠定基础。
1.1.2 泊松积分的定义及其相关定理在介绍泊松积分的基本知识时，我们主要给出一个基本定理，这个基本定理掌握了，泊松积分也就随之理解了。
定理1：设B是以A为中心，R为半径的球，B是它的边界，f（p）在B连续，则泊松积分：u（P）=f（P）H（P，P）d，
其中H（P，P）=，
是问题：u=0，（在B内）u=f， （在B上）的解。
这个定理说明了当边值上的函数值f满足一定的条件时，我们所得出的球，半空间等区域上的泊松积分，就是相应的调和方程第一边值问题上的解。我们在这里只是对球的泊松积分，说明了定理的正确性，当然我们也可以验证其他区域上的泊松积分也是一些方程在某些问题的解，在接下来的例题中，我们会对它进行详细的证明，并给出了一些实用的解题方法。
1.2 波动方程的简介
1.2.1 波动方程的定义
（1）形如=a的方程，就是弦的自由振动方程，通常称为弦振动方程或一维波动方程。而在弦受外力作用，即=a+f（x，t），称为弦的强迫振动方程，非齐次弦振动方程或一维非齐次波动方程。
（2）形如=a++f（x，y，t）（其中a为常数，f（x，y，t）是与外界有关的已知函数），这样的方程就称为膜的强迫振动方程。特别的，当f=0时，则得=a+，称为膜的自由振动方程或二维齐次波动方程。
1.2.2 波动方程的初值问题
波动方程的初值问题要讨论的主要是它的求解的过程与方法。对于不同的方程或同一类方程，由于维数的不同，定解条件的不同，它的定解问题的求解方法也是不同的。我们对波动方程的初值问题的求解方法简单总结如下：①无界弦自由振动的初值问题，可以用达朗贝尔方法求解，带入达朗贝尔公式直接求解即可。②三维齐次波动方程的初值问题，用求平均法求解，得到解的表达式即泊松公式，然后将泊松公式带入求解即可。③二维齐次波动方程的初值问题，可以用降维法，先降维再带入泊松公式求解即可。
当然，波动方程初值问题的求解方法远不止这些，在这里我只是把几个比较简单的常用的并且在我们需要讨论的泊松积分在波动方程中的应用涉及到的方法列举出来了。
2泊松积分在波动方程中的应用
在泊松积分在波动方程中的应用这部分，我们主要分两部分内容进行讨论，一是泊松积分在波动方程中的形式及其推广，二是通过几个具体的例题演示一下泊松积分在波动方程中的应用。
2.1 波动方程中的形式及其推广
一维弦振动方程的初值问题
u-au=f（x，t），x∈R′，T>0，U（x，0）=φ（x），u（x，0）=φ（x），x∈R′，
的求解公式是
u（x，t）=[φ（x+at）+φ（x-at）]+φ（ξ）dξ+
f（ξ，r）dξdτ
当f0时，该公式就是一维齐次弦振动方程d′Alembert的公式。
三维波动方程的初值问题
u-au=f（x，t）， x∈R，t>0，u（x，0）=φ（x），u（x，0）=φ（x），x∈R，
的求解公式是
u（x，t）=ds+ds+dv
利用球坐标，该公式又可以写成
u（x，t）=tφ（x+atsinθcosφ，x+atcosθ）sinθdθdφ
+tφ（x+atsinθcosφ，x+atcosθ，x+atcosθ）sinθdθdφ
+fx+rsinθcosφ，x+rsinθsinφ，x+rcosθ，t-rsinθdθdφdr
当f0时，称该解公式为齐次方程的Poisson公式或Kirchhoff公式。
2.2 例题
例1 ：求解下列初值问题
u-a（u+u+u=0）， （x，x，x）∈R，t>0，u =x+xx，u=0， （x，x，x）∈R，
解：直接代入三维波动方程的Poisson公式的球坐标形式得
u（x，t）=t[（x+at sinθcosφ）+（x+atsinθsinφ）
（x+atcosθ）]sinθdθdφ
=[t（4πx+4πatx+4πxx+43πatx）]
=x+3atx+xx+atx
所以该初值问题的解为
u（x，t）=x+3atx+xx+atx。
例2 ：求二维波动方程的初值问题：
u-a（u+u）=0，u =φ（r），u=ψ（Γ），r=，
的轴对称解u=u（r，t）。
解：二维波动方程的Poisson公式的极坐标形式为
u（x，x，t）=ρdθdρ
+ρdθdρ。
当解为轴对称时，
u（x1，x2，t）=u（r，t），
于是u（r，t）=ρdθdρ+ρdθdρ
=ρdθdρ+ρdθdρ
综上，所求的u（x，t）就是我们所要求的二维波动方程的初值问题下的解。
例 3 ：求解下列初值问题：
u-a（u+u）=cu，（x，x）∈R，t>0，u =φ（x，x），u=ψ（x，x）， （x，x）∈R，
其中c为常数。
解：令v（x1，x2，x3，t）=eu（x1，x2，t），则v满足定解问题：
v-a（v+v+v）=0，（x，x，x）∈R，t>0，v eφ（x，x），v=eψ（x，x），（x，x，x）∈R，
由三维波动方程的泊松公式得：
v（x，x，x，t）=te×φ（x+atsinθcos，x+atsinθsin）sinθdθd+te×ψ（x+atsinθcos，x+atsinθsin）sinθdθd 。
因此，
u（x，x，t）=ev（x1，x2，x3，t）
=teφ（x+atsinθcos，x+atsinθsin）sinθdθd
+teψ（x+atsinθcos，x+atsinθsin）sinθdθd 。
例 4：设Ω是正方形x1，x1，u（x，t）是初值问题
u-4u=0， x∈R，t>0，u（x，0）=φ（x），u（x，0）=ψ（x）， x∈R，
的解，其中
φ（x），ψ（x）=0， x∈Ω，t>0，>0， x∈R，
试指出当t>0时，u（x，t）0的区域。
解：直接利用二维齐次方程初值问题的Poisson公式可以看出，当t时，在正方形区域x1-2t，x1-2t内u（x，t）0;当t>时u（x，t）>0，在R上成立。
例 5 ：证明电报方程的初值问题
u=au-aλu，x∈R′，t>0， u（x，0）=0，u（x，0）=ψ（x），x∈R′，
的解是
u（x，t）=J（λ，δ）ψ（ξ）dξ，
其中s=，J是由下式定义的零阶Bessel函数：
J（ξ）=cos（z sinθ）dθ 。
证明：令v（x，y，t）=u（x，t）cosλy，则v满足初值问题
v=a（v+v）， （x，y）∈R，t>0，v（x，y，0）=0，v（x，y，0）=ψ（x）cosλy，（x，y）∈R，
有二维波动方程的Poisson公式得
v（x，y，t）=dy
=ψ（ξ）dξdy，
这里，∑是ξ-η平面上以p（x，y）为圆心，以at为半径的圆面。做变换
η-y=sinθ，
则有
dη
=cosλy+λsinθdθ
=2cosλycosλsinθdθ
=πJ（λs）cosλy 。
于是v（x，y，t）J（λs）ρ（ξ）dξ，
这样就得到原问题的解为u（x，t）=J（λs）ρ（ξ）dξ 。
我们从而就证明了电报方程的初值问题：
u=au-aλu，x∈R′，t>0，u（x，0）=0，u（x，0）=ψ（x），x∈R′，
的解是u（x，t）=J（λ，δ）ψ（ξ）dξ。
通过以上几个例题，我们对于泊松积分在波动方程中的应用有了更深的理解，并且我们理论联系实际给出了证明，使其印象更加深刻。
3小结
从本文可以看出，作为数学物理方程学科中的重要的两种方程――泊松方程和波动方程在我们生活中有着重要的作用，他们相互间也有着密切的联系，并不是我们想象中的难以理解，而且还有一定的规律可循，只要我们对其认真地研究，对培养数学物理方程思维能力和数学物理方程思想方法将有很大的帮助，同时也可以从这一课题的研究中学习到许多积分计算的技巧。泊松积分的应用非常广泛，不仅在波动方程中，而且在普通物理、遥感技术、生物医学、军事技术都有极其重要的作用，在物理方面的贡献尤为突出，因此，对泊松积分在波动方程中的应用的研究是很有意义的。
参考文献：
[1]张渭滨编著.数学物理方程[M].北京：清华大学出版社， 2007，8.
[2]张慧清，吴小吟，杨小军编著.数学物理方程与特殊函数[M].西北工业大学出版社，2005，8.
[3]陈才生主编.数学物理方程[M].南京：东南大学出版社，2006，8.
[4]庄万，萧礼，洪良辰，刘中元编著.数学物理方程[M].济南：山东科学技术出版社，2005，11．
[5]刘经红，朱起定.Poisson 积分方程的解的存在唯一性.湖南大学学报，2004，1：9-12.
[6]魏志道.Poisson 定理的推广.四川师范大学学报，1996，3：69-72.
【关键词】波动方程；混合积分变换法；行波法；积分变换法
【中图分类号】O175.2【文献标识码】A
本项目得到中国石油大学青年教师教改项目（QN201304）和研究生学位点建设项目（XWS13012）的资助.
在科学与工程技术的实际问题中，经常会遇到大量的偏微分方程，如描述电磁场中电场和磁场强度的麦克斯韦方程等.有效的求解方法对于研究这类偏微分方程所体现的科学意义具有重要的理论和实际作用.在《数学物理方程》类课程教学中，通常我们会介绍四种经典的求解方法，即分离变量法、行波法、积分变换法和格林函数法等.事实上，在某些情形下这些方法可加以综合运用.下面以一阶波动方程为例加以说明.
例utt=uxx+tsinx（x∈R，t>0），u|t=0=0（x∈R），ut|t=0=sinx（x∈R）.
解法1基尔霍夫公式
u（x，t）=12∫x+tx-tsinξdξ+12∫t0dτ∫x+t-τx-t+ττsinξdξ=tsinx.
解法2傅氏积分变换法
首先，方程及初始条件两端分别关于空间变量x取傅氏变换，并注意到
[sinx]=iπ[δ（w+1）-δ（w-1）]，
可得
d2Udt2（ω，t）+ω2U（ω，t）=iπ[δ（w+1）-δ（w-1）]t，
U|t=0=0，Ut|t=0=iπ[δ（w+1）-δ（w-1）].
接下来，求解上述二阶常系数常微分方程得
U（ω，t）=iπ[δ（w+1）-δ（w-1）]ω2・t+ωsinωt-sinωtω.
最后，上式关于w取傅氏逆变换，由定义及δ函数性质，得u（x，t）=12π∫+∞-∞U（ω，p）eiωxdω=tsinx.
解法3拉氏积分变换法
首先，方程两端关于时间变量t取拉氏变换，并考虑初始条件及[t]=1p2，得
d2Udx2-p2U=-1+1p2sinx.
接下来，由特征方程法求解上述二阶常微分方程，得
U（x，p）=C1epx+C2e-px+sinxp2.
注意到弦上各点的位移有界，特别是无穷远处.从而必然有C1=C2=0，即U（x，p）=sinxp2.
最后，上式中关于p取拉氏逆变换，得
u（x，t）=-1sinxp2=tsinx.
解法4混合积分变换法
第一步，首先方程两端关于时间变量t取拉氏变换，并考虑初始条件，得
p2U（x，p）-sinx=d2Udx2（x，p）+sinxp2.
其次，上式进一步关于空间变量x取傅氏变换，可得
p2U（ω，p）-iπ[δ（w+1）-δ（w-1）]=-ω2U（ω，p）+iπp2[δ（w+1）-δ（w-1）].
第二步，整理上式求，得
U（ω，p）=p2+1p2（p2+ω2）・iπ[δ（w+1）-δ（w-1）].
第三步，首先关于上式中w取傅氏逆变换，由定义及δ函数性质，得
U（x，p）=12π∫+∞-∞U（ω，p）eiωxdω=sinxp2.
其次，上式进一步关于p取拉氏逆变换，得
u（x，t）=sinx・-11p2=tsinx.
注1混合积分变换法的基本思想是先关于某一自变量进行一次拉氏积分变换，将含两个变量的偏微分方程转化为含一个参量的常微分方程，然后再关于另一变量取傅氏积分变换，得到易于求解的含两个参量的代数方程，最后再依次取相应的逆变换，即可求得原问题的解.
注2上述求解过程给我们展示了如何综合运用多种方法求解定解问题，但读者不能盲目的任意叠加这些方法.在求解之前，一定要分析清楚所给问题的条件到底适合哪些方法.
【参考文献】
[1]王元明.数学物理方程与特殊函数[M].北京：高等教育出版社，2003.
[2]梁昆淼.数学物理方法[M].第3版.北京：高等教育出版社，1998.
关键词: 瀑布型多重网格方法；矢量有限元；外推；矢量波动方程
中图分类号:O 24182
文献标志码:A文章编号:1672-8513(2011)04-0267-05
A New Cascadic Multi-Grid Method for the Vector Wave Equation
LU Kangmei，LI Chengliang，CAO Yanbin
(School of Mathematics and Computational Science，Guilin University of Electronics Technology，Guilin 541004，China )
Abstract: The vector finite element method，which can effectively avoid the spurious modes，is widely used in the analysis of electromagnetic problems，as is the case in the discrete calculating of the vector wave equation．The finite element extrapolation technique for the eigenvalue problem, which is proposed by Yang Yi-duo，is extended and applied to eigenvalue problems of vector wave equation．In this paper，a new cascadic multi-grid method based on the vector field eigenvalue problem is proposed，based on the extrapolation technique．The numerical results show that the new method is a very accurate and efficient method．
Key words: cascadic multi-grid method；vector finite element；extrapolation；vector wave equation
在微波理论和技术中，谐振腔本征值问题是最基本的问题之一.很多微波部件和系统的分析与最优化设计又往往以该问题的求解为基础．矢量有限元方法是近10年来在电磁场计算中应用比较广泛的一种方法．只要选择了合适的矢量基，所考虑结构的内部和外部边界都能够从数学上自然满足，就能够很好地解决伪解问题；又因为有限元方法的网格划分能很好地模拟实际结构，因此我们选择了矢量有限元对谐振腔进行离散计算．
多重网格方法，对于求解由微分方程离散化得到的方程组来说，是目前最快速高效的方法之一，它的求解的工作量可降为O(n)或O(nlnn)．因此它在计算电磁场问题上得到了广泛的关注，最近一些学者发展了基于棱边元的多重网格算法[1-6]．瀑布型多重网格无需粗网格校正，它比一般的多重网格方法计算量减少了，从而提高了计算效率，然而基于棱边元上的瀑布型多重网格方法的研究还是比较少的，特别是矢量场的本征问题，目前研究的也较少．
外推算法是由林群， 陈传淼等[7-8]引入到有限元求解偏微分方程，杨一都在文献［9］中引进了本征有限元外推的一个新技术， 李永明等[10]将有限元外推应用到了波导本征问题．本文在其基础上， 将其推广应用到矢量波动本征问题， 具体的算例表明其可行性和高精度性．
本文结合外推技术，提出了一种基于矢量场本征问题的瀑布型多重网格方法．数值算例结果说明该方法的有效性和实用性．
1 数学基础
对于一个填充相对介电常数为εr和相对磁导率为μr介质的封闭谐振腔体，对应的矢量波动方程为
其中ε为精度控制参数． 以上迭代过程简记为(λj,uj)=NC｛(j,j)｝．
2.3 外推技术
定理1 设剖分是强正规的， 又设λ为对应于波导本征问题的简单本征值， 相应的规格化本征函数u∈C4(Ω)， 则本征值外推估计为:
λ=13(4λh/2－λh)+ο(h4).
详细证明见文献［9］．
类似的，设某矩形谐振腔在粗矩形剖分单元下计算得到的本征值为λh， 随后在此粗剖分基础上加密剖分1次， 再计算得到的本征值为λh/2 ，则根据定理1 ，此矩形谐振腔本征值可由外推λhw=13(4λh/2－λh)（记为λhw=E(λh/2,λh)），得到更精确的解．
2.4 矢量波动本征问题的瀑布型多重网格算法
对于矢量有限元离散产生的一系列本征方程组Aiui=λiBiui(i=0,1,2,…,l)， 若要求解最细网格层l层上本征方程Alul=λlBlul， 我们结合外推技术给出如下求解该本征问题的新瀑布型多重网格算法．
由图2，我们也可以看出随着网格的加密计算的精度越来越高，从而验证了算法的快速收敛性，且和一般的瀑布型多重网格方法比较计算精度也提高了很多．
4 结语
从上述的算法公式，以及矩形谐振腔的谐振频率计算结果的分析可以看到，利用瀑布型多重网格方法大大减少了计算量，将外推推广应用到谐振腔本征问题求解精度得到了很大的提高，因此在比较少的单元剖分下便可以较大的提高计算精度．计算过程也比较简单，易于编程，因此该算法在现代数值计算中是一种十分实用、简单、高效的新方法．
参考文献:
［1］WATANABE K ,IGARASHI H ．Robustness of nested multigrid method for edge-based finite element analysis ［J］． IEEE Transactions on Magnetics， 2009，45(3)：1088-1091.
［2］TSAI C L ,WANG W S ．An improved multigrid tcchnique for quasi-TEM analysis of a microstrip embedded in an inhomogeneous anisotropic medium［J］．IEEE Trans． Microwave Theory Tech，1997，45(5)： 678-686．
［3］WEISS B ,BIRO O .Edge element multigrid solution of nonlinear magnetostatic problems［J］. Comple,2001，20(2)：357-365．
［4］SCHINNERL M ，SHOBERL J ,KALTENBACHER M． Nested multigrid methods for the fast numerical computation of 3D magnetic fields［J］． IEEE Trans Magn，2000，36：1539-1542．
［5］WATANABE K ,IGARASHI H ，HINMA T．Comparison of geometric and algebraic multigrid methods in edge-based finite-element analysis［J］． IEEE Transactions on Magnetics,2005,41(5)：1672-1675．
［6］金建铭，王建国，葛德彪．电磁场有限元方法［M］．西安：西安电子科技大学出版社，2001:165-186.
［7］陈传淼，黄云清．有限元高精度理论［M］．长沙：湖南科学技术出版社，1995：451-492.
［8］李郴良，陈传淼，许学军．基于超收敛和外推方法的一类新的瀑布型多重网格方法［J］．计算数学，2007，37(9)：1083-1098．
［9］杨一都．本征值有限元外推的一个新技术［J］．贵州大学学报：自然科学版，1989，6(3)：6-11．
［10］李永明，俞集辉．有限元外推法在波导本征值问题中的应用［J］．重庆大学学报:自然科学版，1999，22(4)：78-81．
【关键词】非线性波动方程 初边值问题 整体解 衰减估计
一、 引言及主要结论
本文讨论如下初边值问题：
的整体广义解的存在性及衰减性，其中是空间中具有光滑边界的有界域。
我们用Galerkin方法证明问题(1)―(4)的整体广义解的存在性和唯一性，用扰动能量法证明解的衰减性，主要结论为：
定理 假定
(A4)h是非负有界的二次连续可微的实值函数，满足且对某个，当时，有和.其中为正常数。
在F1上满足相容性条件.其中
则对任意的T>0，问题(1)―(4)存在至少一个整体广义解
若令(A2)中的p=1且和(A1)中的β充分小，则存在正常数c和γ，使得
此外，若f1Lipschitz连续，则解是唯一的。
二、定理的证明
设是空间的一组标准正交基，使得.作问题(1)―(4)的近似解
据Galerkin方法，满足下列常微分方程组的初值问题
据常微分方程的一般理论，问题(5)―(6)存在唯一的局部解随后进行的第一个先验估计将说明uk(t)能被整体延拓到[0，+∞)上.
第一个估计
将方程(5)中的wj换成uk(t)，并在(0，t)上积分，得到
由假设(A3)，Cauchy 不等式和Sobolev迹嵌入定理，得到
合并(7)，(8)，选η充分小，利用Gronwall不等式即得第一个估计
其中L1是不依赖于k∈N和t∈[0，T]的正常数。
第二个估计
首先估计ukttn(0)的L2范数，易得uktt≤L2.其中L2是一个不依赖于k的正常数.接下来，对(5)式两边关于t求导一次，而后将其中的wj换为uktt得到
对上式左边第一项进行估计，可得
在(0，t)上积分(9)式，，得到
类似(8)的做法，，可以得到
合并(10)、(11)，选η充分小，利用Gronwall引理得到第二个估计
其中，L3是一个不依赖于k∈N和t∈[0，T]的正常数。
非线性项的分析
利用Lions引理和Sobolev迹嵌入定理易得和.
由第一个估计和假设(A1)可知，存在函数，使得
利用广义格林公式，由(5)得到
由于Δu∈L2())，从而有
再利用广义格林公式，由(13)和(14)得到
由于
接下来，把(5)中的wj换为 uk并在(0，T)上积分，得到
对上式两边取极限得到
合并(14)、(15)和(16)，利用广义格林公式得到
由假设(A1)得到
利用Lions引理即得(12).
惟一性
设u1和u2是问题(1)-(4)的两个解，则z=u1-u2满足
在(17)中令w=zt(t)，由假设(A2)、(A3)可以得到
其中，λ来自于不等式
另一方面，由f1lipschitz连续和假设(A1)可知，存在常数C使得
在(0，t)上积分(18)得到
利用Gronwall引理，由上式即得.
惟一性得证.
能量的一致衰减
由(1)得到
(19)
若令(A2)中的P=1，则存在正常数δ1和δ2，使得(20)
考虑到假设(A1)，由(19)，(20)得到
(21)
注意到h(0)=0简单的计算易知
(22)
其中
定义修正能量
(23)
假定
(24)
其中λ>0来自于不等式则由(23)，(24)可得
(25)
因此，E(t)的衰减是e(t)衰减的直接结果.
接下来，由(21)，(22)，(23)，(25)及假设(A4)得到
定义扰动能量.其中.不难证明，存在正常数C1>2，C2>0使得(26)
和(27)
综合(26)，(27)可得
其中，C，λ为正常数.证毕.
参考文献：
[1]Yang Zhijian.Global existence， asymptotic behavior and blow up of solutions for a class of nonlinear wave equations with dissipative term. J. Differential Equation，2003，187：520-540.
[2]Tokio Matsuyama.On global solutions and energy decay for the wave equations of kirchhofftype with nonlinear damping term.J.Math.Anal.Appl.1996，204：729-753.
[3]M.M.Cavalcanti，V.N.Domingos Cavalcanti，J.S.Prates Filho and J.A.Soriano.Existence and uniform decay of solutions of a degenerate equation with nonlinear boundary damping and boundary memory source term.Nonlinear Analysis T.M.A.，1999，38：281-294.
[4]M.M.Cavalcanti，V.N.Domingos Cavalcanti.Existence and uniform decay of the wave equation with nonlinear boundary damping and boundary memory source term.Calculus of Variations，2002，15：155-180.
[5]M.Nakao.Asymptotic Stability of the Bounded or Almost Periodic Solutions of the Wave Equation with Nonlinear Dissipative Term，J.Math.Anal.Appl.1997，58，336-343.
基金项目：河南省自然科学基金资助项目，编号0211010500.
(作者单位：河南财经学院信息学院)
[关键字] 三维声波方程 有限差分正演 GPU并行计算 CUDA
[中图分类号]P628+.3[文献码] B [文章编号] 1000-405X（2013）-1-284-1
0前言
三维地震勘探技术是一种信息量大、精度高的石油地球物理勘探方法，并已经在实际生产中大规模应用。在复杂构造地区，三维资料的处理与解释需要与三维正演理论模型进行对比验证，国内外学者做了大量的研究工作并提出了许多高精度的三维正演算法。一般来说，正演方法在网格规模快速增长的时候，计算量急速增长，使得普通的桌面级计算机无法承受运算耗时，必须通过工作站和大型机完成运算。GPU并行运算可以有效降低运算耗时，使大数据量三维正演可以在桌面级计算机上完成。
1 三维有限差分算法的GPU算法
三维GPU并行算法的访存模式比较复杂，因为要使用二维线程块逐一计算4个场分量，同一Block的不同线程不能重复利用处于不同面（即 不同）网格内的波场值，所以这里的数据不能载入共享存储器。在TFL模式中，在沿i方向逐步运算的过程中，当前线程会两次利用y和z场分量，因此只能使用寄存器存储当前和上一层的y和z场分量，当前层的波场值在计算过程中逐步转为上一层的波场值。用 方向场值举例递推TFL模式算法步骤如下：
（1）存储当前层y、z方向场值，设置其为上一层场值；
（2）从全局存储器读取当前层y、z方向场值，将其存储到寄存器；
（3）对当前Block中所有线程进行同步；
（4）将y、z场值从寄存器读出，将其载入共享存储器；
（5）同步当前Block内的所有线程；
（6）更新x并写入global存储器；
（7）将当前层号加1，返回步骤（l）重复执行上述步骤，直到最后一层；
2 数值实验
所选测试模型为三维均匀速度模型，网格模型大小为109*109*109，xyz三个方向的网格步长均为1m，时间步长为0.2ms，介质速度为250m/s，有限差分精度为8阶。
截取200ms时波场快照图，并记录各个算法的运行时间，由于程序设计的原因，程序计时可能会与算法实际的运行时间的产生一定的误差，但该误差值较小不会对实验数据产生影响，其结果由图1所示。
由于受到很多因素的限制，GPU算法的加速比很难达到理论峰值。即便如此，从上图我们也可以看出，GPU对算法的加速效果非常明显，GPU1的加速效果是CPU的60倍，GPU2的加速效果是CPU的88倍。
为了更好的GPU的加速效果和加速能力，本文分别针对CPU、GPU1、GPU2进行了进一步的模拟运算，加大网格规模，测试GPU和CPU的运算能力和加速能力，如图1所示。从图中我们可以看出，随着网格规模的增大CPU运算耗时也越来越大，且增速很快，大规模的三维运算在CPU上运行时难以接受的。而GPU的运算耗时虽然也随网格规模在增大，但是增速不明显，说明GPU并行运算适用于波动方程正演。虽然GPU算法的加速效果很好，但是也有极限，图2展示了加速极限。
3 结论
波动方程正演与射线追踪正演相比有很高的精确性。但是，在实际工作中，波动方程并不会是简单的声波方程，一般都会是弹性波和粘弹性波，并且模型一般都较为复杂，这些都会大大增加计算量，如果只用CPU进行并行和串行运算，计算成本十分高昂。我们可以针对专业的GPU设备设计相适应的算法，可以大幅提高计算效率，令使用GPU算法的程序效率达到CPU的60倍左右。然而目前GPU的理论仍有待开发，程序的实际加速比并未达到理论峰值，相信随着技术的逐步成熟，GPU并行运算大有可为。
参考文献
[1]赵改善.地球物理高性能计算的新选择：GPU计算技术[J].勘探地球物理进展，2007，30（5）：399-404.
[2]金胜汶，陈必远，马在田.三维波动方程有限差分正演方法[J].地球物理学报，1994，37（6）.804-810.
[3]赵改善，李剑峰，王于静等.网格计算技术及其在石油勘探开发中的应用前景[J].石油物探，2005，44（5）：413-420.
关键词：声波方程；隐式差分；先验估计；代数多重网格
1 引言
波动方程是一种主要描述声波，水波和光波的波动现象的双曲型偏微分方程，到目前为止，已有不少人对一维和二维波动方程做过系统的研究。但对三维波动方程的研究的人相对较少，其主要采用的差分方法直接计算存在计算量大收敛速度慢的缺陷。文 中对三维声波方程的采用有限差分法进行差分离散，再采用数值频散法对地震波进行叠前和叠后数值模拟，文 中的精细积分法；文 的采用三维频率空间域有限差分波场延拓算子以“逐步-累加”的方式实现了三维波动方程基准面校正 ；文 采用延拓法实现三维波动方程正演，这些方法受到差分精度及网格剖分步长的严重影响，要想得到比较准确的数值结果须极小的剖分步长和非常高精度的差分，从而导致计算时迭代缓慢，计算量巨大。
偏微分方程数值解一般有两种离散方法，差分离散方法和有限元逼近方法。差分方法有多种格式，三点格式、五点格式、紧差分格式 等。有限元法则是电子计算机时代的产物，有多种专门的应用软件 ，在应用范围方面比差分方法更优越。以上两种方法离散和所得的方程，通常计算量巨大，即使用计算机处理也不是很简单，有些方程用巨型计算机处理都会造成单片CPU内存不足的情况。
多重网格法是求解偏微分方程的高效快速算法，具有收敛速度快，精度高等优点，其收敛速度和效率并不受差分步长的影响，也就是说当离散更精细时，其收敛速度并不减慢，因此多重网格法可大大的减小计算量及占用的CPU内存，加快运算速度。文献 用多重网格法对二维及三维理想化的波动方程做过研究，但不能解决一般的三维波动方程 ，因此笔者在本文中用多重网格法讨论诸如文 中的波动方程，即先采用紧差分格式将上述波动方程的解离散，得到的差分方程组每个时间层上求解类似于求解三维椭圆问题，再采用类似于文 中的多重网格法求解。并对三维波动方程所采用的差分方法进行了讨论，根据文 中一维，二维差分格式解的误差先验估计式，推导了三维差分格式及K维差分格式解的误差先验估计式，并对其进行了证明，利用三维格式格式解的误差先验估计式，证明三维差分格式解的存在性和收敛性，及稳定性。
1、声波方程
在非均匀介质中，三维声波方程为
2、差分方法
以 为时间步长，空间取等距离的立体网格，其边长为 ，网格点为 ， ， ， ， ， ， ，
这里对 项暂不讨论，具体参见文 ，在文 中 ，则 。
可依次迭代求出式（10）各项的值。
3、差分格式解的先验估计式
首先引进如下记号
以上前4个分别为无穷范数（一致范数），2范数，差商范数，和 范数。
定理1 设 为差分方程组
其中
证明 用 乘以（14）的两边，并对 求和，得到
根据文 中一维，二维二阶双曲型方程隐式差分格式解的误差先验估计式及定理1我们可推导出K维二阶双曲型方程隐式差分格式解的误差先验估计式。
定理2 对于K维的二阶线性双曲型方程
的隐式差分格式，设 为差分方程组
的解，则对任意步长比 ，有
其中
4差分格式解的存在性，收敛性和稳定性
存在性
5、多重网格法
多重网格法有V循环和W循环，其中常用的是V循环。通常包括3个要素，松算子、限制算子和插值算子。
令
5.1多重网格算法的实现
（1）利用定理1，和（25）式，估计在所要求的误差范围内 和 的大致取值。
（2）对于第 个时间层，把得出 的值，作为细网格层， 作为粗网格层，由粗网格层开始进行差分。
（3）计算 时刻的误差残量，使用完全加权算子将其限制在相邻的粗网格上。
（4）在粗网格层上用G-S迭代进行一次误差磨光，具体算法参见 。
（5）将粗网格层上的误差校正结果用三线性插值算子返回到细网格层上。
（6）在细网格层上用G-S迭代进行一次误差磨光，具体算法参见 。
（7）令 ，即转入下一个时刻的计算。
（8）重复（2）-（7），直到计算需要达到的时刻。
6、 数值算例
以上问题的精确解为
通过对上述算例的计算，结果表明：，多重网格法比传统的迭代法计算速度快，随着 的减小，越明显。
参 考 文 献
【1】何兵寿，魏修成，刘洋。三维波动方程的数值频散关系及其叠前和叠后数值模拟【J】。石油大学学报（自然科学版） ，2001，25（01） ：067-071。
【2】王润秋。精细积分法三维地震波正演模拟【J】。勘探地球物理进展，2006-12，29（06）：0394-0397。
【3】杨锴。三维波动方程基准面校正方法的应用研究【J】。地球物理学报，2007-7，50（04）：1232-1240。
【4】熊高君。三维混合延拓一步法波动方程正演模拟【J】。成都理工学院学报，1999-1，26（01）631-634。
【5】熊高君。三维波动方程正演及模型应用研究【J】。石油勘探，2005-11，44（06）0554-0556。
【6】孙志忠。偏微分方程数值解法【M】。北京：科学出版社，2007。
【7】马爱军，周传月，王旭。Patran和Nastran有限元分析。北京：清华大学出版社，2005。
【8】葛永斌，吴文权，卢 曦。求解二维扩散方程的加权平均隐式多重网格方法【J】内蒙古大学学报（自然科学版）。2003-9，34（05）0490-0494。
【9】葛永斌，吴文权，田振夫。二维波动方程的高精度隐格式及多重网格算法【J】厦门大学学报。2003-11，42（06） 0691-06996。
【10】葛永斌，田振夫，吴文权。三维波动方程的隐式多重网格方法【J】华中师范大学学报。2005-3，39（01）0009-0012
【11】舒适。几类基于几何和分析信息的代数多重网格法及其应用【D】（湘潭大学）博士论文。2004-7。
【12】谢学斌。代数多重网格法在岩体力学有限元分析中的应用【J】。工程力学，2005-11，22（05）0165-0170
【13】谢学斌。岩体力学三维有限元分析的代数多重网格求解法【J】。岩体力学与工程学报，2006-11，25（11）2358-2363。
【14】肖映雄。 代数多重网格算法研究及其在固体力学计算中的应用【D】。（湘潭大学）博士论文。2007-4。
【15】孙志忠，袁慰平，闻震初。数值分析【M】。南京：东南大学出版社，2006。
【16】胡健伟，汤怀民。微分数值解法【M】。北京：科学出版社，2007。
关键词： 格林函数 冲量定理 推导过程
格林函数又称为点源函数，它表示一个点源在一定的定解条件下所产生的场.由于任意分布的源所产生的场可看成许许多多点源产生的场的叠加，因此格林函数一旦求出，就可算出任意源的场.在一般情况下，应用格林函数法解题的关键在于求格林函数，对于波动与输运这类含时间的定解问题，其格林函数的求解可借助于冲量定理法求得.
本文具体分析了如何利用冲量定理法求解一维波动方程和一维输运方程的格林函数，希望本文的讨论有助于理工类本科生更好地掌握数学物理方法等相关课程.
一、用冲量定理法求解一维波动方程的格林函数
1.一维无界空间中的受迫振动问题
2.一维有界空间的波动定解问题
二、用冲量定理法求解一维输运方程的格林函数
1.一维无界空间的输运定解问题
2.一维有界空间的输运定解问题
3.一维半无界空间的输运定解问题求解
三、结语
不同的定解条件下，一维波动方程和一维输运方程所对应的格林函数不同，但均可用冲量定理法求出相应的格林函数，尔后将格林函数以一定的积分形式表示所讨论的数理方程的定解问题对应的通解.格林函数的冲量定理法，还能用于求解非齐次方程，齐次方程兼带非齐次定解条件的问题.
参考文献：
[1]梁昆淼.数学物理方法[M].北京：人民教育出版社，1979.
[2]陈恕行，秦铁虎，周忆.数学物理方程[M].上海：复旦大学出版社，2003.
[3]梁昆淼.数学物理方法[M].北京：高等教育出版.
[4]梁昆淼.数学物理方法[M].高等教育出版社，1978.6.
[5]姚端正.数学物理方法[M].武汉大学出版社，1991.10.
[6]马文金.计算物理学[M].科学出版社，2005.5.
[7]复旦大学数学系.数学物理方程[M].人民教育出版社，1979.5.
[8]郭玉翠.数学物理方法学习指导[M].清华大学出版社，2006.2.
[9]程建春.数学物理方程及其近似方法[M].科学出版社，2004.8.
[10]李泉凤.电磁场数值计算与电磁铁设计[M].清华大学出版社，2002.10.
[11]李荣华，冯果忱.微分方程数值解法（第三版）[M].北京：高等教育出版社，1996.
[12]李立康，於崇华，朱政华.微分方程数值解法[M].上海：复旦大学出版社，1999.
[13]陆金甫，顾丽珍，陈景良.偏微分方程差分方法[M].高等教育出版社，1988.
[14]林群.微分方程数值解法基础教程（第二版）[M].科学出版社，2004.2.
[15]吕英华.计算电磁学的数值方法――高等院校与通信工程系列教材[M].清华大学出版社，2006.
安徽省高校省级自然科学项目（KJ2014A236）
[关键词]无界域 行波法 达朗贝尔公式 积分变换法 发散思维
在《数学物理方程与特殊函数》这门课中，对于无界域内波动方程的定解问题，我们通常采用行波法求解，从而得出了一维波动方程的达朗贝尔公式（4）。
无界长弦做自由横振动，满足下列定解问题
作变换 ，得方程（1）通解u（x，t）=f1（x+at）+f2（x-at），
其中f1、f2是任意两次连续可微的函数。
在定解条件下的解为
（4）
于是学生认定，无界域波动方程定解问题只能用行波法，不认为傅立叶变换方法、拉普拉斯变换方法是可行的，教材束缚了他们发散性思维和创造性思维。该问题是适定的，即问题的解存在唯一并且是稳定的，因此采用其他解法，都能得到与（4）式相同的解，下面给出此问题的积分变换解法，以培养学生的创造性思维。
一、应用傅里叶变换求解
设U（ω，t）、Φ（ω）、ψ（ω）分别为 的象函数
对（1）式取关于x的傅里叶变换
对（2）（3）式取同样变换
由（5）得通解U（ω，t）=c1cosaωt+c2sinaωt
（6）（7）代入上式得
二、应用拉普拉斯变换法求解
设 分别为 的象函数
对方程（1）取关于t的拉普拉斯变换，得
.
这是一个关于象函数 的二阶非齐次常微分方程，其通解为
.
考虑到 不应为无限大，常数A取为零，同时，
也不应为无限大，常数B也取为零。为保证积分收敛，第一个积分的下限应取为+∞，第二个积分的下限取为-∞，于是
对上式两端进行拉普拉斯逆变换
联立上述各式易得解（4）.
对于无界弦振动的柯西问题，上述解法也是可行的。积分变换法不受方程类型的限制，对于波动方程、拉普拉斯方程、热传导方程的定解问题，都可求解. 在教学过程中，教师要采用提问式和启发式教学，支持学生从多个角度以不同方式得到解决问题的多种方案或最优方案，培养他们的发散性思维和创造性思维。
一、因子面板数据随机波动模型
对金融资产收益率的变动特征刻画主要基于波动率和随机干扰项的共同作用，随机波动模型的研究主要体现在随机波动项的设定上。鉴于随机波动率的不可测性以及正定性，除了采用因子模型以外，还可以采取矩阵指数变换、Cholesky分解以及Wishart模型等来反映随机波动项的协方差成分。这些设定一方面是为了体现随机波动对收益率的影响模式，另一方面也是为了在模型估计、诊断检验和模型比较时获得好的效果。由于金融资产价格变动影响因素的复杂性，很难得出哪种方法更为有效的一致性结论，所以随机波动项的设定需要根据具体的数据类型来确定。因子面板数据随机波动模型由于包含协变量、不可观测因子以及随机误差项，其随机波动性既可以体现在公因子部分，也可以体现在随机误差项中。本文主要考虑统计因子和基础因子两种类型的公因子，其中统计因子包含了某些不可观测的影响因素。通过这样的设定，能够有效地简化随机误差项的结构。其中，rit(i=1，2，…，N，t=1，…，T)表示第i项金融资产在第t时期的对数收益率。假设资产收益率不仅受以往波动率的影响，还受其他一些可观测与不可观测因素的影响。与一元随机波动模型和多元随机波动模型类似，面板随机波动模型同样包括两个部分:均值方程和波动方程。一般情形下，我们可以考虑采用变系数随机效应面板数据模型来反映可观测和不可观测的因子对金融资产收益的影响。这样，面板数据随机波动模型均值部分设定为:对于均值方程式(1)和式(3)，后者采用因子分解来表示随机效应和固定效应，在简化模型随机误差项的结构的同时增强了模型的解释意义。由于因子模型的公因子代表了市场外部的共同冲击，在金融资产分析中表示资产价格的变动受某些不可观测因素的共同影响。在解释变量中选定某些可观测的影响因素，比较直接地体现了影响资产价格变动的某些客观因素。波动方程式(2)和式(4)分别反映随机波动项的条件自回归异方差性和共同影响因素的滞后效应。所以式(1)－式(4)能够对金融资产的平均收益和条件波动做出合理解释。
二、面板数据因子模型的贝叶斯估计
按照分块抽样的思路，协变量系数βi可以作为一个单独的块进行讨论。由于βi反映了解释变量对被解释变量的影响程度和影响方向，在因子面板数据随机波动模型中，βi的估计结果和因子分解一起决定了波动方程和动态因子方程中潜变量的取值和模型的估计结果。λi和ft作为面板随机波动模型的因子载荷和公因子，可以同时进行估计。相比较而言，λi的估计更加复杂，其自由参数较多。如果将ft看作潜在解释变量，在ft取值给定的情况下，λi的估计过程和解释变量系数βi类似，除了参数有所区别外，其后验分布的类型相同。在实际应用中，潜变量hit和qjt虽然来自不同的参数变换，但是二者的变换形式相同，仅仅只有分层分布参数不同，所以估计过程的设定相似。由于潜变量和波动方程的结构类似，我们将采取相同的方法对两组不同的波动方程的系数(αi0，αi1)和(φj0，φj1)进行估计。(一)βi后验分布推断不失一般性，此处我们仅仅考虑个体随机系数面板数据模型。面板数据模型按照解释变量的系数是否随个体变化分为固定系数模型和随机系数模型。随机系数模型是指各个解释变量的系数依个体而变动，这样，各个个体无论是模型的截距项还是斜率项都不同。由于我们在此考虑的解释变量主要是影响资产价格变动的可观测因素，不同的个体之间这些观测因素差别较大，选择固定系数不能很好地体现该差异，所以选择随机系数模型。同时，与其他线性回归模型的贝叶斯估计类似，面板数据随机系数模型的估计也可以采用非分层贝叶斯方法和分层贝叶斯方法。在面板数据随机波动模型中，采用分层贝叶斯方法能够充分利用相关的数据对计算结果不断进行检查，便于找出合适的条件后验分布。所以此处随机系数模型解释变量的系数βi设定为与个体有关，并且采用分层先验的形式。
式(3)的误差项与波动方程有关，所以其似然函数非常复杂。为了便于说明问题的方便，我们首先假设βi(i=1，2，…，N)不存在个体之间的交互效应。其分层先验是从正态分布中独立提取出来，这样根据Koop(2003)，本文假设。由前面对误差项以及因子项的分析可知，随机系数βi的后验分布参数不仅与随机误差项有关，而且与因子分解的结果有关。为了便于对相关参数的算法进行设计，我们总可以采用某种基础变换，使得个体随机系数模型的误差项均值为0N，方差—协方差矩阵为M－1IN，下标N表示矩阵的阶数，M为误差精度。这样，在给定观测数据和误差精度的情况下，就可以设定βi的后验分布参数，为采用Gibbs抽样对其进行估计做准备。由于βi的后验分布不仅与先验信息和实际数据有关，还受模型的误差精度M影响，所以模型剩余部分的设定将直接影响相关参数的估计。与一般面板数据随机效应模型相比，因子面板数据随机波动模型的误差精度不仅与随机误差项有关，而且受因子分解结果的影响。(二)λi后验分布设定为了便于对因子载荷矩阵λi进行识别，需要添加一定的识别条件。在此，假设λi为下三角矩阵，进而其元素λij=0(i＜j，i=1，…，N，j=1，…，p)。对于因子面板数据随机波动模型，λij的先验分布可以设定为:其中，1(•)为示性函数。对面板随机波动模型的误差成分进行因子分解是为了分析不可观测因素对个体和时间的影响。在面板数据固定效应部分β''''ixit已经设定的条件下，因子分解可以看做是对已知信息集的扩展，从中得出公共因子以及相应的特定系数。因子载荷矩阵的非零元素的后验分布仍然服从正态分布。由于公因子个数p一般远小于原始变量的个数N，因子载荷阵的组成元素的后验分布与其所处位置有关。波动项ht的共轭后验分布参数由此得出，其中波动方程的系数主要由波动项的生成过程确定。设定好波动项分层先验参数和后验参数后，随机波动项的系数可以采用MCMC结合粒子滤波算法给出。
在因子随机波动模型中，动态线性分层模型(DLHM)式(35)－(39)代表了随机波动过程的演进步骤。由于对数误差波动部分采用了分块的方法，通过构造块hT实现对N个独立的随机波动过程的系数αi进行联合估计，这样能够有效提高运算效率得出合理估计结果。对数误差波动项ht与联合系数αi的后验分布设定过程的主要区别在于分块移动能改变ht的后验分布的类型。式(35)已经给出了采用独立的单个随机波动项的分布，即正态分布。如果对块hT进行抽样，其分布类型将不再满足原有条件。对一元随机波动模型，Chib等人(2002)［12］采用Kim等人(1998)提出的7个成分正态分布作为对数χ2分布的近似。对于多元随机波动模型和面板随机波动模型，由以上所假设的独立条件，在将多元波动方程进行分解后，N个变量分解成了N个独立的一元随机波动过程，在均值方程的随机误差项满足对数χ2分布的条件下，同样可以采用7成分正态分布对随机波动项的分布予以近似。此时，随机误差项所服从的对数χ2分布可以表示为:
三、金融资产收益的影响因素分析
在一元和多元随机波动模型的基础上，面板随机波动模型从数据类型方面进行了扩展，将模型的研究对象从一元时间序列数据和多元时间序列数据推广到了面板数据，并考察了更多影响资产价格和收益率变动的因素。其中，不可观测因素用潜在因子表示。整个模型由三部分组成，包括面板均值方程、波动方程和因子方程。均值方程体现了资产的平均收益，波动方程主要是为了刻画金融资产的波动特征，因子方程反映不可观测的影响因素。与传统模型相比，面板因子随机波动模型的构建主要是在多元随机波动模型的建模思路上引进面板数据模型的分析方法，进一步在因子随机波动模型的基础上研究潜在波动性。由前文的分析可知，与多元随机波动模型相比，尤其是与多元因子随机波动模型相比，因子面板数据随机波动模型不仅能够分析被解释变量受某些不可观测因素影响，还能够具体度量可观测因素的影响程度大小。由于随机效应模型系数随个体的变化而变化，通过对个体影响因素的分析，能够合理实现资源配置，优化投资组合，以及控制金融资产的风险。所构建的模型虽然形式上比较复杂，但是各组成部分具有明确的现实意义，分别反映了金融资产收益率的波动特征和主要影响因素。为了具体分析金融资产配置的各种影响成分，尤其是不可观测的影响因素，本文采用中国股票市场的数据，利用因子面板数据随机波动模型研究金融资产收益变动的几个影响因素。此处随机选取了中国股票市场上海证券交易所上证50指数的20只成分股，涵盖银行、证券、钢铁、地产、通信、贸易、建筑、汽车、能源和采矿等10大行业，具体股票代码见表1和表2。为了便于显示，省略了每只股票代码前面共有的三个数字“600”，例如“000”代表“600000”浦发银行。数据来源于国泰安数据中心的CSMAR数据库，原始数据交易日期从2010年4月8日至2012年12月31日，时间跨度为1000天。考虑这些日期有的不是交易日，有的交易日部分股票因种种原因临时停牌，为了研究的需要，最终只保留了股票均进行交易并且具有可比性的499个交易日的数据。
我们主要根据以上设定的因子面板数据随机波动模型分析影响多元资产收益率变动的某些可观测和不可观测因素。由于场外因素的复杂性，在此仅考虑场内因素。被解释变量为每只股票的非预期收益率(考虑分红派现等因素)，解释变量包括每日交易量(股，记为“Vol”)、交易金额(元，记为“Amo”)以及流通市值(千元，记为“Val”)的变动，即一阶差分变换。最后进行标准化，得到三个解释变量Var1、Var2、Var3，该变换过程可以表示为。其中，下标t表示时间，表示个体的下标i省略，μ和σ分别表示下标对应变量的均值和标准差。标准化的目的是为了对比可观测因素对被解释变量的影响强度。均值方程设定为不包含常数项的个体随机系数面板数据因子模型，具体形式同式(1)－(4)。采用联合估计得出因子载荷和模型系数的估计结果。此处提取了三个公因子，公因子个数的选择方法根据Bai和Ng(2002)。三个因子载荷的估计结果见表1，个体随机系数的估计结果见表2。从表1可以看出，各个个体的因子载荷的分布并不均匀。与经典的多元因子分析不同，在因子模型中，我们并没有将因子载荷进行正交旋转，因为此处的目的是为了体现资本收益率受某些共同因素的作用结果。这三个公因子虽然没有具体含义，但是代表了三个共同的潜在影响因素。对所研究的20只股票进行分析，因子“1”上载荷值较大的股票包括五矿发展(058)、包钢稀土(111)、阳泉煤业(348)和江西铜业(362)等。因子“2”上载荷较大的股票有华夏银行(015)、民生银行(016)和招商银行(036)等。因子“3”载荷值总体较小，上汽集团(104)和兖州煤业(188)等股票受因子“3”影响较大。各个公因子在不同股票上的因子载荷大小体现了这些股票的冲击程度。如果能确定这些公因子的来源，在构建投资组合时，可以考虑采用这些不可观测信息针对各个因素的来源增加投资收益。
协变量系数估计虽然与潜在因子和因子载荷的估计同时进行，但是从模型的结构可以理解为提取出公共因素以后再分析各个个体受可观测因素影响的程度。在对原始变量进行标准化处理之后，能更加明显地观察三个解释变量作用于被解释变量的方向和程度。从表2可知，每日流通市值的变动(var3)相对于交易量(var1)和交易额(var2)的变动对收益率的影响更大，其系数估计结果也更为显著。由表2，每个变量系数估计结果的下方为对应变量的估计标准误差，选取5%的显著性水平，除了交易额(var2)在华夏银行(015)的估计结果不显著外，其他各变量的系数估计均通过了5%的显著性检验。从表2估计结果可知，除了包钢股份(010)和江西铜业(362)交易额(var2)的变动与股票收益率为负相关关系，其余变量在所有股票上系数估计结果为正数，都表现为正相关关系。可以解释为当交易量或者流通市值增加时，所有股票收益率呈上升趋势。当交易额增加时，除了包钢股份(010)和西铜业(362)收益率呈下降趋势外，其余股票均呈上升趋势。表2还可以进一步观测每只股票的收益率与各个变量的变动的具体影响大小。在金融资产配置的分析中，与多元因子模型一样，面板因子模型缺乏较好的预测效果。因为因子随机波动模型包含许多可观测因素，对被解释变量进行预测时，首先要能获取这些可观测因素的信息。正如一个硬币有两面一样，因子面板数据随机波动模型在资产配置构建和风险管理中，可以直观地分析这些可观测因素对构建投资组合和风险的影响，并且能够分析不同金融资产之间的差异，如上例中包钢股份(010)和江西铜业(362)所表现出的特点。选择可观测因素时，还可以选择场外因素，尤其是一些通过大众渠道可以获取的场外消息，经过量化以后，可以加入模型中，这样就能全面考虑股票市场风险控制中市场风险和系统性风险。应用例子进一步证实面板数据因子随机波动模型能较好拟合股票市场金融资产收益率受可观测和不可观测因素共同影响的时变波动性和异方差特征。
四、结论
因子面板数据随机波动模型与一般因子随机波动模型相同之处在于存在两个波动方程，这两个波动方程分别代表了资产的随机波动和共同冲击。虽然估计时我们将两个波动方程合并成了一个，然而由于因子面板数据随机波动模型加入了可观测因素，依赖于内部数据的多元MCMC算法并不能直接采用被解释变量或者金融资产收益率数据，因此必须采取某种联合估计方法。同时，由于模型中可观测和不可观测部分的估计方法不一致，因此面板随机效应模型与滤波和抽样算法需要整体考虑，本文采用的是前向滤波倒向抽样整体估计算法。模型中可观测因素可以包括资本市场内部因素和外部因素，在金融资产影响因素的应用分析中，仅仅考虑了市场内部因素。前面已经说明，之所以没有考虑市场外部因素是由于其较难量化，另外一个原因是我们认为采用面板数据因子模型可以分别通过公因子来反映共同冲击，由于共同冲击既有对整个市场的冲击，也有针对某个行业或某只股票的冲击，所以在进行部分行业随机波动比较时，可以采用行业因子代替潜在因子，此时模型结构类似于Fama和French(1992)的基础因子模型。
分层贝叶斯估计中，各个参数的先验分布和后验分布参数(超参数)的设定对MCMC算法的收敛速度有一定影响。进行应用分析时，在参数分布形态已经确定的情况下，根据经验数据来选择先验参数仍然很重要。由于需要设定的参数较多，并且在因子面板数据模型中需要考虑协变量的经验分布，因此先验参数的设定不仅与解释变量有关，还与被解释变量有关。在实行分块抽样时，充分考虑到了各部分之间的相互联系。由于金融资产收益率同时受可观测的市场因素和不可观测的潜在因素的影响，不可观测因素受某些随机因素的驱动，因此可以在多元随机波动模型的基础上，考虑可观测因素和潜在因子对收益率波动特征的刻画。高维面板数据随机波动模型包括大量参数和潜在因子需要估计。本文引入了基于FFBS的联合估计方法对模型进行估计，潜在波动性的设定在给出某些提议的先验信息基础上通过研究其平稳后验分布实行。应用例子进一步证实了面板数据因子随机波动模型能够较好地拟合股票市场金融资产收益率的时变波动性和异方差特征，可以用于构建投资组合和风险管理。
关键词：叠后逆时偏移 优化递推算法 实际资料处理
1 引言
在1982年的SEG年会上，Whitmore提出了一种全新的波动理论偏移方法，这就是后来迅速发展起来的逆时偏移。
2 优化递推算法方程原理
逆时偏移与传统的偏移方法有很大的不同。逆时偏移是在时间轴上实现外推。传统的波动理论偏移方法基于单程波方程，而逆时偏移则基于全波方程。逆时偏移实际上可以看作是反时间方向的正演模拟过程。
叠后逆时偏移的成像原理是基于爆炸反射面模型。根据爆炸反射面模型，将介质速度取为实际速度的一半，这样可以把地下反射面上的每一个点当作是一个二次震源，所有的二次震源在零时刻同时起爆，在最大时刻到达接收点，接收点接收到的信息可以认为是来自地下各个二次震源的直达波的信息。采用的基本流程如（图2）：
4 模型数据及实际资料试处理
针对SEG年会上的岩丘模型（Bird Model），图中蓝色部分为岩丘形状，为高速异常体，本研究采用此模型检测逆时偏移成像方法对复杂介质的成像效果，其中重点关注岩丘的边界以及岩丘下的构造成像，通过声波波动方程高阶有限差分方法对模型进行正演模拟299炮，道距25米，满覆盖次数75次，炮间距50米。图3代表所建速度深度模型，图4表示基于模型数据处理后的全波动方程叠后深度偏移叠加面，它的处理效果和成像精度（框中部分）均逊于优化递推算法叠后逆时偏移（图5）。
本研究选取2009年新疆二维实际资料做了测试，针对实际资料新疆区块-line355线的叠后时间偏移（图6）与优化递推算法叠后逆时偏移效果对比可以看到，后者无论在各地层展布特征、厚度变化规律，还是构造格架等都具有着更高的分辨率（框中部分），同相轴连续性、波组接触关系，内幕构造更加清晰、断裂发育特征更加清楚，取得了比较理想的效果。
5 结论与认识
（1）基于模型正演数据的优化递推算法叠后逆时偏移，取得了很好的成像效果，它较好的消除或减弱了层间反射，比全声波波动方程法具有更高的成像精度。
（2）优化递推算法叠后逆时偏移应用于实际资料试处理也取得了一定成效，显示该方法具有比较好的应用前景，值得进一步深入研究。
参考文献
董良国.复杂地表条件下地震波传播数值模拟，勘探地球物理进展， 2005， 28， 631～445