Modelling Orebody Structures:
Block Merging Algorithms and Block Model Spatial Restructuring Strategies Given Mesh Surfaces of Geological Boundaries

In this paper, the term ‘surface’ encompasses both 2.5D and 3D surfaces. The former refers to open surfaces or warped 2D planes; these generally include mineralisation, hydration and stratigraphic surfaces in a mining context. The latter refers to closed surfaces that envelop a volume in 3D space. Examples include compact 2-manifolds that are topologically equivalent to a sphere. A simply-connected polyhedron surface would satisfy this requirement and have an Euler-Poincaré characteristic of 2. Closed surfaces may represent regions with local enrichment in an ore deposit or pockets with high level of contaminants.

In our envisaged application, the resultant block structure facilitates volumetric estimation of geochemical or geometallurgical properties within a mine. A block-based representation provides spatial localisation while block labelling provides differentiation between geological domains. These properties are important for two reasons. First, it enables spatially varying signals and non-stationary geostatistics to be learned and captured at an appropriate granularity for mining. Individual blocks can be populated with local estimates of the geochemistry, material type composition or other physical attributes. Second, it prevents sample averaging from being applied across discontinuities (different domains) which can lead to incorrect probabilistic inferences and skewed predictions if applied unknowingly across a boundary with a steep attribute gradient. Although polyhedral, topological and surface-oriented models offer some interesting alternatives, block models have been well established and widely deployed within the resource industry. Introducing a radically different approach will likely cause significant disruption and require changes, acceptance and adaptation by various stakeholders. The proposed block model restructuring strategy is designed to work seamlessly with existing models. It makes changes which are transparent and compatible with subsequent processes such as mining exacavation which also operates on a block scale.

The goal is to identify blocks in the input model which intersect with the given surface(s). These represent areas where model refinement is needed in order to minimise the surface approximation error. To establish a sense of scale, the input blocks (also called “parent blocks”) typically measure $25\times 25\times 5$m in an axis-aligned local frame, “axis-aligned” means the edges of each block are parallel to one of the (x, y or z) axes.

Block decomposition is performed on surface-intersecting parent blocks to improve spatial localisation. This process divides each block into smaller blocks of some minimum size (e.g., $6.25\times 6.25\times 1.25$m). These sub-blocks are disjoint and together, they span the whole parent block. Geometry tests are applied to determine which sub-block inside each surface-intersecting parent block actually intersects with a given surface. This process places blocks into one of three categories: (a) parent blocks that never intersect with any surface; or else sub-blocks inside surface-intersecting parent blocks that (b) intersect with a surface or (c) do not intersect with any surface.

This component coalesces sub-blocks to produce larger blocks in order to prevent over-segmentation. Recall that one of the objective is to minimise the total number of blocks in the output model. Therefore, if an array of $3\times 4\times 2$ sub-blocks could be merged into a perfect rectangular prism, it would be more efficient to represent them collectively by a single block with a new centroid and their combined dimensions.

In the event that the sub-blocks (within a surface-intersecting parent block) intersect with multiple different surfaces (say $s_{1}$ and $s_{2}$) then merging considerations will be applied separately to the sub-blocks intersecting with $s_{1}$ and $s_{2}$. It is important to observe that block consolidation is applied to both category “b” and category “c” sub-blocks as defined above. Also, merging of sub-blocks across parent block boundaries is not permitted. In Sec. 3.3, we describe in detail the proposed block merging approach which draws inspiration from the coordinate-ascent algorithm.

As a general principle, block tagging assigns abstract block labels to differentiate blocks located above and below a given surface. The objective is to support tagging with respect to multiple surfaces and two labelling policies. The first policy distinguishes surface-intersecting blocks from non-intersecting blocks. The second policy forces a binary decision and labels surface-intersecting blocks as strictly above or below a surface. The scheme also offers the flexibility of using a surface to limit the scope of an update, thereby leaving previously assigned labels intact above / below a surface. The processes described throughout (Sec. 2.1 – Sec. 2.4) are deterministic and applicable in an iterative setting. This completes our brief overview of the system.

The method we employed for determining if an axis-aligned rectangular prism intersects with a triangular patch from a surface is described by Akenine-Möller in [18]. This approach applies the Separating Axis Theorem (SAT) which states that two convex polyhedra, A and B, are disjoint if they can be separated along either an axis parallel to a normal of a face of either A or B, or along an orthogonal axis computed from the cross product of an edge from A with an edge from B. The technical details can be found in Appendix A.

In essence, block-surface intersection assessment consists of a series of “block versus triangle” comparisons where the triangles considered for each block are selected based on spatial proximity. The block-triangle overlap assessment operates on one basic principle: a “no intersection” decision is reached as soon as one of the tests returns false. A block-triangle intersection is found only when all 13 tests return true, when it failed to find any separation. These tests are used to identify the surface triangles (if any) that intersect with each parent block.

To maximise computation efficiency, each block is tested against a subset of triangles on each mesh surface rather than the entire mesh surface. A kD-tree accelerator (a variation of binary space partitioning tree) [19] is constructed using the 3D bounding boxes of each triangle. The subset includes only triangles whose bounding box overlaps with the block; only these candidates can intersect with the block. This pruning step limits the number of comparisons and speeds up computation considerably. The indices harvested here can subsequently reduce the test effort in the block structure decomposition stage.

The key premise is to decompose surface-intersecting blocks into smaller blocks to improve spatial localisation. The basic intuition is that the surface discretisation error decreases as spatial resolution increases. Precision increases when smaller blocks are used to approximate the surface curvature where the blocks meet the surface.

Block structure decomposition entails the following. For each block ($b$) that intersects with the surface, divide it volumetrically into multiple sub-blocks using the specified minimum block dimensions $(\Delta^{\text{block}}_{x,\min},\Delta^{\text{block}}_{y,\min},\Delta^{\text{block}}_{z,\min})$. The main constraints are that sub-blocks cannot overlap and they must be completely contained by the parent block whose volume is the union of all associated sub-blocks.¹¹1Although we typically require parent blocks $(\Delta^{\text{block}}_{x}[b],\Delta^{\text{block}}_{y}[b],\Delta^{\text{block}}_{z}[b])$ to be divisible by the minimum block dimensions, the method works fine even if fractional blocks emerge during the division, i.e., the last block toward the end is smaller than the minimum size; the ratios in each dimension $n^{\text{block}}_{x},n^{\text{block}}_{y},n^{\text{block}}_{z}\notin\mathbb{Z}$ can be non-integers. This essentially means there is no fundamental limit on the minimum spatial resolution. Within each surface-intersecting parent block, we also identify all sub-blocks that intersect with a surface and which surface they intersect with. This is accomplished using the associative mapping [obtained during block-surface overlap detection] which limits the relevant surface triangles to a small subset for each surface-intersecting parent block. The relevant attributes captured in the output include a list of surface-intersecting parent blocks, and attributes for each sub-block: viz., its centroid, dimensions, parent block index, intersecting surface and position within the parent block.

The consolidation component focuses on merging sub-blocks inside surface-intersecting parent blocks to reduce block fragmentation. Its basic objective is to minimise the block count, although other mining or geologically relevant criteria such as block aspect ratio can be optimised for. These options are considered later in Sec. 6.1 and Sec. 7.1. For now, these high resolution sub-blocks may themselves intersect or not intersect with any surface. This is indicated in the BlockStructureDecomposition result which serves as input. This component returns the SubBlockConsolidation result which describes the consolidated block structure. This encompasses all parent blocks processed, including those which do not intersect with any surface, as well as sub-blocks or super-blocks that constitute the surface-intersecting parent blocks.

The proposed merging algorithm is inspired by coordinate-ascent optimisation and may be summarised as follows.

• •

  The algorithm is inspired by “coordinate ascent” where the search proceeds along successive coordinate directions in each iteration. The goal is to grow each block (a rectangular prism) from a single cell and find the maximum extent of spatial expansion, $\mathbf{k}=(k_{x},k_{y},k_{z})$, without infringing other blocks or cells that belong to a different class.

• •

  For merging purpose, each parent block is partitioned uniformly down to the minimum block size. The smallest unit (with minimum block size) is referred as a cell. Block dimensions are expressed in terms of the number of cells that span in the x, y and z directions. Accordingly, if a block $b$ with cell dimensions $\mathbf{k}=(k_{x},k_{y},k_{z})$ is anchored at position $\mathbf{c}^{(b)}\in\mathbb{R}^{3}$, its bounding box would stretch from $\mathbf{c}^{(b)}-\frac{1}{2}\boldsymbol{\Delta}^{\text{block}}_{\min}$ to $\mathbf{c}^{(b)}+(\mathbf{k}-\frac{1}{2})\boldsymbol{\Delta}^{\text{block}}_{\min}$ for $k_{x},k_{y},k_{z}\in\mathbb{Z}\geq 1$.

• •

  Merging states are managed using a binary occupancy map (boolean 3D array) with cell dimensions identical to the parent block. Given a block anchor position $\mathbf{c}^{(b)}$ with $\mathbf{k}$ initialised to $(1,1,1)$, expansion steps are considered in each direction $\boldsymbol{\delta}\in\mathbb{Z}^{3}$ which must alternate through the sequence $\boldsymbol{\delta}_{0}=(1,0,0)$, $\boldsymbol{\delta}_{1}=(0,1,0)$ and $\boldsymbol{\delta}_{2}=(0,0,1)$.

• •

  A step in the direction $\boldsymbol{\delta}$ is feasible if all of the cells within the bounding box ($\mathbf{c}^{(b)}-\frac{1}{2}\boldsymbol{\Delta}^{\text{block}}_{\min}$, $\mathbf{c}^{(b)}+(\mathbf{k}+\boldsymbol{\delta}-\frac{1}{2})\boldsymbol{\Delta}^{\text{block}}_{\min}$) are 1 (active). In this case, the spatial extent is updated via $\mathbf{k}\leftarrow\mathbf{k}+\boldsymbol{\delta}$. Each iteration steps through $\boldsymbol{\delta}_{0}$, $\boldsymbol{\delta}_{1}$, $\boldsymbol{\delta}_{2}$ in turn. This continues until no expansion is possible in any direction, at which point the centroid and dimensions of the merged block are computed and the corresponding cells in the occupancy map are set to 0 (marked as inactive).

• •

  Sub-block merging terminates for a parent block when all cells in the occupancy map are set to 0.

• •

  The volume of a consolidated block is effectively the outer product $[0,k_{x})\otimes[0,k_{y})\otimes[0,k_{z})$ where $k_{x}$, $k_{y}$, $k_{z}$ each represents some integer multiple of the minimum block size, $(\Delta^{\text{block}}_{x,\min},\Delta^{\text{block}}_{y,\min},\Delta^{\text{block}}_{z,\min})$, with respect to x, y and z.

• •

  Different coordinate directions are used cyclically during the procedure. At all times, it must ensure the expansion does not include alien blocks in the encompassing cube, e.g., an L-shape block within a $2\times 2\times 1$ cube is not allowed.

The algorithm is explained by way of an example in Fig. 3. To illustrate the algorithm, let us refer to blocks that belong to the class under consideration as the active cells. In the present context, this refers to either surface-intersecting sub-blocks, $\mathcal{S}_{\text{intersect}}$, or the non-intersecting sub-blocks, $\mathcal{S}_{\text{non-intersect}}$, inside a parent block. The inactive cells (labelled 0) refer to the complementary set to the active cells (labelled 1).

In Fig. 3, suppose the active sub-blocks (cells) refer to $\mathcal{S}_{\text{intersect}}$. The initial state shows the decomposition of a parent block in terms of three x/y cross-sections. There are $5\times 3\times 3$ sub-blocks (cells) within the parent block and 31 are considered “active”. The white square cells all intersect with the same surface. It so happens the first cell encountered (with index $(i_{Z}\cdot n_{Y}+i_{Y})\cdot n_{X}+i_{X}=0$) is the first active block. The algorithm considers expansion along each axis (x, y and z) in turn. In cycle 1, an expansion step in the x-direction is possible. The progressive expansion of the coalesced block is represented by black square cells. In cycle 2, further growth in the y-direction is not possible due to impediment by one or more “inactive” blocks (denoted by $\mathord{\scalerel*{\cdot}{gX}}$), however, the x and z dimensions each allow one step expansion. In cycle 3, expansion continues in the x direction, resulting in the formation of a $4\times 2\times 3$ super-block. At this point, a coalesced block is extracted as no further growth in any direction is now possible. Subsequently, the cells which have just been merged are marked out-of-bounds. The process repeats itself, starting with the next active cell it encounters in the raster-scan order.

Once all the sub-blocks (cells) are coalesced within the surface-intersecting parent block, we are left with 3 merged blocks (coloured in different shades of red in Fig. 5). These have cell dimensions²²2These represent multiplying factors relative to the minimum block size. ($k_{x},k_{y},k_{z}$) of $(4,2,3)$, $(1,1,3)$ and $(2,1,2)$ and relative centroids of $(\frac{4}{10},\frac{2}{6},\frac{3}{6})$, $(\frac{9}{10},\frac{1}{6},\frac{3}{6})$ and $(\frac{6}{10},\frac{5}{6},\frac{4}{6})$, respectively, with respect to the parent block’s minimum vertex and parent block dimensions.

• •

  When multiple surfaces are considered, the coalesce function will next be applied to other sub-blocks (within the same parent block) that intersect with other surfaces. This seldom happens but it is possible if the surfaces are close.

• •

  Finally, the coalesce function is applied to inactive sub-blocks (within the same parent block) which do not intersect with any surface. The result for inactive blocks (coloured in different shades of blue) are shown in Fig. 5.

To summarise, coalescing is applied only to surface-intersecting parent blocks which are subject to decomposition. The sub-blocks contained within may intersect with one or more surfaces, or not intersect with any surface at all. After consolidation, the output contains both the original and merged blocks. Based on the sub-block and surface intersection status, the consolidated blocks may be arranged into two separate sets, $\mathcal{B}_{\text{intersect}}$ and $\mathcal{B}_{\text{non-intersect}}$, as shown in Fig. 5. Surface intersecting merged blocks ${\color[rgb]{0.6,0,0.3}\blacksquare}{\color[rgb]{0.8,0,0.4}\blacksquare}{\color[rgb]{1.,0.5,0.5}\blacksquare}$ (within surface-intersecting parent blocks) are placed in $\mathcal{B}_{\text{intersect}}$. Non surface-intersecting merged blocks ${\color[rgb]{.75,.75,1}\blacksquare}{\color[rgb]{0,.4,1}\blacksquare}{\color[rgb]{0,0,.4}\blacksquare}$ (within surface-intersecting parent blocks) are placed in $\mathcal{B}_{\text{non-intersect}}$. Non surface-intersecting parent blocks are also appended to this. The coordinate ascent algorithm provides a method for merging cells (sub-blocks of minimum size) within the confines of a parent block. Coalesced sub-blocks must share the same label and form a rectangular prism.

The main points to appreciate in Fig. 6 are that

• •

  Each triangular patch of the surface has a normal vector associated with it. For example, the normal $\mathbf{n}_{t_{k}}$ associated with triangle $t_{k}$ points in the upward direction. This normal is computed by taking the cross product between two of its edges, for instance, $(\mathbf{v}_{k,1}-\mathbf{v}_{k,0})\times(\mathbf{v}_{k,2}-\mathbf{v}_{k,0})$. This arrow would point in the opposite direction (rotate by 180^o) if the edge traversal direction is reversed; for instance, by swapping two vertices $\mathbf{v}_{k,1}$ and $\mathbf{v}_{k,2}$ in the triangle.

• •

  To determine “which side of the surface” a block is on, it suffices to consider the triangular patch located closest to the block. After establishing the positive side as the space ‘above’ the surface, one can say block $b_{1}$ is located above the surface and has a positive sign ($\sigma_{b_{1}}=+1$) because the inner product $(\mathbf{c}_{t_{k}}-\mathbf{c}_{b_{1}})\cdot\mathbf{n}_{t_{k}}$ is negative. Here, $\mathbf{n}_{t_{k}}$ and $\mathbf{c}_{t_{k}}$ represent the normal and centroid of triangle $t_{k}$, likewise $\mathbf{c}_{b_{1}}$ represents the centroid of block $b_{1}$. Conversely, block $b_{0}$ is below the surface and has a negative sign ($\sigma_{b_{0}}=-1$) because the inner product $(\mathbf{c}_{t_{j}}-\mathbf{c}_{b_{0}})\cdot\mathbf{n}_{t_{j}}$ is positive. This ‘projection onto normal’ approach provides the first method for block tagging.

A volumetric block region (one comprising blocks of varying size) may be constructed in multiple passes and labelled sequentially by a single or a combination of surfaces using logical and/or operators based on the plane partitioning property of open surfaces, or expressed as the intersection/union of multiple closed surfaces. As an example, a fault surface may partition a space into two disjoint sets, say, blocks to the left (respectively, right) of the divide. Input blocks in contact with the fault line may be partitioned (subblocked) to conform to the shape of the fault surface. Subsequently, block tagging is applied independently to these two sets. The bounding surfaces that correspond across the discontinuity will provide similar block labelling instructions despite their physical separation. This ensures blocks within the same layer will receive the same label notwithstanding the vertical offset introduced by the fault.

The running example has thus far taken only a regular block structure as input. In this section, we demonstrate that the framework can also modify the spatial structure of a model with irregular (non-uniform) block dimensions. Of particular interest is that tilted surfaces are used to model a hypothetical dyke channel running through bedded layers. This highlights two significant features: (a) ability to iteratively improve the spatial structure of an existing block model whilst preserving the labels for horizontal strata which have been previously assigned; (b) ability to work with oblique surfaces and produce correct result when the surface orientation is ambiguous (positive may not point upward), thus user has to specify precisely what is meant by the positive direction in relation to the supplied tagging instructions. An example of this is shown in Fig. 11.

Fig. 12 shows the recut block model follows the curvature of the tilted surfaces. In (k), blocks within the dyke are removed for clarity. Pre-existing labels (outside the dyke) remain intact. The spatial structure is only modified around the tilted surfaces. In Fig. 12(l)–(m), blocks within the dyke and those located on the east / west interfaces are shown in isolation. Fig. 12(n)–(p) show the existing labels for blocks in A, B, C and D have been perfectly preserved. As expected, new labels — {8, 10} and E respectively — have only been assigned to blocks that intersect with and sandwiched between the tilted surfaces.

In regard to block tagging and ‘which-side-of-the-surface’ determination, the ‘projection-onto-normal’ method described in Section 2.4 works well when the surface is smooth and triangle mesh is dense and uniform. Potential issues arise when the surface exhibits local jitters and the triangles are sparse. For instance, the mesh resolution is low relative to the parent block size when triangle patches stretch over distances of up to one kilometer in certain areas. As an illustration, Fig. 13 shows a block associated with triangle $t_{j}$ (the nearest patch based on block-triangle centroid-to-centroid distance). This association yields the wrong result, a negative sign with respect to the normal $\mathbf{n}_{t_{j}}$ is obtained (according to the plane partitioning test) even though it lies above the surface. We observe the projection of the block lies outside the support interval of the referenced triangle $t_{j}$ and the same comparison with $t_{k}$ which is further away would produce the right result (a positive sign). This problem can be remedied by upsampling the mesh surface to increase its density. Fig. 14 shows an example where triangles are recursively split along the longest edge until the maximum patch area and length of all edges fall below the thresholds of 1250m² and 100m. A better solution, however, is ray-tracing [23]. In general terms, ray-tracing determines whether a block is above or below an open surface (resp., inside or outside a closed surface) by counting the number of intersections between the surface and a ray casted from the block. This recommended approach is described in Appendix B. The key advantage is that ray tracing is not susceptible to variation in surface mesh density (the sign inversion problem due to folding), furthermore, it does not require dense surfaces or consistent (e.g. clockwise) ordering of triangle vertices. Therefore, ray-tracing will be the method of choice for block tagging moving forward. It replaces the projection-onto-normal approach in our proposal.

The block consolidation component, as it currently stands, considers the sub-blocks that belong to the surface-intersecting set ($\mathcal{B}_{\text{intersect}}$) and non surface-intersecting set ($\mathcal{B}_{\text{non-intersect}}$) independently. As Fig. 5 has shown, the sub-blocks (cells) within each respective set are merged separately to form larger rectangular prisms within the confines of the parent blocks. While this split is useful for extracting surface-intersecting sub-blocks, it has two drawbacks. In terms of boundary localisation accuracy, Fig. 15(a) shows that surface-intersecting sub-blocks — with centroids located on different sides of the surface — are merged together irrespective of whether it is predominantly above or below the surface. In terms of compaction, merging potential is limited because adjacent cells from $\mathcal{B}_{\text{intersect}}$ and $\mathcal{B}_{\text{non-intersect}}$ cannot be coalesced even if they lie on the same side of the surface.

To reinforce the first point, the two cells marked with “?” in Fig. 15(a) ought to be labelled as above the surface. However, under the current regime, they are considered jointly with the three cells immediately to the right that also intersect with the surface, thus they are treated collectively as a 1-by-5 merged block. Since the centroid (black dot) lies marginally below the surface, the merged block will also be labelled as such. This distorts the boundary as it introduces a vertical bias of around $\Delta_{\text{y,min}}^{\text{(block)}}/2$ to the two left-most cells.

To accurately localise the boundary and achieve the result shown in Fig. 15(b), ray tracing is used to determine the location of cells with respect to each relevant surface that intersects the parent block. Given $S$ surfaces, the {0=above (or no intersection), 1=below} decisions naturally produce up to $2^{S}$ categories (or states) which would be treated separately during sub-blocks consolidation in lieu of $\mathcal{B}_{\text{intersect}}$ and $\mathcal{B}_{\text{non-intersect}}$. In practice, however, we suggest encoding the location with respect to each surface $s$ using 3 bits $b_{3s+2}b_{3s+1}b_{3s}$, where the mutually exclusive bits are set to 1 to denote the outcomes of above, below and untested, respectively. The rationale is that when $b_{3s}=0$, domain identification needs not be attempted during block tagging with respect to surface $s$, since the decision has already been made here prior to sub-blocks consolidation (either $b_{3s+2}=1$ or $b_{3s+1}=1$) and this information is passed on. In fact, applying ray-casting at the cellular-level (highest resolution) yields more accurate results near the surface than applying to merged blocks, especially for undulating surfaces with high local curvature.

To summarise, administering ray-casting before sub-blocks consolidation helps divide cells along surface boundaries; this increases boundary localisation accuracy and ensures merging is performed within the right domains with maximum potential. Table 1 provides a comparison of the techniques discussed.

The first application we preview relates to block modelling for a typical mine site located in Western Australia. The block model spatial restructuring results shown in Fig. 16 illustrate a block-wise partitioning of the mine site into different domains. This is achieved using surfaces which were created to separate geological domains (mineralised, hydrated and waste) within the Brockman Iron Formation which contains members of interbedded BIF and shale bands in the Hamersley Basin Iron Province [24]. Although ore-genesis theories vary depending on the minerals or commodity, the ability to model formations and features such as igneous intrusions in ore deposits is of general interest in areas not limited to mining, but also in further understanding the structural geology of mineral deposits. Using open and closed surfaces to represent structures of varying complexity — this may encompass volumes with exceptional geochemical or geophysical attributes — it is possible to extract waste pockets with high concentration of trace elements, or regions with magnetic / gravity anomaly [25]. Equally, if the surfaces represent the boundary of aquifers separated by clay and lignite seams [26], the process may serve as a basis for creating a structural hydrogeological model to study hydraulic and transport conditions in geotechnical or environmental risk assessment. The techniques developed for shaping a 3D block model can be used potentially in a variety of contexts, including surface buffer analysis (in GIS and structural modelling) for triangle mesh 3D boundary representation of localised objects [27]. In the next section, we focus on a specific application of block merging to reduce spatial fragmentation in a block model.

The dissolved sub-block boundaries convention can be useful for healing a fractured block model, for instance, over a region where a false geological boundary was given in a previous surface update. Under the dissolved sub-block boundaries convention, coordinate-ascent can start from a clean slate and merge sub-blocks in a fragmented area back to the fullest extent in cases where individual sub-block dimensions or internal boundary alignments are otherwise incompatible. It is not bound by the consequences of prior model restructuring decisions. The technical details are given in Appendix D.7. For readers simply looking for a basic explanation, Fig. 19 reinforces these points and illustrates what healing really means in practice.

Fig. 19(a) shows an existing (poorly estimated) boundary created using limited data at an earlier point in time. Fig. 19(b) shows the existing block model structure created using this surface which unbeknown to the user is not faithful to the actual boundary shown in gray. At a later point in time (perhaps after several months have elapsed), more assay measurements have been taken at previously unsampled locations. These denser observations help improve the estimated boundary contact points and a new (more accurate) surface is produced as a result; see bottom curve in (a). A second iteration of the block model spatial restructuring process is applied using the new surface. This workflow significantly improves the boundary localisation property of the new block model in (c), however it leaves behind remnant sub-blocks from the previous iteration. By relaxing the internal constraints essentially by dissolving the sub-block boundaries, block merging under the dissolved sub-block boundaries convention helps promote healing in fractured areas. As evident in (d), fragmented blocks are coalesced into larger blocks in areas where the misplaced boundary (blue surface) once occupied.

The final block merging strategy (see Algorithm 3 in Appendix E) can optimise results with respect to different objectives. For instance, it can minimise the block count or avoid extremely elongated blocks by optimising for the block aspect ratio. It also exploits symmetry by following different scan patterns to alleviate ordering effects associated with using a fixed starting point. These issues are further discussed in Appendixes D.9 and D.10.

The block merging algorithms are formally described in Appendix E which constitutes part of the Supplementary Material.

Two experiments were performed targeting two different usage scenarios. The first experiment described in Sec. 7 uses only the block merging algorithms to reduce block fragmentation in an existing block model where mesh surfaces are not available and the minimum block size has been fixed by the supplier of the model a priori. This restricted setting requires only the use of the sub-blocks consolidation component depicted in the system diagram (see Fig. 2). The main objective is to analyse its performance in terms of block compaction and multi-threaded execution.

The second experiment described in Sec. 8 provides an end-to-end evaluation of the entire system (from block model creation to block tagging) which applies sub-block consolidation and ray-tracing in an iterative setting. This involves over 80 surfaces and creates approximately 30 different domains. The goal is to compare the proposed strategies with two hierarchical subblocking techniques based on octrees, to highlight the constraints of dyadic decomposition and the importance of inter-scale block merging.

Fig. 20 illustrates the block merging performance of Algorithm 2 vs the baseline with emphasis on block aspect ratios. In this plot, performance is measured in terms of block aspect ratio. The volume-weighted aspect ratio is computed for each parent block and sorted in ascending order. Large aspect ratios correspond to long narrow blocks whereas small aspect ratios indicate more balanced, less extreme block dimensions which are more preferred in a mining context. In this plot, the lower the curve, the more desirable it is. The main observation is that the proposed block merge strategy is superior to the baseline. Even when the coordinate-ascent procedure follows only the standard scan pattern, it produces higher quality merged blocks (dotted line is below the dash line). The performance margin increases significantly when multiple scans are deployed in Algorithm 3 to minimise the aspect ratio explicitly (note: solid line is consistently below the dash line).

Fig. 21 provides an alternative view of the same result from a spatial perspective. A birds eye view and two cross-sections from the proposed method are shown. Blocks are coloured to highlight differences in aspect ratio between the baseline and proposed method. The darker the red, the more superior is the proposed method. Conversely, blue blocks show the proposed method is inferior. The dominance of the red blocks show the proposed block merge method is able to consistently produce blocks with less extreme aspect ratios across the site. A magnified view for two regions of interest are shown in Fig. 22.

To verify these results, we compare the baseline and proposed method for parent blocks with a high contrast in $\log\left(\frac{AspectRatio(\text{baseline})}{AspectRatio(\text{proposed})}\right)$. Fig. 23 shows the sub-block structure within a parent block of interest, it shows from left to right: the input blocks, and merge results from the baseline and proposed methods. The key observation is the presence of thin narrow blocks in the baseline result; these disappeared under the proposed method. Algorithm 3 discovers block merging opportunities missed by the baseline method. It is able to minimise the block aspect ratio by conducting multiple scans during coordinate-ascent.

Fig. 24 shows a second example. In the top panel, blocks are coloured by block volume such that bigger blocks are coloured in warm colours (yellow, for instance). This shows the proposed method is more effective at consolidating smaller blocks into larger blocks. The bottom panel shows the merge result for the same parent block from a different vantage point. Here, blocks are coloured by domain labels. This illustrates an interesting case where there is a limit to how much merging can be achieved because input blocks from different classes are processed independently even when they share the same parent.

Table 2 reports the number of output blocks produced by the baseline and proposed methods. The block count column shows the proposed method is more efficient at coalescing blocks than the baseline method, as the block count under “persistent block memory” is lower. The block count reported for the proposed method under “dissolved sub-block boundaries” is lower still, this is consistent with our expectation and the reasoning given in Section D.7 based on constraint relaxation.

Table 3 reports the processing time for block merging and demonstrates the scalable nature of the proposed algorithm.⁵⁵5Experiments were conducted on a linux machine with 32 CPUs, 8 cores (Intel Xeon® CPU E5-2680 @ 2.70GHz), 64GB RAM and 2G swap memory The parallel nature of the algorithm is a consequence of processing all input blocks that belong to the same parent block together independent of other parent blocks. This allows processing to be compartmentalised. From the description in Sec. 3.3, it is clear that block consolidation and block tagging are applied locally to spatially disjoint regions; this is exploited in a multi-threaded implementation.

Merging under the “dissolved sub-block boundaries” convention given 697,097 blocks (of varying sizes) from the input model, is equivalent to merging 29,154,116 blocks (with single cell dimensions) directly under the “persistent block memory” convention. The problem size of “dissolved”${\ddagger}$ grows by a factor of 41.8, however, the computation time increased sub-linearly by a factor of only 9.5, this is partly due to the removal of sub-block alignment / compatibility constraints when internal sub-block boundaries are dissolved. Of course, the run time still increases relative to “persistent”, ultimately this may be viewed as a trade-off between speed and compaction when Table 3 is viewed along side Table 2. The CPU utilisation and memory use columns are included for completeness, but these figures are otherwise unremarkable. Running 16 threads, the physical memory footprint amounts to $1\%$ of the RAM available.

In Table 4, it can be seen that performing octree decomposition without merging produces a large number of blocks (4241289 for $D\!=\!3$ levels). With octree merging, the block count is reduced to 43.167%. The proposed strategies have a block count of about 25% relative to standard octree decomposition. It is worth mentioning that despite this block count disparity, each block model represents the modelled region with the same fidelity and the total block volume for each domain classification remains the same irrespective of $D$ and the block count. In one sense, the ‘Octree+Merge’ method provides an empirical bound on block model efficiency (as measured by block count) that other strategies can surpass. The impressive efficiency gains achieved by the proposed strategies may be attributed to inter-scale block merging without dyadic constraints. For instance, it is possible to combine two adjacent blocks $\mathcal{B}_{i}$ and $\mathcal{B}_{j}$ in neighbouring octants if their labels are the same and their cell-dimensions, say, $\mathbf{k}^{(i)}=(1,2,1)$ and $\mathbf{k}^{(j)}=(2,2,1)$ are compatible. This would not be possible using the naïve ‘Octree+Merge’ method as it involves blocks from different octants.

From Table 4, it is evident the block count roughly grows by a factor of 4 with each additional decomposition level. This is made clear in Table 5 and expressed as a ratio $B_{d}/B_{d^{\prime}}$ between levels $d$ and $d^{\prime}$. Intuitively, the growth rate relates to the surface contact area with the blocks which scales with $O(n^{2})$ where $n$ represents the number of cells measured in one dimension. Interestingly, the proposed coordinate-ascent block merging strategies have a lower growth rate compared to ‘Octree+Merge’ — moving from 3 to 4 levels, the block count increases by a factor of 3.6408 as opposed to 4.0721; from 3 to 5 levels, the growth rate is 13.8152 as opposed to 16.9611. In this respect, the proposed strategies scale better and they are more efficient at representing the same domain information. We caution however that block modelling is not ultimately about an endless pursuit for increasing precision. Block granularity (which may be specified in terms of the number of decompositions, $D$) is often limited by practical considerations. In our application, individual blocks in the model serve as containers for grade estimation. The blockwise spatial distribution of chemical attributes in turn informs what, where and how ore material are excavated. The heavy equipments used to extract and transport such material operate with physical or kinematic constraints (such as bucket size and turn radius), so the required spatial resolution in the case of mining is probably limited to 1.5m to 6m at best. This naturally imposes a lower bound on the spatial fidelity needed from a model. It informs the choice of $D$ in our evaluation, as $D\!=\!3$ and $D\!=\!5$ yield a resolution of 1.5625m and 6.25m, respectively, in a horizontal plane.

For grade estimation at a given level of spatial fidelity, a lower block count is desirable not so much for storage but from a computational perspective. A compact model contains fewer blocks that require the inference engine to estimate into, and by extension, it has fewer values to update when new information becomes available. For instance, the difference between ‘Octree+Merge’ and ‘Proposed’ can be as much as 50% (based on $D\!=\!5$ in Table 4). A secondary consideration is an optimisation preference for more uniform blocks (as opposed to thin, narrow blocks) that have lower aspect ratios. In a pure octree decomposition scheme, the blocks always maintain the same aspect ratio (in our experiment, this value is 50/20 or 2.5). However, this is obtained at significant cost, resulting in a substantial increase in block count. In Table 6, the main observation is that the weighted aspect ratios in the ‘Proposed-D’ column are a marked improvement over ‘Proposed-P’ when merging constraints are relaxed (when input block boundaries are no longer persistent in the sense defined in Sec. 6.1). It also increases slower than ‘Proposed-P’ as the resolution (number of levels $D$) increases, and tends toward the values in the ‘Octree+Merge’ column.

The diversity of blocks within each model is examined in Fig. 26. The graph shows the cumulative distribution of block dimensions for $D\!=\!3$ decomposition levels, with block volumes arranged in ascending order (from left to right) and sorted by aspect ratio in the event of a tie. The staircase nature of the orange curve for ‘Octree+Merge’ shows that there is a much smaller pool of possible block size permutations produced via octree decomposition with merging. To simplify the analysis, we measure blocks in multiples of $\boldsymbol{\Delta}_{\text{min}}^{\text{D\!=\!3}}=(6.25,6.25,2.5)$, for instance, using $\mathbf{k}=(2,1,3)$ to denote a block with real dimensions $(12.5,6.25,7.5)$. Without merging, octree decomposition alone produces an overwhelming 80.44% of blocks with cell dimensions $\mathbf{k}=(1,1,1)$, while 96.86% and 98.93% of blocks have dimensions at/below $(2,2,2)$ and $(4,4,4)$, respectively. The red curve demonstrates that the blocks in the ‘Octree’ model are indeed highly fragmented. With intra-scale merging, the orange curve shows there is a larger spread of blocks and the proportion occupied by the smallest blocks has substantially reduced. Blocks with cell dimensions $(1,1,1)$, $(1,1,2)$, $(1,2,1)$ and $(2,1,1)$ now occupy the 24.09^th, 31.26^st, 39.98^th and 52.08^nd pecentiles, respectively. Because this strategy does not permit inter-scale merging (outside of an octant), it cannot produce blocks with cell dimensions $(1,1,3)$, $(1,3,1)$ or $(3,1,1)$ and indeed $(1,1,4)$. These lost opportunities are capitalised by the proposed strategies. However, the ‘Octree+Merge’ strategy can still produce blocks of size $(1,2,2)$, $(2,1,2)$ and $(2,2,1)$; the so called 4-connected cells occupy the 58.45^th, 68.14^th and 78.64^th percentiles, respectively.

Moving up the spatial hierarchy, the $(2,2,2)$, $(2,2,4)$, $(2,4,2)$ and $(4,2,2)$ blocks correspond to the 84.01^th, 85.50^th, 87.36^th and 90.13^th percentiles; as opposed to 96.86% at the $(2,2,2)$ stage for octree decomposition without merging. With intra-scale merging, the $(2,4,4)$, $(4,2,4)$ and $(4,4,2)$ blocks occupy the 91.07^st, 92.99^nd and 95.25^th percentiles. Again, blocks with intermediate cell dimensions, such as $(2,3,4)$ and $(3,2,5)$, can only be synthesized by the proposed techniques through inter-scale merging and relaxation of the dyadic constraints. Similar observations hold for more levels.

The octree comparisons with the proposed methods highlight the constraints imposed by hierarchical decomposition and the importance of block merging. In cases where the block structure is subsequently consumed by a block estimation process, minimising the block count (while maintaining the same level of spatial fidelity and domain classification granularity) can substantially reduce the computational effort associated with probabilistic inference of block attributes. The analysis has shown, for instance, when a minimum block size of $(1.5625,1.5625,0.625)$ and 5 decomposition levels are applied to the test site, the block count of the proposed model and ‘Octree+Merge’ model are 20.52% and 41.92% respectively relative to standard octree decomposition without block merging. The ‘Octree+Merge’ approach is probably more indicative of the flexible subblocking approach used in the Datamine software [1] which adjusts the split depending on the angle of intersection of a particular block with the surface controlling the sub-division. It is however unclear if a split extends all the way down within a parent block, whether different permutations and aspect ratio optimisation are considered in Datamine. Returning to the block count comparison, the factor of 2 reduction may be attributed to inter-scale merging and relaxation of dyadic constraints as illustrated in Sec. 8.3. The proposed coordinate-ascent block merging approach also scales better. For instance, the block count growth rate from 3 to 5 levels is kept at 13.8152 as opposed to 16.9611 (see Table 4). From the cumulative distribution in Fig. 26, we observed the block models obtained using the proposed methods are spatially less fragmented.

Although the performance metrics all point in favour of the proposed methods, the real strength of the proposal is its flexibility. First, the block merging algorithm was used on its own to consolidate and heal a spatially fragmented block model given only the parent block size and minimum block dimensions (see Sec. 7). Second, the same algorithm was used as part of a coherent block model spatial restructuring strategy to coalesce sub-blocks partitioned by intersecting surfaces (see Sec. 8). Finally, it provides a successive refinement framework for updating the spatial structure of a block model; this takes a dynamic (non-static) view and emphasizes the evolving nature of the model. In contrast, some commercial softwares [1] produce block models that cannot be altered or easily manipulated. The iterative setting implies the input blocks can have different dimensions. The standard octree paradigm is not equipped to deal with input blocks with varying sizes. As a case in point, given an input block with cell dimensions $(4,2,8)$, a three level dyadic decomposition would need to be applied asymmetrically to each axis, via $d_{1}(z),d_{2}(x,z),d_{3}(x,y,z)$, to not violate the minimum block size (lower bound and integer multiple) conditions. Furthermore, an input block with cell dimensions $(3,5,2)$ would leave the octree approach with no viable sub-blocking options unless the block is divided non-uniformly which would set it on a similar path as the proposed strategy notwithstanding the hierarchical embedding (rigid subblock alignment) constraints.

A crucial part of ‘iterative refinement’ is that large portions of an existing model should remain intact, useful information in unaffected parts of the model are not invalidated through the update process, and changes affect only localised regions that interact with new surfaces. This allows creation of a new block model from scratch given a set of surfaces (representing geological boundaries or other types of delineation), as well as automatic correction for misplaced boundaries within an existing model. In [31], the current proposal is portrayed as a system within a system, where new information captured by sampled assay data are used to manipulate the shape of relevant surfaces to maximise agreement with the latest observations. These are used in turn to update the block model to improve the accuracy of its domain structure by maximising the resolution and positional integrity of the blocks near those surfaces. This has been shown in [31] to have a positive effect on processes downstream. In particular, grade estimation performance improves as boundary (smearing) artefacts are reduced during inferencing.

The overall contributions of this paper may be summarised as follows.

1. 1.

   Developing a flexible framework capable of altering the spatial structure of a block model given boundary constraints in the form of triangle mesh surfaces and a set of tagging instructions.

2. 2.

   Describing computational techniques for all system components

   • •

     From block-surface intersection detection (identifying areas requiring an update), block structure decomposition (improving boundary localisation), block consolidation (increasing compactness of representation) to block tagging (encoding the location of blocks relative to the given surfaces).

3. 3.

   Devising a new algorithm (coordinate-ascent inspired block merging) that is amenable to multi-threading.

4. 4.

   Demonstrating real applications: The algorithm was used firstly to reduce block fragmentation in an existing mine resource estimation block model where parameters are fixed and no surfaces are given; secondly as an integrated component within the block model creation/augmentation workflow to coalesce subblocks partitioned by intersecting surfaces.

5. 5.

   Providing insights through quantitative analysis, highlighting the constraints of dyadic hierarchical decomposition and importance of inter-scale merging for achieving a compact block model; discussing why this matters for subsequent block-based estimation processes (see Sec. 8.1–8.4).

The first point includes correcting the position of previously misplaced boundaries in an existing model (see Fig. 19). Philosophically, it challenges the perception that voxel-based models are hard to modify once they have been constructed. It demonstrates that local refinement (including consolidation of fragmented blocks from previous updates) is both possible and can be done efficiently in an iterative setting; with changes to the model confined to areas that interact with new surfaces. The second point is about greater transparency. In 3D geological modelling papers, models are often generated using commercial software where the fundamental techniques employed are unclear. This paper sheds light on this and discusses pit falls and robust solutions from an engineering perspective. The third point emphasizes flexibility. Depending on user requirement on whether internal subblock boundaries are to be respected, it can operate in two modes (see Fig. 18): with persistent block memory to ensure an input block always maps uniquely to a merged output block, or with dissolved subblock boundaries to achieve higher block compaction (a lower block count). The specific objective in the fourth point was to capture geological structures of an ore body in a block model at a specified resolution with a minimal block count. However, the method can potentially be applied elsewhere.