1. Introduction
In geological exploration, seismic tomography imaging is extensively employed due to its high accuracy and resolution capabilities [
1,
2]. The propagation velocity of seismic waves is influenced by the physical properties of rock, such as density and elastic modulus, as well as the structural features of rock masses. This technique leverages the propagation characteristics of seismic waves through geological strata to analyze and infer the physical properties and structural attributes of formations, including karst terrains, faults, and fractured rock masses [
3,
4,
5,
6].
However, the results of seismic wave CT are typically scattered and consist of independent two-dimensional profile data, making it challenging to analyze data correlation and connectivity across profiles. Consequently, it fails to intuitively convey the distribution, direction, spatial location, and morphology of underground anomalies, thus lacking in intuitive comprehension, particularly for non-geologists. This limitation complicates the establishment of three-dimensional spatial spreading scenarios, adversely affecting the visualization of detection results [
7]. Understanding the anomalous distribution, spatial positioning, and morphological features of three-dimensional strata structures based solely on 2D profiles is a complex process.
To address these limitations, researchers have explored three-dimensional visualization techniques for seismic tomography CT. Li Zuneng [
8] employed Voxler 3D visualization software to establish a three-dimensional numerical model of seismic wave velocities. This model was processed through slicing, volume rendering, and surface rendering to analyze and delineate seismic wave velocity ranges of various lithological bodies, stratigraphic layers, and boundaries of karst development zones. Hu Junjie et al. [
9] utilized cross-hole electromagnetic CT methods to detect karst development in bridge foundation areas, demonstrating the regional distribution of karst development zones through three-dimensional modeling. Yu Bo [
10] applied quasi-three-dimensional cross-hole seismic CT detection technology and Voxler software for three-dimensional data gridding, achieving visualization of geological structures, bedrock surfaces, and the three-dimensional spatial morphology of karst in the survey area. Liang Juan [
11], using Blender software and Python programming language, rapidly constructed dynamic three-dimensional visualization models of seismic wave propagation, dynamically simulating and analyzing wave effects using the aforementioned methods.
There have been significant international efforts in the development of three-dimensional geological visualization software. For instance, the Vulcan software developed by MAPTEK in Australia [
12] visualizes three-dimensional geological models with exploration data and provides modules for road design and quarry modeling, which are crucial for open-pit mining planning. Similarly, the Mining Visualization System (MvS) series developed by CTECH in the United States [
12] enables visualization of geological models and related spatial analyses. This software supports true three-dimensional geological modeling, statistical analysis, structural analysis, and geological engineering calculations.
With advancements in computer hardware and software, three-dimensional visualization techniques for seismic tomography CT have matured, offering richer display effects [
13]. However, the current predominant approach involves presenting and analyzing geological data within specialized visualization software, which requires complex operational steps and professional technical support. These software packages often necessitate high-performance computers and dedicated graphics processors, limiting their accessibility and widespread adoption [
14]. In contrast, significant advancements have been made in WebGL-based visualization for medical CT, demonstrating that user-friendly web interfaces can provide precise and efficient algorithmic services [
15,
16,
17]. WebGL-based visualization addresses the limitations of poor mobility performance on PCs and offers significant advantages in data transmission speed and display smoothness. In the field of underground engineering, researchers have also explored WebGL-based visualization of stratigraphic models. For example, Xinshan Xu [
18] developed a convex hull fit for the body side of stratigraphy and TINs for top and bottom surfaces, achieving visualization of geological body models. Despite these advances, the application of WebGL technology for the visualization of exploration results on the web remains underdeveloped.
This paper proposes a three-dimensional visualization method for seismic tomography CT results based on WebGL technology, which is independent of specific software and offers high visualization efficiency. This method utilizes wave velocity reference points to create grids and applies vertex and surface coloring to the model. By processing model profiles with color adjustments based on CT detection data, it achieves model coloration and visualization display. The visualization results are showcased on an intelligent construction platform with response times controlled within 2 s. The proposed visualization method was validated and applied to underground chamber projects at a hydropower station. Comparative validation with geological three-dimensional design software, independently developed by the China Northwest Institute of Civil Engineering and Architecture, demonstrated the effectiveness and practicality of the method. The WebGL-based visualization results can reveal weak geological layers, display the distribution of surrounding rock types in underground structures, and provide intuitive references for support schemes of underground structures, thus promoting the digital management process during engineering construction periods.
2. Research Method
There are many methods for processing discrete data by difference, and the most widely used methods at present include Kriging interpolation, inverse square distance interpolation, and DSI interpolation.
The Kriging method is a geostatistical approach used for spatial data interpolation. It is widely employed in geological exploration, environmental science, and geographic information systems [
19]. This method estimates the values at unknown locations based on known discrete point data to achieve numerical predictions for the entire region. The advantage of the Kriging interpolation method is that it considers the correlation between spatial data, providing accurate estimates for unknown points and assessing the confidence of the interpolation results. However, its drawbacks include the need for extensive computations on large-scale data and the complexity of parameter tuning and model selection [
20,
21].
The IDW method is a commonly used interpolation technique that calculates the value of an unknown point based on the weight of the value of known points, with the weight being inversely proportional to the distance from the unknown point [
22,
23]. The advantage of this method is its simplicity and ease of use, allowing it to reasonably weight neighboring points according to their distance, so that points closer to the unknown point contribute more to the estimated value. However, it also has disadvantages, such as the potential for uneven distribution of dense and sparse points, which may lead to discontinuous interpolation results. Therefore, in practical application, it needs to be balanced and adjusted according to the specific situation [
24].
The DSI interpolation method was first proposed by Professor Mallet of the University of Nancy in France [
25]. The basic idea of DSI interpolation is to establish an interconnected network among discrete data points. If the known node values on the network satisfy certain constraint conditions, the unknown node values can be obtained by solving linear equations. The DSI difference method only relies on the topological relationships between nodes and is an interpolation algorithm that is not limited by dimension. The interpolation calculation of unknown points only depends on the topological grid, not on the plane or spatial coordinates, and the calculation process is generally iterative, which is conducive to dynamically approaching sampling points and is suitable for interpolation problems of multi-valued surfaces. The DSI interpolation method can consistently maintain the geometric detail features of the optimized grid and the original grid, impose discrete control point constraints, ensure consistency with the original grid during optimization, and can greatly improve the grid quality, thereby significantly improving the data visualization effect. Since the DSI interpolation method has high interpolation accuracy and calculation efficiency, it is widely used in geological modeling and earth physical data processing.
WebGL (Web Graphics Library) is a JavaScript API that presents interactive 3D and 2D graphics in web browsers. WebGL can achieve efficient real-time graphic rendering and is more suitable for large-scale data visualization than traditional 2D graphic technology. It can provide smoother rendering effects and faster response speeds and can enable graphic viewing at any time on the web end.
In this paper, the DSI interpolation method is used to perform differential processing on the discrete seismic wave CT result data. At the same time, in order to improve the visualization degree and convenience of browsing, the seismic wave CT three-dimensional results are realized in the web-based platform through WebGL technology, and the three-dimensional visualization display is achieved.
It is important to note that although the DSI interpolation method performs well in handling noisy or incomplete seismic data, its effectiveness is also influenced by factors such as data quality, the extent of missing data, and interpolation parameters. Therefore, when processing noisy seismic data with the DSI interpolation method, it is often combined with other denoising techniques, such as filtering methods, to identify and filter out noise in the seismic data. This approach utilizes the selectivity of effective signals and noise in specific domains to eliminate the noise. When dealing with incomplete seismic data, the DSI interpolation method takes into account the spatial correlation of the seismic data, using the spatial distribution and correlation of known data points to predict the values of unknown data points, thereby filling in the missing data.
3. Visualization Process
3.1. Creation of Wave Velocity Reference Point Grid
The color grid defined in this paper consists of multiple vertices and the corresponding wave velocity values. By setting the bounds, the wave velocity at each point in the model can be quickly calculated. Finally, an appropriate interpolation method is chosen to determine the wave velocity at each point within the mesh.
3.1.1. Determine the Bounds
The bounds should encompass the entire model that needs to be rendered. The bounds form the smallest hexahedron with edges parallel to the coordinate axes and vertices with three-dimensional integer coordinate values. The minimum values of the BIM model coordinate points in the X, Y, and Z directions are
and
. The maximum values are
and
, respectively. The bounds can be represented as follows.
The coordinates
x,
y, and
z of all points at the boundary and interior of the bounding box satisfy the following equation.
3.1.2. Grid Distribution Density
In geological 3D modeling systems, grid technology requires that each unit can store information and perform data interpolation operations. Compared to finite element models, the number of grids used for geological analysis services is 2–3 orders of magnitude higher. Generally, the number of mesh elements in a three-dimensional finite element model is within one million, with some models reaching several million elements. However, the quality grading and parameter values of geological rock masses often involve grids numbering in the millions, typically tens of millions, and in some cases, even billions.
Generate a point set Ω based on the set step size λ and store it in an array. The size of
is given by
. Any point in the point set can be represented by
and can be expressed as the following equation.
where
Considering computational efficiency, the step size λ is primarily determined by the complexity of the actual scene and model. In practical applications, it is recommended to generate fewer than 500,000 grid reference points to ensure real-time computational efficiency.
3.1.3. Interpolation and Processing of Wave Velocity of Mesh
Table 1 shows a sample of the detection data. The set of all points from the probe data is denoted by Ψ. The wave speeds of all points in Ω need to be obtained through interpolation based on the wave speed values of all points in Ψ.
DSI uses discrete sampling points to estimate values between these points through a smooth interpolation method, generating smooth surfaces or curves. This feature makes DSI suitable for applications in image processing, computer vision, and 3D reconstruction. The primary goal of the DSI algorithm is to estimate values between a set of discrete samples by employing a smooth interpolation method to produce smooth surfaces or curves. DSI establishes an objective function R*(ϕ) = R(ϕ) + ρ(ϕ) for calculating the optimal solution on mesh nodes, where R(ϕ) is the global roughness function and ρ(ϕ) is the linear constraint violation function. By minimizing this objective function, two main objectives are achieved. First, minimizing R*(ϕ) ensures that the value of the function at any given node approximates the average of the values at the surrounding nodes. Second, minimizing ρ(ϕ) ensures that the raw sampled data are converted into linear constraints on the nodes, making the node values as close as possible to the sampled data by maximizing compliance with the linear constraints [
26,
27].
Table 1 provides an example of the detection data. The DSI method is represented by G. For any point p∈Ω, the computation of its velocity v can be expressed by the following equation.
3.2. Color Rendering of Model Surface Based on WebGL Technology
Model surface coloring requires determining the corresponding color for each point in the model based on the calculated wave speed values, followed by the rendering process. The rendering described in this paper is based on WebGL technology. WebGL allows developers to directly use the GPU to render 3D graphics within a browser, incorporating rich visual effects into web design. Color rendering through WebGL relies on shaders [
28,
29,
30]. Shaders enable customization of the graphics card rendering algorithm, allowing the image to achieve the desired visual effect. There are two types of shaders: vertex shaders and fragment shaders, which handle data at the vertex and pixel levels, respectively.
3.2.1. Preprocessing of Vertex Colors
Colors are uniformly encoded using the RGB color mode, which employs three variables to represent the intensity of red, green, and blue [
31,
32]. This mode provides a unique numerical code for each color. Subsequently, color bands are customized to preprocess the colors of points with different wave speeds. The color band consists of the following components.
① minNum (minimum value).
② maxNum (maximum value).
③ count (5 ≤ count ≤ 30, Color quantity).
④ minColor (starting color RGB, minimum value).
⑤ maxColor (termination color RGB, maximum value).
The set of color bands is denoted by color
. Arranging the colors in RGB from small to large, each color band
and the corresponding wave velocity value
can be expressed, respectively.
3.2.2. Drawing of Model Surface Colors
In this paper, the surface of the model is colored using shaders. During the graphics rendering process, the graphics rendering pipeline visualizes the vertex data of the model. When the graphics card processes the triangle drawing task, it generates additional pixel fragments during the rasterization stage to ensure accurate rendering of the triangle [
33,
34,
35].
The positions of the fragments are determined based on their relative positions within the triangle. The graphics card interpolates the input variables of the fragment shader using this location information. Subsequently, the fragment shader automatically interpolates the color of each pixel within the triangle. To achieve the effect of seismic wave CT contour cloud image rendering, the color of each pixel must be accurately calculated in the fragment shader. This allows the graphics card to present a colored contour cloud map on the surface or section of the model after rendering is complete, providing an intuitive visual presentation for scientific research and data analysis.
The specific operation steps are as follows.
(1) Coordinate information input: input model vertices, projection matrix, observation matrix, model matrix , and the texture buffer containing reference point set data into the shader.
(2) Information transmission: transfer the vertex coordinates of the model to the tile shader in the vertex shader.
(3) Coordinate calculation: calculate the world coordinates of the current point in the tile shader.
(4) Reference point query: obtain the coordinates of the 8 closest reference points to the current point based on the world coordinate
, and calculate the wave velocity values of these 8 reference points in the texture buffer. The position relationships of these points are illustrated in
Figure 1. Assuming that the world coordinates of the current point are
, the relevant points can be calculated using the step size λ. The eight reference points are as follows.
where,
.
(5) Calculation of wave velocity values: perform trilinear interpolation on the wave velocity values of the 8 reference points to obtain the wave velocity v of vertex p0.
(6) Visualization result: obtain the final color of the pixel based on the wave velocity value and color band division. This 3D model can display wave velocity information through surface color.
By following these steps, the graphics card can accurately render the model surface with colors that represent the seismic wave velocities, providing a clear visual representation for analysis.
3.3. Color Determination of Model Profiles
To improve the update speed of real-time display results in this study, the pixel color of points outside the reference point set is set to semi-transparent white, and the calculation and processing of profile contours are omitted. The steps for drawing a profile are as follows:
(1) Generate the four vertices of the profile based on the location information that the user intends to view.
(2) Pass the profile vertices, projection matrix, observation matrix, model matrix , and texture buffer containing reference point set data into the shader.
(3) Transfer the vertex coordinates of the profile to the tile shader in the vertex shader.
(4) Calculate the world coordinates of the current point in the tile shader using the formula.
(5) If the coordinate point lies within the range of the reference point set, retrieve the coordinates of the eight closest reference points based on the world coordinate Vworld, and obtain the wave velocity values of these eight reference points from the texture buffer. If it is outside the range of the reference point set, set the pixel color of the point to semi-transparent white.
(6) Calculate the wave velocity value of pixels within the reference point set range.
(7) Using the wave velocity values and color band division results, generate a profile that represents the wave velocity information using surface colors.
5. Applications in Engineering
5.1. Project Introduction
The hydropower station features a concrete face rockfill dam, an underground water diversion and power generation system on the left bank, and a spillway and flood discharge tunnel on the right bank. The project includes multiple underground caverns, with diversion and flood discharge tunnels situated on the right bank and water diversion and power generation systems located on the left bank. Key structures include the underground powerhouse, main transformer cavern, and tailrace gate operation room, which are arranged in parallel with their axes oriented directly north–south, forming an angle of approximately 30° to 40° with the rock strata. The underground powerhouse has dimensions of 183.5 m × 25.3 m × 64.25 m, the main transformer cavern measures 117.2 m × 16.5 m × 30.8 m, and the tailrace gate operation room is 115.6 m × 10.2 m × 59.1 m. The cavern spans are 25.3 m, 16.5 m, and 10.2 m, respectively, with distances of 39.7 m, 19.8 m, and 19.8 m between these structures.
The scale of the underground cavern group is substantial, with the largest span being 25.3 m for the powerhouse. The caverns are interspersed with thin rock layers between some, and the geological conditions are complex. Preliminary surveys indicate that the rock formations primarily consist of T3Z2 (3) to T3Z2 (5). The main powerhouse, main transformer cavern, and tailrace gate operation room are surrounded by type IV2 rock. Several unfavorable geological features, such as the steeply dipping f31 fault and PH strips, pose challenges to the stability of the surrounding rock.
By combining seismic wave CT visualization methods with investigation data and underground engineering excavation, we achieve a visual display of geological properties and surrounding rock classifications.
5.2. Comparison and Validation
Utilizing Web-GL technology, the physical exploration results of this underground cavern are visualized and compared with the visualization results obtained from geological 3D design software (V1.0). This geological 3D design software, developed independently by the Northwest Institute of China Electric Construction, employs the DSI interpolation algorithm as its core modeling substrate. It integrates a geological database and uses survey data or manual judgment as modeling constraints, achieving step-by-step approximation and fitting to ensure geological analysis and rationality. This software has been widely promoted and applied in the survey and design of many hydropower station projects, such as Jinchuan and Mar block. It boasts extensive practical experience in the professional application of rock mass grading, geotechnical engineering analysis, and design within the 3D geological model.
Figure 3 presents the visualization results generated by the geological 3D design software, while
Figure 4 displays the visualization results based on Web-GL technology (Version: 1.0.1). The two visualizations are consistent with each other, and the CT wave velocity calculations yield identical results. The 3D visualization method of seismic wave CT results, as implemented on the Web terminal in this paper, can accurately represent the physical exploration findings. This method offers a convenient display technique for the digital management of construction projects.
5.3. Results and Analyses
Based on the intelligent construction platform described in
Section 3, the geophysical visualization results of the underground cavity group can be obtained.
Figure 5 illustrates the enclosing mesh body that completely contains the underground caverns.
Figure 6 presents the attribute data and rock mass quality classification data within the enclosure mesh. The blue area denotes the weak zone, which requires special attention during the excavation and support process of the underground powerhouse cavern group, while the red area represents areas with good rock mass quality.
By assigning the geophysical wave velocity data on the cubic network to the surface of the underground engineering structure model, the seismic wave velocity distribution of the underground powerhouse cavern group is depicted in
Figure 7. The non-uniformity of wave velocity distribution reflects the differences in the physical properties of underground rocks, providing an important basis for evaluating the quality of surrounding rocks.
Figure 8 shows the result of assigning the surrounding rock classification data from
Figure 6 to the underground engineering model. The distribution of surrounding rock quality is uneven, with multiple weak zones and unstable areas. Corresponding reinforcement measures need to be implemented during the design and construction process in these areas to ensure the stability and safety of the underground cavity group.
The web-based 3D visualization technology of seismic waves intuitively displays the geological conditions and surrounding rock mass distribution, providing robust support for engineering design and construction. Additionally, web-based visualization offers high efficiency and significantly enhances digital management during the construction period of engineering projects. The 3D visualization technology of seismic waves on the web can intuitively display the geological conditions and rock mass distribution of underground powerhouse caverns, providing essential support for engineering design and construction. Additionally, web-based visualization updates are highly efficient, promoting effective digital management throughout the construction period of engineering projects.
The three-dimensional visualization method of seismic wave CT results has broad scalability. This method relies on advanced algorithms, and as these algorithms are continuously optimized and upgraded, the visualization effects and computational efficiency can be improved. With the ongoing enhancement of computer hardware performance, such as the acceleration capabilities of GPUs and increased memory capacity, the seismic wave CT three-dimensional visualization method can better utilize these hardware resources to achieve higher resolutions and visualize more complex scenes. This method is not only applicable to the field of seismic exploration but can also be extended to other areas such as geological exploration, geotechnical engineering, and environmental engineering. For instance, when conducting engineering construction in karst-prone areas, it is often necessary to grout karst caves. The precise acquisition of the location and three-dimensional morphology of underground karst formations is crucial for assessing the degree of karst development and calculating the volume of grouting required. Through three-dimensional visualization technology of seismic wave CT results, the morphology of karst caves can be finely detected, enabling more accurate calculations of grouting volume and effectively reducing material waste in construction projects. Mineral reserves are an important basis for the evaluation of mineral resources, serving as a critical reference for assessing the industrial significance of ore deposits, determining the production scale of mining enterprises, investment scale, and service life. The three-dimensional visualization technology of seismic wave CT results enables accurate estimation of underground mineral reserves, which is significant for the exploration and production design of mines, production planning, reserve management, and production economic management. In addition, this technology can also be applied to the fine detection of adverse geological bodies underground and the delineation of groundwater bodies and lens profiles.
6. Conclusions
Aiming at the problems that the post-processing of CT detection results in geological specialties is restricted by specific software, low processing efficiency, cumbersome processing procedures, and is unable to be displayed conveniently, this paper proposes a three-dimensional visualization method of seismic wave CT results based on WebGL technology. The specific conclusion is as follows.
1. Firstly, generate an integer enclosing box based on the target 3D model, and calculate the wave velocity values of all reference points using the DSI interpolation calculation method. Secondly, calculate the wave velocity values of all the points to be rendered using the wave velocity values of the reference points. Finally, determine the color values of the points to be rendered and render the model and profile in real-time using shaders.
2. The visualization results of the research method in this paper were displayed on the intelligent building platform with a response time controlled within 2 s. The comparative validation of the visualization results between the proposed method and the geological 3D design software (independently developed by the Northwest China Institute of Civil Engineering and Architecture) proves the effectiveness, efficiency, and practicality of the proposed method.
3. The method was used to represent the geological information and classification of the enclosing rocks of an underground cavern, and the visualization results provide an intuitive demonstration of weak strata and unstable areas. The results present a visual reference for the support scheme of the underground structure and facilitate the digital management process during the construction of the project.