CN111859253B - 一种基于fpga的高阶波动方程求解方法 - Google Patents

一种基于fpga的高阶波动方程求解方法 Download PDF

Info

Publication number
CN111859253B
CN111859253B CN202010653228.6A CN202010653228A CN111859253B CN 111859253 B CN111859253 B CN 111859253B CN 202010653228 A CN202010653228 A CN 202010653228A CN 111859253 B CN111859253 B CN 111859253B
Authority
CN
China
Prior art keywords
plus
fpga
equals
wave equation
solving
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010653228.6A
Other languages
English (en)
Other versions
CN111859253A (zh
Inventor
高雨豪
田志明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Xuehu Technology Co ltd
Original Assignee
Shanghai Xuehu Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Xuehu Technology Co ltd filed Critical Shanghai Xuehu Technology Co ltd
Priority to CN202010653228.6A priority Critical patent/CN111859253B/zh
Publication of CN111859253A publication Critical patent/CN111859253A/zh
Application granted granted Critical
Publication of CN111859253B publication Critical patent/CN111859253B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)
  • Complex Calculations (AREA)

Abstract

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

Description

一种基于FPGA的高阶波动方程求解方法
技术领域
本发明涉及石油勘探、地质成像领域，具体是一种在地质成像过程中基于FPGA的高阶波动方程求解方法。
背景技术
随着油气田勘探程度的不断深化,普通易探明油气田几乎没有剩余，但油气资源又是现代工业所必须产品，故只能向结构复杂地域、海洋进行探索。现有呈像方法多使用克希霍夫射线法来进行呈像处理，此方法对盐下构造、复杂地质构造探照度不够，从而使用逆时偏移(RTM)方法来对地下地质结构进行呈像。这是一种基于波动方程在时间-空间域用高阶差分方方程求解声波偏微分方程的呈像方法，此方法可以真实的模拟波在地下的传播过程，无倾角限制，适合海洋及复杂地质结构，在复杂三维结构呈像方面有明显优势，呈像精度较射线法更高。
在油气地质勘探过程中使用逆时偏移方法进行呈像涉及到对波动方程的求解，此过程计算量巨大，现有技术方案使用的为GPU加速卡进行多核多线程，相较于传统使用CPU计算有较大的提升。
但在现今阶段GPU加速卡价格昂贵且功耗高，现有方案均是通过调用NVIDIA公司发布CUDA平台对高阶差分方法波动方程进行求解，由于CUDA平台为GPU通用性方案，且GPU硬件结构的限制，故没有针对高阶差分方程求解波动方程方法进行优化，从而造成GPU在求解波动方程时有一定局限性，在保留计算精度的情况下计算花费时间比较长。
本发明利用了FPGA硬件特性，针对高阶差分方程求解波动方程做了硬件上的特异化的优化，对方程进行了等价拆分、组合，从而达到流水的方式对方程进行计算求解,充分利用硬件资源，较现有GPU方案同精度下更快，功耗更低。
发明内容
本发明的目的在于提供一种基于FPGA的高阶波动方程求解方法，以解决上述背景技术中提出的问题。
为实现上述目的，本发明提供如下技术方案：
一种基于FPGA的高阶波动方程求解方法，将三维双程波动方程降为三个一维方程进行求解即分解成x,y,z三个分量，求解主体为波动方程的差分方程求解方式，所述基于FPGA的高阶波动方程求解方法包括：
S101：将需要计算的点坐标数据从硬盘或其他存储介质中传入FPGA内存中，FPGA芯片将数据从内存中读取到芯片内部；
S102：将传至芯片内部的数据按照一定的顺序进行组合计算，按照差分方程计算公式，需读取P(i,j,k,n)，p(i,j,k,n-1)，p(i-1,j,k,n)、p(i+1,j,k,n)，p(i-2,j,k,n)、p(i+2,j,k,n),...,p(i-m,j,k,n)、p(i+m,j,k,n)，p(i,j-1,k,n)、p(i,j+1,k,n)，p(i,j-2,k,n)、p(i,j+2,k,n),...,p(i,j-m,k,n)、p(i,j+m,k,n)，p(i,j,k-1,n)、p(i,j,k+1,n)，p(i,j,k-2,n)、p(i,j,k+2,n),...,p(i,j,k-m,n)、p(i,j,k+m,n)和对应的v(i,j,k)处的速度场值，m为所计算阶数M/2，M均为偶数；
S103：将对应值乘以对应差分系数进行相加，a等于w0乘以P(i,j,k,n)加上w1乘以p(i-1,j,k,n)加上p(i+1,j,k,n)的和，b等于w2乘以p(i-2,j,k,n)加上p(i+2,j,k,n)的和，...,c等于Wm乘以p(i-m,j,k,n)加上p(i+m,j,k,n)的和,d等于w0乘以P(i,j,k,n)加上w1乘以p(i,j-1,k,n)加上p(i,j+1,k,n)的和，e等于w2乘以p(i,j-2,k,n)加上p(i,j+2,k,n)的和，...,f等于Wm乘以p(i,j-m,k,n)加上p(i,j+m,k,n)的和,g等于w0乘以P(i,j,k,n)加上w1乘以p(i,j,k-1,n)加上p(i,j,k+1,n)的和，h等于w2乘以p(i,j,k-2,n)加上p(i,j,k+2,n)的和，...,l等于Wm乘以p(i,j,k-m,n)加上p(i,j,k+m,n)的和；
S104：对上诉计算得出的结果分别乘以对应速度场在其所对应分量上的系数，
XX乘以a加上b加上，...，c的和，再加上YY乘以d加上e加上，...，f的和,再加上ZZ乘以g加上h加上l的和，最后再加上2倍的p(i,j,k,n)，减去p(i,j,k,n-1)，即可得出最终p(i,j,k,n+1)的值，此点下一时刻的值；
S105：重复步骤二、步骤三、步骤四，直到将所有当前时刻点遍历完毕，得到下一时刻波场数值P(n+1)；
S106：基于P(n+1)波场值，重复步骤四，直到遍历完成所有时刻的值，得到最终P(nmax)值。
作为本发明进一步的方案：所述波动方程为：
其中，p(x,y,z,t)为波场，v(x,y,z)为介质速度。
作为本发明再进一步的方案：所述差分方程为：
其中，M为阶数。
与现有技术相比，本发明的有益效果是：
1、本发明针对FPGA设备应用，与以往使用CPU，GPU等冯诺依曼架构计算单元相比，FPGA具有更高的可定制性，以及更高的并发性。
2、FPGA器件属于专用集成电路中的一种半定制电路，是可编程的逻辑列阵，能够有效的解决原有的器件门电路数不足的问题。FPGA的基本结构包括可编程输入输出单元，可配置逻辑块，数字时钟管理模块，嵌入式块RAM，布线资源，内嵌专用硬核，底层内嵌功能单元。由于FPGA具有布线资源丰富，可重复编程和集成度高。
3、基于FPGA的可定制，可编程的特性，我们针对算法求解过程做了等价变形，使得在计算时充分发挥FPGA芯片的并发性能，使得计算指令无需等待，以流水的方式实现无缝连接的计算，省去了CPU、GPU等传统冯诺依曼架构芯片的指令调度、指令执行等待、数据搬运时间。从而大大加速了计算过程。据统计比同价位GPU快5倍以上，较CPU快百倍以上。
4、同时，市场现存的GPU有较大的内存限制，导致数据无法装载完全，从而限制GPU运算速度，或者针对某些规模数据无法计算。而当前FPGA实现方式通过流水来提高并发量，有效缓解了内存的压力，使得在FPGA上可实现大规模数据的计算
附图说明
图1为一种基于FPGA的高阶波动方程求解方法的原理图。
图2为一种基于FPGA的高阶波动方程求解方法的流程图。
图中：Pn：当前时间切片波长数据，Pn-1：前一时间切片波长数据，V：速度参数场，wm：差分系数，Pn+1：后一时间切片波长数据
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
请参阅图1～2，本发明实施例中，一种基于FPGA的高阶波动方程求解方法，将三维双程波动方程降为三个一维方程进行求解即分解成x,y,z三个分量，求解主体为波动方程的差分方程求解方式，所述基于FPGA的高阶波动方程求解方法包括：
S101：将需要计算的点坐标数据从硬盘或其他存储介质中传入FPGA内存中，FPGA芯片将数据从内存中读取到芯片内部；
S102：将传至芯片内部的数据按照一定的顺序进行组合计算，按照差分方程计算公式，需读取P(i,j,k,n)，p(i,j,k,n-1)，p(i-1,j,k,n)、p(i+1,j,k,n)，p(i-2,j,k,n)、p(i+2,j,k,n),...,p(i-m,j,k,n)、p(i+m,j,k,n)，p(i,j-1,k,n)、p(i,j+1,k,n)，p(i,j-2,k,n)、p(i,j+2,k,n),...,p(i,j-m,k,n)、p(i,j+m,k,n)，p(i,j,k-1,n)、p(i,j,k+1,n)，p(i,j,k-2,n)、p(i,j,k+2,n),...,p(i,j,k-m,n)、p(i,j,k+m,n)和对应的v(i,j,k)处的速度场值，m为所计算阶数M/2，M均为偶数；
S103：将对应值乘以对应差分系数进行相加，a等于w0乘以P(i,j,k,n)加上w1乘以p(i-1,j,k,n)加上p(i+1,j,k,n)的和，b等于w2乘以p(i-2,j,k,n)加上p(i+2,j,k,n)的和，...,c等于Wm乘以p(i-m,j,k,n)加上p(i+m,j,k,n)的和,d等于w0乘以P(i,j,k,n)加上w1乘以p(i,j-1,k,n)加上p(i,j+1,k,n)的和，e等于w2乘以p(i,j-2,k,n)加上p(i,j+2,k,n)的和，...,f等于Wm乘以p(i,j-m,k,n)加上p(i,j+m,k,n)的和,g等于w0乘以P(i,j,k,n)加上w1乘以p(i,j,k-1,n)加上p(i,j,k+1,n)的和，h等于w2乘以p(i,j,k-2,n)加上p(i,j,k+2,n)的和，...,l等于Wm乘以p(i,j,k-m,n)加上p(i,j,k+m,n)的和；
S104：对上诉计算得出的结果分别乘以对应速度场在其所对应分量上的系数，
XX乘以a加上b加上，...，c的和，再加上YY乘以d加上e加上，...，f的和,再加上ZZ乘以g加上h加上l的和，最后再加上2倍的p(i,j,k,n)，减去p(i,j,k,n-1)，即可得出最终p(i,j,k,n+1)的值，此点下一时刻的值；
S105：重复步骤二、步骤三、步骤四，直到将所有当前时刻点遍历完毕，得到下一时刻波场数值P(n+1)；
S106：基于P(n+1)波场值，重复步骤四，直到遍历完成所有时刻的值，得到最终P(nmax)值。
本发明的一具体实施例中，所述波动方程为：
其中，p(x,y,z,t)为波场，v(x,y,z)为介质速度；
所述差分方程为：
其中，M为阶数。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (3)

1.一种基于FPGA的高阶波动方程求解方法，其特征在于，将三维双程波动方程降为三个一维方程进行求解即分解成x,y,z三个分量，求解主体为波动方程的差分方程求解方式，所述基于FPGA的高阶波动方程求解方法包括：
S101：将需要计算的点坐标数据从硬盘或其他存储介质中传入FPGA内存中，FPGA芯片将数据从内存中读取到芯片内部；
S102：将传至芯片内部的数据按照顺序进行组合计算，按照差分方程计算公式，需读取P(i,j,k,n)，p(i,j,k,n-1)，p(i-1,j,k,n)、p(i+1,j,k,n)，p(i-2,j,k,n)、p(i+2,j,k,n),...,p(i-m,j,k,n)、p(i+m,j,k,n)，p(i,j-1,k,n)、p(i,j+1,k,n)，p(i,j-2,k,n)、p(i,j+2,k,n),...,p(i,j-m,k,n)、p(i,j+m,k,n)，p(i,j,k-1,n)、p(i,j,k+1,n)，p(i,j,k-2,n)、p(i,j,k+2,n),...,p(i,j,k-m,n)、p(i,j,k+m,n)和对应的v(i,j,k)处的速度场值，m为所计算阶数M/2，M均为偶数；
S103：将对应值乘以对应差分系数进行相加，a等于w0乘以P(i,j,k,n)加上w1乘以p(i-1,j,k,n)加上p(i+1,j,k,n)的和，b等于w2乘以p(i-2,j,k,n)加上p(i+2,j,k,n)的和，...,c等于Wm乘以p(i-m,j,k,n)加上p(i+m,j,k,n)的和,d等于w0乘以P(i,j,k,n)加上w1乘以p(i,j-1,k,n)加上p(i,j+1,k,n)的和，e等于w2乘以p(i,j-2,k,n)加上p(i,j+2,k,n)的和，...,f等于Wm乘以p(i,j-m,k,n)加上p(i,j+m,k,n)的和,g等于w0乘以P(i,j,k,n)加上w1乘以p(i,j,k-1,n)加上p(i,j,k+1,n)的和，h等于w2乘以p(i,j,k-2,n)加上p(i,j,k+2,n)的和，...,l等于Wm乘以p(i,j,k-m,n)加上p(i,j,k+m,n)的和；
S104：对上述计算得出的结果分别乘以对应速度场在其所对应分量上的系数；
XX乘以a加上b加上，...，c的和，再加上YY乘以d加上e加上，...，f的和,再加上ZZ乘以g加上h加上l的和，最后再加上2倍的p(i,j,k,n)，减去p(i,j,k,n-1)，即可得出最终p(i,j,k,n+1)的值，此点下一时刻的值；
S105：重复步骤S102、步骤S103、步骤S104，直到将所有当前时刻点遍历完毕，得到下一时刻波场数值P(n+1)；
S106：基于P(n+1)波场值，重复步骤S104，直到遍历完成所有时刻的值，得到最终P(nmax)值。
2.根据权利要求1所述的一种基于FPGA的高阶波动方程求解方法，其特征在于，所述波动方程为：
其中，p(x,y,z,t)为波场，v(x,y,z)为介质速度。
3.根据权利要求1所述的一种基于FPGA的高阶波动方程求解方法，其特征在于，所述差分方程为：
其中，M为阶数。
CN202010653228.6A 2020-07-08 2020-07-08 一种基于fpga的高阶波动方程求解方法 Active CN111859253B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010653228.6A CN111859253B (zh) 2020-07-08 2020-07-08 一种基于fpga的高阶波动方程求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010653228.6A CN111859253B (zh) 2020-07-08 2020-07-08 一种基于fpga的高阶波动方程求解方法

Publications (2)

Publication Number Publication Date
CN111859253A CN111859253A (zh) 2020-10-30
CN111859253B true CN111859253B (zh) 2023-09-22

Family

ID=73153234

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010653228.6A Active CN111859253B (zh) 2020-07-08 2020-07-08 一种基于fpga的高阶波动方程求解方法

Country Status (1)

Country Link
CN (1) CN111859253B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1271783A2 (en) * 2001-06-29 2003-01-02 STMicroelectronics Ltd. FPGA with a simplified interface between the program memory and the programmable logic blocks
US7573770B1 (en) * 2007-07-16 2009-08-11 Lattice Semiconductor Corporation Distributed front-end FIFO for source-synchronized interfaces with non-continuous clocks
CA2797434A1 (en) * 2010-05-12 2011-11-17 Shell Internationale Research Maatschappij B.V. Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
WO2015042386A1 (en) * 2013-09-20 2015-03-26 Westerngeco Llc Eikonal solver for quasi p-waves in anisotropic media
CN108566205A (zh) * 2018-04-19 2018-09-21 南通大学 一种基于fpga实现的d/a转换器

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101172506B1 (ko) * 2010-07-30 2012-08-10 서울대학교산학협력단 균일 분포된 가상 송신원을 이용한 지하 구조의 영상화 방법
US8898480B2 (en) * 2012-06-20 2014-11-25 Microsoft Corporation Managing use of a field programmable gate array with reprogammable cryptographic operations
US9230091B2 (en) * 2012-06-20 2016-01-05 Microsoft Technology Licensing, Llc Managing use of a field programmable gate array with isolated components

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1271783A2 (en) * 2001-06-29 2003-01-02 STMicroelectronics Ltd. FPGA with a simplified interface between the program memory and the programmable logic blocks
US7573770B1 (en) * 2007-07-16 2009-08-11 Lattice Semiconductor Corporation Distributed front-end FIFO for source-synchronized interfaces with non-continuous clocks
CA2797434A1 (en) * 2010-05-12 2011-11-17 Shell Internationale Research Maatschappij B.V. Seismic p-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
WO2015042386A1 (en) * 2013-09-20 2015-03-26 Westerngeco Llc Eikonal solver for quasi p-waves in anisotropic media
CN108566205A (zh) * 2018-04-19 2018-09-21 南通大学 一种基于fpga实现的d/a转换器

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
E. Motuk 等.Implementation of finite difference schemes for the wave equation on FPGA.《 Proceedings. (ICASSP '05). IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005.》.2005,237-240. *
叠前逆时深度偏移成像的实现与应用;刘定进;杨勤勇;方伍宝;孔祥宁;张慧宇;;石油物探(第06期);545-549+527 *
基于CPU/GPU协同加速叠前逆时偏移方法研究;高新成;李春生;;陕西理工学院学报(自然科学版)(第01期);44-49 *
基于GPU并行加速的逆时偏移成像方法;张慧;蔡其新;秦广胜;高爱荣;;石油地球物理勘探(第05期);707-710+670+854 *
基于时空域自适应高阶有限差分的声波叠前逆时偏移;严红勇;刘洋;;地球物理学报(第03期);971-984 *
格林函数法求解含有高阶时间导数的波动方程的初边值问题;阚洪弟;《河南教育学院学报(自然科学版)》;第28卷(第01期);47-49+58 *

Also Published As

Publication number Publication date
CN111859253A (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
Guo et al. A survey of FPGA-based neural network accelerator
Li et al. A method of using GPU to accelerate seismic pre‐stack time migration
WO2011147075A1 (en) Optimiztion-based simulated annealing for integrated circuit placement
CN111814332B (zh) 一种基于fpga的pml边界三维地震波传播模拟方法
CN111859253B (zh) 一种基于fpga的高阶波动方程求解方法
CN111563582A (zh) 一种在fpga上实现及优化加速卷积神经网络的方法
Ceyhan et al. Towards a faster and improved ADCIRC (ADvanced Multi-Dimensional CIRCulation) model
CN108038262A (zh) 一种考虑sssi效应的楼层反应谱简化计算方法
CN109725346B (zh) 一种时间-空间高阶vti矩形网格有限差分方法及装置
CN109583006B (zh) 一种基于循环切割和重排的现场可编程门阵列卷积层的动态优化方法
CN115201896B (zh) 吸收衰减介质逆时偏移方法、装置、成像方法、介质
CN110765715A (zh) 一种面向gpu芯片渲染输出单元性能仿真方法及平台
CN102314215A (zh) 集成电路系统中小数乘法器的低功耗优化方法
Li et al. A power-efficient optimizing framework FPGA accelerator for YOLO
CN105893326B (zh) 基于fpga实现65536点fft的装置和方法
CN112596105A (zh) 一种基于fpga的地震波正向延拓方法
CN112862080B (zh) EfficientNet的注意力机制的硬件计算方法
CN115293978A (zh) 卷积运算电路和方法、图像处理设备
CN104156606B (zh) 基于混合计算架构的用于核素迁移控制方程的求解方法
Zhang et al. Design of a Convolutional Neural Network Accelerator based on PYNQ
Meng et al. GPU parallel acceleration of transient simulations of open channel and pipe combined flows
CN113791445A (zh) 一种基于fpga全波形反演的频率成分波场信息提取方法
Chen et al. eSSpMV: An Embedded-FPGA-based Hardware Accelerator for Symmetric Sparse Matrix-Vector Multiplication
Zhou et al. Time domain boundary element method for semi-infinite domain problems using CSR storage method
Xu et al. Towards Energy-Efficient Asynchronous Circuit Design with Flip-Flop-to-Latch Replacement

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant