** 应用 Python 在组件式 GIS 中集成地质统计学 **

** ——以 Python 和 SuperMap Object 进行趋势面分析为例 **

马维峰李林王晓蕊

** 摘 ** ** 要 ** ：随着 GIS 应用在地学应用领域的拓展和加深，有必要为 GIS 工具引入和整合地质统计学等分析工具。应用 COM 组件技术和 Python 及其科学计算扩展模块 SciPy 可以迅速高效的开发出具有 COM 接口的地质统计学分析模块，并整合到组件式 GIS 工具中。

** 关键词 ** ：地质统计学，趋势面， SuperMap Object ，组件式 GIS ， Python ， SciPy

1. 引言

地质统计学是数学地质领域中一门发展迅速且应用广泛的新兴学科，广泛应用在异常评价、找矿勘探、矿体圈定、储量计算等地学生产和科研等领域（赵鹏大， 1994 ，侯景儒 1996 ， 1998 ）。趋势分析或趋势面分析，作为地质统计学中的一种，就是要把研究的地质现象进行分解，表示为数学方程为：

_ Z = Z 0 + A + R _ _ （ 1 _ _ ） _

_ Z _ 代表地质现象， _ Z 0 _ 代表区域性因素（背景）， _ 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 ）根据已知的观测数据（ _ X i _ _ ， Y i _ _ ， Z i _ ），计算方程组的求解矩阵；（ 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],