Increasing Resolution of Seismic Hazard Mapping on the Example of the North of Middle Russian Highland

Abstract

The standard problem of engineering geophysics, solved for road and house building and other construction types, is in the localization of areas with increased mobility in the upper part of a geological cross-section and in the parameterization of this mobility in terms of seismic intensity. There is a standard approach, according to which researchers assess the elastic strength characteristics of the core to a depth of about 30 m, implement the accumulation of seismogram observations, simulate accelerograms for particular conditions and, taking into account the data of complex geophysical methods, calculate the increment of seismic intensity as one of the parameters of a seismic hazard. The final result of this approach has the form of a seismogenic hazard map and a set of recommendations including the consideration of identified areas with a significant increasing seismic intensity increment, due to the peculiarities of the geological structure of polygons. This result is reliable, but very expensive, and requires the development of primary estimations of the rock massif with reduced resistance to external loads, which would optimize the efforts in engineering drilling and in field geophysical measurements in order to densify their spatial grid in the vicinity of a priori known positions with an increased seismogenic hazard. In addition, relatively sparse grids of wells, as well as local geophysical profiles laid under conditions of a complicated landscape, do not accurately localize risky areas in order to focus the attention of builders on strengthening the specific part of raised constructions. Following the wishes of our customers and relying on long-term testing of our interpretational developments, we formed an approach to primary hazard forecasting based on remote sensing data and digital elevation models, which can be classified as data with relatively free access. This article presents the results of research which was based on these free-of-charge data and which was developed in the field of construction of ground engineering structures for agricultural purposes, where one of the factors of mobility in the upper part of a cross-section is intensive karstification. Basically, the construction area according to the general seismic zoning maps is seismologically passive, though the relatively fast dynamics of karst determines the relevance of the detailed seismic zoning. The results of our interpretations are verified by deep geological and structural reconstructions based on wave analogies. The representativeness of the final forecast was confirmed by subsequent seismic assessments, which is related to the scientific novelty of the presented article. The authors’ technology for the qualitative and quantitative interpretation of remote sensing data and digital elevation models with high resolution provides the opportunity to increase the spatial resolution of seismic microzonation forecasts, implemented by standard geophysical methods, and it determines the practical significance of completed research.

1. Introduction

Strictly speaking, the problem, which we solve, has a narrow engineering nature; it is reduced to parametric mapping of the elements of fracture tectonics and their ranking. Nevertheless, linking our work to the general research cycle, we involuntarily touch upon a wide-ranging area of engineering and geological surveys associated with the assessment of primary (endogenous) and secondary seismicity. Without departing from the essence of our work, we should note that these assessments are traditionally divided into the macroseismic method and seismic microzonation. The basis of macroseismic assessments is formed by small-scale maps of a general seismic hazard, specified on the basis of systematic seismological observations, updating the geological maps and population survey. In general, the structure of maps of a macroseismic (regional seismic) hazard is determined by active continental margins, young folded structures and the position of regional suture zones. The development of macroseismic assessments has been systematically implemented since the 1960s; it is associated, first of all, with the publications by American geophysicists (B. Bender (1972), C.A. Cornell (1968), etc.), through whose efforts the first versions of computer programs for seismic hazard calculation were developed. From the other side, in the USSR and Russian Federation, seismic hazard assessments had a systematic, methodological and physico-mathematical character too; this way of research is associated with V.I. Ulomov. (1999), S.V. Medvedev (1960) and others. Seismic microzonation is based on macroseismic determinations and parametrically takes into account the influence of soil conditions on seismic activity in the form of integral parameters of the seismic intensity increment. This increment has the form of correction to the initial seismic intensity, determined at the stage of macroseismic studies; this correction is varied in the range from −1 to +1 and calculated taking into account the set of geological and structural factors. In most countries, seismic microzonation is carried out in the parameters of seismic ground motions, which can be used immediately for further engineering computations [1,2]. In the Russian Federation, this microzonation is traditionally performed in units of intensity with accompanying assessments of resonance frequencies and simulation of ground vibration accelerograms [3,4]. In addition to obvious geological and geodynamic anomalies that determine the growth of a seismic hazard, the processes and objects affecting the growth of the seismic intensity increment include areas of intense gradient in the heights of the Earth’s relief, shallow depths of the groundwater table and significant dynamics of exogenous processes (karst and gully formation). The genetic relationship between geological and geodynamic anomalies, as well as the mentioned objects and processes, on the one hand, and the elements of fracture tectonics that form geodynamic zones, on the other hand, is apparent.

The idea of mapping geodynamic zones in the form of extended geostructural elements based on the analysis of remote sensing images is tested when compiling general seismic hazard maps for different probabilities of a seismogenic event [5]. However, the use of the distance basis in this case requires some additional level of a priori information concerning the position and intensity of seismogenic disturbances relative to the geodynamic zone, as well as the characteristics of the geodynamic zone itself (fracture width, surface fracture length, etc.). Such assessments are feasible on a regional scale, while with the transition to detailed seismic zoning and further to seismic microzoning, the level of a priori information required for representative conclusions on the detailing of seismogenic hazardous zones noticeably increases. Considering the high cost of field work, which determines this level of a priori information, we can talk about a natural tendency to assess the seismic hazard based on a family of indirect signs related to individual engineering structures [6]. At the stage when these structures are only being planned, attempts to find a correlation between the characteristics of the landscape observed in remote sensing images and the characteristics of seismicity recorded in some databases can be called obvious. The methodological basis for detecting this correlation is rather modest, since it relies on particular cases of landscape anomalies and does not contain a fundamental concept. For example, the method of localizing seismogenic hazardous areas based on mapping the lineament density maxima [7] can be called quite popular. As another example, we can point to the study of the empirical relationship of individual faults, manifested in the system of lineament structures on remote sensing images, on the one hand, and earthquake epicenters, on the other hand [8]. There is an interesting work on the assessment of the variability of the directional rose of lineaments under conditions of high-intensity seismogenic disturbances [9]. Summarizing publications and works similar to this, it can be argued that there is a direct functional relationship between the characteristics of seismicity and the characteristics of landscape variability: our earlier work confirms this [10] and allows us to apply quantitative estimates of this variability in a systemic form in this study.

At the initial stages of geophysical engineering and ecological engineering support of (road and house) building, when the known data are limited only by coordinates of the area of interests, reliable mapping of geodynamic zones and their quantitative parameterization using free-of-charge data can form the essence of the primary forecast of seismically hazardous areas for their subsequent verification with standard and more expensive approaches. Actually, the processing of remote sensing data and digital elevation models in engineering problems is the standard work for geomorphologists. The advantage of our approach in comparison to classical decoding of remote sensing data is in using quantitative and quantitative interpretation, similar to assessments applied in geophysics, while geomorphology operates with empirical classifications and visual expert estimations.

Turning to the choice of area of interest, one can note from the point of view of macroseismic assessment and seismic microzonation that this area was not anything remarkable. It was a commonplace district with a relatively smoothed relief of the Earth’s surface, intended for construction of greenhouses and located far from any significant endogenous seismic manifestations—in the central part of ancient Russian plate, in the north of the Central Russian Upland. Thus, the requirements of regulatory documents for geophysical engineering support of building would be purely formal ones, if one does not take into account the peculiarities of geological structures both at small-scale and large-scale levels. The investigated area is located on the joint of two structures—the northern slope of the Voronezh massif and the southern slope of the Moscow depression (Central Russian syneclise). In addition, it is mapped within the periphery of a large circular geological structure being 250 km in diameter. This circular structure is observed in the territory of the Tula and Oryol regions (Tula dome-ring structure) and coincides with a significant part of Central Russian Upland (Figure 1a) [11]. Hypothetically, dome-shaped formations are the result of generation of a hot spot or, in other words, the impact of endogenous plumes. On the Earth’s surface, similar vertical motion of plumes is often correlated with elevations, and along their periphery, they are limited by ring and concentric faults, including deep ones. These processes are complicated by the formation of a set of different order linear tectonic faults, which form, in particular, a fairly branched network in the Tula region. The most famous of them is the regional linear fracture with a latitudinal direction, crossing through the entire Oka basin [12], which includes the settlements Tula and Shchekino (Figure 1b). There are few kilometers from Shchekino to the area of planned construction. In the eastern part of the Tula region, there is one more regional fault with a northeastern direction—it is propagated approximately along the line of cities of Efremov, Uzlovaya, Novomoskovsk and Venev, i.e., it is placed approximately 35–40 km from the object of investigation. This fault probably originates from the Mediterranean coast and is traced to the territory of the Russian Federation in the northeastern direction up to the Yamal Peninsula.

Despite the location of the Tula region being in the zone of low-intensity seismic responses and rather far from the epicenters of earthquakes registered in the Alpine belt, the seismic activity is not negligible here. This is indirectly indicated by the morphology of the map of general seismic zoning (2015-C: the period of shaking is around 5000 years with a 1% probability of exceeding the calculated seismic intensity within 50 years) [14,15,16]. Although, in the vicinity of the investigated object, this map has a monotonic and low-intensity structure, there is a gradient zone from scores 5 to 6–7 of the Richter scale at the distance of about 130–140 km. Seeing the strike of anomalous seismic zones in sublatitudinal or northeastern directions, one can suppose their genetic relationship with the regional disjunctives mentioned above.

The vertical cross-section of the region [17] is characterized by monoclinal bedding of structural and compositional complexes with nonconformity at the contact of the Vendian, Archean-Proterozoic formations and Devonian deposits. In the relief of contacts of the stratified mining massif, the synform and Vendian layer thinning in the vicinity of investigated area can be observed. This feature of the geological structure is supposed to be the marker of the presence of regional deep faults which directly cross the estimated territory and are inherited by structural anomalies. Moving up through the cross-section, thick Devonian deposits as well as Carboniferous formations can be noted, including widespread limestones, which caused karst formation and, as the result, the danger of vertical subsidence during interaction of near-surface structural and compositional complexes with elements of fracture tectonics.

Karst phenomena observed north and south of Tula city [18,19,20,21] are mainly caused by interaction between groundwater and Mississippian limestones as well as Devonian gypsum (Figure 1c). The karst in Upper Devonian gypsum is especially intense. In the investigation region, the karst is observed in various forms: sinkholes, hollows, gullies, karst lakes, disappearing rivers, coastal depressions, karst depressions, niches and underground voids. The intensity of karst development is estimated by the area of its manifestation and by the specific volume of mapped karst cavities. This intensity is higher in the Tula region in comparison with the neighboring territories—this is due to the greater fragmentation and relatively high fracturing of carbonate strata. The dependence of karst formation both on dynamics of endogenous fracture tectonics and on exogenous factors [13,22,23] leads to continuous dynamics of karst development and to a nonzero hazard of the new vertical subsidence under the fast healing of karst formations. This implies the relevance of geodynamic evaluation of the area [20,24,25], especially under conditions of extensive mine workings, which are observed within the investigated object. Previously, the authors carried out the investigation of a seismogenic hazard in the vicinity of geodynamically active geoblocks because of activation of fault tectonics under conditions of hydraulic fracturing [10], with application of non-potential geofields. This article presents the extension of a systematic approach to the analysis of these geofields within Earth’s geodynamically stable crust.

2. Data and Methods

The key element of the database is remote sensing data combined with digital models of the Earth’s surface elevations, including, in particular, aerophotos and satellite images of the terrain. At the middle scale (1:50,000–1:100,000), remote sensing data were downloaded from open access USGS website EarthExplorer; at the large scale, remote sensing images were produced by an unmanned aerial vehicle survey based on the GeoScan technology. A detailed digital elevation model was developed in two ways: the first one (main) included digitizing originally analogous topographic maps available through a licensed Android application, ATLOGIS Geoinformatics; the second way was used for the mentioned digitizing update—it was based on stereophotography and laser sensing during the unmanned aerial vehicle survey.

Remote sensing data combined with digital elevation models are defined as scaled-down models of the terrain, where geological and structural anomalies are encoded in combinations of landscape components, and these combinations are encoded in the set of half-tone photo anomalies. Geological heterogeneities, which are not explicitly observed at the surface due to masking by soil and vegetation cover, intense weathering processes and anthropogenic influences, are clearly manifested in remote sensing data. Physically, this is explained by the evolutionary process of landscape formation within the considered territory, where all landscape components functionally depend on geological heterogeneities. At the stage of expert decoding [26], the lineament field is reconstructed both with remote sensing data and DME for reliable verification of results. The lineament as a linearized (straightened) landscape element is traced through heterogeneous forms of half-tone anomalies in remote sensing data and heterogeneous forms of relief in order to map the elements of fracture tectonics, obtaining the latest activation [25,27]. Expert decoding is based on visual estimations and manual, little-automated operations with one or the set of most contrasting spectrum bands. Here, the contrast is traditionally defined as parameter

$$K=\left(B_2-B_1\right)/B_2,0<K<1,$$

computed by the values of the optical density field between pairs of neighboring objects and averaged over the entire considered area.

Scalar fields of parameters derived from the initial scalar field by means of application of differential operators and estimation of exponential parameters are a hint in monitoring the results of expert decoding. The calculation of shadow forms of the geofield surface $f\left(x,y\right)$, when this surface is illuminated by a homocentric source [28], can be defined as the simplest transform of this type. If $\beta$ is the angle formed by the incident ray and the horizon plane, $\gamma$ is the angle between the conditional direction to the north and the projection of the incident ray onto the horizontal plane, counted clockwise, and

$$\varphi =arctg\left({\left(\raisebox{1ex}{$\partial f$}\!\left/ \!\raisebox{-1ex}{$\partial x$}\right.\right)}^2+{\left(\raisebox{1ex}{$\partial f$}\!\left/ \!\raisebox{-1ex}{$\partial y$}\right.\right)}^2\right)$$

is the slope angle of the virtual surface of the geofield, and

$$\phi ={\left(arctg\left(\raisebox{1ex}{$\partial f$}\!\left/ \!\raisebox{-1ex}{$\partial y$}\right.-\raisebox{1ex}{$\partial f$}\!\left/ \!\raisebox{-1ex}{$\partial x$}\right.\right)\right)}^2$$

is the angular sector of the view, where

$$\frac{\partial f}{\partial x}=\frac{\left(f_{31}+2f_{32}+f_{33}\right)-\left(f_{11}+2f_{12}+f_{13}\right)}{8\cdot K_S},$$

$$\frac{\partial f}{\partial y}=\frac{\left(f_{13}+2f_{23}+f_{33}\right)-\left(f_{11}+2f_{21}+f_{31}\right)}{8\cdot K_S},$$

$K_S=9$ is the number of matrix elements composed by values of the geofield in the sliding window:

$$3 \times 3 \left(\begin{array}{ccc}{f}_{11}& f_{12}& f_{13}\\ f_{21}& f_{22}& f_{23}\\ f_{31}& f_{32}& f_{33}\end{array}\right)$$

Then, the shadow forms of the geofield surface are represented by the parameter

$$S_R=255\cdot (\cos \left(90-\alpha \right)\cdot \cos \left(\varphi \right)+\sin \left(90-\alpha \right)\cdot \sin \left(\varphi \right)\cdot \cos \left(\gamma -\phi \right)).$$

Contrasting of the local half-tone structures is implemented by preliminary application of differential operators, for example, the Sobel operator [29,30] used in the convolution mode in the 3 × 3 window with the following transfer functions:

$${\Omega }_x=\left(\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right)\mathrm{and}{\Omega }_y=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right)$$

The result of recalculation is assigned to the center of the sliding window.

The left transfer function contrasts the geostructural lines of the sublatitudinal strike, and the right transfer function underlines geostructures of the submeridional strike (Figure 2a,b). Another example is borrowed from thermodynamics (Figure 2c) and reflects the ratio of ordered components of the geofield against disordered (chaotic) components in the structure of the scalar field (Figure 2c) [31,32]:

$$H=-{\displaystyle \sum _{i=1}^n\left(p_i\cdot {\log }_2{p}_i\right)}$$

entropy, where $p_i$ is the probability of the $i$-th state of the physical system (in practice, it is the frequency of realization of the $i$-th state, which is described as the interval of field values between two adjacent isolines).

The key procedure in the applied system of transforms is developed under this project in automated lineament analysis which accompanies and verifies expert lineament decoding. The algorithm includes searching for extremum points in the system of values of the initial scalar field as well as in the system of values of the calculated field of the horizontal gradient modulus. The next procedure is rotation of elementary lineaments around each of the mentioned points based on a standard rotation matrix and local spline interpolation. The final strike azimuth of a particular lineament is chosen by the criterion of the minimum of the dispersion functional. The result of such recalculation is represented by the set of coaxial structures derived by the generalization procedure. The latter is represented, first of all, by a statistical sorting of one-oriented structures with elimination of quasiperiodic lineament nets. Three standard procedures form the basis of presorting: determination of the direction rose of isolines of the scalar field based on the structure of the autocorrelation function; estimation of variance of strike azimuths of local lineaments; elimination of the set of elementary lineaments with the highest strike azimuth variances. The second component of the generalization procedure includes combining the tips of the nearest coaxial lineaments into macroscopic straightened structures (geostructural lines) with the subsequent smoothing of these structures. The algorithm also includes logical procedures, such as elimination of “blind” intersections of lineament structures of different azimuths as well as elimination of the priority of extended linear elements in comparison with arch concentric elements. The generalization of results of lineament decoding in the form of a map of the lineament density distribution is an important operation for final interpretation. Since the methodological foundations of seismic hazard maps are based on a direct analogy between lineaments and fracture tectonics [14,33], the increasing lineament density can be correctly interpreted as increasing the fragmentation in the upper part of the cross-section. Numerically, the procedure is based on counting the number of tips of selected lineaments per unit area of the integration interval (Figure 2d).

Morphotectonic analysis is an independent operation that verifies the morphostructural reconstructions described above and information about the position of areas of reduced stability in upper part of the cross-section. Here, the basic element includes the selection of geoblock boundaries, in which, firstly, the expert analysis of anomalies of the spatial stationarity parameter (entropy parameter) and, secondly, the optimal separation of this signal into long-wave and short-wave components are applied. The expert analysis of the scalar entropy field is performed visually with involvement of contour lines and half-tone images of entropy anomalies. The principal elements of morphostructural zoning include the existence of isometric and/or elongated anomalies, fixation of their linear size, the dominant strike azimuth of axes of anomalies and/or local gradient zones, existence of zones of sharp space-related gradients, variation in general “image” of isolines, variation in space-related stationarity of the signal and the existence of a significant jog in the structure of individual isolines or a set of isolines. Separation of the signal into components is initially associated with a methodological element for identifying the base and erosion-active layers in the relief structure. In the considered case, it is enough to select the threshold wavelength (space-related frequency) of the signal, separating the spectrum of this signal into two areas with fundamentally different structures of harmonics. The distribution of the threshold wavelength along the coordinate axes can differ in remarkable anisotropy. This is marked by the azimuthal variation in the autocorrelation radius approximated by the elliptical contour. As a consequence, the signal separation by the wavelength is reduced to averaging in the sliding elliptical window, with semiaxes being equal to the autocorrelation radius in two mutually orthogonal directions. Finally, extended elements are traced in the structure of components with different wavelengths $\lambda$, correlated with a priori known manifestations of fracture tectonics, manifested at different depths.

The complex of areal estimations described above is completed by deep geostructural reconstructions based on the analysis of landforms of the Earth’s surface. The physical foundations of the method and their correspondence to geological and engineering formulation of the problem are considered in [10]. The main idea is in the interpretation of the presence of periodic and quasiperiodic components in the structure of the spatial geological-geophysical signal by the standing wave, shaping

$$K\left(x\right)=A\cos \left(kx\right)\cos \left(\omega t-\phi \right)$$

in the spatial structure of nonequilibrium geological media [34,35]. Its wave dynamics are caused by periodic transgressive–regressive processes and regular tectonic-magmatic activations, covering both the entire Earth’s crust and the particular small vicinity of the point within the considered geoblock. Since the geological environment is a spatially distributed system, the documented vibrations of its individual elements are able to propagate in the volume of rock complexes, forming the waves. In the volume of the layered mining massif, they are recorded, in particular, as an alternation of synforms and antiforms with a definite step in the relief of the interface between any pair of structural–constitutional complexes. In other words, the wave image of the structuring of the mining environment is an extremely real approximation of phenomena associated with a nonequilibrium distribution of heterogeneities.

Analytical continuation of potential and nonpotential geosignals into the area of geological or technogenic heterogeneity (the source of a geofield anomaly) under such a model is equivalent to the transition to areas with increasing surface tension force, leading to the suppression of short-wave and low-amplitude components. If the geofield is given by a spatial distribution of absolute heights of the Earth’s surface relief, then under conditions of its noticeable differentiation and small anthropogenic change, individual elements of this geofield are hydrostatically (isostatically) compensated in the volume of the geological environment. Regarding general considerations, the process of hydrostatic compensation of positive and negative landforms is confirmed on the geochronological scale, when the rocks are characterized by elastic-plastic dynamics similar to the dynamics of a viscous fluid [36,37]. As the result, according to the skin effect, the larger shape of the Earth’s relief corresponds to the deeper hydrostatic compensation, which conforms to the ideas mentioned above about the growth of surface tension forces as a function of the depth.

3. The Results of Geostructural Reconstructions

The standard method of processing data of different scales, taking into account significant anthropogenic influences on the landscape and masking the terrain due to Quaternary structural–constitutional complexes and vegetation cover, was applied for analysis. The smaller scale of remote sensing images corresponds to the lower amplitude of anthropogenic impact on the modern landscape, and to the more contrasting responses of elements of fracture tectonics in the optical density field of remote sensing data and digital elevation models. Thus, the satellite image of the regional scale, being around 1:200,000 (Figure 3a), reflects the elements of planetary fractures, masked directly at the surface and obtaining definite peculiarities. One can see three systems of regional geodynamic zones: sublatitudinal, submeridional, with a northeastern strike, which, firstly, have a quasiperiodic character, and, secondly, a small curvature for each linear structure. The structure of the northwestern strike azimuth with a high contrast of manifestation on the remote image has aperiodicity and a significant curvature that define its confinement to the suture geoblock zone. As Figure 3c shows, this zone may reflect the response in the modern landscape of the main fault, which complicates the southwestern flank of the Pachelm Trough [38].

With regard to the submeridional geodynamic zone, there is the effect of shear displacement of elements of the geodynamic zone of the northwestern strike within the investigated area. Taking this into account, one can identify the northeastern strike of the axis along which the normal stresses are applied to the mining massif within the investigated area (Figure 3b). Then, applying the concept of conjugate fractures and the associated model of ellipsoid of deformations (Figure 3a), the generalized field of normal and shear stresses in projection onto the cartographical plane of the investigated territory is derived.

In mid-scale detailing of the geostructural image, trying to identify the generalized contour of the stressed rock massif, gradient operators and the generalization procedure discussed above were used. As a result, the image of the elliptical structure was mapped, fragmented along the set of sublatitudinal geodynamic zones, and along the zones of the northeastern strike and, to a lesser extent, the northwestern strike (Figure 4a). One can note the elliptical structure in Figure 4a surprisingly inherits the deformation ellipsoid shown in Figure 3a: there are coincidences both in the strike of the axes of the ellipses and in the distribution of the reconstructed normal and shear displacements which verify the mapped structural elements.

The contour of the elliptical structure is supposed to be the boundary of a certain weakened zone [41], which allows applying the analogy with the model of development of stress concentrators for localization of areas with an increased seismic hazard of the mining massif destruction. The mentioned model is related to the simplest mechanical formation, on the surface of which normal stresses σ are applied. It causes stress concentrators appearing in the bulk of the formation and in the small vicinity of the cavity (weakened zone), and the propagation of the zone of material plastic flow, developing into a crack. At the final stage of geostructural estimations in the cartographic plane, the large-scale data (approximately 1:14,000) were considered, reflected in Figure 4b. The smoothed morphology of the Earth’s relief implies the use of multiseason remote sensing images and digital elevation models in order to contrast the latest vertical subsidences within the investigated area.

The schematic map of geodynamic zones in Figure 5b is generalized up to the level of a lineament density map (Figure 3c) and verified at the stage of morphotectonic analysis (Figure 5a).

First of all, one can see the manifestation of a ring structure in the map of geoblock structuring (Figure 5a), which verifies the middle-scale reconstructions (Figure 4a). The area of contact of more than two geoblocks is considered as a forecast marker with regard to the localization of areas of “Possible Earthquake Sources” (PES), where there is significant variation in theof entropy parameter (parameter of spatial stationarity of the Earth’s relief) (Figure 5b). At final stage of analysis, the ranked PES areas were compared with the lineament density map, where the maxima detected the areas of increased fragmentation of the geological environment. The coincidences of PES areas of the third and the fourth rank with the located maxima of the lineament density field were indicated by triangular markers and drawn in the final map of prognostic areas of the seismogenic hazard (Figure 5c). It is particularly remarkable that most of the selected PES areas were attracted to a priori specified zones of subsidence.

Considering the independence of decoding criteria, the combination of heterogeneous data (remote sensing images and digital elevation model) at different scale levels and the application of mutually verifying interpretation methods, the forecasting map in Figure 6c can be considered as the completed result requiring confirmation by geophysical estimations. Strictly speaking, this result is not the standard indicator adopted in the technology of seismic microzonation: using a combination of several parameters and recalculation methods, we just localize the region of reduced stability of the rock massif with no quantitative assessment of its stresses, displacements and increment in seismic intensity. Nevertheless, as we noted above, the problem is precisely in the localization of areas of a high seismic hazard for their subsequent quantitative description by more expensive methods of seismic microzonation.

4. Discussion

Taking into account the last thesis about geophysical detailing, the interpretation procedure was extended by the solution to an inverse problem based on the recalculation of digital elevation models of the Earth’s relief in the volume of the geological environment. The final aim of this procedure was geostructural reconstruction in the plane of the vertical cross-section (Figure 6). The algorithm of recalculation used the wave analogies described above.

The profile line is drawn across the dominant strikes, reconstructed in Figure 5a,b. The area of interest is located in the center of this profile and occupies its fifth part: seeing a small step of discretization in the digital elevation model (it is about 0.5 m), this should allow both reconstruction of geostructural features in the investigated area and determining the presence of geostructural features of different ranks in the vicinity of this area.

In contrast to the reconstructions based on geophysical (seismic and electrical) data, in the considered case, the inverse problem is solved by the less rigorous model, in which the depth $h$ is computed as the relative parameter with application of empirical proportion $h=\lambda /\sqrt{2}$ between the linear size $\lambda$ of landscape anomalies and the depth of their endogenous roots. The effect of multilevel isostatic compensation widely known in gravity prospecting is supposed to be the physical model: the larger the size of the positive or negative form of the Earth’s relief, the deeper the form hydrostatically compensated (see the Pratt and Airy model). The profile for deep reconstruction and, accordingly, for the selection of absolute heights of the Earth’s surface relief is chosen in the sublatitudinal direction, crossing the investigated area in its middle part (Figure 6a). The signal detail is the first meters, which allows speaking about the representativeness of the geostructural image in small parts of the cross-section. The selected PES areas are attracted to the zones of regional fault manifestation, to structural–constitutional complex thinning (sharp lateral variability of their composition) and to the elements of increased permeability zones. The deep reconstruction of the cross-section in Figure 6b is considered to be an independent way of presorting PES areas localized at the stage of reconstruction in the cartographic plane (Figure 5b). According to industrial experience, this approach is popular for estimation of objects for which there are no engineering and geophysical data due to limited access and a complicated landscape.

Upon achieving the leading forecast, the authors implemented seismic measurements within the industrial order (Figure 7), organized according to the “common depth point” (CDP) technique [42,43].

The measurement system was central, with no offset. The spacing between geophones was equal to 2.5 m, and the step of “VibroPoints” (VP) was equal to 5 m. A 48-channel array composed by SGD-AD ground geophones with an eigen frequency of 10 Hz was used. The increase in the CDP order at the junction of geophone arrays led to the application of an additional 12 VP along the flank observation system with variable offset from 2.5 to 57.5 m. Thus, the CDP order changed from 12 to 24 because of the low relative intensity of the target wave in the fractured rocks. Office processing of field data was carried out by RadExPro software, designed for integrated processing of data of surface engineering seismics as well as for quality control of field seismic data [44]. In total, during the processing, 21 sets of seismic determinations were identified, which fitted into three types of structuring of the velocity and depth propagation of the wave extremum: one—according to CDP, and two—according to the “refraction correlation method” (RCM). Following RCM interpretation, the layer which had the roof at the depth of 20 m and an increased velocity distribution within the range of 789–1317 m/s was selected. This layer had a widespread occurrence within the area of interest. The comparison of locations of the different ranked PES areas (see Figure 5c) and the negative anomalies of the P-wave velocity (zones of increased watering and deconsolidation, see Figure 7) demonstrates a significant spatial correlation. This correlation especially concerns PES zones of the third and the fourth rank as well as the result of their presorting according to the lineament density map. Thus, the objectivity of the result of the method of PES areas’ leading forecast based on qualitative and quantitative analysis of remote sensing data and digital elevation models of the Earth’s relief can be noted.

5. Conclusions

Investigations are completed based on a transition from small to large scales: from the location of the polygon position at the contact of Precambrian structures up to deriving the particular detections of regional faults in the vicinity of the area of interest. Despite the seismologically stable structural-tectonic position, this area is placed in the zone of widespread limestones, which causes karst development and a final reduction in the stability of the upper part of the geological cross-section. The karst development intensity and final high-amplitude subsidence of the Earth’s surface are associated with both the activation of tectonic faults and exogenous processes: against the background of fast-treated karstogenic formations, the appropriate erosion processes and watering initiate continuous dynamics of karst. These processes reflect the importance of the morphostructural and geodynamic analysis of remote sensing data and digital elevation models of the Earth’s surface, focused on mapping the inherited manifestation of elements of fracture tectonics and the geoblock structure. This analysis is especially effective under conditions of an initially smoothed relief and developed sea sediments widespread in the area of interest. The presented way of investigation includes qualitative and quantitative approaches to interpretation, where the first one is focused on lateral tracing of heterogeneities, whereas the second one is reduced to deep reconstruction. In the case of quantitative interpretation, a digital elevation model reflecting the functional dependence between the landscape and pre-Quaternary geology is considered within the concept of wave structuring of the Earth’s crust. The importance of the considered approach is explained by the data type—remote sensing data and digital elevation models are related to broad access materials for any continental part of the Earth’s crust. The mentioned lateral tracing is based on parametric decoding, solving the problem of tracing geodynamic zones and reconstruction of the geoblock structure as well as patterns in the stress field of the mining massif. The positions of possible focuses of secondary earthquakes are verified both by decoding at different scale levels and by applying functionally independent parametric estimations. Deep reconstructions are based on the wave model developed in publications by Petrov O.V. and Movchan I.B. [35], in which the jump-like changes in standing wavelengths are associated with the law of their dispersion and reflect the position of the interfaces in the stratified geological environment. The subvertical position of these interfaces is interpreted as the response of through zones of increased geological environment permeability. The correlation of possible sources of secondary earthquakes with the zones of increased permeability allows predicting areas of increased seismological magnitude with high reliability. Subsequent industrial activity by the shallow seismic technique confirms the reliability of PES areas localized at the stage of remote sensing data processing in combination with a digital elevation model.