-
摘要: 利用不同尺度下的位场数据进行边界检测可获得相应的多尺度边界信息，其在区域地质构造单元划分、地质填图、矿产资源圈定和结晶基底研究等方面具有广泛应用。通过联合二维经验模态分解（bi-dimensional empirical mode decomposition, BEMD）和小波模极大值（wavelet modulus maximum, WMM）方法实现位场数据多尺度分解及各尺度异常边界检测，首先，证明了BEMD方法能够实现对位场的多尺度分解，且分解各尺度异常具有实际地质含义；然后，对构建的模型重力异常进行边界检测，结果显示, WMM方法具有良好的抗噪声干扰能力，能够清晰、准确地检测出组合模型体的边界位置、走向及分布范围；最后，将该联合方法用于三峡地区布格重力异常的多尺度分解和各级固有模态函数的边界检测，结合径向对数功率谱曲线对分解各尺度异常进行场源深度估计，得到了三峡地区不同深度层位上的构造或岩体线性分布特征，并对各尺度异常边界检测结果与区域地质构造特征进行了比较分析和解释，验证了该方法的优越性和有效性。Abstract:Objectives Multi-scale boundary information can be obtained using potential field data at different scales, which having an important practical value on account of the tectonic division, geological mapping orebody delineation and the research of crystalline basement.Methods We use the bi-dimensional empirical mode decomposition (BEMD) and the wavelet modulus maximum (WMM) methods to obtain the multi-scale boundary detection of potential field data.Firstly, it is shown that the BEMD method can realize the multi-scale decomposition of potential field data, and the anomalies at different scales possess practical geological significance. Secondly, the model test established by this paper shows that WMM method has the ability of resisting the noise interference so as to detect the boundary of the combined model about the distribution, location and direction clearly and accurately using the gravity anomalies. Finally, this paper also achieves the results of multi-scale decomposition of the Bouguer gravity anomalies in the Three Gorges region by BEMD algorithm, and estimates the approximate depth of the buried source by the radial logarithmic power spectrum of each BIMF (bi-dimensional intrinsic mode function)and residue, and we deal with each BIMF using WMM method so as to obtain the distribution characteristic of the structure and rock mass in different depths.Results and Conclusions The results of boundary detection at all levels are analyzed and explained according to the regional geology structure feature, the advantage and effectiveness of this method is verified in the end.
-
全球离散格网系统(discrete global grid system，DGGS)是一类以地球作为组织框架的空间参考系统。相较于传统空间数据模型，其更适合解决大尺度的问题[1]，且单元运算可完全借助编码实现，有助于提高数据处理效率。六边形具有拓扑关系一致、采样效率高等理想几何特性，有助于实现空间分析算法[2-3]及数据的可视化表达[4]，但其不具备叠合一致性，无法直接建立多分辨率格网层次关系。Gibson等[5]提出一种多维数据结构，其二维特例可通过广义平衡三进制(generalized balanced ternary，GBT)运算实现不同层次六边形单元索引。然而，GBT采用高分辨率单元聚合生成低分辨率单元，在球面上会导致格网破碎或重叠。
基于多面体剖分建立六边形全球离散格网系统仍被认为是目前最有效的方法。正二十面体面数最多，投影到球面后变形最小[6]，因而被广泛使用。Sahr[7]认为，正二十面体等积三孔六边形格网系统(icosahedral snyder equal area aperture 3 hexagonal grid, ISEA3H)具有较好的几何属性，他借助与GBT类似的码元分布，建立改进广义平衡三进制(modified generalized balanced ternary，MGBT)编码方案，但其编码不唯一，仅适用于表达矢量数据。之后，Sahr[8]又将MGBT分解为两种子单元模块的组合，建立三孔六边形树(aperture 3 hexagon tree，A3HT)唯一编码方案，并扩展至正二十面体。Vince[9]建立了正二十面体三孔六边形格网的进制索引，并实现了傅里叶变换。童晓冲等[10]建立了四孔六边形格网系统的四元平衡结构(hexagonal quad balanced structure，HQBS)，并实现正二十面体格网编码运算。
数学家已从理论上证明，封闭球面拓扑等价于封闭正二十面体表面，全球离散格网系统编码运算拓扑等价于正二十面体格网系统编码运算。然而无限平面格网的编码运算规则不适用于正二十面体封闭表面，如何实现正二十面体上的高效率编码运算是目前的研究难点。A3HT通过指定正二十面体上首层单元的递归剖分规则建立正二十面体格网系统，该方案不是直接利用编码运算实现单元操作，而是先借助二维整数坐标系计算结果，再将之转换为编码[8]。HQBS采用分面编码到全局编码的映射解决上述问题，但只能保证编码运算在单个三角面上连续，跨面操作仍需借助笛卡尔坐标运算实现。上述两种方案的计算效率均不高，理论上也不严密。Vince[9]通过对三孔六边形格网层次结构的裁剪拼接建立正二十面体格网模型，由于区分奇(偶)层编码，运算规律较复杂，且相邻层次单元方向会交替变化，不利于格网层次分析。
本文结合四孔六边形格网系统在正二十面体表面的分布特点，基于六边形格点四叉树定义顶点瓦片与面瓦片结构，提出正二十面体四孔六边形格网系统编码方案，归纳编码运算规律，通过对比实验验证方案的有效性与优越性。
1 平面四孔六边形格网编码运算
四孔六边形格网单元剖分过程如图 1所示，为保证编码的唯一性，本文的格网剖分选择图 1阴影部分所示的子单元模块，它们关于父单元对称，几何属性更优。
1.1 编码方案
为便于讨论，定义平面格网系统单元中心为格点(lattice point)等效替代单元。首层格网单元可自定义。令Ln表示第n层格点集合(n≥1)，根据编码运算需要，本文定义的L1在复数平面上的位置如图 2所示，图中$ R $轴和$ I $轴分别表示复数坐标系的实轴和虚轴。根据图 1生成的L2与L1层次关系见图 3，为方便表示，图 3只选取L1中心7个格点为例，每个格点在L2各生成4个子格点，中心子格点与其父单元格点重合。按照上述思路，可依次生成Ln($ n>2$)。
根据文献[11]，在复数平面上，令$ \theta =\frac{1}{2}+\frac{\sqrt{3}}{2}\mathrm{i} $，L1用以下5个集合$ {D}_{j}\prime (1\ll j\ll 5)$表示：$ {D}_{1}\prime =\left\{0, {\theta }^{e}\left|\mathrm{ }1\le e\le 6\right.\right\} $，$ {D}_{2}\prime =\left\{2{\theta }^{e}\left|\mathrm{ }1\le e\le 6\right.\right\}$, $ {{D}_{3}}^{\prime }=\left\{{\theta }^{e}+{\theta }^{(e+1)\mathrm{\%}6}\left|1\le e\le 6\right.\right\} $，$ {D}_{4}\prime =\left\{2\alpha \left|\alpha \in {{D}_{3}}^{\prime }\right.\right\} $，$ {D}_{5}\prime =\left\{2{\theta }^{{e}_{1}}+{\theta }^{{e}_{2}}\left|\mathrm{ }\left|{e}_{1}-{e}_{2}\right|\right.=\mathrm{1, 1}\le {e}_{1}, {e}_{2}\le 6\mathrm{或}{e}_{1}=1, {e}_{2}=6\mathrm{或}{e}_{2}=1, {e}_{1}=6\right\}$。其中，%表示取余；$ e\mathrm{、}{e}_{1}\mathrm{、}{e}_{2}\in Z\mathrm{ }$。令$ \omega =-\frac{1}{2}+\frac{\sqrt{3}}{2}\mathrm{i}\mathrm{ } $，$ D=\left\{0, \omega , {\omega }^{1}, {\omega }^{2}, {\omega }^{3}\right\} $，经过严格证明，四孔六边形格网系统的格点集合可唯一表示为:
$$ \left\{\begin{array}{l}{L}_{n}={L}_{1}+\stackrel{n}{\sum \limits_{j=2}}\frac{1}{{2}^{j}}D\\ {L}_{1}=\frac{1}{2}D\prime \end{array}\right. $$ (1) 式中，$ D\prime ={D}_{1}\prime \bigcup {D}_{2}\prime \bigcup {D}_{3}\prime \bigcup {D}_{4}\prime \bigcup {D}_{5}\prime $。在$ {L}_{1}$中任选一个元素，对于后n-1层(即$ 2\le j\le n$)，每层在集合D中任选一个元素参与计算，$ n $层相加结果为$ {L}_{n} $中一个元素。按上述方法，所有元素组合的计算结果集合为$ n$层格点集合$ {L}_{n}$。如图 3所示，$ {L}_{1}$中心7个格点集合可表示为$ \left\{0, \frac{1}{2}{\theta }^{e}\left|\mathrm{ }1\le e\le 6\right.\right\} $, 假设从中任意取一格点，坐标用$ \alpha $表示，则其在第二层的4个子格点集合可表示为:$ \alpha +\frac{1}{{2}^{2}}D=\alpha +\frac{1}{4}\{0, \omega , {\omega }^{2}, {\omega }^{3}\} $。
式(1)表明，四孔六边形格点系统本质上是一种定位计数系统，格点与实数域上的二进制、十进制等数字的本质完全相同，是复数域上一种特定形式的“数”。该定位计数系统的进制为2，数字位集合为$ D\prime $和D。对于任意一格点复数坐标$ x\in {L}_{n} $，可用数字位和进制的加乘表示，即：
$$ x=\frac{1}{2}{d}_{1}+\frac{1}{4}{d}_{2}+\cdots +\frac{1}{{2}^{n}}{d}_{n} $$ (2) 式中，d1∈$ D\prime $；d2, d3…dn∈D。这与十进制数类似，因此，省略幂指数的底，式(2)又可记为：
$$ x={d}_{1}{d}_{2}\cdots {d}_{n} $$ (3) 式中，复数数字位di($ 1\le i\le n$)可用相应的整数替换，由此得到该格网系统的整数编码集合。对数字位d1做如下替换，当$ {d}_{1}\in {D}_{1}\prime $时，将d1替换为坐标的幂次$ e$，即$ {d}_{1}\to e\left(1\le e\le 6\right) $; 若$ {d}_{1}=0 $，则${d}_{1}\to 0 $；当$ {d}_{1}\in {D}_{2}\prime $时，$ {d}_{1}\to 100e\left(1\le e\le 6\right) $，即100与e相乘；当$ {d}_{1}\in {D}_{3}\prime $时，$ {d}_{1}\to 10\left((e+1)\mathrm{\%}6\right)\left(1\le e\le 6\right)$；当$ {d}_{1}\in {D}_{4}\prime $时，$ {d}_{1}\to 100{d}_{1}\prime $，其中$ {d}_{1}\prime $为$ \alpha $替换后码元；当$ {d}_{1}\in {D}_{5}\prime $，$ {d}_{1}\to 100{e}_{1}+{e}_{2}$。对于式(3)中$ {d}_{j}(2\le j\le n) $，若$ {d}_{j}={\omega }^{e} $，则${d}_{j}\to e(1\le e\le 3)$；若$ {d}_{j}=0 $，则$ {d}_{j}\to 0$。
根据上述替换规则，数字位集合D'被替换为整数码元集合$ E\prime $$ =\left\{\mathrm{0, 1}, \mathrm{2, 3}, \mathrm{4, 5}, \mathrm{6, 10, 20, 30}, \right.$40 $ {\rm{50}},{\rm{60}},{\rm{100}},{\rm{200}},{\rm{300}},{\rm{400}},{\rm{500}},{\rm{600}},{\rm{106}},{\rm{601}},{\rm{201}},{\rm{102}},{\rm{302}},{\rm{203}},{\rm{403}},{\rm{304}},{\rm{504}},{\rm{405}},{\rm{604}},{\rm{506}},{\rm{1}}000,2{\rm{000}},{\rm{3000}},{\rm{4000}},{\rm{5000}},{\rm{6}}\left. {000} \right\}, $换为码元集合$ E=\left\{\mathrm{0, 1}, \mathrm{2, 3}\right\}$，由此得到任一格点x的唯一编码标识:
$$ x={e}_{1}{e}_{2}\cdots {e}_{n} $$ (4) 式中，码元$ {e}_{1}\in E\prime $，$ {e}_{j}\in E(2\le j\le n) $。为方便区分，实际表示时在$ {e}_{1} $后加逗号分隔，例如编码$ 20\mathrm{、}213 $为第四层编码。该编码方案是一种典型的四叉树结构，首层每个格点都分别生长为一棵独立的四叉树，因此将这一结构命名为六边形格点四叉树(hexagonal lattice quad tree，HLQT)。
图 4为首层格点$ {D}_{1}\prime $在L3中的子格点编码示意图，由于$ {D}_{1}\prime $对应首层码元均为1位数字，图中不用逗号分隔。
![]()
图 4 格网第三层编码及运算示意图Figure 4. Diagram of Codes in the Third Level and Example of Code Arithmetic1.2 编码运算
编码运算采用更适合计算机处理的一维码元序列代替传统二维实数坐标，从而提高计算效率。
设Hn为HLQT第n层编码集合(n≥1)，编码$ \alpha \mathrm{、}\beta \in {H}_{n}$，符号“$ \oplus $”表示编码加法运算。HLQT编码加法运算满足向量相加的平行四边形法则。如图 4所示，$ \alpha \oplus \beta =123\oplus 010=233$。经过归纳、验证，无限二维平面上的HLQT编码运算满足交换群的性质[11]。
设$ n $层编码$ \alpha ={a}_{1}{a}_{2}\cdots {a}_{n} $，$ \beta ={b}_{1}{b}_{2}\cdots {b}_{n} $，编码加法计算方式与十进制加法类似，从右侧末位码元(即$ {a}_{n}$或$ {b}_{n} $)开始，向左侧首位码元(即$ {a}_{1}$或$ {b}_{1} $)逐位计算，但进位机制略有差异：HLQT加法运算当前位计算结果可能在其左侧多个加法位上产生进位码元。编码加法运算可分解为不同加法位对应编码相加之和：
$$ \begin{array}{c}\alpha \oplus \beta ={a}_{1}{a}_{2}\cdots {a}_{n}\oplus {b}_{1}{b}_{2}\cdots {b}_{n}=({\varepsilon }_{n}^{{a}_{n}}\oplus {\varepsilon }_{n}^{{b}_{n}})\oplus \\ ({\varepsilon }_{n-1}^{{a}_{n-1}}\oplus {\varepsilon }_{n-1}^{{b}_{n-1}})\oplus \cdots \oplus ({\varepsilon }_{2}^{{a}_{2}}\oplus {\varepsilon }_{2}^{{b}_{2}})\oplus ({\varepsilon }_{1}^{{a}_{1}}\oplus {\varepsilon }_{1}^{{b}_{1}})\end{array} $$ (5) 式中，$ {\varepsilon }_{j}^{e}=\underset{j-1\mathrm{ }\mathrm{个}}{\underset{⏟}{0\cdots 0}}e\underset{n-j\mathrm{ }\mathrm{个}}{\underset{⏟}{0\cdots 0}} $，显然，$ {\varepsilon }_{j}^{e}\in {H}_{n}(1\le j\le n) $，j-1与n-j表示码元0的个数。
对于非首位加法位上的码元相加，即$ {\varepsilon }_{j}^{{e}_{1}}\oplus {\varepsilon }_{j}^{{e}_{2}}(1 <j\le n) $，根据上文编码方案可得：$ {e}_{1}, {e}_{2}\in E $，具体计算结果以及进位规则如下：
1) 若$ {e}_{1}={e}_{2}$，$ {\varepsilon }_{j}^{{e}_{1}}\oplus {\varepsilon }_{j}^{{e}_{2}}=\left\{\begin{array}{l}\underset{j-2\mathrm{ }\mathrm{个}}{\underset{⏟}{0\cdots 0}}{e}_{1}0，j>2\\ \left(2\right.{e}_{1})0，j=2\end{array}\right.$；
2) 若$ {e}_{1}\ne {e}_{2}$，且$ {e}_{1}, {e}_{2} $均不为0，$ {\varepsilon }_{j}^{{e}_{1}}\oplus {\varepsilon }_{j}^{{e}_{2}}=\left(\frac{{e}_{1}+{e}_{2}}{{4}^{\left|{e}_{1}-{e}_{2}\right|-1}}\right)\underset{j-1\mathrm{ }\mathrm{个}}{\underset{⏟}{e\cdots e}} $，其中$ e$为$ \{0, {e}_{1}, {e}_{2}\} $对$ E $的补集中的元素；
3) 若$ {e}_{1}\ne {e}_{2} $，且$ {e}_{1}, {e}_{2}$中有一个为0，设$ {e}_{1}$为码元0，则$ {\varepsilon }_{j}^{0}\oplus {\varepsilon }_{j}^{{e}_{2}}=\underset{j-1\mathrm{ }\mathrm{个}}{\underset{⏟}{0\cdots 0}}{e}_{2} $。
计算首位码元，包括$ {a}_{1}$、$ {b}_{1}$及非首位码元相加对首位产生的进位码元。具体运算规律以表 1加法查找表的形式给出，表 1中的首行及左列分别为相加码元。
表 1 加法运算查找表Table 1. Addition Table of Seven Code Cells相加码元 相加码元 0 1 2 3 4 5 6 10 20 30 40 50 60 100 200 300 400 500 600 0 0 1 2 3 4 5 6 10 20 30 40 50 60 100 200 300 400 500 600 1 1 100 20 2 0 6 10 106 102 200 3 5 600 201 30 4 60 601 2 2 20 200 30 3 0 1 100 201 203 300 4 6 102 302 40 5 10 3 3 2 30 300 40 4 0 1 200 302 304 400 5 20 203 403 50 6 4 4 0 3 40 400 50 5 6 2 300 403 405 500 1 30 304 504 60 5 5 6 0 4 50 500 60 600 1 3 400 504 506 10 2 40 405 605 6 6 10 1 0 5 60 600 601 100 2 4 500 605 106 20 3 50 506 2 编码运算在正二十面体上的扩展
HLQT是平面上的编码结构，正二十面体表面为封闭三角面集合，HLQT编码运算规则不能完全适用，为了避免复杂、低效的跨三角面运算，需要将HLQT编码扩展到正二十面体上, 并设计简洁、高效的编码运算方案。
2.1 正二十面体顶点结构
正二十面体共有20个三角面及12个顶点，每个顶点周围都连续排列有5个共顶点的三角面(图 5(a))，将该顶点结构作为正二十面体表面基础结构。沿任意一条边$ {V}_{1}{V}_{i}(i=\mathrm{2, 3}, \mathrm{4, 5}, 6)$将其剪开，展开至平面，得到平面上围绕公共顶点连续排列的5个三角形(图 5(b)中T1~T5)。图 5(b)等效于从围绕公共顶点连续排列、完整覆盖平面的6个三角形中去除其中1个三角形(从图 5(c)中去除三角形T6)。这一变化过程为将平面格网扩展至顶点结构，建立正二十面体上的格点数学模型提供了思路。
2.2 面瓦片与顶点瓦片
以图 5(c)作为研究对象，作三角面内接六边形，则还存在一个以V1为中心、与6个内接六边形为相邻的顶点六边形。定义三角面内接六边形内格点集合为$ {P}_{n}^{i}(1\le i\le 6)$，其中顶点六边形内格点集合为$ {P}_{n}^{0} $，相邻六边形边界格点集合为$ {E}_{n} $，令以V1为中心的格点集合$ {R}_{n}\prime ={P}_{n}^{0}\bigcup {E}_{n} $。当n=1时，以不同彩色单元区分$ {P}_{1}^{i}(1\le i\le 6) $，如图 6(a)所示，$ {P}_{1}^{0} $与$ {E}_{1} $均用黑色单元表示，$ {R}_{1}\prime ={P}_{1}^{0}\bigcup {E}_{1} $。以图 6(a)中格点\$、$ {e}_{1} $示例$ {e}_{2} $的分布特点，$ {E}_{n} $为$ {e}_{1}$与$ {P}_{1}^{0} $外接六边形相交边界上的格点，$ {P}_{1}^{3}$为$ {e}_{2} $与$ {P}_{1}^{3} $外接六边形相交边界上的格点。按照上述规则定义的格点集合$ {P}_{1}^{4}$与$ {P}_{n}^{i}(1\le i\le 6) $(黑色单元)所覆盖的平面部分，单元之间不重叠且不遗漏。
为建立顶点结构上的格点集合，借助图 5所示的正二十面体顶点结构变化，去除一个$ {R}_{n}\prime {P}_{n}^{i}(1\le i\le 6) $以及${R}_{n}\prime $在Ti内格点，拼接处保留一条边界上的格点。以图 6(a)为例，去除$ {R}_{1}\prime {P}_{1}^{6}$以及\$在T6内格点并保留$ {V}_{1}{V}_{2} $边界格点，结果如图 6(b)所示。将删减后的$ {R}_{n}\prime $内格点集合定义为$ {R}_{n} $，如图 6(c)所示拼合边界，得顶点结构格点集合，如图 6(d)所示。需要注意的是，图 6(d)中正二十面体顶点处六边形单元在拼接后变为五边形。将$ {P}_{n}$与$ {R}_{n} $分别命名为面瓦片与顶点瓦片。
记正二十面体表面的四孔六边形离散格网系统第n层对应的格点集合为$ {G}_{n}(n\ge 1) $，$ {G}_{n} $可分为两个部分：第一部分为20个面瓦片$ {P}_{n}^{j}(j=\mathrm{0, 1}, 2\cdots 19) $，j为面瓦片中心所属三角面编号；第二部分为12个顶点瓦片$ {R}_{n}^{k}(k=\mathrm{0, 1}, 2\cdots 11)$，k为顶点瓦片中心对应顶点编号。即：
$$ {G}_{n}=\stackrel{19}{\bigcup \limits_{j=0}}{P}_{n}^{j}\bigcup \stackrel{11}{\bigcup \limits_{k=0}}{R}_{n}^{k} $$ (6) 式中，$ \bigcup $表示各格点集合的并集。$ {P}_{n}^{j}$之间除了中心不同之外完全等价，只需在格点编码前添加三角面索引编号加以区别；同理，各$ {R}_{n}^{k}$只需在格点编码前添加顶点编号加以区别。顶点瓦片与面瓦片编码方案与HLQT编码方案相同，顶点瓦片只是对面瓦片单元的稍加删减与扩充。
2.3 编码运算规律
面瓦片位于单个三角面内，满足HLQT编码运算规律。顶点瓦片位于5个相连三角面上，运算涉及跨三角面单元的处理。由于顶点瓦片由平面六边形裁剪拼接而得，在运算不经过拼接边界时，仍满足HLQT编码运算规律；当运算经过拼接边界，由于裁剪了一个三角面的编码单元，平面编码运算规律不再适用，运算需跨越一个三角面。
如图 6所示，边界$ {V}_{1}{V}_{2}$上格点在拼接后与$ {V}_{1}{V}_{7}$上对应格点一一重合，例如$ \alpha $与$ \alpha \prime $、$ \beta $与$ \beta \prime $。当运算经过格点$ \alpha $(图 6(a))且向$ {V}_{1}{V}_{7}$方向时，可将格点$ \alpha $转换至$ {V}_{1}{V}_{7} $上对应格点$ \alpha \prime $编码再运算。经验证，两条边界上格点编码对应转换关系可完全通过编码运算描述，设边界上对应格点编码分别为$ \alpha ={b}_{1}{b}_{2}\cdots {b}_{n} $和$ \alpha \prime ={c}_{1}{c}_{2}\cdots {c}_{n} $，码元转换关系如下：
$$ {c}_{1}=\left\{\begin{array}{l}\begin{array}{l}{b}_{1}, {b}_{1}\in \left\{\mathrm{0, 2}\right\}\\ 201, {b}_{1}=203\end{array}\\ \begin{array}{l}20, {b}_{1}=30\\ 2 000, {b}_{1}=3 000\end{array}\end{array}\right. $$ (7) $$ {c}_{j}=\left\{\begin{array}{l}{b}_{j}, {b}_{1}\in \left\{\mathrm{0, 1}\right\}\\ 2, {b}_{j}=3\\ 3, {b}_{j}=2\end{array}\right. $$ (8) 式中，$ 2\le j\le n $。由此可实现顶点瓦片上纯编码跨面单元运算。
以接边处跨面邻近单元搜索为例，图 7为$ {R}_{2}$编码分布，由于首位码元可能不止一位数字，单元内编码用逗号分隔首位与第2位码元，为了便于示意接边编码转换，灰色填充单元表示被剪掉的边界编码。邻近单元跨面的单元有两类：(1)格点位于边界上的单元。例如$ \alpha =\mathrm{2, 2}$，则$ \alpha \prime =\mathrm{2, 3}$，$ \alpha $的邻近单元中有3个(蓝色单元)在拼接后实际为$ \alpha \prime $的邻近单元，根据编码加法运算和式(7)、式(8)求顶点瓦片上$ \alpha $的6个邻近单元为：$ \alpha \oplus \mathrm{5, 1}=\mathrm{3, 3} $；$ \alpha \oplus \mathrm{0, 2}=\mathrm{3, 0} $；$ \alpha \oplus \mathrm{3, 3}=\mathrm{3, 1} $；$ \alpha \prime \oplus \mathrm{5, 1}=\mathrm{1, 2} $；$ \alpha \prime \oplus \mathrm{0, 3}=\mathrm{1, 0} $；$ \alpha \prime \oplus \mathrm{1, 2}=\mathrm{1, 1} $。(2)一条边与边界重合的单元，例如$ \beta =\mathrm{3, 1} $，只有一个邻近单元$ \mathrm{1, 1}$跨面，该邻近单元计算过程如下：$ \beta \oplus \mathrm{0, 3}=\mathrm{2, 2} $；$ \mathrm{2, 2}\stackrel{\alpha \to \alpha \prime }{\to }\mathrm{2, 3}$；$ \mathrm{2, 3}\oplus \mathrm{1, 2}=\mathrm{1, 1} $。其余5个邻近单元可通过正常编码加法运算得到。
此外，还需注意顶点单元为五边形，只有5个邻近单元，加上5个方向单位向量的对应编码即可。
3 格网生成与跨面编码运算效率对比实验
本文在Visual C++ 2012平台上实现了基于正二十面体的四孔六边形格网构造、显示算法。根据用户输入的格网层次$ n(n\ge 1) $，首先生成编码集合$ {P}_{n}$、$ {R}_{n}\prime $，然后根据§2.2的规则去除$ {R}_{n}\prime $中部分编码得到$ {R}_{n} $，接着分别复制$ {P}_{n} $及$ {R}_{n}$到正二十面体20个三角面及12个顶点上得到$ {P}_{n}^{j}(j=\mathrm{0, 1}, 2\cdots 19) $和$ {R}_{n}^{k}(k=\mathrm{0, 1}, 2\cdots 11) $，通过施奈德等积多面体投影生成球面格网，最后显示正二十面体格网与球面格网的三维图形。整个算法流程如图 8所示，运行结果如图 9、图 10所示，$ {P}_{n}^{j} $及$ {R}_{n}^{k}$分别用红色及白色填充。
为验证本文提出的顶点结构在编码跨面运算中的优势，将本文方案与HQBS对比，设计了跨面邻近单元搜索实验，即§2.3所述两类单元邻近搜索方法，并统计各层符合条件的编码个数。全部程序编译为Release版本，并在一台兼容机(软件配置：Windows 7 x64旗舰版+SP1；硬件配置：Intel Core i5-6500 CPU@3.20 GHz，8 GB RAM，KingSton 120 GB SSD)上进行测试。统计两种方案在6~19层的运算效率，即在单位时间内查找邻近单元的平均个数。单位时间设置为1 ms，效率比=$ \frac{\mathrm{本}\mathrm{文}\mathrm{方}\mathrm{案}\mathrm{效}\mathrm{率}}{\mathrm{H}\mathrm{Q}\mathrm{B}\mathrm{S}\mathrm{效}\mathrm{率}}, $实验结果见表 2，效率对比曲线、HQBS效率曲线分别如图 11、图 12所示。
表 2 邻近单元搜索对比实验结果Table 2. Comparative Result of Search for Adjacent Cell层次 单元个数 HQBS效率 本文方案效率 效率比 6 84 168.2 8 886.6 52.8 7 170 318.1 7 746.4 24.4 8 340 329.3 6 297.5 19.7 9 682 324.6 6 027.7 18.6 10 1 364 292.3 5 437.3 18.6 11 2 730 276.6 4 871.3 17.6 12 5 460 255.9 4 500.1 17.6 13 10 922 249.2 3 683.5 14.8 14 21 844 242.8 3 399.0 14.0 15 43 690 232.7 3 346.5 14.4 16 87 382 224.6 3 044.8 13.6 17 174 764 213.6 2 828.0 13.2 18 349 526 203.9 2 556.2 12.5 19 699 052 190.3 2 383.4 12.5 ![]()
图 11 邻近搜索效率对比曲线Figure 11. Comparison of Searching Efficiency Curves for Adjacent Cell分析上述实验结果，得出以下结论：(1)本文方案的跨面邻近单元搜索效率约为HQBS的19.6倍。这主要因为HQBS跨面运算需要借助效率不高的浮点数运算来实现三角面间的旋转、平移，而本文提出的纯编码跨面运算规律只涉及整数编码的直接转换，简单高效，更适合计算机处理。(2)两种方案的编码加法运算效率均随层次升高呈下降趋势。主要原因是：编码长度随层次升高线性增加，编码加法运算需逐个码元处理，编码越长，计算量越大。由于本文方案的绝对效率是HQBS的几十倍，HQBS效率的数值变化在图 11上难以明显看出，图 12表明HQBS的运算效率也呈下降趋势。(3)本文方案效率足以满足绝大多数应用的需求。以对编码运算效率要求较高的实时地形三维可视化为例，对于视点较远的大区域场景，通常采用低分辨率层次格网承载数据，尽管格网单元数目较多，但是效率极高的低层次编码运算可保证处理效率。对于视点较近的局部场景，通常采用高分辨率层次格网承载数据，经过可视区域裁剪后，有效单元数目较少且基本稳定，效率较高的高层次编码运算亦可保证处理效率。
4 结语
本文结合四孔六边形格网系统在正二十面体上的分布特点，利用顶点瓦片与面瓦片，巧妙地将无限平面上的四孔六边形格网编码运算方案退化到封闭正二十面体表面。与现有方案相比，本文方案不仅形式简洁、理论严密，而且具有较高的跨面编码运算效率。本文研究思路可扩展到其他规则格网系统上，为编码运算方案的设计、分析和优化提供理论支撑。
-
![]()
图 3 组合模型空间分布示意图
Figure 3. Schematic Diagram of the Space Distribution of Combined Model
![]()
图 8 研究区域布格重力异常的BEMD分解结果
Figure 8. Decomposition Results of Bouguer Gravity Anomalies by BEMD Method in Study Area
![]()
图 9 BEMD分解各分量的径向对数功率谱曲线
Figure 9. Radial Logarithmic Power Spectrum Curve of Bi-dimensional Intrinsic Mode Function at All Levels by BEMD Method
![]()
图 10 研究区域布格重力异常BEMD分解重构异常
Figure 10. Reconstruction Bouguer Anomalies in Study Area by BEMD Method
![]()
图 11 BEMD分解各级固有模态函数的小波模极大值边界识别结果
Figure 11. Results of Edge Detection of Bi-dimensional Intrinsic Mode Function at All Levels by WMM Method
![]()
图 12 三峡地区布格重力异常BEMD分解重构异常小波模极大值边界识别结果
Figure 12. Results of Edge Detection of the Reconstruction Bouguer Anomalies of BEMD by WMM Method
表 1 组合模型相关参数
Table 1 Related Parameters of the Combined Models
模型编号 （长, 宽, 高）/m 上顶埋深/m 下底埋深/m 剩余密度/(g·cm-3) 模型1 (600, 300, 300) 150 450 0.8 模型2 (650, 300, 250) 100 350 0.6 模型3 (500, 350, 400) 200 600 1.0 
下载: 导出CSV
-
[1] 张壹, 张双喜, 梁青, 等. 重磁边界识别方法在西准噶尔地区三维地质填图中的应用[J]. 地球科学, 2015, 40(3): 431-440 https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201503005.htm Zhang Yi, Zhang Shuangxi, Liang Qing, et al. Application of Boundary Identifying Technologies Using Gravity and Magnetic Maps in Three-Dimensional Geological Mapping of Western Junggar Area[J]. Earth Science, 2015, 40(3): 431-440 https://www.cnki.com.cn/Article/CJFDTOTAL-DQKX201503005.htm
[2] 张双喜, 张青杉, 陈超, 等. 利用磁异常模量进行磁性体边界检测[J]. 地质与勘探, 2015, 51(6): 1025-1032 https://www.cnki.com.cn/Article/CJFDTOTAL-DZKT201506003.htm Zhang Shuangxi, Zhang Qingshan, Chen Chao, et al. Edge Detection of Magnetic Bodies Using the Modulus of Magnetic Anomalies[J]. Geology and Exploration, 2015, 51(6): 1025-1032 https://www.cnki.com.cn/Article/CJFDTOTAL-DZKT201506003.htm
[3] Mallat S, Zhong S. Characterization of Signals from Multiscale Edges[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(7): 710-732 doi: 10.1109/34.142909
[4] 刘金钊, 王同庆, 陈兆辉, 等. 利用滑动窗口的多项式拟合算法进行重力位场区域-剩余异常分离[J]. 武汉大学学报·信息科学版, 2018, 43(10): 1478-1482 doi: 10.13203/j.whugis20170027 Liu Jinzhao, Wang Tongqing, Chen Zhaohui, et al. Regional and Residual Anomaly Separation Based on Robust 2D Polynomial Fitting with Moving Window[J]. Geomatics and Information Science of Wuhan University, 2018, 43(10): 1478-1482 doi: 10.13203/j.whugis20170027
[5] Liu J Z, Li H L. Separation and Interpretation of Gravity Field Data Based on Two Dimensional Normal Space-Scale Transform (NSST2D) Algorithm: A Case Study of Kauring Airborne Gravity Test Site, Western Australia[J]. Pure and Applied Geophysics, 2019, 176(6): 2513-2528 doi: 10.1007/s00024-019-02131-5
[6] Nunes J C, Bouaoune Y, Delechelle E, et al. Image Analysis by Bidimensional Empirical Mode Decomposition[J]. Image and Vision Computing, 2003, 21(12): 1019-1026 doi: 10.1016/S0262-8856(03)00094-5
[7] Nunes J C, Guyot S, Deléchelle E. Texture Analysis Based on Local Analysis of the Bidimensional Empirical Mode Decomposition[J]. Machine Vision and Applications, 2005, 16(3): 177-188 doi: 10.1007/s00138-004-0170-5
[8] 张双喜, 陈超, 王林松, 等. 基于BEMD的小波模极大值方法在位场多尺度边界检测中的应用[C]// 第一届中国地球科学联合学术年会, 珠海, 广东, 2014 Zhang Shuangxi, Chen Chao, Wang Linsong, et al. The Application of Multi-scale Boundary Detection in Potential Field by Wavelet Modulus Maximum Method Based on BEMD[C]//The 1st Annual Meeting of Chinese Geoscience Union (CGU), Zhuhai, Guangdong, China, 2014
[9] 张双喜, 陈超, 王林松, 等. 二维经验模态分解及其在位场去噪和分离中的应用[J]. 地球物理学进展, 2015, 30(6): 2855-2862 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201506054.htm Zhang Shuangxi, Chen Chao, Wang Linsong, et al. The Bidimensional Empirical Mode Decomposition and Its Applications to Denoising and Separation of Potential Field[J]. Progress in Geophysics, 2015, 30(6): 2855-2862 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWJ201506054.htm
[10] 张双喜, 陈兆辉, 王同庆, 等. 利用二维经验模态分解提取川滇地区流动重力异常特征[J]. 大地测量与地球动力学, 2018, 38(4): 407-413 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201804016.htm Zhang Shuangxi, Chen Zhaohui, Wang Tongqing, et al. The Extraction of Abnormal Feature of Mobile Gravity in the Sichuan-Yunnan Region Using BEMD Method[J]. Journal of Geodesy and Geodynamics, 2018, 38(4): 407-413 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201804016.htm
[11] Huang N E, Shen Z, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[J]. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903-995 doi: 10.1098/rspa.1998.0193
[12] Huang N E, Shen Z, Long S R. A New View of Nonlinear Water Waves: The Hilbert Spectrum[J]. Annual Review of Fluid Mechanics, 1999, 31(1): 417-457 doi: 10.1146/annurev.fluid.31.1.417
[13] Nunes J C, Bouaoune Y, Delechelle E, et al. Image Analysis by Bidimensional Empirical Mode Decomposition[J]. Image and Vision Computing, 2003, 21(12): 1019-1026 doi: 10.1016/S0262-8856(03)00094-5
[14] Nunes J C, Guyot S, Deléchelle E. Texture Analysis Based on Local Analysis of the Bidimensional Empirical Mode Decomposition[J]. Machine Vision and Applications, 2005, 16(3): 177-188 doi: 10.1007/s00138-004-0170-5
[15] Ferreira F J F, de Souza J, de S Bongiolo A, et al. Enhancement of the Total Horizontal Gradient of Magnetic Anomalies Using the Tilt Angle[J]. GEOPHYSICS, 2013, 78(3): J33-J41 doi: 10.1190/geo2011-0441.1
[16] 李强, 赵旭, 蔡晋安, 等. 三峡水库坝址及邻区中上地壳S波速度结构[J]. 地震学报, 2011, 33(1): 39-50 doi: 10.3969/j.issn.0253-3782.2011.01.004 Li Qiang, Zhao Xu, Cai Jin'an, et al. S-Wave Velocity Structure of Upper and Middle Crust Beneath the Three Gorges Reservoir Dam and Adjacent Region[J]. Acta Seismologica Sinica, 2011, 33(1): 39-50 doi: 10.3969/j.issn.0253-3782.2011.01.004
[17] 王孔伟, 张帆, 邱殿明. 三峡库区黄陵背斜形成机理及与滑坡群关系[J]. 吉林大学学报(地球科学版), 2015, 45(4): 1142-1154 https://www.cnki.com.cn/Article/CJFDTOTAL-CCDZ201504017.htm Wang Kongwei, Zhang Fan, Qiu Dianming. Relation of Huangling Anticline and Landslide Group in the Three Gorges Reservoir Area[J]. Journal of Jilin University (Earth Science Edition), 2015, 45(4): 1142-1154 https://www.cnki.com.cn/Article/CJFDTOTAL-CCDZ201504017.htm
[18] 吴海波, 申学林, 王杰, 等. 三峡库区上地壳三维速度结构的双差层析成像研究[J]. 地球物理学报, 2018, 61(7): 2802-2814 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201807013.htm Wu Haibo, Shen Xuelin, Wang Jie, et al. Three-Dimensional Velocity Structure of Upper Crust in the Three Gorges Reservoir Area Derived from Double-Difference Tomography[J]. Chinese Journal of Geophysics, 2018, 61(7): 2802-2814 https://www.cnki.com.cn/Article/CJFDTOTAL-DQWX201807013.htm
[19] Spector A, Grant F S. Statistical Models for Interpreting Aeromagnetic Data[J]. GEOPHYSICS, 1970, 35(2): 293-302 doi: 10.1190/1.1440092
[20] Bhattacharyya B K, Leu L K. Spectral Analysis of Gravity and Magnetic Anomalies Due to Two-Dimensional Structures[J]. Geophysics, 1975, 40(6): 993-1013 doi: 10.1190/1.1440593
[21] 张毅, 陈超, 梁青, 等. 三峡地区中上地壳密度结果[J]. 地球科学, 2012, 37(增刊): 213-222 Zhang Yi, Chen Chao, Liang Qing, et al. Density Structure of Upper and Middle Crust in Three Gorges Reservoir Area[J]. Earth Science, 2012, 37(s1): 213-222
-
期刊类型引用(5)
1. 王蕊，贲进，梁晓宇，梁启爽，涂祖瑞. 菱形三十面体格网系统构建的等积投影方法研究. 武汉大学学报(信息科学版). 2025(03): 579-586 . 
百度学术
2. 姜博辉，周为峰. GeoHash、Google S2和Uber H33种全球地理格网编码方法对比分析. 地理与地理信息科学. 2024(02): 19-28 . 百度学术
3. 赵龙，李国庆，姚晓闯，马跃. 正二十面体六边形全球离散格网编码运算. 地球信息科学学报. 2023(02): 239-251 . 百度学术
4. 周建彬，贲进，王蕊，郑明阳. 四孔六边形全球离散格网一致瓦片层次结构编码运算. 武汉大学学报(信息科学版). 2023(04): 639-646 . 百度学术
5. 郑明阳，贲进，周建彬，王蕊. 局部区域多孔径六边形格网系统快速生成算法. 武汉大学学报(信息科学版). 2022(09): 1376-1382 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 895
- HTML全文浏览量: 212
- PDF下载量: 53
- 被引次数: 8
