1 引言

随着信息的高速发展, 合理的、高效率的数据挖掘对于信息的分析和研究至关重要.庞杂的地貌信息作为一种亟待人们探索和研究的空间数据信息对地球物理学的研究有着重要的意义^[1-3].近些年来, 数据挖掘广泛应用于石油、矿物等的预测和勘探领域, 对地貌信息的数据进行勘探和挖掘已经成为信息科学研究的难点和热点^[4-8].

在分析和研究地貌特征中, 地貌信息的挖掘和获取途径有多种, 其中GPS和探地雷达方法的应用较为常见.张勤等^[9]利用GPS和InSAR从3个时间段分析了西安地裂缝沉降的不均匀性, 探究了西安地面沉降与其时空演化的关系.郑文俊等^[10]、张清志等^[11]、殷海涛等^[12]利用GPS分别研究了青藏高原断裂带的变形情况和地表滑坡的特征, 为地貌的形成和抗震减灾提供了数据支持.何瑞珍等^[13]和王春辉等^[14]则利用探地雷达的方法对地表进行了测量, 准确实现了土质的检测和地表含水量的模拟.此外, 利用GIS软件对地貌信息进行模拟和仿真也较为常见, 文献[15-18]利用MAPGIS对地貌、地质信息进行了科学的统计和分析, 实现了对矿产资源、环境等信息的监测和预测, 为地貌的研究提供了可靠而详实的信息.近年来, 信息熵方法在地貌信息研究中的应用也日益广泛.王红等^[19]、李雄伟等^[20]提出了用信息熵作为地貌数据信息量大小的度量方法, 利用信息熵分析了地貌高程信息, 解释了熵与地貌高程信息分布与匹配的联系, 实现了对地形数据库中要素层的全面评价.朱春红等^[21]、解河海等^[22]、姜栋等^[23]利用高精度的数字高程模型数据从坡度, 坡向等专题研究和分析了地貌特征和地形信息因子, 为实现地貌特征的综合评价提供了有效的依据.

本文基于双重逆极限空间提出了一种新的关于地貌的不完备信息的数据挖掘方法, 将未知的地貌信息数据映射到双重逆极限空间中去, 通过研究两类空间之间的拓扑传递性完成对未知地貌信息数据集合的挖掘与分析.该方法把繁杂的特征地貌信息整体视为拓扑空间中的集合, 建立了适合于不完备地貌信息的拓扑空间模型, 以此为信息源, 可以提取所需的各种特征参数进行复合分析, 得到不同专题的综合评价, 在一定程度上提高了地貌信息整体特征的描述和未知区域预测的精度.

2 基于地貌信息的不完备信息空间的建立

根据不完备信息的系统四元组, S=〈U, A, V, f〉称为一个信息系统^[24].其中U={U_i|i=1, 2, …, n}表示对象的集合, A={A_k|k=1, 2, …, m}表示属性的集合, V_a表示a的值域, f表示U×A →V的一个信息函数, 它为每个对象在每个属性上赋予了一个不同的信息值.如果∃a∈A, 使得对于某个对象而言, V_a取值为空, 则该信息系统称之为不完备信息系统.地貌就是地壳表层内外共同作用到某一发展阶段所形成的三维地表形态, 是内外营力、介质性质和时间的函数^[25].其表达式为

(1)

其中, M表示地貌形态; F表示内外营力对地表的作用; x, y, z为空间坐标; m表示构成地貌的介质性质, 包括岩性与构造两个方面, t表示作用时间.运用子基描述特征原始信息子空间, 进而把地貌信息系统的原始信息空间描述为一个拓扑空间.设地貌信息空间为X, 将信息对象全体描述为

(2)

其中, x^u(u=1, 2, …, N)为X的对象集信息空间, 是地貌单元集合.

由于各类地貌单元不同的属性对应不同的值域, 所以一定存在两组子基{A_i}、{V_i}, 使X=.每一类子空间x^u中包含A个属性, 并且每一类属性A_i中包含V个值域, 则:

(3)

(4)

其中, A和V为整数, x₁^u, x₂^u, …, x_A^u为第u个地貌单元对象子空间x^u的属性子集. A为第u个地貌单元属性子空间的维数. x_{A, 1}^u, x_{A, 2}^u, …, x_{A, V}^u为第u个地貌单元对象子空间属性子集x_A^u的值域信息子空间. V为其值域子空间的维数.由于各个子空间的属性、值域的数量A、V并不相同, 为了保证子空间的维数相等, A、V取子空间最高维数的值, 其余低于A、V的子空间维数的特征量用零补齐.

由上述可知, X是一个集合, x=(x^u|u=1, 2, …, N)是X的一个对象子集族, 由拓扑空间定义可知, (X, f)是一个拓扑空间, 构成了不完备信息系统空间.

3 地貌信息的双重逆极限空间挖掘算法

对于不完备的地貌信息系统(X, f), 由于f不一定是同胚映射, 因而不具备可逆性, 导致系统(X, f)拓扑空间X上只是一个半动力系统, 不能满足拓扑传递性, 从而无法保证进行准确的数据挖掘.而逆极限系统中的诱导映射是同胚映射, 所以逆极限系统是一个拓扑动力系统, 从而可以利用逆极限系统的一些动力性质来揭示系统(X, f)相应的动力性质.将繁多的地貌信息描述成拓扑空间点集, 利用拓扑传递性实现这些数据信息相互表达, 仅仅依靠单一的映射, 根本无法将这些相对独立的数据信息完全且有效的涵盖, 所以引入紧致空间上两个交换的自映射并且通过这两个自映射任意次数的迭代, 实现对不完备地貌信息系统中的所有数据信息的表达.下面引入双重逆极限空间X作为目标信息提取空间.

3.1 双重逆极限空间的建立

设f: X→X, g: X→X是两个连续映射且满足fg=gf.考虑到不完备信息填充是以对象集合为基本单位进行替补, 记集合X={=(x_{i, j}^u)_{i, j=-∞}^{+ ∞}| f^mgⁿ(x_{i, j}^u)=x_{i-m, j-n}^u, ∀m, n≥0, i, j∈Z}为X的双重逆极限集合, 在X上引进度量

(5)

其中, x, y ∈X, d是X空间的度量, x_{i, j}^u=fⁱg^j(x^u)表示对象集合x^u经过两个连续映射f: X→ X、g: X→ X分别迭代i次和j次后得到的集合, 则空间(X, d)称为双重序列{X, f^~g}的双重逆极限空间.并且记为X=, 则X亦为紧致度量空间.

根据上述定义, (X, f^~g)为不完备信息空间, (X, f^~g)是与其对应的双重逆极限空间, 其中, f, g是逆极限空间的映射, 将双重逆极限的子空间描述成有限元素组合的表示形式, 即可分为

(6)

其中, (u=1, 2, …, N)为X的对象信息空间.为了保证空间(X, f^~g)与(X, f^~g)同胚, 应用局部映射原理建立原始空间和双重逆极限空间之间对象集与对象集之间的变换关系---拓扑反变映射, 利用反变映射挖掘原始子空间中的未知数据.建立拟线性映射: F_{i, j}使得F_{i, j}(x^u)=_{i, j}, 则F_{i, j}: X→X为一一对应的连续映射, 并且下列式子成立:

(7)

(8)

对于卡氏积空间X, X⊂x_{i. j}, 由引理命题得X中的基开集具有以下形式: U=F^-1_{im, jm} (Um), 其中n∈N, U_m⊂X_{im, jm}为开集.取i₀=max{i_m|m=1, 2, …, N}, j₀=max{j_m|m=1, 2, …, N}, 令U=, 则U为=X中的开子集.

因此, 引用式子(7)、(8)有如下推导:

(9)

所以F_{i, j}^-1(U)为X的一组基.由此可知∀∈X的开子集, 使得F^-1=x^j, x^j是X中的一个对象子集, 由上述证明可知F_{i, j}是同胚映射, 从而保证空间(X, f^~g)与空间(X, f^~g)同胚, 可以实现拓扑同构的两个空间在F_{i, j}作用下进行完备准确的相互转化和表达.

另外, 由前文不完备信息空间定义可知, (X, f^~g)是可分空间.可分性是同胚映射的拓扑不变性, 所以由同胚映射F_{i, j}得到的空间(X, f^~g)也为可分空间, 因此式子(6)的表示是正确的.

3.2 挖掘算法的提出

设f、g为紧致度量空间X的连续映射, 且fg=gf, 称f^~g为拓扑传递, 如果对于X中的任意开子集U和V, 都存在m, n∈Z, 使得f^mgⁿ(U)∩V≠ φ.根据式子(9), 此定义的式子可变为

(10)

其中F_{i. j}为X→X的同胚映射, f^mgⁿ(F_{i0, j0}^-1 (U))∩V=k, k为正整数表示集合的基数, V和U分别是X和X的对象集.

为简化计算, 取同胚映射F_{m, n}: X→ X为原双重空间系统的连续映射f^mgⁿ, m, n∈ Z, 令V=f^mgⁿ(F_{i, j}^-1(U)), 根据文中提出的双重逆极限空间的定义可得:

(11)

所以

(12)

其中, V是不完备地貌信息系统中存在未知数据的数据集合, x^u是X的对象集.

根据第一节所述, 将不完备地貌信息集合构建成关于对象集U、属性集A、值域集V三个坐标轴的立体空间模型, 使每个数据信息在此立体空间中都会被精确的描述出来.

地貌形态中的对象集地貌单元为

(13)

地貌单元属性集为

(14)

地貌单元属性的值域集为

(15)

根据双重逆极限空间的定义:

(16)

引出空间X对应的双重逆极限空间X作为目标特征信息空间, 对X空间中的对象集(x^u|u=1, 2, …, N)和X中的对象集(|u=1, 2, …, N)有:

(17)

根据文中前述证明, 空间X和其对应的双重逆极限空间X皆为紧致空间, 并且F_{m, n}为同胚映射, 因此两空间具有拓扑相似关系, 利用拓扑传递性分析空间X和X的拓扑相似度.下面将式(12)中的同胚映射函数明晰化, 为了计算方便, 取ξ=1, 则: V=F_{m, n}(xⁱ), 考虑到同胚映射F_{m, n}, 令:

(18)

上式为地貌单元对象集x^u+1到对象集x^u的有向信息传递指数, 与式(17)相符.其中H(x^u)是x^u的函数, H(x^u(u+1))是x^u与x^u+1的联合函数, 分别表示x^u的内在复杂性特征和x^u、x^u+1的联合复杂性特征, 即它们的变化规律和演变特征.将所有的地貌单元视为简单的信息源, 并分别对具有相等的高程值的地貌单元做频率统计, 从而得到所有的地貌单元的信息量, 其中,

(19)

(20)

这里将地貌中的全部地貌单元作为全体对象u, v是高程值的取值种类, 即地貌单元的属性值域集的取值种类.对具有不同高程值x_{A, V}^u的地貌单元分别进行统计, 求出概率, 即, 再根据式(19)即可计算出该地貌关于地貌单元的高程值的信息量大小, 所以它可以作为一个度量来描述地形的特性, 令:

(21)

其中, U为不完备地貌信息数据集合中有数据缺损的对象集合.基于概率相似思想, 当k取值最大时, U与V两个对象集合结构最为相似, 且V是由拓扑传递性映射得到, 因此可以认为U与V拓扑等价, 将式(19)、(20)代入(18)式中, 最后得出:

(22)

其中, , , 分别表示对象集x^u、x^u+1、x^u和x^u+1中各个元素信息在X集合中出现的概率.至此, 得到了不完备地貌信息系统关于地貌单元信息量的大小的结果.

另外, 地貌的地形状态特征由地貌高程差、粗糙度、相关长度、坡度方差等描述.作为地形参数, 这些地貌特征因子都可以反映出地貌信息大小的程度.随着地形粗糙度, 坡度方差的增加, 地形信息会更加丰富, 相关长度越小则地形数据的地理性越强, 地形信息也会越丰富.并且, 公式(19)的计算结果与地貌单元的信息量呈负相关, 即也与高程差等地貌特征呈负相关, 因此, 通过公式(22)计算得出的地貌单元的坡度、坡向之间的比较, 可以相对得知与已知高程的地貌单元的高程差, 从而最终确定未知的地貌单元的大致高程, 最后, 利用公式(22)的计算结果与上述地形参数的信息含量的丰富度进行匹配和对比, 通过GIS软件进行模拟后便可以得出大致的地貌形态.

4 应用实例

首先, 利用GIS软件基于规则网络模型, 用含有未知数据的地貌信息系统模拟一个地貌表面, 将这个地貌表面按照对象集、属性集、属性值域集将其区域空间划分为规则的网格单元, 每个网格单元都代表地貌信息中的一个对象集, 每个地貌单元中都有其特有的内外营力等属性, 每个地貌单元都对应一个高程值作为其信息的值域集合.在随机模拟的过程中, 有意将那些未知的缺损数据信息平均分散开来, 使模拟出的地貌形态不至于由于局部的地形信息的匮乏致使地表过分离散.同时, 地貌单元信息的缺损会使地貌单元在数字高程模型中无法表达, 只能由相邻的地貌单元表示出其周围的高程和起伏, 因此不完备的地貌信息系统中的地貌单元在数字高程模型中的数量相对稀疏, 因此地貌整体表现为起伏度不大, 地表较为光滑, 坡度、高程差、坡向等得不到良好的表达.见图 1.

由于在不完备地貌信息系统数据挖掘当中, 地貌随时间的演化是一个相当漫长的过程, 因此, 在这里地貌信息系统中的时间t、内外作用力F等属性集合不做过多考虑, 将地貌演化时间定格在一瞬, 因此相对来说地貌现有的高程、坡度、坡向等属性参考价值较大.考虑到这些, 运用GIS软件, 针对上述不完备地貌信息基于其Grd模型对所有地貌单元和点集的高程、坡度、坡向进行提取, 提炼出地貌信息中地貌单元的有效高程值、坡向和坡度, 并且将未知的高程信息的地貌单元信息按照3: 1的比例均匀的分散到已知的地貌单元信息当中并且制成信息表格. 表 1是部分地貌高程信息列表.

表 2、表 3分别是同一块地表的坡度和坡向的数据信息:

由以上两个列表可知在一块极微小的地表及其周围的地貌单元, 其地貌的坡度、坡向的变化几乎微乎其微, 根据前文所述步骤分别将有未知数据的地貌单元集合通过逆极限映射映射到与其拓扑相似的双重逆极限空间中.如:横向第24个地貌单元对应的高程值对象集合x²⁴(*, 59. 977, 59. 966, *, 59. 943)对应的双重逆极限空间的拓扑相似集合V²⁴为

(23)

其中, p²⁴、p²⁵分别表示地貌单元集合x²⁴、x²⁵中每一个元素在整个信息系统中出现的概率.经计算可得: V²⁴=(59. 978, *, 59. 966, 59. 954, *), 由于对象集x²⁴中属性Y₄和Y₇对应的值域分别为59. 977和59. 943, 因此将V²⁴中属性Y₄、Y₇的值域*替换为59. 977和59. 943, 可以得到V²⁴的数据集合为V²⁴=(59. 978, 59. 977, 59. 966, 59. 954, 59. 943), 运用相同的方法对表 1数据进行处理, 见表 4所示.

根据以上计算方法分别对所有的具有未知数据的地貌单元进行数据处理, 最终得到那些未知地貌单元的高程值、坡度、坡向.运用GIS软件的Grd模块对原具有数据缺损的地貌单元逐一进行高程值的编辑, 最终得到图 2.

由图 2可知将拥有空缺数据的地貌信息系统输入GIS数据库模拟得到的是一个模糊的轮廓不清晰的地貌形态, 利用本文提出的方法, 将计算得到的地貌单元数据集合重新输入GIS数据库中, 模拟后得到的地貌形态清晰, 轮廓明显, 地貌起伏度得到了良好的表达, 能更切合实际的接近或"还原"原有的地貌形态.

5 结论

本文从拓扑学的角度将不完备的地貌信息系统数据集合描述为一个拓扑空间, 根据逆极限定义, 利用映射诱导法建立与原地貌信息空间拓扑同构的双重逆极限空间作为目标信息空间, 然后利用拓扑传递性分析原地貌信息空间与目标信息空间的拓扑相似度, 实现不完备的地貌信息空间数据挖掘, 最后利用文中提出的理论方法对地貌信息系统进行数据挖掘, 为了方便计算和验证, 用GIS软件进行前后模拟对比, 实验结果表明上述理论方法是可行且准确的.相较于传统的地貌数据挖掘方法, 本文提出的理论方法是基于拓扑学角度实现的地貌信息数据挖掘, 为了使地貌信息系统中繁杂没有联系的数据更加直观、便于分析, 本文将地貌单元、地貌特征因子等参量划分为对象集、属性集、属性值域集, 从而避开了与地貌信息挖掘无关的数据计算, 使算法得到了简化.