应用Python在组件式GIS中集成地质统计学
——以Python和SuperMap Object进行趋势面分析为例
马维峰 李林 王晓蕊
摘 要：随着GIS应用在地学应用领域的拓展和加深，有必要为GIS工具引入和整合地质统计学等分析工具。应用COM组件技术和Python及其科学计算扩展模块SciPy可以迅速高效的开发出具有COM接口的地质统计学分析模块，并整合到组件式GIS工具中。
关键词：地质统计学，趋势面，SuperMap Object，组件式GIS，Python，SciPy
1. 引言
地质统计学是数学地质领域中一门发展迅速且应用广泛的新兴学科，广泛应用在异常评价、找矿勘探、矿体圈定、储量计算等地学生产和科研等领域（赵鹏大，1994，侯景儒1996，1998）。趋势分析或趋势面分析，作为地质统计学中的一种，就是要把研究的地质现象进行分解，表示为数学方程为：
Z = Z0 + A + R （1）
Z代表地质现象，Z0代表区域性因素（背景），A代表局部性因素（异常），R表示随机因素（干扰值）。地学研究中，问题往往都归结为寻找区域上的异常值的问题。趋势面分析的方法为将地学信息的某个指标或某些指标的综合表示为一个空间曲面，然后人工根据数学方法构造一个曲面代表区域性因素（背景值），同时应用统计方法消除随机因素，以达到突出和圈定异常值的目的（於崇文，1980，赵旭东，1992）。
传统的地质统计学分析或趋势面分析是应用地学专业软件包或其他具有回归分析能力的数学、统计学软件包来进行的，然后结合地图等工具进行空间和区域上的解释。目前，随着GIS在地学实际工作和研究工作中的引入，在很多结合地质统计学的地学研究中也引入了GIS方法和工具，但现有GIS工具并不支持地质统计学分析，因此空间信息和指标在必须在GIS工具和地质统计学工具之间相互导入导出，影响和限制了工作效率的提高，造成了不必要的数据信息精度的损失，有时甚至由于数据文件格式的封闭性，造成工作不能开展或中断。
随着GIS应用在地学应用领域的拓展和加深，为GIS工具引入和整合地质统计学等分析工具势在必行，因此有必要开展地质统计学和GIS整合技术和方法的研究。本文将以趋势面分析为例，探讨应用COM组件技术和Python及其科学计算扩展模块SciPy开发地质统计学分析工具的方法，以及如何在组件式GIS工具中整合地质统计学扩展模块。
2. 地质统计学与GIS工具的集成方式目前在GIS基础软件与空间分析模型和地学模型常用的集成方式主要有以下几种（宋关福，1998，方裕，2001，冯克忠，2003）： （1）通过数据文件在GIS平台软件和空间分析模型或地学模型之间建立联系，属于松散集成方式。其优势是GIS平台的开发和空间分析模型的开发都可以独立进行，只要遵循一致的数据交换格式就可以进行；其缺点是通过文件方式交换数据，不适合大量且频繁的数据交换的情况；（2）直接使用GIS软件提供的二次开发语言编制应用分析模型，例如Mapinfo，Arcinfo都有其二次开发语言。这种方式在开发一些小系统时比较合适，但由于GIS所提供的二次开发语言大都相对简单，难以开发相对复杂的分析模型；（3）基于组件技术的扩展方式。组件技术的促进了组件式地理信息系统的产生，大大简化了GIS整合空间分析模型或地学模型的难度。基于组件的接口调用方式，一方面保证了GIS系统开发和分析模型开发之间的独立性，另一方面又实现了二者的无缝集成。因此可以基于组件技术，在组件式地理信息系统中，如ArcGIS 8.0及以上版本、SuperMap Object等的开发或扩展中，为GIS工具引入和整合地质统计学分析工具。
作为动态脚本语言的Python具有语法简单，易于学习，解释性、无需编译，易于部署和维护等优点，与直接使用静态编程语言C、C++、VB等语言相比，使用Python开发效率更高，开发和维护难度较低，作为系统扩展和快速开发语言具有一般静态编程语言不具有的优势（Hammond，2000，Lutz，2001，http://www.python.org）；另一方面，Python具有高层的内建数据结构，优秀的科学计算扩展库SciPy（包括科学计算、统计分析、可视化等模块，http://www.scipy.org），因此在ArcGIS、SuperMap Object等基于COM的组件式地理信息系统平台中，可以通过Python来实现和整合地质统计学分析工具。
3. 系统设计与实现 3.1. 系统结构地质统计学分析工具的设计需要达到以下目的：（1）地质统计学分析工具将以一个COM组件库的形式暴露给调用者，组件可以通过COM接口被支持COM扩展或集成的GIS工具调用；（2）组件的内部实现、结构对调用者隐藏，调用者只通过COM接口与组件进行交互，这样，组件的部署，维护对组件调用者完全无关，组件本身可以在任何时候修改、升级而不影响调用者的使用；（3）组件的数据结构、数据文件不与任何特定GIS平台绑定，以保证组件可以在多个GIS工具中使用或集成。
我们将使用Python及其科学计算扩展模块，应用面向对象的方法来开发地质统计学模块，并使用PythonCOM（Hammond ，2000，http://www.activestate.com）扩展模块来编写可以在COM中调用的接口[①]（图 1）。因此，系统主要由两大部分组成：（1）Python编写的地质统计学分析计算组件库，本部分使用Python实现，可以在任何可以运行Python的环境中运行；（2）应用PythonCOM实现的地质统计学分析计算组件库的COM接口，通过此接口，任何可以调用COM组件的编程环境和语言都可以调用此组件库，PythonCOM必须在Windows环境下运行。任何可以通过COM集成或扩展的GIS工具都可以通过PythonCOM实现的接口，使用此地质统计学分析计算组件库。
GIS工具
组件调用者
（ArcGIS）
（SuperMap obeject）
…
Python编写的地质统计学组件库
COM
调用
接口
趋势面分析
回归分析
聚类分析
……
地质统计学组件库
Python执行环境
标准库
SciPy扩展库
PythonCOM扩展
Windows API + 其他运行库 + …
图 1 地质统计学组件库系统结构图
3.2. 系统实现我们首先使用Python，借助其强大的科学计算扩展模块SciPy编写地质统计学的有关模块。SciPy内置了大量优秀的科学计算工具包，例如矩阵运算、线性代数[②]、统计分析、FFT等等，可以在程序编写时直接使用；SciPy中内置的各类科学计算算法，大多都是使用C语言或者Fortran语言应用成熟的算法或代码库实现，而后集成在在Python中的。因此，使用SciPy实现地质统计学的有关算法，一方面可以借助Python良好的编程性，迅速高效的实现；另一方面，由于SciPy本身是使用C语言或者Fortran语言实现的基础算法，速度和使用C语言或者Fortran语言编写同类算法差别不大，所以由SciPy实现的地质统计学运算速度将不会因为Python是脚本语言而打折扣。
下面我们通过使用SciPy实现趋势面分析为例来介绍整个系统的实现。趋势面一般使用2种数学函数表示，一种是多项式函数，另一种是傅立叶级数，其中多项式趋势面最为常用。趋势面求解过程包括以下几个步骤：（1）根据需要求解的趋势面的数学表达式（多项式及其次数或傅立叶级数），使用最小二乘法建立求解方程组，这一部可以参考有关参考文献（於崇文，1980，赵旭东，1992），直接使用已经演算好的方程组；（2）根据已知的观测数据（Xi，Yi，Zi），计算方程组的求解矩阵；（3）应用高斯消元法求解方程组；（4）通过求解结果，计算趋势面的拟合度和残差（即公式1中的A）。这4步中，第三步是编程实现较难实现的一步，但SciPy提供了矩阵等抽象数据结构，并实现线性方程组求解等算法，可以直接使用，使用示例见图2[③]，对于使用矩阵运算还是线性方程求解函数，结果和效率是一致的，对于使用者，如同使用Matlab等工具，可以不考虑底层算法，因为SciPy已经优化过了，不需要使用者操心。对于实际应用，首先计算系数矩阵A和b、然后直接就可以求解待定系数矩阵，得到趋势面函数的表达式。
# 导入SciPy模块
>>> from scipy import *
# 建立矩阵A和b
>>> A = mat('[1 3 5; 2 5 1; 2 3 8]')
>>> b = mat('[10;8;3]')
# 通过矩阵运算求解方程组
>>> A.I*b
Matrix([[-9.28],
[ 5.16],
[ 0.76]])
# 通过调用线性方程求解函数求解方程
>>> linalg.solve(A,b)
array([[-9.28],
[ 5.16],
[ 0.76]])
图 2 在Python命令行环境下求解线性方程组得两种方法
完成Python代码编写后，下一步就是将此代码封装成一个COM对象，可以定义一个新的类，也可以直接将这个类封装为COM对象，实现需要的属性和方法，使用PythonCOM扩展，实现ProgID, CLSID，public_methods等属性，在注册表里注册这个Python类，使其成为一个可以被调用的COM对象①。
这样，对于地学分析，就可以调用此对象，在GIS等可视化工具中，将原始数据、趋势面、残差表示为等高线或曲面，进行地学分析。
3.3. 组件使用下面将以在SuperMap Object中使用和调用Python编写的趋势面分析组件为例来说明如何将此组件集成到GIS工具中[④]。
（1）数据组织：通过空间分析或查询，将需要分析的数据集应用“Query”方法生成一个“Recordset”中，通过“soRecordset.GetFieldValue”方法获取需要分析的Z值，通过“soRecordset.GetGeometry()”方法获得需要分析的几何对象，一般是点（soGeoPoint），通过其获得X，Y坐标，写入临时文件或数组，循环得到所有的X，Y，Z。
（2）趋势面分析方法的调用：我们编写的应用Python求解趋势面分析组件使用文本文件方式传递数据，其优势是可以传递大数据量的数据，GIS平台无关等，劣势是效率不如内存方式高，但对于地学问题，由于数据量一般都较大，因此使用文本文件方式较为合适。调用趋势面分析需要通过以下方法创建对象并调用（图3），strDataFile是输入的数据文件路径和文件名，strResultPath是输出的数据文件的路径，组件将结果、拟合度、残差数据文件分别写入此目录。
Dim objTrend As Object
Set objTrend = CreateObject("GeoStatistics.Trend")
'调用2次趋势面分析
objTrend.Poly2(strDataFile, strResultPath)
图 3 在VB下调用趋势面分析
（3）在SuperMap Object中生成需要的数据集：根据Python的输出，在SuperMap Object中可以新建数据集，将这些数据读入。新建一个数据集“objDataset”，应用“objDataset.Query("", true)”，生成一个“Recordset”，然后从数据文件中读出（原始数据，残差）或计算出（趋势面）X，Y，Z，根据X，Y生成“soGeoPoint”对象，利用“Recordset”对象的“AddNew”方法添加此对象到“objDataset”，并将原始数据、趋势面、残差的值（Z）赋给相应得“Field”，即可动态生成需要的数据集。
（4）在SuperMap Object中进行可视化分析：利用前面生成的数据集，我们可以使用SuperMap Desktop等工具，由这些数据集生成Grid数据，然后进行3维可视化显示和分析，也可以直接应用编程方法，通过“soGridAnalyst.Point2DToGrid”等类似方法生成需要的数据集。当生成了Grid数据集后，即可通过可视化方法，结合GIS的空间分析方法，如叠加分析、缓冲区分析，进行地学机理的分析。
4. 结论
应用Python及其扩展模块SciPy可以高效的实现地学统计分析，使用PythonCOM扩展模块可以将其封装为COM组件，在组件式GIS中调用，为GIS工具引入地质统计学分析功能。相对于其他开发语言，Python优秀的科学计算、统计分析、可视化扩展模块SciPy，以及Python动态语言的优势，可以提高地质统计学组件的开发效率和开发质量，缩短开发周期。基于COM的组件式开发和调用一方面保证了GIS系统的开发和Python扩展的开发的独立性，对于Python对象开发的组件，可以在任意时刻随意修改（包括编译打包后），不会影响GIS组件对其的调用；另一方面又实现了Python扩展组件和GIS工具的无缝集成。
参考文献
1. 赵鹏大，李紫金，胡旺亮，矿床统计预测（M）．北京：地质出版，1994，10-195
2. 侯景儒，尹镇南，李维明等，实用地质统计学（M）．北京：地质出版社，1998
3. 侯景儒，中国地质统计学(空间信息统计学)发展的回顾及前景（J），地质与勘探，1996，32（1）：20-25
4. 赵旭东，石油数学地质概论，北京：石油工业出版社，1992，65-81
5. 於崇文，数学地质的方法与应用，北京：冶金工业出版社1980，
6. 冯克忠，万庆，励惠国，基于组件技术的GIS广义空间分析，地球信息科学，2003，1：62-66
7. 宋关福，钟耳顺；组件式地理信息系统研究与开发，中国图象图形学报，1998，3（4）：314-317
8. 方裕，周成虎等，第四代GIS软件研究，中国图形图像学报，2001，6（9）：817-823
9. Hammond M., Robinson A., Python Programming on Win32, O'Reilly, 2000
10. Lutz M., Programming Python, 2nd Edition, O'Reilly, 2001
11. Rossum G., Extending and Embedding the Python, http://www.python.org, 2003
12. http://www.python.org
13. http://www.activestate.com
14. http://www.scipy.org
[①]PythonCOM是Python在Windows环境下将Python模块封装成COM对象的一个Python扩展，详细介绍见马维峰等，在组件式GIS开发中集成Python方法研究，未发表
[②] SciPy的矩阵运算和线性代数模块构建于ATLAS LAPACK和BLAS，是非常优秀的线性代数代码库
[③] 由于本问题实际求解代码太长，不便列入文中，但思路和求解方法完全一致
[④] 本部分所使用的一些名词、对象及其方法都是SuperMap Object的对象或方法，不熟悉者可以参考其有关文档
本文地址：http://com.8s8s.com/it/it42646.htm