Multi-Scalar Oblique Photogrammetry-Supported 3D webGIS Approach to Preventive Mining-Induced Deformation Analysis

Abstract

Underground coal mining will inevitably cause serious ground deformation, and therefore, preventive mining-induced deformation analysis (MIDA) is of great importance in assisting mining planning and decision-making. Current web-based Geographic Information System (webGIS)-based applications usually use 2D GIS data and lack a holistic framework. This study presents a multi-scalar oblique photogrammetry-supported unified 3D webGIS framework for MIDA applications to fill this gap. The developed web platform uses multiple open-source JavaScript libraries, and the prototype system provides user-friendly interfaces for GIS data collecting and corresponding database establishment, geo-visualization and query, dynamic prediction, and spatial overlapping analysis within the same framework. The proposed framework was tested and evaluated in the Qianyingzi mining area in eastern China. The results demonstrated that multi-scalar oblique photogrammetry balances data quality and acquisition efficiency and provides a good source of GIS datasets, and the web-based platform has a good absolute and relative spatial accuracy verified by two types of validation data. Practical application results proved the feasibility and reliability of the system. The developed web-based MIDA prototype system attains an appealing performance and can be easily extended to similar geoscience applications.

1. Introduction

Long-term and intensive underground mining activity in China creates space, lowering the groundwater level of the area and contributing to the surface subsidence of the mining area. Mining subsidence, a common mining-induced hazard, can lead to geological and ecological environment issues, such as farmland degradation and house or infrastructure destruction [1,2,3]. On the other hand, similar to landslide and earthquake hazards, mining-induced disasters are difficult to monitor and control in real time. Consequently, mining-induced deformation analysis (MIDA), which provides early warning information about these damages, is becoming an extremely important aspect of mining activities. MIDA aims to observe the impacted range intuitively and to produce quantitative assessment results of the ground objects caused by underground coal mining. To this end, MIDA-related applications based on Geographic Information Systems (GISs) is commonly used for this work [4]. A GIS is a comprehensive framework that is used to capture, store, query, analyze, and display spatial data in various applications [4]. Several studies have been reported to assess mine subsidence damages from different perspectives. For instance, Cai et al. [5] developed an integrated system based on an ArcGIS platform to predict subsidence due to underground mining. Malinowska and Hejmanowski [6] studied building damage assessment in mining surfaces based on building strength and terrain deformation predicted by the subsidence computation model (Polish model) in a GIS environment. Duarte et al. [7] took a self-burning coal waste pile as a research object to monitor the temperature variations by producing dynamic maps based on a GIS open-source environment—GIS-ECOAL. However, the above MIDA systems were generally used as standalone and desktop-based applications and thus are limited in data-sharing and multi-source data integration and data processing.

In order to ensure high efficiency, the web-based GIS technology (webGIS) can be used to analyze the involved data effectively. The construction of the webGIS breaks through the visualization and analysis means based on data and the web and helps people to plan, make decisions, and complete many kinds of collaborative work [8]. Specifically, the webGIS uses the Internet to express GIS-related functions reasonably and overcomes the limitation of the GIS functions related to local area networks. Therefore, the method has a wide range of applications. In recent years, webGISs have developed rapidly, and many representative studies have emerged, such as Bing and Google Maps [9]. The system represents the description of geographic information in graphic form and has developed into a network computer with good ability and good performance. For instance, Duarte et al. [10] created a susceptibility map of contamination in the areas surrounding the self-burning coal waste pile using an Analytical Hierarchy Process (AHP) and Fuzzy AHP approaches and further developed a webGIS application incorporating all the needed information in a GIS environment. Cao et al. [11] constructed the Trace Elements in Chinese Coals Database Management System (TECC) and analyzed the spatial distribution of harmful trace elements by applying spatial data algorithms and the visual technology of the webGIS.

The existing webGIS-based MIDA systems discussed above usually exploited 2D data, i.e., raster data with time series remote sensing images or various observed computer-aided design (CAD)-based vector data, to express spatial information. However, because MIDA results feature multi-variable, high-dimensional, and time-varying elements, it is a great challenge for non-professional personnel to analyze and apply these 2D-based MIDA system results. Actually, MIDA-related GIS data or analysis results should be represented in a three-dimensional (3D) visual environment and hence, MIDA-related data visualization and analysis naturally requires a 3D or stereoscopic display environment. Owing to the rapid development of the Web Graphics Library (WebGL), the 3D visualization function of Web browsers has been greatly expanded. Due to no plug-in requirements, cross-browser compatibility (Chrome, Firefox, Internet Explorer, etc.), graphical processing unit (GPU) program drawing, and other features, the WebGL is very suitable for network 3D display. 3D visualization technology based on the WebGL has been developed rapidly [12]. The “virtual earth” technology based on the WebGL has opened up a new way for the display and processing of GIS data. For example, 3D virtual visualization technologies (such as Google Earth, NASA World Wind, and Cesium), which integrate geo-spatial data from all over the world, enable researchers to interact in real 3D virtual scenes and it has proven to be a very ideal visual environment [9]. Many studies have pointed out the role of 3D visualization technology in the fields of environmental protection [9,12], hydrology [13], and oceanographic science [14]. These studies portray a different hazardous event’s occurrence, development, and evolution by a professional numerical simulation model at different scale levels. Although the flowcharts of these studies are similar to the MIDA applications, they do not fully cover GIS data collecting, storage, and spatial analysis, and thus, GIS functionality is not deeply involved in these studies. It should be noted that, recently, the 3D webGIS-based approach to preventive conservation of heritage buildings has become a hot topic [15,16,17,18]. These studies took full use of GIS and WebGL visualization techniques by integrating multiple sources, i.e., building information modeling (BIM), 3D point clouds, real-scene 3D models, 360 panoramas, or the data of the Internet of things (IoT), to carry out building information management and health monitoring.

Good resources of 3D GIS data are necessary to construct a 3D webGIS application. 3D modeling based on real scenes is an important data source in GIS research, which mostly comes from oblique images and laser scanning. The joint application of the structure from motion (SfM) technique and unmanned aerial vehicles (UAVs) makes it possible to reduce cost and time. As a new remote sensing method, it has become one of the hotspots in geoscience research [17]. It includes some computer vision algorithms, such as feature extraction and bundle adjustment (BA) [19,20]. So far, SfM-based 3D reconstruction has been widely used in many fields, such as geological disaster monitoring [21,22], open-pit mining planning [23], cultural heritage protection [24,25], and the structural health monitoring of infrastructures (such as dams) [26]. On the other hand, light detection and ranging (LiDAR) scanning is among the most accurate remote sensing methods available to acquire spatial data, ensuring measurement accuracy at a few millimeters or centimeters [27]. Laser scanning techniques can be divided into two types: space-based and ground-based. Ground-based LiDAR is a commonly used method to measure the three-dimensional spatial structure or geometric shape of objects, which includes the monitoring of building component deformation and surface deformation [28]. However, ground-based or airborne-based LiDAR techniques are not always feasible for the large-scale scenes investigated in this work due to high investment or leasing costs and inevitable occlusions that restrict their broader applications. In some cases, ground-based LiDAR is often adopted only as a complementary method of oblique photogrammetry to acquire more detailed information about the occluded object [28]. A few studies have been conducted on web-based systems to integrate, process, and disseminate LiDAR data in a 3D webGIS framework [29,30]. Further studies concentrated on cultural resource protection by importing more data types (such as 360° panoramas, remote sensing, and the Internet of things) into the 3D webGIS application [16] functionalities.

The development of 3D data collecting methods and 3D webGIS applications in various fields has made great progress, all of which inspired us to construct a 3D webGIS-based MIDA application. To our best knowledge, few researchers try to develop a unified framework that couples data collecting, storage, visualization, dynamic prediction of mine surface deformation, and spatial analysis into a 3D webGIS-based application for a MIDA application. Similar research was carried out by Gábor et al. [8]; they established a 2D and 3D webGIS application in the service of the open-pit mining industry by integrating photogrammetric models, laser-scanned point clouds, and geological strata data, in which the basic GIS functionalities were provided, such as measuring geographical data, drawing 2D profiles, and querying spatial databases. Their webGIS system only used small-scale open-pit mines and lacked a MIDA application. Hence, this study provides early warning information caused by mining surface deformation. There is no custom network geographic information system, only a few different, difficult to learn, and separate software can achieve such requirements. It is because of this shortcoming that we have formulated the proposed system.

The remaining contents of this study include the following: Section 2 describes the system design and components of the proposed methodology and the related technologies of this web-based framework; Section 3 evaluates the proposed system on case-based application tests to demonstrate its capabilities, feasibility, and reliability; Section 4 summarizes and discusses our findings; and Section 5 summarizes our work.

2. Materials and Methods

2.1. Study Area

The Qianyingzi coal mining area in Suzhou, Anhui province of eastern China (Figure 1) was selected as the case study. Its geographical location is 116°20′ E–116°57′ E, 32°43′ N–32°56′ E with a north-south span of 10.8 km and an east-west span of 9 km, which results in a total area of about 97.2 km² within the mining boundary. Most importantly, as a part of the Huaibei plain of China, the mining area is a densely populated and important agricultural region, making it an important grain production base in China. The coal mine was built in the 1980s and has been mining for more than 40 years. A large depression appeared on the ground because of the long-term mining, directly leading to serious damage to the ground scene.

The concerned mining face is named W3233 and is depicted in Figure 1c. Specifically, Figure 2 gives the actual position of it from underground view mode, which can be accessed and viewed from our 3D webGIS platform. The features of this working face are as follows: the strike length and dip length are about 3000 m and 360 m, respectively, and the average dip angle of the coal seam α is 2°. The mining depth H₀ is approximately 460 m, and the average mining thickness m is 3.32 m. This mining working face was planned and started on 11 July 2022, and will be stopped on 15 October 2025.

Table 1. Prediction parameters related to mining working face of W3233.

q	tanβ	θ0 (°)	S1 (m)	S2 (m)	S3 (m)	S4 (m)
0.81	2.01	89.42	19.41	81.25	16.53	−11.58

presents the detailed prediction parameters related to this mining working field, comprehensively deduced from previous field observation data and used for surface deformation estimation. In Table 1, q, tanβ, and θ₀ stand for subsidence coefficient, the tangent of the main influence angle, and the mining transference angle, respectively, and (S₁, S₂, S₃, S₄) is the left, right, upper, and lower deviation of the inflection point [1].

2.2. Architecture of the Proposed 3D webGIS-Based MIDA Framework

There are three ways of development for network applications: local architecture, remote architecture, and hybrid architecture [9]. The hybrid architecture combines the advantages of local architecture and remote architecture, so it can take advantage of the more powerful hardware provided by the server to improve the performance of network applications and enable users to interact, visualize, and perform spatial analysis in the user’s browser. In view of the large amount of operation and computation of the current network MIDA system, a hybrid structure method is adopted to realize it. Figure 3 outlines the overall architectural components and the relationship between them, which are similar to the solutions of most network applications. The model is divided into four layers: data layer, model layer, application server layer, and presentation layer (client layer).

On the server side, the data from different sources, different types, different formats, and different scales are managed. In My Structured Query Language (MySQL) (Oracle Corporation, Austin, TX, USA) [31], the general open-source relational database management system (RDBMS) is used to store and manage the GIS entity data. Relational database (RDB) is a database comprising relation variables, which is straightforward-to-represent data in tables [32]. Access to data occurs through database access interfaces (e.g., data source adapters), web service interfaces, and open standards, such as web feature service (WFS) or web map service (WMS), that are set forth by the Open Geospatial Consortium (OGC) [9]. Large files, such as three-dimensional entities and entity data, can be stored directly in the file system. The input and output files of the application, depending on case, were arranged and encapsulated into the project directory specified on the server. It centralizes all kinds of materials used in MIDA system into the same folder (or subfolder), which can not only improve the efficiency of file management, but also facilitate users reloading files. The module level (MIDA engine) was configured on the server side, providing computing services, including the probability integral steady-state prediction model (PIM) and GIS spatial analysis model. Web application middleware was introduced to the server level to host a variety of dedicated web services for Hypertext Transfer Protocol (HTTP) requests and replies. This project was implemented through representational state transfer (REST) application programming interface (API) and event-based execution environment. Based on Microsoft Internet Information Service (IIS) 10.0, Node, JavaScript provides users with custom network services [28]. In addition, interactive data languages (IDL) and GeoServer [29] were configured at this level to realize automatic download and publication of data on the Internet via standard OGC web services. The presentation layer was a web graphical user interface (GUI) designed by HTML5, JavaScript, and third-party libraries in a Visual Studio Code environment [31,32,33,34,35]. All the technologies used in this study are open source or free development and do not require any additional programs to be installed on the user’s computer which mainly include Cesium adopted as webGIS API [33], Turf.js used for calculation and analysis [35], MySQL used for database management [31], and GeoServer served as image service management [34].

2.3. Data Source Overview

The proposed framework mainly involves two types of GIS datasets: model and entity datasets. The former is further classified into two kinds, i.e., the realistic 3D model of the study area obtained by oblique photogrammetry and the virtual 3D model reconstructed with third-party software supported by photogrammetric data. As for the latter, in the GIS domain, entity data were first classified into three geographic types according to their geometry shape (i.e., point, line, or polygon type). Further, each entity type in our case study included several object types and the conventional 2D vector-based data model was adopted, opposed to the raster-based model for remotely sensed imagery or triangle-based 3D models. Additionally, this study uses a geospatial information data exchange format based on JavaScript Object Notation (GeoJSON) [36] format for data transmission and storage. Object data is divided into primary key, geometry data, and attribute data (

Table 2. Data table format for storing a certain type of GIS entity data.

GIS Entity Type (Point, Line, or Polygon)
ID (Primary key)	Geometry Data (GeoJSON)	Attribute Data (GeoJSON)

), where IDentifier (ID) is a database self-growing sequence. The geometric data consists of geometric types (points, lines, polygons) and coordinates in GeoJSON format, while the attribute data consist of key value pairs for names and values and thus, it can be easily stored in a database for interchange between documents and queries. Figure 4 gives an example of GeoJSON format for a point feature structure.

Our study considers 2D-shaped geographical features. The real 3D features in the Cesium environment mainly correspond to 3D photo-realistic and 3D virtual models, respectively. On the other hand, for each GIS dataset type, its data collecting method and the raw outputted data format are heterogeneous. Thereby, disseminating these geospatial data with online map services initially requires the creation of GIS data layers and relies on the map servers to render them in web-friendly data formats (such as 3D Tiles [37]) and then publish the data as WMS services for rendering and visualization in the web application’s map view. Therefore, data format conversion must always be addressed for different GIS dataset types. A brief overview of the dataset involved in our case study is tabulated in

Table 3. Brief overview of data source involved in our case study.

Dataset Type	Data Layer	Geometry Type	Attribute Data	Collecting Method	Data Format	Storage Method
3D model	3D photo-realistic model (for the study region)	3D triangle	—	Multi-scalar oblique photogrammetry	3D Tiles	File system
	3D virtual model (for power towers)	3D triangle	—	Photogrammetric data-supported 3D virtual modeling	glTF
	3D virtual model (for power lines)	3D vector	—	Measured from 3D virtual model manually	GeoJSON
GIS entity	Power tower	2D point	Administrator, voltage, etc.	Monomer modeling from a realistic 3D model based on the Cesium virtual globe	GeoJSON	MySQL database
	Road	2D line	Name, level, width, material, etc.
	River/ditch		Name, length, width, etc.
	House	2D polygon	Ownership, area, structure, etc.
	Farmland		Ownership, area, usage, etc.
	Mining working face		Coal thickness, prediction parameters, etc.	Extracted from planned CAD drawings	GeoJSON	MySQL database

. Detailed information on them is introduced below.

2.4. Multi-Scalar Realistic 3D Modeling and Data Format Conversion

2.4.1. Multi-Scalar Oblique Photogrammetry

The high precision, high fidelity, and measurability of the real-scene 3D model produced by multi-view oblique images is an important progression in GIS research. In recent years, this method has been widely used in intelligent city management, emergency command and rescue, land surveying and mapping, and many other aspects. In the implementation of aerial survey, it is more convenient to obtain real-scene images, which greatly reduces the cost and working time. Currently, image-collecting devices used for oblique photogrammetry can be coarsely classified into two types: fixed-wing UAVs, usually equipped with an oblique camera of five tilt lenses, and multi-rotor rotatory UAVs with an on-board single or oblique camera. Fixed-wing UAVs usually have higher image acquisition efficiency due to their high flying speed and larger spatial coverage of the target area [38]. However, the exposure delay when a camera takes a picture requires a relatively long time (~5 ms) from focusing photosensitizing to digital storage. Therefore, the value of the ground sample distance (GSD) with fixed UAV is relatively high (usually higher than 4 cm) due to the camera’s exposure time limitation interval under higher flying speed, degrading the quality of the resulting 3D real-scene model. The multi-rotor UAV is contrary to that of the fixed-wing UAV.

Considering the factors discussed above and in balancing the efficiency, accuracy, and texture quality of the produced 3D realistic model, this paper proposes the following multi-scalar oblique photogrammetric scheme:

• Coarse-scale: For the entire study domain with a large area, we use a fixed-wing UAV equipped with an oblique camera with a five-tilt lens to collect the oblique images. This hardware setup enhances data acquisition efficiency for large-scale areas without losing the geometric accuracy required by practical applications. The coarse-scalar photogrammetry in our case study is mainly used to reconstruct the topographic scene with surface continuity and no large terrain undulations.
• Meso-scale: The mining area’s scattered villages and industrial sites are crucially concerned by mining enterprises and local government due to their social and economic importance. These regions of interest (ROIs), labeled with yellow frames in Figure 1, mainly include houses and industrial equipment above the ground surface possibly affected by mining-induced surface deformation. For this scale, a heavy-lift multi-rotor UAV equipped with a camera of five lenses can be used to obtain accurate 3D models to meet the subsequent GIS requirements.
• Fine-scale: For some objects (e.g., tall buildings), obtaining their complete 3D models with coarse-scale or meso-scale or even by combining the above two due to flight height limitations or occlusions is challenging. Hence, the lightweight and low-cost rotatory UAV equipped with a single camera of adjusted shooting angles should be adopted to capture more multi-view photos of this type of object, enhancing the completeness and quality of the corresponding 3D photo-realistic model.

The above three scales complement each other, and a high-quality 3D model can be obtained by fusing multi-scalar and multi-view images. Figure 5 depicts an optional scale combination strategy including the specifications of the equipment, flight plan parameters including flight altitude (H), GSD, and the overlap in the forward direction (O_F) and side direction (O_S) for each scale, which was also chosen as a feasible solution to multi-scalar photogrammetry in our case study.

Concerning the combinative photogrammetry scheme presented in Figure 5, for the coarse-scale and meso-scale, both fixed-wing UAVs (iFly V5) and four-rotor UAVs (M300 RTK) are equipped with an oblique camera of five tilt lenses, which were used to capture the oblique images at a fixed height and five fixed shooting angles. This configuration takes a zigzag flight route to realize oblique photogrammetry (Figure 6a). As for the quadcopter platform (Mavic M3E) mounted with a single camera used in the fine-scale, it can also take multi-view oblique images by adjusting the viewing angles every time we conduct the flight mission (e.g., 90° with nadir viewpoint, 45° with slanted view) according to the same predefined flight path. In this context, the photos were taken in parallel stripes arranged transversely and longitudinally to the boundaries of the target area, creating a grid of stripes parallel and perpendicular to each other (Figure 6b). However, apart from the cross-flight path, an important feature of Mavic M3E is that it can adopt a circular flight path around the target object at multiple horizontal levels to capture the demanded oblique images with a changeable shooting angle automatically or manually (Figure 6c). This flight path can yield a dominant resolution at the millimeter level and create a high-quality 3D model [39]. On a fine scale, the Mavic M3E can capture images of tall buildings with high image resolution.

2.4.2. SfM-Based Photogrammetric Processing

The images were processed using the SfM-based photogrammetry algorithm [20]. The SfM was applied to obtain camera parameters and reconstruct the 3D model from a set of feature matches simultaneously. Indeed, if the calibration matrix K is known, the rotation matrix R and translation vector t can be deduced. Given m perspective images of n fixed points, the 2D image projection is related to the 3D point in the scene and can be calculated by triangulation with the calibration, rotation, and translation matrices, as follows [20]:

$${\stackrel{\mathrm{~}}{x}}_{ij}={ M}_{ij}{\stackrel{\mathrm{~}}{X}}_j, i=1,\dots ,m; j =1,\dots ,n$$

(1)

where ${\stackrel{\mathrm{~}}{x}}_{ij}$ is the image point coordinate of the 3D point ${\stackrel{\mathrm{~}}{X}}_j$ projected onto the ith camera image plane and ${\stackrel{\mathrm{~}}{x}}_{ij}~{\left[u_{ij}, v_{ij}, 1\right]}^{\mathrm{T}}$ and ${\stackrel{\mathrm{~}}{X}}_j ~{ \left[X_j, Y_j, Z_j, 1\right]}^{\mathrm{T}}$ are homogeneous coordinate versions, respectively. The symbol$~$means equality up to scale. M is the projection matrix, expressed as:

$$M =\left[\begin{array}{cc}K& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^{\mathrm{T}}& 1\end{array}\right]$$

(2)

To minimize the distance between the projected and observed image points, an optimization method known as the BA algorithm is used over the collective set of camera and scene parameters, which is expressed as:

$$g\left(C, \stackrel{\mathrm{~}}{X}\right)=\sum _{i=1}^n\sum _{j=1}^m{\omega }_{ij}{‖{\stackrel{\mathrm{~}}{x}}_{ij}-P\left(C_i, {\stackrel{\mathrm{~}}{X}}_j\right)‖}^2$$

(3)

where ${\omega }_{ij}$ is an indicator variable: ${\omega }_{ij}=1$ if camera $i$ observes point j, otherwise, ${\omega }_{ij}=0$. C is the collection of camera parameters for a single camera with $C_i=\left\{M, R, t, K\right\}$, P is the projection equation, and $P\left(C_i, {\stackrel{\mathrm{~}}{X}}_j\right)$ is the projected image point.

The objective function g is the sum of the squared projection errors. BA is used to calculate the camera calibration parameters and posture of images that minimize this objective function. This procedure is repeated until the parameters are suitable for all photos, and the results can be used to generate a relatively sparse set of 3D points. After BA, a dense point cloud is constructed using multi-view stereo (MVS) algorithms that create a 3D mesh with triangles and texture into a 3D photo-realistic model.

On the other hand, before conducting the flight mission, a specific number of ground control points (GCPs) must be deployed on the target region as a constraint to refine the model reconstruction in the SfM processing step and to generate the UAV-based model with the geolocation and scale. At the same time, a certain amount of GCPs should be selected as checkpoints (CPs) to assess the accuracy of the final 3D model. The 3D location (X, Y, H) of GCPs is measured by real-time kinematic (RTK) mode based on global navigation satellite system (GNSS) in the predefined coordinate system, and the observation time at each point is about 15 s. This topographic and geodetic instrument measures the coordinates of the center of the GCPs with a horizontal precision between 0.014 m to 0.030 m, and vertical accuracy between 0.030 m to 0.050 m at each point. To ensure this level of accuracy, RTK corrections are carried out by using the closest national permanent GPS reference stations around the research site [40]. In this study, ($X, Y$) is referred to as Gauss planar rectangular coordinate in the China Geodetic Coordinate System 2000 (CGCS2000) framework [41]. H is the leveling height of a ground point relative to the geoid (sea level) that is defined based on the 1985 National Height Datum in China [41]. For clarity, hereafter, the 3D spatial reference coordinates represented by (X, Y, H) are called surveying coordinates (SC). The resulting photogrammetric data products have the same coordinate system as the GCPs.

2.4.3. 3D Tiles Generation of the Photo-Realistic 3D Model

Oblique photogrammetry modeling can generate multiple formats of models. The commonly used 3D realistic model is the OpenSceneGraph binary (OSGB) format. This data format cannot be directly loaded in the Cesium globe. Therefore, to convert 3D model data from the OSGB format to the 3D Tiles standard, the osgbTo3d-tiles tool is used. The converted results were placed in the appropriate subfolder of the 3D model dataset and uploaded to a web server, from which any application compatible with this standard may access it.

2.5. 3D Virtual Reconstruction of the Power Corridors

Another type of 3D model data relates to power corridors, i.e., high-voltage power towers and associated power lines listed in Table 3, which are vulnerable to mining-induced surface deformations. However, although the multi-scalar oblique photogrammetry technique obtains a high-resolution photo-realistic 3D model, it cannot effectively cope with objects having elongated and poor texture structures in the aerial images. Therefore, such objects will greatly reduce the 3D scene quality and visual effects. Recently, an active remote sensing technique exploiting airborne or mobile laser scanning with dense point clouds and more detailed information has gained much interest as a complementary method to photogrammetry. This method is promising for power corridor mapping and 3D modeling. However, the automatic reconstruction of power corridors from a raw, dense point cloud is challenging due to its incompleteness, noise, and anisotropy in density [42,43,44,45,46]. On the other hand, for many power corridors distributed in our case study region, an extra laser scanning campaign solely used to reconstruct the power corridors is not a cost- and time-effective method. To address these shortcomings, we present a simple yet effective approach to reconstructing the power corridors supported by photogrammetric data. Our method reconstructs the power towers and power lines separately based on different technical routes. Taking the power towers typically presented with four tower foundations as examples, Figure 5 introduces our method’s workflow diagram. The details of the workflow are outlined below.

2.5.1. Collecting the Description Data of Power Corridors

We used a third-party 3D modeling software, 3ds Max (Version 2021), to reconstruct the 3D virtual models of the power corridors. The required size data were obtained using the linear triangulation (LT) method typically employed in photogrammetry and computer vision applications [47]. Then, the parameters needed by Cesium for rendering and 3D visualization, i.e., 3D locations, posture angles, and hanging points of the power towers, were precisely collected. The implementation procedure is discussed below.

• Determining the 3D location of a tower in Cesium environment: The 3D location of a tower is defined as the geometric center (red dot in Figure 7a) of the bottom profile determined by its four-tower bases. Although the reconstructed 3D tower structure elements are not complete or clear, the positions of the tower bases can be recognized and measured in their corresponding 3D realistic models. We loaded and displayed the 3D realistic model into Cesium globe on the client side and derived the 3D location of a tower by averaging the manually measured 3D positions of its four-tower bases. Cesium allows the client to provide 3D visualization of a real 3D model on the Earth, which is represented by the World Geodetic System 1984 (WGS84) ellipsoid. Therefore, the computed 3D point location of the tower is represented in WGS84 coordinates and referred to as T(lon_T, lat_T, alt_T) (Figure 7a), which contains two parts: (1) geographic coordinates (lon_T, lat_T) that refer to the longitude and latitude values, and (2) alt_T that refers to the geodetic height relative to the ellipsoid surface along the normal direction.
• 3D virtual modeling of a tower supported by photogrammetric data: This step includes three parts: (1) collecting the required size data of a tower based on the LT method, (2) 3D virtual modeling of a tower, and (3) data conversion and storage. First, the LT method that uses the two equations for perspective projection to solve for the 3D object point was used to derive the data size of a tower. To exploit the LT method, we search the images in which a tower occurred. This can be realized by reprojecting the 3D location of a tower to the image space using a known projective matrix of each image given by the AT processing step earlier, according to Equations (1) and (2). At this step, we used the photos captured with coarse-scale photogrammetry as candidate images, and the number of retrieved target images was denoted as p. Further, to improve positioning accuracy, a well-conditioned model should satisfy the constraint of the intersection angle γ, which is the angle between the two adjacent images mapped by the same points. This angle is calculated using Equation (4). The smaller the value of γ, the greater the coordinate error of the object point. The range of γ is set to 30° to 90° in this paper. Note that the 3D location coordinates of a tower used for photogrammetric reprojection should use the SC system. Therefore, a coordinate transformation process is required from WGS84 to the SC system. In this study, the geographic coordinates of the tower location (${lon}_{\mathrm{T}},{ lat}_{\mathrm{T}}$) are converted to the 2D Gauss planar coordinate ($X_{\mathrm{T}}, Y_{\mathrm{T}}$) using the cartographic projections and coordinate transformations library (PROJ), which is a generic library for transferring coordinates among various geospatial reference systems [48]. Regarding the 1D leveling height H_T, we used the known fitting model concerning the case region between geoid and ellipsoidal height for height transformation from ${alt}_{\mathrm{T}}$ to H_T. The converted 3D point location of this tower is represented as ${\stackrel{\mathrm{~}}{X}}_{\mathrm{T}}={\left[X_{\mathrm{T}},Y_{\mathrm{T}},H_{\mathrm{T}}\right]}^{\mathrm{T}}$ to fit the need of Equation (4). After optimal view selection, the actual size data was derived using the LT method. For example, in Figure 7b, we assumed that only two optimal views were retrieved for the tower depicted in Figure 7a. Then, we manually measured the projected image coordinates of the two endpoints (red crosses) of this tower’s side trestle (green lines in Figure 7b) in each optimal image. After that, their corresponding 3D spatial coordinates in the SC system were computed based on the LT method according to Equation (1), using the corresponding homogeneous coordinates of the image and 3D points. This side’s trestle length (i.e., spatial distance labeled with 21.218 m) is further calculated. Others demanded the same way to obtain size data. Second, once the required data size of a tower was known, we used the third-party 3D modeling software to complete the virtual 3D reconstruction of the tower in the actual model scale. This study employs the 3ds Max software, which is commonly used in 3D modeling, and the result was exported in OBJect (OBJ) format, widely used in exchanging 3D software models [49]. Third, the OBJ standard was converted to Graphics Library Transmission Format (glTF). using the obj2gltf tool to load and display the model in the Cesium virtual globe. We chose glTF because it enjoys broad support and free licensing, can be saved as text or binary, and supports the entire scene with assemblies, colors, and materials. It should be noted that the structure type of the towers in the same power corridor may differ, meaning that one power corridor may comprise several tower types. In our study, we discriminated them with a model filename, e.g., i.glTF, where i is an identifier of the ith tower type, ensuring that the same type of tower has a unique identifier. The storage and access methods are the same as that of the realistic 3D model as discussed above.

  $$\gamma =\arccos \frac{\left(t_i-{\stackrel{\mathrm{~}}{X}}_{\mathrm{T}}\right)·\left(t_j-{\stackrel{\mathrm{~}}{X}}_{\mathrm{T}}\right)}{‖t_i-{\stackrel{\mathrm{~}}{X}}_{\mathrm{T}}‖·‖t_j-{\stackrel{\mathrm{~}}{X}}_{\mathrm{T}}‖} (i, j=1, \dots , p) \mathrm{with} i\ne j$$

  (4)
• Determining the orientation parameters of a tower in Cesium virtual globe: We obtained the glTF-based 3D virtual models of the towers. However, when rendered in a Cesium virtual globe, the reconstructed 3D virtual tower model utilizes a relative coordinate system with six-degree-of-freedom (DOF). As we obtained the 3D location of the power T(lon_T, lat_T, alt_T) in step (1), there are three other orientation parameters defined by three angles (i.e., heading, pitch, and roll) relative to the Cartesian coordinate system that is a three-dimensional, earth-centered reference system used in Cesium virtual globe, in which locations are identified by their x (East), y (North), and z (Zenith) values. When the Cartesian coordinate origin was placed at the location as T(lon_T, lat_T, alt_T), the corresponding geometric meanings of the above three angles are depicted in Figure 7a. The pitch and roll angles were set to zero in our study. The heading angle, corresponding to the orientation angle defined in the 2D Gauss planar coordinate system commonly used in surveying and mapping fields, was computed using the line-segment vector formed by the two sequentially selected tower bases (blue dots in Figure 7a) in the Cesium environment. This value should be within 0° to 360°. We have developed a tool contained in “GIS tools” (Figure 3) to compute the value of the heading angle. The detailed processes of using this tool are indicated in Figure 8.
• Collecting the description data of all power towers: Repeat steps (1) to (3) to obtain the needed description data required by Cesium. Note that step (2) can be omitted for the same tower type.
• Measuring the hanging points of a power tower in a Cesium environment: As an indispensable part of the power corridors, power lines are linked to the components of a tower (such as insulators and cross-arms) by hanging points. Therefore, collecting the hanging points of all towers becomes an essential step in the 3D reconstruction of power lines. Since we manually built the virtual power tower models to replace the original ones, if we still collect the hanging points using the LT method, the gap between power lines and towers would inevitably appear due to the collecting errors caused by using a different data source. Especially when the difference in data fineness between levels is large, the model’s discontinuity is evident during the switching of the viewing angles. To overcome this problem, the hanging points should be measured from the 3D virtual models of the towers. However, we cannot directly measure the hanging points from its raw 3D virtual model since it uses a relative object coordinate system. To this end, we developed a “power line plotting” tool on the client side (Figure 3). As depicted in Figure 7c, this tool aims to manually collect the hanging points (yellow crosses) in the Cesium virtual globe after loading and rendering the 3D virtual tower models using the glTF standard. The 3D locations based on the WGS84 coordinates of the hanging points were stored in a GeoJSON format in a predefined directory.

2.5.2. 3D Reconstruction of the Power Corridors

We complete all power corridors’ description data following the above steps. For each power corridor, these description data should comprise all towers’ geo-located 3D positions, heading angles, and 3D virtual models, identified with their filename and hanging points. Based on these data or parameters, we can easily load and render the power towers in the Cesium virtual globe. The power lines between two adjacent towers (i.e., over a span) have the shape of a catenary curve and thus can be modeled by fitting catenary curves to the measured hanging points [50]. As presented in Algorithm 1, Cesium-based rendering and visualization methods were used for the 3D reconstruction of power corridors. Figure 7d shows the partially reconstructed results of one power corridor.

Algorithm 1: 3D reconstruction (rendering) of power corridors in Cesium virtual globe
Input: n_pc: the number of power corridors, n_tower: the number of towers in a power corridor, n_TowerType: The number of tower types, model data: M = {n_type, url}: description data of a power corridor: T = {n_ModelType, (lat, lon, alt), (heading, pitch, roll)};
begin
1:	for i in n_pc
2:		for j in n_tower
3:			for k in n_towertype
4:			// Draw power-towers
5:				if T[j].nModelType = M[k].n_type
6:					Add a 3D model (url[k], T[j].(lon, lat, alt), T[j].(heading, pitch = 0, roll = 0));
7:				end if
8:			end for
9:			// Draw powerlines
10:			Get hanging points of the adjacent towers (T[j].(lon, lat, alt), T[j+1].(lon, lat, alt));
11:			Calculate the catenary curve;
12:			Render powerlines between the two adjacent towers;
13:		end for
14:	end for
end

2.6. GIS Entity Data Collecting and Storage

Although the oblique photogrammetric 3D modeling technique reconstructs a 3D scene, the model comprises a point cloud and triangulated irregular network (TIN) and is represented by the “skin”, which means that all objects (GIS entities) are physically grouped together. This avoids being clicked, queried, and implemented to other basic GIS functions, inhibiting its capacity to meet sophisticated GIS needs. Therefore, generating monomer models from TIN-based 3D models is mandatory.

Some solutions have used the real 3D model to produce monomer 3D models, mainly including slicing, ID, and dynamic rendering-based methods. The slicing-based method directly cuts the target objects on the TIN mesh surface of the 3D real model to derive the corresponding individual models. However, this operation destroyed the Level of Detail (LOD) and model texture, leading to a sharp decline in 3D viewing efficiency and visual quality when we rendered them in the Cesium virtual globe. The ID-based method realized monomer modeling by searching the TIN points belonging to the same object that must be isolated and then assigning the same ID value for these searched TIN points. This method’s difficulty is determining the TIN points of the same object. Unlike the above two methods, the third method adopted in our study manually collects the feature points of an object on the bottom profile. After that, the Cesium-based rendering mode of “fitting to ground” depicts the outer geometry of an object in the form of a “skin”. This method does not destroy the structure of the raw data and completely represents an object, with the monomerized results typically suitable for GIS operations and applications. Figure 9 depicts the monomerized results of the building using the above three methods.

After completing an object’s monomer operation, the attribute information needed by the GIS application should be linked to it. As mentioned above, the dynamic rendering-based monomer method should be conducted on the Cesium Earth. Therefore, to effectively perform monomerizing and establish the corresponding GIS databases, we developed a graphics interactive platform in the Cesium environment by coupling the two sequential steps initiated from the developed functional module of “data management” (see Section 3.2 for details). Figure 10 presents the process of GIS entity data collection and storage.

We take the building rendered pink in Figure 10a as an example. After receiving the real-scene 3D model data in 3D Tiles format from the server, an editing mechanism enables users to add or edit this building in the scene. Based on the pop-up interface, we interactively measure the geometry data depicted by the building corners on the bottom profile and store them as geographic coordinates in the format of longitude and latitude. At the same time, attribute information can be added and linked to this building in the pop-up interface. After that, a data record related to this building was stored in a predefined MySQL database for a house entity and can be reloaded later, if necessary.

Figure 10b presents four data records corresponding to four buildings. Apart from mining working-face data, other GIS entity data collection methods (Table 3) can used, like for the building according to the rules we set. Regarding the mining working-face data corresponding to the polygon entity, its geometry and property data can be acquired from the planned CAD drawings or thematic maps. Such data in China are generally represented in a Gauss planar coordinate system. They should be converted to geographic coordinates and transferred to JSON format to make the data more readable and writable in a web-friendly format. The collected data associated with attribute information of mining-work face are stored in its GIS database. On the other hand, based on the established GIS database, the query function acts similarly to a bidirectional interaction mechanism initiated by clicking an object or a record in the GIS data table. Figure 10b illustrates the queried information linked to the object highlighted in red and displayed in a pop-up panel.

2.7. Dynamic Deformation Prediction and GIS Overlapping Analysis Models

The two main functions were developed related to MIDA, all of which are similar to the actual requirements of the mining engineer/enterprise. One is the PIM, and the other is the GIS spatial overlapping analysis model. The former uses the dynamic prediction method (DPM) of probability integral, which considers the correlation between the process of surface subsidence and the spatial position and time. The DPM is the final state of surface subsidence caused by underground coal mining, namely PIM, multiplied by the time effect function. The PIM is based on the random medium theory and has been extensively used in China. The PIM can be used to illustrate the prediction model of an arbitrary point on the surface, and the formula of unit mining surface subsidence basin is determined according to Equation (5):

$$W_e\left(x\right)=\frac{1}{r}{e}^{\pi \frac{x^2}{r^2}}$$

(5)

Considering the 3D characteristics in Figure 11, the surface subsidence of any point A(x, y) due to the unit mining B(s, t) can be expressed as Equation (6):

$$W_e\left(x, y\right)=\frac{1}{r^2}{e}^{\pi \frac{{\left(x-s\right)}^2+{\left(y-t\right)}^2}{r^2}}$$

(6)

When the mining is O₁CDE, the subsidence prediction formula of point A caused by the whole mining is as follows:

$$W\left(x, y\right)=W_0{\int }_0^l{\int }_0^L\frac{1}{r^2}{e}^{\pi \frac{{\left(x-s\right)}^2+{\left(y-t\right)}^2}{r^2}} \mathrm{d}t\mathrm{d}s$$

(7)

where $W_0=mq\cos \alpha$ is the maximum subsidence at supercritical mining, $r =H_0/\tan \beta$ is the main influence radius, $l=D_3\text{-}{S}_3\text{-}{S}_4$ is the strike length, and $L =\left(D_1\text{-}{S}_1\text{-}{S}_2\right)·\sin \left({\theta }_0+\alpha \right)/\sin {\theta }_0$ is the dip length.

To use this model, it is essential to determine the prediction parameters $q$, $\tan \beta$, $S_1$, $S_2$, $S_3$, $S_4$, and ${\theta }_0$. The sample data for our case study are listed in Table 1. Based on the subsidence formula, aside from the surface subsidence (referred to as w from now on), other deformation values of any point A(x, y) comprising slanting, curvature, horizontal movement, and displacement along the strike direction (90°) and dip direction (0°) for the surface (X, Y) can be deduced, hereinafter referred to as (i90, i0), (k90, k0), (u90, u0) and (s90, s0), respectively. After that, the deformation isolines for each surface deformation can be obtained using a certain interpolation calculation method based on its corresponding scattered data points. Interested readers can refer to [1] for more detailed PIM information. When we accomplish the dynamic deformation prediction calculation represented by deformation isolines, a subsequent GIS overlapping spatial analysis based on clipping operation can be conducted, which requires the intersection between the damage areas depicted by closed deformation contours and the GIS entities, including point, line, or polygon geometry. Based on the intersected results, the new GIS entity data layers can be produced hereafter. It should be noted that the two main functions used Visual Studio 2017 C# and were encapsulated in the server-client as a dynamic link library (DLL) format, which was called when the Cesium web-based framework triggered the MIDA engine.

3. Results

3.1. Results of Realistic 3D Model and Quality Assessment

According to the multi-scalar oblique photogrammetric data acquisition scheme, the total number of captured images for each scale (i.e., coarse-scale, meso-scale, and fine-scale) was about 330,000, 89,000, and 6000, respectively. The 3D coordinates of the GCPs were well-distributed in the study area and measured by a Hi-Target GNSS RTK (Guangzhou, China), which can accept full-band satellite signal, affording a horizontal precision between 1.4 cm and 3.0 cm and vertical accuracy between 3.0 cm and 5.0 cm at each point. To evaluate the real accuracy of the resulting 3D photo-realistic model, some GCPs were selected as CPs, and the ratio of GCPs/CPs is 522/209.

All captured images were processed using the SfM-based digital photogrammetry algorithm implemented by Context Capture (CC) [51], which is suitable for oblique photogrammetric image processing and uses the SfM-MVS technique. The resulting 3D photo-realistic model was exported into an OSGB format. Given that the computing task of AT and the 3D reconstruction is enormous, and the improvement of photogrammetry accuracy and the expansion of the measurement range, the total amount of data is massive. Thus, due to the lower efficiency when using a single workstation, it is necessary to perform parallel and distributed calculations on the photogrammetric processing of AT and 3D reconstruction. The hardware setup in this work involved high-frequency computing nodes (CPU (Central Processing Unit): Core i7-8700K with 3.70 GHz, operating system: Win10) with one NVIDIA GeForce RTX 2080 Ti GPU). Related information about the photogrammetric process is summarized in

Table 4. Related information about the photogrammetric process.

	The Total Number of Images	GSD	The Number of Compute Nodes	Time of Processing
AT	425,000	0.5–4.6 cm	20	94.3 h
3D modeling				220.6 h

.

After the 3D model reconstruction is accomplished, the effect and accuracy of the resulting 3D model should be checked. The precision indexes to evaluate the real accuracy of the 3D model rely on the root mean squared error (RMSE), including horizontal (RMSE_XY) and vertical (RMSE_H) directions calculated by Equations (8) and (9).

$${RMSE}_{\mathrm{XY}}=\sqrt{{\sum }_{i=1}^N\left(∆X_i^2+∆Y_i^2\right)/N}$$

(8)

$${RMSE}_H=\sqrt{{\sum }_{i=1}^N\left(∆H_i^2\right)/N}$$

(9)

where N is the number of CPs, and i represents the ith CP, and $∆X_i$, $∆Y_i$, and $∆H_i$ are the deviations between field survey data and the measured data from a 3D realistic model.

Table 5. Statistical results of RMSEs with different combinations of scales.

No.	Combination of Scales	RMSEXY/cm	RMSEH/cm
1	coarse-scale	8.2	9.1
2	coarse-scale + meso-scale	4.2	4.9
3	coarse-scale + meso-scale + fine-scale	3.9	4.3

reports the RMSEs in the horizontal and vertical positions of the CPs for different scale combinations, highlighting that, compared to the first combination, the second and the third combination significantly improve the geometric accuracy by about 50% in horizontal and vertical directions. In our case study, the accuracies of UAV-based photogrammetric surveys are expected to be 5–10 cm in the horizontal (X, Y) position and 5–15 cm in the vertical (H) position. However, the accuracy values obtained in different scales or combinations were all less than 10 cm, indicating that the results are reliable and usable.

Figure 12 depicts the 3D model output, which is about 260 gigabytes. Aiming at visual quality inspection, Figure 13 illustrates the photogrammetric surveying results with the graphics format of the two ROIs, including AT processing, 3D realistic models, corresponding local models, and details by merging different photogrammetric scales. By comparing the local 3D models of the same object in the first ROI displayed in Figure 13c (only with coarse-scale) and f (with the combination of coarse- and meso-scale) corresponding to the two local parts labeled with red rectangles in Figure 13b,e, respectively, we observe that multi-scalar oblique photogrammetry can better restore the actual shape for the complex object of the industrial equipment. Similarly, for the second ROI, when compared to the local 3D model of a tall building (~60 m) by fusing coarse- and meso-scale as depicted in Figure 13i, the corresponding result with all three scales indicates that the holes are repaired, the wall is flatter, the corners and edges are more clear, and the housing sector can be accurately collected, all of which can be comparatively observed from their corresponding local details labeled with a red ellipse. Figure 13 and Table 5 reveal that the hybrid photogrammetry scheme effectively balances the geometric precision and visual quality of the produced 3D model.

3.2. Practical Application and Spatial Accuracy Evaluation of the Prototype System

A prototype system was constructed and implemented, serving as a MIDA application testbed. The system adopts the browser/server architecture, i.e., the server side was developed with Visual Studio C#, the client side was developed with Vue.js, and the 3D engine on the client side is Cesium. At present, the developed webGIS system was only deployed and accessed in a local area network (LAN) concerning data security. The platform was designed as an account-based application. Thus, user registration was required to access the full functions of the system. After logging into the website, certain registered users were authorized to perform the system’s function. On the other hand, a well-designed interface could resolve some of the complexities of the proposed system, making it more useable and, thereby, widely accepted and adopted by the users. Therefore, this study focused on designing interactive and user-friendly GUIs that allow experts and non-experts to operate the system easily. Different modules structure the GUIs, making the system easy to use and expandable. Figure 14 presents the main GUIs of the web-based MIDA platform. In particular, Figure 15 illustrates the enlarged local 3D real model of an industrial square provided by a platform with a bird-view mode.

The GUIs comprised a map viewer, various GIS entities and analysis results visualization, and GIS analysis tools, mainly including some basic GIS tools, a coordinate conversion tool, and the developed MIDA-related functional modules. The map viewer was a Cesium-based 3D virtual globe that provided a WGS84 spatial reference and enabled intuitive geo-visualization of the analysis results and multi-source geospatial data on a default base map (e.g., satellite images) obtained from the OGC’s WMS. The developed functional components were docked at the upper-left of the map viewer. The model datasets, including the 3D photo-realistic model, 3D virtual model of power corridors, and other base maps, were contained in the “Data Layers” component and used as the background data layers. Some basic GIS tools, including point, distance, area, and orientation angle measurement, could facilitate the manifold and practical applications. It should be noted that the tool of coordinate conversion frequently used in the system was also developed on the client side to realize bidirectional conversion between WGS84 and the SC system.

The two important functional modules, namely “Data Management” and “MIDA”, were user interfaces where users could interact with the system. The “Data Management” module collects the required GIS entity data and stores it in the MySQL database described in Section 2.6. After all the required GIS datasets were collected and saved, in the pop-up window named “Data Management”, shown in the left part of the main interface, we can review these data in a data grid table and query them by a predefined filter condition. Also, we can further edit them (modify or delete) to ensure the reliability of the collected GIS entity data.

After completing the data management operation, the next operation of MIDA can be performed. The steps are as follows.

• Selecting a deformation type: By pushing the “MIDA” menu, a pop-up sub-interface named “MIDA” would be displayed in the central part of the main interface. In this subwindow, we should first select a deformation type from the combo list, which consists of different deformation types, as mentioned in Section 2.7.
• Computing the deformation contours: Once a deformation type is confirmed, the system can post a PIM execution HTTP request to the server. When the server receives this client-side request, the MIDA engine is triggered to perform surface deformation calculations. After that, the corresponding deformation isolines will be computed on the serving side and returned to the client side to be rendered and visualized on the background layers. These deformation isolines represent different impact areas in the shape of 2D polygon contours.
• Selecting the interested GIS entity data and performing overlapping analysis: The GIS vector data are organized with down-left data layer lists in the “MIDA” sub-interface according to entity types (i.e., point, line, and polygon), and for each entity type, it was also classified into several subtypes by a set of expandable nodes. We can choose the interested data type in this data layers list for GIS analysis. Subsequently, we can manually click and select a deformation isoline to trigger the GIS overlapping analysis model engine for performing deformation analysis, i.e., we used the selected 2D deformation contour and the interested data layers chosen above as a background polygon and target geometries, respectively, to perform an intersecting calculation for interested data layers. The computational results are pushed to the client-side application for visual rendering in real time. After this process, the analysis results will be produced and stored in the GIS database. At the same time, the selected interested data layers list will be updated and displayed in the down-right data grid table of the “MIDA” sub-interface.

Take two deformation types of subsidence and horizontal movement as examples to illustrate the results of the MIDA presented in Figure 16, in which the deformation contours and the selected isolines used as background polygons for intersection are rendered with green and red colors, respectively. Note that for clarity, the sparse representative sample data of the above two types of deformation isolines were extracted and used for GIS analysis. Figure 16a,b illustrate the damage analysis results using subsidence isoline (w) with a value of 10 mm and horizontal movement along the strike direction (u90) with a value of −10 mm, respectively. After the MIDA, only the intersected results of different types of objects inside the impacted region are rendered and displayed, and one can also click to select an arbitrary object inside this region to retrieve its attribute information from the results database. Once an object with a certain entity type is chosen, it will be highlighted with red, and its corresponding attribute information panel instantly appears. Figure 16a,b illustrate this process for a road (red line) and farmland (red polygon region).

We tested the spatial accuracy of the proposed 3D webGIS system with two types of validation data. The first type is the geo-referenced ground points with distinguished features such as road or house corners and GCPs; the other is the spatial overlapping results data. The 3D locations of the former were measured in the raw, realistic 3D model and represented in the SC system by CC Acute Viewer software (Version 10.20) [51]; these were used as the true values and denoted as $X_{\mathrm{rel}}$,$Y_{\mathrm{rel}}$, and $H_{\mathrm{rel}}$, respectively. Their corresponding coordinates measured in the web platform are denoted as $X_{\mathrm{web}}$,$Y_{\mathrm{web}}$ and$H_{\mathrm{web}}$, respectively. This type of validation data is mainly used to assess the absolute accuracy of the system. As for the latter, overlapping analysis results of two kinds of entities represented by line and polygon geometry were used to evaluate the relative accuracy. In particular, only the clipped parts of lines and polygons inside the background polygon are considered. The main reason is that it can be used to verify the relative precision, correctness, and reliability of the GIS overlapping analysis algorithm developed by the authors. The true values (length and area) were obtained by the spatial overlapping analysis tool in ArcGIS software (Version 10.2) and denoted as D_Arc (length) and S_Arc (area), respectively. The corresponding measured results automatically computed and returned from the system are denoted as D_web and S_web, respectively. The absolute differences ($∆X$, $∆Y$, ΔH) and relative differences ($∆D$, $∆S$) were calculated between the true and observed values. The absolute precision utilizes the RMSE_S as indexes, and the calculation method can refer to Equations (8) and (9). Regarding the relative accuracy of length and area, a deviation ratio (DR) is first defined to evaluate the length and area errors. It can be calculated by ΔD/D_Arc and ΔS/S_Arc, respectively, after the average DR (ADR) is computed and used as a relative accuracy index. For both types of validation data, 31 pieces of data were collected, and the accuracy results are tabulated in

Table 6. Assessment results of absolute accuracies of the web-based application system.

No.	ΔX (cm)	ΔY (cm)	ΔH (cm)
1	1.9	7.5	8.3
2	3.1	0.8	6.9
3	6.6	0.9	10.2
4	1.6	0.5	4.8
5	0.7	5.8	7.9
6	1.2	1.6	2.6
7	2.2	0.8	8.5
8	1.8	1.2	7.9
9	3.2	1.3	10.5
10	0.9	1.2	8.3
11	2.2	1.5	6.9
12	2.7	1.4	7.2
13	1.6	0.8	5.1
…	…	…
31	2.6	1.4	2.5
RMSEs	2.3		7.6

and

Table 7. Assessment results of relative accuracies of the web-based application system.

No.	DArc (m)	Dweb (m)	ΔD (cm)	SArc (m2)	Sweb (m2)	ΔS (m2)
1	907.832	907.811	2.1	21,440.744	21,443.263	2.519
2	953.583	953.552	3.1	53,470.971	53,470.645	0.326
3	359.964	359.963	0.1	5218.522	5218.501	0.021
4	1397.391	1397.442	5.1	101,955.341	101,949.511	5.830
5	1110.600	1110.432	16.8	75,583.731	75,583.951	0.220
6	1052.295	1052.341	4.6	68,907.812	68,907.633	0.179
7	779.522	779.475	4.7	29,032.383	29,029.294	3.089
8	574.484	574.522	3.8	2378.262	2378.951	0.689
9	527.241	527.191	5.0	8265.741	8264.173	1.568
10	1073.310	1073.382	7.2	43,165.614	43,165.503	0.111
11	644.373	644.373	0.0	24,570.173	24,570.112	0.061
12	1059.894	1060.001	10.7	67,984.801	67,986.371	1.570
13	414.961	414.923	3.8	9132.491	9131.052	1.439
…	…	…	…	…	…	…
31	311.771	311.825	5.4	5851.232	5852.943	1.711
ADRs	0.0059%				0.0074%

; Table 6 indicates that the absolute accuracies of plane and height coordinates represented by RMSEs are 2.3 cm and 7.6 cm, respectively. As for the relative accuracies defined with ADR, Table 7 demonstrated that the accuracies of both length and area are very high; the corresponding values of ADR are 0.0059% and 0.0074%, respectively. These results exhibit much smaller deviations between the true and the observed data and demonstrate the reliability of the GIS overlapping analysis algorithm and the web platform and can meet the practical engineering requirements.

4. Discussion

Owing to the important social effect of mining activities, there is an urgent need to obtain early alarming risk information. The main contributions of this paper are to propose a unified multi-scalar oblique photogrammetry-supported web-based framework including 3D realistic and virtual modeling, GIS vector data acquisition and storing, and the corresponding web-based prototype system development. The proposed framework allows non-professionals to perform complicated MIDA tasks and use a 3D webGIS in the mining industry to a large extent. The proposed framework is still worth discussing as follows.

Compared to traditional fieldwork, UAVs are suitable for surveying large areas of the Earth’s surface because they generate low-cost oblique photogrammetry, contribute to constructing realistic 3D models, and allow for recording and measuring the areas of interest. The resulting data products can be used as an important data source for GIS applications. Despite this, in practice, there is always a conflict between model quality and flight efficiency. Aiming at this problem, a technical pathway of multi-scalar photogrammetry was proposed by fully considering the shooting performance of different UAVs to balance the data acquisition efficiency and model quality, which is particularly essential for the large-scale regions examined in this work. The visual inspection and precision assessment results verify that multi-scalar photogrammetric 3D modeling ensures the model’s effectiveness and achieves centimeter-level accuracy, which meets practical engineering needs.

Taking the oblique photogrammetric data as an important data source, a technical workflow for GIS dataset production and a storage method were proposed. Regarding modeling 3D data, aside from the 3D realistic model, the 3D virtual model of power corridors distributed in the mining area was constructed using third-party 3D modeling software based on the actual size derived by the photogrammetric LT method. The above two 3D model datasets in our 3D webGIS system are the background data layers to enhance the 3D virtual environment. Concerning GIS vector data, many involved GIS data are collected from the 3D realistic model in the Cesium environment and managed by the GIS database. In our case study, the main consideration of selecting photogrammetric data as the main data source is to reduce the computational burden and complexity of multi-source data processing usually adopted by existing studies.

Moreover, a web-based graphics platform was developed to visualize geographic data and web-based MIDA tasks. The proposed framework reduces the mining industry personnel’s knowledge requirements for mining deformation. The users do not need to pay special attention to the details of MIDA. Additionally, each key process in the framework has its data interface operation commands, which reduces the difficulty for mining designers to master and understand expert knowledge quickly. On the other hand, during the implementation of this prototype, open-source technologies and tools were used to provide a more open, transparent, and flexible solution. Multiple interactive data visualization servers and terrain servers use WebGL to store, distribute, and display 3D geospatial data more conveniently in a web browser environment without any plugin. A practical application test in the coal mine is carried out in the case study. It demonstrates that each process within the framework can be manipulated and performed seamlessly and smoothly. The spatial accuracy results of the web-based system by two types of validation data further indicated that the application accuracies have a good agreement with the true values.

5. Conclusions

A preventive MIDA application is of significant importance to the mining industry. However, there are still gaps in the current research and practical applications. The major challenges include the difficulty of obtaining the precise 3D GIS data, the lack of an intuitive visualization environment, and the lack of a unified MIDA framework. To fill these gaps, an oblique photogrammetry-supported unified workflow for MIDA applications is proposed in this paper. In this framework, a multi-scalar oblique photogrammetric scheme is first designed and used to collect and process essential data for a 3D webGIS. The visualization module realizes the 3D visualization of GIS entity data and MIDA results based on WebGL. The proposed framework is applied to a large-scale underground coal mine area, and the results prove the feasibility and reliability of the system. Despite these contributions, some limitations need to be addressed. First, the important step of GIS data collecting and storing in our framework is carried out manually. This process should be further optimized to reduce human involvement. Furthermore, although the GIS entities in our case study present 3D shapes by dynamic rendering (e.g., houses) or manual 3D modeling (e.g., power towers), the corresponding geometric components are dependent on each other. Future research can be carried out by integrating BIM and webGIS environments, as GIS describes the large-scale macro environment. In contrast, BIM describes the micro-information of object geometry and semantic attributes inside an object, complementing each other. This feature is important because when an interesting object is destroyed in a mining area, we can repair or reconstruct it according to the recorded documentation.