TFInterpy: 高性能空间插值Python包

软件X 20（2022）101229原始软件出版物TFInterpy：高性能空间插值Python包陈志文，钟宝荣长江大学计算机科学学院，荆州，434000ar t i cl e i nf o文章历史记录：收到2022年收到修订版，2022年7月15日接受，2022年MSC：65D0586A32关键词：PythonTensorFlow插值Kriginga b st ra ct插值算法是空间分析的重要工具。克里格法满足最佳线性无偏估计，广泛应用于需要高精度的场景。但是，当使用克里格方法处理大数据集时，程序运行时间可能会长得令人无法接受。为了解决这个问题，我们基于TensorFlow框架开发了TFInterpy这个Python包为插值算法提供了一个开源、跨平台、易于使用的API，并在应用于大规模任务时实现了显著的加速版权所有©2022作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-22-00084可复制胶囊的永久链接GNU通用公共许可证，第3版（GPL v3）使用git的代码版本控制系统使用Python 3的软件代码语言、工具和服务编译要求，操作环境依赖Python3，TensorFlow，Keras，NumPy，Pandas，Scipy，VTK，PyQT 5如果可用开发人员文档/手册链接问题支持电子邮件orchenz@qq.com1. 动机和意义在大多数情况下，空间分析只能收集少量的样本点。为了建立目标区域的视觉模型，有必要使用空间插值技术来估计未采样位置的相关属性值反距离加权法（IDW）和克里格法（Kriging）在环境科学[1]、地质勘探[2]、水文地质[3]等领域有着广泛的应用IDW方法是最简单的插值方法，只考虑样本点与目标点之间的距离，而样本点的权值对目标点的贡献与远处IDW也可以看作是一种双曲核径向基函数方法。Kriging方法是*通讯作者。电子邮件地址：barry@yangtzeu.edu.cn（钟宝荣）.https://doi.org/10.1016/j.softx.2022.101229作为最优无偏估计，它可以使方差最小化[4]。与IDW相比，Kriging方法考虑了更复杂的空间相关性。它使用变差函数来拟合距离和差异之间的相关性与其他算法相比，该算法具有更高的插值精度。克里格法包含一系列衍生自普通克里格法（OK）的简单克里金法（SK）、指示克里金法（IK）和协同克里金法（CK）。通常，这些算法需要求解一组线性方程，称为克里格方程系统。求解克里格方程系统涉及矩阵乘法和矩阵求逆运算，并且在此之前计算变差函数的值。因此，当要估计的目标点的数量很大时，该程序可能非常耗时，特别是在3D建模中，目标点的数量可能达到数百万或更多。因此，提高Kriging算法的执行速度一直是Kriging算法研究的热点。2352-7110/©2022作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx陈志文、钟宝荣软件X 20（2022）1012292∑我--⎦⎣1∑好⎦⎢⎣⎦..关于提高插补过程的速度，已经有许多研究提出了解决方案。Ferreira 等 人 [5] 使 用 OpenMP 编 程 开 发 了 具 有 可 并 行 计 算 的Kriging。Hu和Shu [6]建立了一个基于MPI的分布式存储系统，并实现了普通克里格（OK）的并行处理。Cheng [7]，Gutiérrez de Ravé等人[8]，Liu [9]相继提出了基于CUDA编程的Kriging用于GPU加速。Misra等人。[10]实现了一个基于Spark框架的分布式运行OK。所有上述方案都导致了运行速度的显著提高，但遗憾的是，这些研究都没有将其源代码开源。此外，上述解决方案需要很高的编程技能，这与克里格算法本身的复杂性一起，使得再现这些解决方案非常困难。Python语言有一个成熟的软件ecosys-TEM和许多其他优点，使其越来越多地使用在科学计算和工程应用中。Krig- ing算法是空间分析的重要工具，有Python包为其提供API，包括Geostat- sPy [11]和PyKrige [12]。GeostatsPy的目标是使用Python语言实现GSLIB [13目前，GeostatsPy仅包含二维空间的Kriging算法，数值运算依赖于GSLIB提供的可执行文件。PyKrige提供各种Kriging算法，包括用于3D插值的接口，但不提供用于3D可视化的工具。为了解决上述问题，我们开发了TFInterpy（ 基 于 TensorFlow 的 插 值 Python 包 ） 。 TFInterpy 包 含 基 于TensorFlow的IDW、OK和SK（简单克里金法）实现，这意味着插值算法可以通过TensorFlow并行实现，并可以使用GPU设备的计算能力。TFInterpy还提供了这些算法的常规版本。TFInterpy是一个建立在许多高度可移植的Python包之上的包，因此TFInterpy还具有高度可移植性，允许用户创建2.1. 软件构架TFInterpy依赖于许多优秀的开源Python包。NumPy [15]用于数据操作; Pandas [16]用于数据导入和导出; SciPy [17]中的KD树用于邻域搜索，最小二乘法用于拟合变差函数; TensorFlow用于实现可并行插值算法; Matplotlib [18]用于2D绘图; VTK [19]用于3D可视化和PyQT 5 [20]用于构建GUI工具。表1列出了TFIn- terpy的主要组件、相应的功能描述以及所需的第三方依赖库。图1显示了TFInterpy中组件之间的依赖关系。图1中的package形状表示Python中的包或模块，包含子项的package形状表示包;否则，它表示模块。图1中的箭头表示依赖关系;例如，指向“utils”的“krige "表示”krige"依赖于“utils”。2.2. 软件功能2.2.1. 插值空间插值可以表示为观测值的加权组合Eq. 其中n表示邻域点的数目（即，邻域大小），Z是目标点的估计结果，Zi是第i个邻域点的属性值，λi表示第i个邻域点（采样点）的影响。IDW的权重计算规则为Eq. （2），其中hi表示第i个采样点与目标点之间的欧几里德距离，指数α表示该距离的影响程度，α值越高将赋予越低的权值。α的值通常取为2。nZ=λiZi（1）i=1 1/hα简单易用是TFInterpy的设计目标之一它提供默认λi=ni=1我（1/hα）（二）用于计算变异函数的函数，这可以降低非专业人员的使用门槛。同时，TFInterpy还提供基本的2D和3D数据可视化解决方案，让用户即时了解插值结果。此外，TFInterpy下的GUI子包提供了一个基于其他核心模块构建的2. 软件描述TFInterpy设计用于仅考虑相邻样本点的局部插值，但也可以通过调整邻域搜索大小来实现全局插值。TFInterpy提供了各种可能用于空间分析过程的功能，包括数据加载、插值计算、数据可视化和数据导出。模块之间的耦合保持在低水平，因此如果用户只需要使用部分函数，则不需要安装所有依赖项。例如，不需要安装TensorFlow和Keras [14]（TensorFlow后端）假设有目标点x和n个邻域点x1，x2，. . .，xn，普通克里格方程系统为Eq. 其中，μ是为了使误差最小化而引入的拉格朗日乘数（即，Kriging方差）[4]，而Ordinary Kriging方差的公式为Eq.（四）、OK方法的加权组合的计算可以参考等式（Eq. （1），和Zi在方程。（1）对应于与邻域点xi对应的属性值。简单克里格方程系统是Eq。（6），简单克里金方差的公式为Eq. 并且OK方法的加权组合的计算可以参考Eq. （五）、在Eq中的m（5）是区域化变量Z（x）的数学期望[4]，m取实际应用中样本点属性的平均值。Eq之间的差异。（3）和（6）是Eq。（3）将权重之和λ约束为1。方程中的协方差C（x， y）（3）、（4）、（6）和（7）是通过计算变差函数得到的。C（x1，x1）· · ·C（x1，xn）1λ1C（x，x1）当仅使用传统的Inter-.⎣. . .... ⎥⎢...=⎥polation算法与PyKrige和GeostatsPy，TFInterpy中的插值函数不要使用结构化网格的坐标和尺寸为了保持最大的灵活性，TFInterpy使用目标点坐标作为参数。分散，结构化C（x n，x1）· · ·C（xn，xn）11···1 0nλnC（x，xn）−µ（三）网格和非结构化网格可以传递目标点坐标到函数界面，得到插补结果。σ2=C（x，x）−∑λiC（x，xi）+μ（4）i=1陈志文、钟宝荣软件X 20（2022）1012293⎦⎣⎣..C（x）x）的· · ·C（x x）∑σ=C（x，x）−λC（x，x） （7）ii.C（xx）表1TFInterpy的主要组成部分组件类型描述依赖关系网格模块用于建立直线网格numpyIDW模块IDW算法麻木的，麻木的Krige模块SK和OK算法麻木的，麻木的变异函数模块用于计算变差函数麻木的，麻木的变差曲线模块用于计算实验变差函数numpy变差图模型模块变差函数模型numpyutils模块包含辅助函数numpyTF分装上面列出的算法的TensorFlow实现numpy，scipy，tensorflowGSLIB分装用于读写GSLIB文件熊猫VTK分装用于基于VTK的数据渲染和模型导出numpy，vtkGUI分装基于上述所有功能的GUI工具matplotlib，numpy，pandas，PyQt5、scipy、vtkFig. 1. TFInterpy中组件之间的相互关系n n意味着接收的数据大小是任意的，其对应于Z=∑λiZi+m（1−∑λi）（5）i=1i=1TensorFlow中的“批量大小”概念。在等式的左侧的矩阵。（6），C（x1，x1）· · ·C（x1，xn）λ1C（x，x1）.右边的Eq。（6），并且采样点。好吧n，. . ..n，n好吧n=、中国（6）变差函数将位置偏移映射到半方差，2SKi=1IDW和Kriging方法都被封装为插值器类，每个插值器类都提供一个名为“execute”的函数调用“execute”函数的最低参数要求最近邻点数、变差函数等参数可选，函数将返回插值结果。基于TensorFlow的类在TFInterpy中以“TF”命名，包括TFIDW、TFSK和TFOK。这些类提供了与传统的实现类一致的函数接口。这些 类 的 实 现 原 理 是 使 用 Keras 的 函 数 API 将 计 算 部 分 构 建 为TensorFlow数据流图，它表示Keras的“Model”对象。例如图2显示结构TFSK的数据流图形状元组中的该公式在Eq. 其中，Z（x）表示以下性质：位置x处的区域化变量，h是样本点之间在欧几里得距离和方向上的间隔[4]，其可以被视为向量参数。在地质统计学中的二阶平稳性假设下，两点之间的协方差仅与距离和方向有关[4]，并且方程中的C（xi，xj）。（3）和（6）可以表示为表示为C（hi，j），并且hi，j表示xi和xj之间的位置偏移。 如等式1所示。在二阶平稳性假设[4]下，协方差函数和变差函数是等价的，即C（hi，j）可以表示为γ（hi，j）。1γ（h）=2var（Z（x）−Z（x+h））（8）γ（h）=C（0）−C（h）（9）变差函数的计算涉及到实验变差函数点对的寻找和理论变差函数模型的拟合，涉及到许多复杂的概念[4]。TFInterpy提供了相应的函数⎦2.2.2. 变差函数计算1λnn陈志文、钟宝荣软件X 20（2022）1012294×图三. 代码片段1.图二. TFSK的数据流图。第2行使用采样点计算嵌套变差函数;第4行代码使用采样点构造OK插值;第6行代码计算插值结果和克里金方差;第8、9行代码使用网格数据构造图中的代码片段。5展示了如何使用TFOK插件。辅助函数可以将嵌套的变差函数转换行中的batch_size10指定为并行操作传入的批处理数。它不应该将batch_size设置得太小;否则，程序将不会得到加速，甚至比常规实现慢。因此batch_size的默认值设置为10000。3. 说明性实例来处理这些步骤。然而，这些功能需要用户设置许多参数，非专家无法理解如何正确设置它们以及这些参数将如何影响插值。因此，TFInterpy还提供了用于计算变差函数的默认函数该过程基于散射的空间分布特性来推断每个参数的值。将表示方向的角度范围划分为多个部分，并使用具有最多对的方向上的距离方差对来拟合变差函数。这种方法不能保证最优建模，但允许用户根据第一次插值的结果来判断如何调整参数2.3. 示例代码段分析图中的代码片段。3展示了如何使用IDW算法。第3行的代码使用样本点构造了一个常规的IDW验证;第5行的代码执行留一交叉验证并打印结果;第7行的代码调用了IDW验证来估计网格点;第11行的代码基于TensorFlow实现构造了一个TFIDW验证，第13行的代码指定了TensorFlow使用GPU设备。在图1所示的代码片段中，如图3所示，变量图中的代码片段。 4显示了如何使用带有默认全方向变差函数的OK函数。中的代码本节显示2D和3D场景插值的示例我们已经在Windows和Linux操作系统上验证了GUI工具示例。硬件环境：AMD锐龙75800H CPU，16 GB RAM，NVIDIA GeForce RTX 3050笔记本GPU。3.1. 2D场景以下示例显示了使用大规模高分辨率数字高程模型的子地图进行2D插值的结果[21]。一个TFSK是通过随机选择100个样本点的原始数据，并构建一个2D网格一致的原始地图大小（100 100）进行插值。图6（a）显示了原始高程数据的子图，黑点是样本点，6（b）显示了估计结果，图。图6（c）显示了估计值和原始数据之间的绝对误差，图6（d）显示了相应的克里格方差。3.2. GUI工具在这个例子中，直线网格3D插值的结果显示了SGeMS [22]教程示例中的数据。界面左侧边栏上的工具允许用户导入、插入、显示和导出数据。用户可以选中“使用tensorflow”复选框，以决定是否使用在TensorFlow中实现的算法。图图7显示了每个方向上变差函数的拟合图。图8示出了使用TFOK执行留一法交叉验证的结果，并且图8示出了使用TFOK执行留一法交叉验证的结果。9示出了使用TFOK执行插值的结果。陈志文、钟宝荣软件X 20（2022）1012295见图4。 代码片段2.图五. 代码片段3.图六、使用TFSK执行2D插值。3.3. 性能比较我们使用不同的网格粒度来比较OK在不同实现中的执行时间，包括Geostat-sPy，PyKrige和TFInterpy。结果如表2所示，在TFInterpy中实现OK。可以看出，随着目标点数量的增加，TFInterpy-TFOK的性能优势变得越来越与GeostatsPy相比，TFOK（使用CPU）可以实现20倍以上的性能提升，与PyKrige相比，它可以实现高达10倍的性能提升由于从计算机内存中复制数据的时间很长，因此TFOK陈志文、钟宝荣软件X 20（2022）1012296见图7。 变异函数拟合结果。图8.第八条。TFO K 的留一法交叉验证结果。在实验环境下，对显卡内存和CPU与GPU的计算4. 影响在以往的研究中，插值算法的加速主要围绕并行执行、利用GPU的计算能力、分布式计算等方面展开，并且经常使用可移植性不高的编程技术。TFInterpy基于Python语言和TensorFlow框架开发，提供完全不同的解决方案。我们可以通过指定设备名称来使用CPU，GPU 甚 至 TPU 设 备 ， 并 且 它 具 有 分 布 式 计 算 的 潜 力 。 此 外 ，TFInterpy具有高度的可移植性和可扩展性。TFInterpy是一个开源的跨平台软件包，与目前可用的类似Python库相比，它具有更高的计算效率和更好的可视化TFInterpy是灵活的，定义了一个默认的算法执行过程，同时保持参数可定制。它为用户提供了一个简洁的API，这大大降低了非专业人士陈志文、钟宝荣软件X 20（2022）1012297见图9。 TFOK表2每个时间间隔的运行时间适用于不同的网格大小（单位：秒）。1× 105≈ 317× 3171× 106= 1000× 1000 2011.351 121.711 82.397 45.737 11.6831× 1073163× 3163 2784.843 1250.980 849.974 452.567 112.331使用.在潜在用户中，地质学家最能欣赏TFInterpy带来的性能提升，因为他们的空间分析任务通常需要处理更大的数据集，这更有可能利用TensorFlow5. 结论空间插值在许多学科中都是必不可少的，而克里格插值方法可以获得很高的插值精度，但计算量很大。虽然以前的许多研究已经描述了如何减少程序运行时间，但仍然没有开源的解决方案。因此，我们提出了TFInterpy，它实现了基于TensorFlow的算法加速。本文介绍了TFInterpy的各种特性，这些特性为空间插值提供了必要和辅助的API支持，并在现有Python包的基础上提供了显著的性能改进在未来的工作中，更多的插值算法将添加到TFInterpy，我们将研究分布式执行的可行性。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作数据可用性数据将根据要求提供。确认本研究得到了国家科技重大专项前期项目的资助，项目编号为2016ZX05055。引用[1]Metahni S，Coudert L，Gloaguen E，Guemiza K，Mercier G，Blais J-F.比较不同的插值方法和序贯高斯模拟来估算土壤中As、Cr、Cu、PCP和二恶英/呋喃的污染量。环境污染2019;252：409-19.http://dx.doi.org/10.1016/j.envpol.2019.05.122，URLhttps://www.sciencedirect.com/science/article/pii/S0269749119306311。[2] Khaba L，Griffiths JA.莱索托北部Muela水库泥沙淤积造成的水库容量损失 计 算 。Int Soil Water Conserve Res 2017;5 （ 2 ） ： 130-40.http://dx.doi.org/10.1016/j的网站。iswcr.2017.05.005，URLhttps://www.sciencedirect.com/science/article/pii/S2095633916300041。[3]Belkhiri L，Tiri A，Mouni L.利用Kriging和Co-Kriging插值的地下水质量空间分布GroundwSustainDev2020;11：100473.http://dx.doi.org/10.1016/j.gsd.2020.100473，网址https：//www.sciencedirect.com/science/article/pii/S2352801X20300461网站。[4] Oliver M，Webster R.地质统计学教程指南：计算和建模变异函数和克里金法。Catena2014;113：56-69.http://dx.doi.org/10.1016/j.catena.2013.09.006，URLhttps://www.sciencedirect.com/science/article/pii/S0341816213002385。[5]Ferreira A，Ponte F，Varella L. GSLIB库的KT3D应用程序使用英特尔至强融核处理器和OpenMP进行并行化。2014年第76届EAGE会议和展览。 2014，（ 1 ） ： 欧 洲相 关 ，地 球 科 学 家 工 程 师 的 作 用 ;2014 ， p 。 1 比 5 。http://dx.doi.org/10。网格大小插值器GeostatsPy-OKPyKrige-OKTFInterpy-OKTFInterpy-TFOK（GPU）TFInterpy-TFOK（CPU）1× 104= 100× 10023.977230.2991.25812.2640.8288.1402.0706.2390.9792.067陈志文、钟宝荣软件X 20（2022）10122983997/2214-4609.20141291，URLhttps://www.earthdoc.org/content/papers/10.3997/2214-4609.20141291。[6] Hu H，Shu H.一种改进的粗粒度并行算法，用于普通Kriging插值的计算加速。ComputGeosci2015;78：44-52.http://dx.doi.org/10.1016/j.cageo.2015.02.011，网址https：//www.sciencedirect.com/science/article/pii/S0098300415000333网站。[7] 郑 T 利 用 支 持 CUDA 的 GPU 加 速 通 用 Kriging 插 值 算 法 Comput Geosci2013;54：178-83.http://dx.doi.org/10.1016/j.cageo.2012.11.013，URLhttps://www.sciencedirect.com/science/article/pii/S0098300412003901。[8] Gutiérrez de Ravé E ， Jiménez-Hornero F ， Ariza-Villaverde A ， Gómez-López J.使用图形处理单元（GPGPU）上的通用计算来加速普通克里金算法。Comput Geosci 2014;64：1-6. http://dx.doi.org/10.1016/j.cageo.2013.11.004，URLwww.sciencedirect.com/science/article/pii/S0098300413002914。[9] 刘毅基于gpu的三维克里格插值算法研究与开发（硕士论文），电子科技大学;2021年。[10]张文，王文，等.一种基于spark框架的分布式克立格算法.北京：计算机科学出版社，2000，24（3）：117 - 118 Int J Data Sci Anal2020;10（3）：249-64.[11]Pyrcz M ， Jo H ， Kupenko A ， Liu W ， Gigliotti A ， Salomaki T ， et al.GeostatsPy. 2021年，https://github.com/GeostatsGuy/GeostatsPy。[12] Murphy B，Müller S，Yurchak R. GeoStat-framework/PyKrige：v1.6.1.Zenodo; 2021年，http://dx.doi.org/10.5281/zenodo.5380342。[13]Deutsch CV，Zynnel AG. Gslib。地统计软件库和用户指南。卷三六九纽约：牛津大学出版社，1998.[14]Chollet F，et al. Keras. 2015年，https://keras.io。[15]Harris CR ， Millman KJ ， Van Der Walt SJ ， Gommers R ， Virtanen P ，Cournapeau D ， et al. Array programming with NumPy. Nature2020;585（7825）：357-62.[16]Pandas：a foundational Python library for data analysisand statistics. In：Python用于高性能和科学计算。 卷十四岁（9）：Seattle; 2011，p.一比九[17]Virtanen P ， Gommers R ， Oliphant T ， Haberland M ， Reddy T ，Cournapeau D ， et al. SciPy 1.0 ： fundamental algorithms for scientificcomputinginPython.自然方法2020;17：1-12.http://dx.doi.org/10.1038/s41592-019-0686-2.[18]亨特JD。Matplotlib：2D图形环境。计算机科学与工程2007;9（03）：90-5.[19]Hanwell MD ， Martin KM ， Chaudhary A ， Avila LS. 可 视 化 工 具 包（ VTK ） ： 为 现 代 图 形 卡 重 写 渲 染 代 码 。 SoftwareX 2015;1-2 ： 9- 1 - 2.http://dx.doi.org/10.1016/j.softx.2015.04.001，URLwww.sciencedirect.com/science/article/pii/S2352711015000035。[20]江 岸 电 脑 有 限 公 司 PyQT5 。 2022 年 ， https ： //www.riverbankcomputing.com/software/pyqt网站。[2022年2月4日]。[21]NOAA国家地球物理数据中心阿拉斯加塔蒂特勒克8/15角秒MHHW海岸数字高程模型。NOAA国家环境信息中心;2011年，https://www.ncei.noaa.gov/metadata/geoportal/rest/metadata/item/gov.noaa.ngdc.mgg.dem：700/html。[2022年1月14日[22]雷 米 ·N S-GeMS ： 斯 坦 福 地 质 统 计 建 模 软 件 ： 新 算 法 开 发 工 具 。 In ：Leuangthong O ， Deutsch CV ， editors. 地 质 统 计 学 班 夫 ， 2004 年 。Dordrecht：SpringerNetherlands;2005，p.865-71.http://dx.doi.org/10.1007/978-1-4020-3610-1_89网站。