一种基于FPGA的PML边界三维地震波传播模拟方法技术

本发明专利技术公开了一种基于FPGA的PML边界三维地震波传播模拟方法，本发明专利技术针对FPGA设备应用，与以往使用CPU，GPU等冯诺依曼架构计算单元相比，FPGA具有更高的可定制性，以及更高的并发性。FPGA器件属于专用集成电路中的一种半定制电路，是可编程的逻辑列阵，能够有效的解决原有的器件门电路数不足的问题。FPGA的基本结构包括可编程输入输出单元，可配置逻辑块，数字时钟管理模块，嵌入式块RAM，布线资源，内嵌专用硬核，底层内嵌功能单元，FPGA具有布线资源丰富，可重复变成，集成度高优势。

【技术实现步骤摘要】
一种基于FPGA的PML边界三维地震波传播模拟方法
本专利技术涉及石油勘探、地质成像领域，具体是一种基于FPGA的PML边界三维地震波传播模拟方法。
技术介绍
随着油气田勘探程度的不断深化,普通易探明油气田几乎没有剩余，但油气资源又是现代工业所必须产品，故只能向结构复杂地域、海洋进行探索。现有呈像方法多使用克希霍夫射线法来进行呈像处理，此方法对盐下构造、复杂地质构造探照度不够，从而使用逆时偏移(RTM)方法来对地下地质结构进行呈像。这是一种基于波动方程在时间-空间域用高阶差分方方程求解声波偏微分方程的呈像方法，此方法可以真实的模拟波在地下的传播过程，无倾角限制，适合海洋及复杂地质结构，在复杂三维结构呈像方面有明显优势，呈像精度较射线法更高。在油气地质勘探过程中使用逆时偏移方法进行呈像涉及到对波动方程的求解，此过程计算量巨大，现有技术方案使用的为GPU加速卡进行多核多线程，相较于传统使用CPU计算有较大的提升。但在现今阶段GPU加速卡价格昂贵且功耗高，现有方案均是通过调用NVIDIA公司发布CUDA平台对高阶差分方法波动方程进行求解，由于CUDA平台为GPU通用性方案，且GPU硬件结构的限制，故没有针对高阶差分方程求解波动方程方法进行优化，从而造成GPU在求解波动方程时有一定局限性，在保留计算精度的情况下计算花费时间比较长。PML边界处理是地震成像过程中使用3DRTM所必须的对波场数据的处理约束，减除反射波，从而使得边界清晰，成像效果准确。本专利技术利用了FPGA硬件特性，针对高阶差分方程求解波动方程做了硬件上的特异化的优化，对方程进行了等价拆分、组合。
技术实现思路
本专利技术的目的在于提供一种基于FPGA的PML边界三维地震波传播模拟方法，以解决上述
技术介绍
中提出的问题。为实现上述目的，本专利技术提供如下技术方案：一种基于FPGA的PML边界三维地震波传播模拟方法，将三维二次波动方程分解成x,y,z三分量，并引入PML边界条件进行求解，求解主体为波动方程的差分方程求解方式，所述基于FPGA的PML边界三维地震波传播模拟方法包括：S101：将所需变量从外部存储介质传入芯片内部；S102：求解差分方程的中间量B，计算对应xyz方向上分量，首先计算X方向上分量，例如：计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；a等于w1乘以p(i-(1/2),j,k,n)加上p(i+(1/2),j,k,n)的和，b等于w2乘以p(i-(3/2),j,k,n)加上p(i+(3/2),j,k,n)的和，c等于w3乘以p(i-(5/2),j,k,n)加上p(i+(5/2),j,k,n)的和，d等于w4乘以p(i-(7/2),j,k,n)加上p(i+(7/2),j,k,n)的和，直到m等于W(M/2)乘以p(i-((M-1)/2)),j,k,n)加上p(i+((M-1)/2),j,k,n))的和，其计算方程为：S103：求解差分方程的中间量B，计算y分量值，例如：计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；a等于w1乘以p(i,j-(1/2),k,n)加上p(i,j+(1/2),k,n)的和，b等于w2乘以p(i,j-(3/2),k,n)加上p(i,j+(3/2),k,n)的和，c等于w3乘以p(i,j-(5/2),k,n)加上p(i,j+(5/2),k,n)的和，d等于w4乘以p(i,j-(7/2),k,n)加上p(i,j+(7/2),k,n)的和，直到m等于W(M/2)乘以p(i,j-((M-1)/2)),k,n)加上p(i,j+((M-1)/2),k,n))的和，其计算方程为：S104：求解差分方程的中间量B，计算z分量值，例如：计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；a等于w1乘以p(i,j,k-(1/2),n)加上p(i,j,k+(1/2),n)的和，b等于w2乘以p(i,j,k-(3/2),n)加上p(i,j,k+(3/2),n)的和，c等于w3乘以p(i,j,k-(5/2),n)加上p(i,j,k+(5/2),n)的和，d等于w4乘以p(i,j,k-(7/2),n)加上p(i,j,k+(7/2),n)的和，直到m等于W(M/2)乘以p(i,j,k-((M-1)/2)),n)加上p(i,j,k+((M-1)/2),n))的和，其计算方程为：S105：求解差分方程的P，a1等于w1乘以B(i-(1/2),j,k,n)加上B(i+(1/2),j,k,n)的和，b1等于w2乘以B(i-(3/2),j,k,n)加上B(i+(3/2),j,k,n)的和，c1等于w3乘以B(i-(5/2),j,k,n)加上B(i+(5/2),j,k,n)的和，d1等于w4乘以B(i-(7/2),j,k,n)加上B(i+(7/2),j,k,n)的和，直到m1等于W(M/2)乘以B(i-((M-1)/2)),j,k,n)加上B(i+((M-1)/2),j,k,n))的和；将上述得出的结果进行求和a1+b1+c1+d1+...+m1得到X，再对X乘以速度相关系数；a2等于w1乘以B(i,j-(1/2),k,n)加上B(i,j+(1/2),k,n)的和，b2等于w2乘以B(i,j-(3/2),k,n)加上B(i,j+(3/2),k,n)的和，c2等于w3乘以B(i,j-(5/2),k,n)加上B(i,j+(5/2),k,n)的和，d2等于w4乘以B(i,j-(7/2),k,n)加上B(i,j+(7/2),k,n)的和，直到m2等于W(M/2)乘以B(i,j-((M-1)/2)),k,n)加上B(i,j+((M-1)/2),k,n))的和；将上述得出的结果进行求和a2+b2+c2+2d+e2+...+m2得到Y，再对Y乘以速度相关系数；a3等于w1乘以B(i,j,k-(1/2),n)加上B(i,j,k+(1/2),n)的和，b3等于w2乘以B(i,j,k-(3/2),n)加上B(i,j,k+(3/2),n)的和，c3等于w3乘以B(i,j,k-(5/2),n)加上B(i,j,k+(5/2),n)的和，d3等于w4乘以B(i,j,k-(7/2),n)加上B(i,j,k+(7/2),n)的和，直到m3等于W(M/2)乘以B(i,j,k-((M-1)/2)),n)加上B(i,j,k+((M-1)/2),n))的和，再加上2p(i,j,k,n)减去p(i,j,k,n-1)得出最终P(i,j,k,n+1)的结果，其计算方程为：S106：重复上述步骤一到步骤五直到得到最终所有点的n+1时间阶输的结果，此时n+1时间阶计算完毕；S107：对基于步骤六的结果，对所有时间阶数进行递增计算，得出最后Nmax最大时间阶数结果，即为当前波动方程差分方程方法求解最终结果。作为本专利技术再进一步的方本文档来自技高网...

【技术保护点】
1.一种基于FPGA的PML边界三维地震波传播模拟方法，其特征在于，将三维二次波动方程分解成x,y,z三分量，并引入PML边界条件进行求解，求解主体为波动方程的差分方程求解方式，所述基于FPGA的PML边界三维地震波传播模拟方法包括：/nS101：将所需变量从外部存储介质传入芯片内部；/nS102：求解差分方程的中间量B，计算对应xyz方向上分量，首先计算X方向上分量，计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；/na等于w1乘以p(i-(1/2),j,k,n)加上p(i+(1/2),j,k,n)的和，b等于w2乘以p(i-(3/2),j,k,n)加上p(i+(3/2),j,k,n)的和，c等于w3乘以p(i-(5/2),j,k,n)加上p(i+(5/2),j,k,n)的和，d等于w4乘以p(i-(7/2),j,k,n)加上p(i+(7/2),j,k,n)的和，直到m等于W(M/2)乘以p(i-((M-1)/2)),j,k,n)加上p(i+((M-1)/2),j,k,n))的和，其计算方程为：/n

【技术特征摘要】
1.一种基于FPGA的PML边界三维地震波传播模拟方法，其特征在于，将三维二次波动方程分解成x,y,z三分量，并引入PML边界条件进行求解，求解主体为波动方程的差分方程求解方式，所述基于FPGA的PML边界三维地震波传播模拟方法包括：
S101：将所需变量从外部存储介质传入芯片内部；
S102：求解差分方程的中间量B，计算对应xyz方向上分量，首先计算X方向上分量，计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；
a等于w1乘以p(i-(1/2),j,k,n)加上p(i+(1/2),j,k,n)的和，b等于w2乘以p(i-(3/2),j,k,n)加上p(i+(3/2),j,k,n)的和，c等于w3乘以p(i-(5/2),j,k,n)加上p(i+(5/2),j,k,n)的和，d等于w4乘以p(i-(7/2),j,k,n)加上p(i+(7/2),j,k,n)的和，直到m等于W(M/2)乘以p(i-((M-1)/2)),j,k,n)加上p(i+((M-1)/2),j,k,n))的和，其计算方程为：



S103：求解差分方程的中间量B，计算y分量值，计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；
a等于w1乘以p(i,j-(1/2),k,n)加上p(i,j+(1/2),k,n)的和，b等于w2乘以p(i,j-(3/2),k,n)加上p(i,j+(3/2),k,n)的和，c等于w3乘以p(i,j-(5/2),k,n)加上p(i,j+(5/2),k,n)的和，d等于w4乘以p(i,j-(7/2),k,n)加上p(i,j+(7/2),k,n)的和，直到m等于W(M/2)乘以p(i,j-((M-1)/2)),k,n)加上p(i,j+((M-1)/2),k,n))的和，其计算方程为：



S104：求解差分方程的中间量B，计算z分量值，计算如下a、b、c、d等共(M/2)+1项的和，项数随M的变化而变，M均为偶数；
a等于w1乘以p(i,j,k-(1/2),n)加上p(i,j,k+(1/2),n)的和，b等于w2乘以p(i,j,k-(3/2),n)加上p(i,j,k+(3/2),n)的和，c等于w3乘以p(i,j,k-(5/2),n)加上p(i,j,k+(5/2),n)的和，d等于w4乘以p(i,j,k-(7/2),n)加上p(i,j,k+(7/2),n)的和，直到m等于W(M/2)乘以p(i,j,k-((M-1)/2)),n)加上p(i,j,k+((M-1)/2),n))的和，其计算方程为：



S105：求解差分方程的P，a1等于w1乘以B(i-(1/2),j,k,n)加上B(i+(1/2),j,k,n)的和，b1等于w2乘以B(i-(3/2),j,k,n)加上B(i+(3/2),j,k,n)的和，c1等于w3乘以B(i-(5/2),j,k,...

【专利技术属性】
技术研发人员：陆敏旻，杨文珺，
申请(专利权)人：上海雪湖科技有限公司，
类型：发明
国别省市：上海;31