On the Issues of Spatial Modeling of Non-Standard Proﬁles by the Example of Electromagnetic Emission Measurement Data

: This paper examines the possibility of the spatial modelling of the Earth’s natural pulsed-electromagnetic-ﬁeld measured values, which form a closed proﬁle without the data inside. This geophysical method allows us to map active tectonic movement which breaches the integrity of pipes. During the experiment, 4.5 km of proﬁles were measured in the Admiralteysky district of St. Petersburg, Russia. Regular electromotive force (EMF) values and anomalous EMF values were obtained, ranging from 0 to 900 µ V and above 900 µ V, respectively. The anomalous values are associated with tectonic faults in the bedrock. The data obtained are characterized by complex spatial anisotropy associated with the development of two groups of tectonic faults of different orientations. The authors have considered the problems of the spatial modeling of the data obtained. The main problems, the solutions to which should allow the obtaining of adequate models, have been identiﬁed. Based on the analysis of the measurement results, geological features of the studied areas, as well as variography, the following possible solutions were proposed: changing the measurement technique; dividing the data array according to the main directions of anisotropy; the need to introduce additional correction coefﬁcients. The problem revealed in this article requires further research on the basis of the obtained results, which will reduce the cost and timing of such studies, and, as a result, give an opportunity to take into account active tectonic disturbances during the construction and scheduled maintenance of underground utilities, which is especially important within the framework of the concept of sustainable development.


Introduction
Ensuring the safety of human life in the urban environment, as well as limiting the negative impact of economic activities on the environment, is an integral part of the concept of the sustainable development of territories [1,2]. Close attention is also paid to the planning of residential areas and related engineering infrastructure-road junctions, electrical and heating networks [3], water supply and drainage systems, sewerage, etc.-as well as the improvement of territories [4] and the creation of social infrastructure [5,6]. These are extremely important and, beyond any doubt, relevant issues. However, in addition to creating new infrastructure, there is a need to optimally integrate existing infrastructures, subjected to rationalization, into this concept [7,8].
Water supply and drainage systems, as well as heating networks in many cities, were built a long time ago and without proper consideration of geological conditions [9,10].
However, geology has a significant impact on the integrity of linear objects due to the intensification of movement along tectonic faults [11], the influence of microbiota in fault zones [12], and erosive processes [13]. If, in new districts, these conditions can be accounted for during construction, in old districts, it seems optimal to partially change the structures by strengthening the sections crossing the fault. Regardless of further technological measures, taking into account active tectonic faults will potentially reduce the scope of repair work and potential adverse effects from pipeline ruptures.
In mining, when designing mine workings, and gas storage collectors, as well as at other special hazardous facilities, in particular those related to the nuclear power industry, tectonic faults and their intersection nodes have long been taken into account as areas of the most probable occurrence of dynamic phenomena and destruction of the massif (collapses, falls, rock bumps, etc.) [14][15][16][17][18][19]. The preliminary geodynamic zoning of territories is carried out for these objects [20].
Geodynamic zoning is a method for studying the block structure and the stress state of a rock mass, and identifying, assessing, and controlling currently active structures, such as active faults, geodynamic zones, and associated risk zones, in order to improve the geodynamic and environmental safety of industrial and civil objects during the development of subsoil and the Earth's surface.
For objects located at great depths, this work is mandatory, but for shallow-buried objects the performance of such work may also be expedient, especially when considering this issue in terms of sustainable development.
Previously, the authors have already compared data on accidents involving underground utilities and information on the position of known tectonic faults. This work made it possible to conclude that there is a high convergence of potentially geodynamically dangerous areas and areas of accidents [11]. In that study, the data were represented in graphs with characteristic features of the profile (electromagnetic (EM) sources, faults, georeferencing, etc.) plotted along the abscissa axis. However, this method does not make it possible to fully assess the nature of the distribution of values over the area, the presence of patterns, local anomalies that have no geological explanation, etc.
A logical continuation of the work already done is a systematization of the accumulated data on measurements of the Earth's natural pulsed electromagnetic field (ENPEMF) and geological data available from open sources, as well as their convenient and informative visualization. Integration of these data into the model will make their accounting, when planning and re-constructing urban infrastructure, simple and understandable even for those not proficient in this field.
St. Petersburg (59 • 57 00" north latitude, 30 • 19 00" east longitude), which in previous works has already served as an experimental area for studies, is well suited for a continuation thereof [11]. The geology of the city has been well studied to great depths, due to the deep burial of tunnels and subway lobbies (about 100 m) and high complexity of construction due to the features of the soil, characterized by significant heterogeneity and the development of swamps [21,22]. The geology of the city will be more thoroughly described in Section 2.1.
Studies of various scientific schools are devoted to the construction of maps based on the ENPEMF parameters. For example, in [23], based on the construction of maps of the electromagnetic field flux density, it is recommended to lay a pipeline, as well as to arrange drainage wells in places with minimum values of the electromagnetic field density. However, the paper does not disclose the method of map construction; therefore, it is difficult to determine its accuracy. In the study of [24], natural neighbor interpolation and inverse distance weight (Kriging) methods are used to estimate the electromagnetic field levels that workers are exposed to in the supply chain of an air-insulated substation. However, just like in the previous work, the data set is a regular grid of measurements. The paper [25] presents the results of magnetotelluric and electromagnetic time-domain soundings on the basis of which a geoelectric model was built, which makes it possible to determine the architectural geometry of the geological environment and variations in the physical properties of various blocks of the Earth's crust limited by large electrical inhomogeneities. Thus, this work also confirms the relationship between the parameters of the electromagnetic field and tectonic disturbances.
Various geostatistical methods are actively used to solve similar problems, in particular for modeling the spatial distribution of environmental data [26] and assessing the homogeneity of data [27]. These works, along with numerous other materials [28,29] containing examples of spatial modeling in geology, ecology, geophysics, etc., are a good basis for conducting this study. Nonetheless, they do not contain a solution to the main problem associated with modeling the data of measurements of the ENPEMF, which acts as a marker of active tectonic faults.
Nonetheless, in all works, the data array is represented either by a regular grid or by a random distribution of data over the entire study area. The task of constructing a spatial model in the conditions of a limited data set was not set. At the same time, it is worth noting that the anisotropy of the experimental data obtained in this study, on the one hand, is not spatially sustained, and, on the other hand, the territory could not be simply and unambiguously zoned due to the ubiquitous spread of the tectonics in three azimuthal directions. In addition, it was difficult to take measurements within the study area for three reasons: the dense and complex urban development of the city center, associated labor costs and time spending, as well as the research technique being based on profiling the area wherein tectonics are crossed with a dense line of measurements.
Thus, the aim of the study is to assess the possibility of constructing a spatial interpolation model in the conditions of random values of experimental data and the spatial distribution of data along the perimeter in the absence of values within the measured area. It should be noted that the authors of [30] claim that kriging is an effective way to estimate random fields with a local concentration. These are all prerequisites and possibilities to achieve the goal of choosing the methods of future interpolation and spatial visualization of the ENPEMF measurement data, but it is necessary to solve a number of related tasks. The main tasks include:

1.
Search for anisotropy of measurement data displaying tectonic faults.

2.
Determination of an effective way of displaying anisotropy when constructing a spatial interpolation model under experimental conditions. 3.
Assessment of the possibility of determining active tectonics during large-scale areal measurement of the ENPEMF parameters.

Study Area
A closed profile in the center of the city, in the Admiralteysky district ( Figure 1), was chosen as the reference one. From the viewpoint of building system, this site is suitable, since there are no industrial facilities, and the height of the residential buildings is insignificant (from 3 to 5 floors) and consistent. From a geological point of view, it is interesting because it is crossed by six tectonic faults at once. St. Petersburg is located in the northwest of the Russian Plate, which is an area of the foundation immersed under the sedimentary cover [33]. The occurrence of rocks is monoclinal, with an orientation to the south. The changing of layers to younger ones can be observed in the same direction.
St. Petersburg is characterized by crystalline rocks of the Archean and Proterozoic groups, mainly represented by granites, gneisses, and diorites. Basement rocks, on average, occur at a depth of 180-240 m [34], while in the south the depth can reach up to 300 m.
Above, are the deposits of the Vendian system: the Redkinskiy and Kotlinskiy horizons. The Redkinskiy horizon is represented by sandstones and siltstones with a thickness of 10-30 m, as well as by claystone-like clays and siltstones with a thickness of at most 15 m. The Kotlin deposits in the lower sub-formation are represented by sandstones and siltstones, and, in the upper sub-formation, by clays with rare sandstone interlayers. The thickness of this layer in some places reaches 150 m, but the surface of the Vendian complex is strongly eroded and varies greatly [35].
Quaternary deposits overlap the territory throughout the city. Their average thickness varies from 20 to 30 m, but in places of buried river valleys it reaches 100-130 m. The sedimentary cover is characterized by frequent lithological variability but is mainly represented by weak sandy-clayey soils.
The profile of the city is shown in Figure 2. The main features of the selected profile will be presented further in Section 2.3. St. Petersburg is located in the northwest of the Russian Plate, which is an area of the foundation immersed under the sedimentary cover [33]. The occurrence of rocks is monoclinal, with an orientation to the south. The changing of layers to younger ones can be observed in the same direction.
St. Petersburg is characterized by crystalline rocks of the Archean and Proterozoic groups, mainly represented by granites, gneisses, and diorites. Basement rocks, on average, occur at a depth of 180-240 m [34], while in the south the depth can reach up to 300 m.
Above, are the deposits of the Vendian system: the Redkinskiy and Kotlinskiy horizons. The Redkinskiy horizon is represented by sandstones and siltstones with a thickness of 10-30 m, as well as by claystone-like clays and siltstones with a thickness of at most 15 m. The Kotlin deposits in the lower sub-formation are represented by sandstones and siltstones, and, in the upper sub-formation, by clays with rare sandstone interlayers. The thickness of this layer in some places reaches 150 m, but the surface of the Vendian complex is strongly eroded and varies greatly [35].
Quaternary deposits overlap the territory throughout the city. Their average thickness varies from 20 to 30 m, but in places of buried river valleys it reaches 100-130 m. The sedimentary cover is characterized by frequent lithological variability but is mainly represented by weak sandy-clayey soils.
The profile of the city is shown in Figure 2. Below, a brief stratigraphic characteristic of the rocks in the investigated part of the city (in Figure 2, the region in the area of the Neva River) is presented.
The electromagnetic properties of soils were given in previous studies, wherein it was found that the use of the ENPEMF registration method is suitable for these soils. The properties of rocks do not prevent the field from passing through the strata, and the depth of the studied sources, which are geodynamically active faults (GDAFs), is insignificant, which also cannot result in attenuation of the electromagnetic wave causing an anomalous value in the measurement zone [11]. From the viewpoint of the fundamental theory, the waveform of the electromagnetic pulse generated by GDAFs, its frequency, period, and other parameters typical for the electromagnetic wave, are studied by various scientific schools, whose representatives are Mulev, S.N. [36], Frid, V. [37], Simpson, F. [38], and Sabaka, T. J. [39]. However, it is worth noting that studies are most often carried out for a specific area with its own geological structure and electromagnetic properties of soils [40,41].
Drainage systems, the most vulnerable from the viewpoint of accidents, represent linear objects of urban communications, which are usually shallow buried (15-90 m [42], in quaternary deposits). The sewers are buried deeper-in the Upper Kotlin clays (in the southern part of the city, this depth corresponds to the Lower Cambrian clays).
From a tectonic point of view, St. Petersburg is situated at the junction of two large structures: the Baltic Shield and the northwestern part of the Russian Plate. In this area, movements of the foundation blocks were recorded even in the modern, quaternary period. In addition, it was found that basement faults are prolonged in sedimentary rocksclays and sandstones-disintegrating them [34]. The Verkhnekotlinsky clays are also fractured, as evidenced by the data on subway tunnels. The city is characterized by fault tectonics of northeastern, northwestern and sublatitudinal strike. In addition, there are paleovalleys in the city.
Previously, scientists from the Mining University conducted an analysis of sections constructed during the engineering of a subway line, which resulted in a map of tectonic faults, shown in Figure 3 [43,44]. This map is the geological basis for research, as it contains the most recent and complete information among open sources.
For large cities it is not possible to cover the entire area with a dense line of measurements, especially at the initial stages of the study when it is fraught with errors in the choice of the measurement pitch, the length of profiles, as well as in the method of data Below, a brief stratigraphic characteristic of the rocks in the investigated part of the city (in Figure 2, the region in the area of the Neva River) is presented.
The electromagnetic properties of soils were given in previous studies, wherein it was found that the use of the ENPEMF registration method is suitable for these soils. The properties of rocks do not prevent the field from passing through the strata, and the depth of the studied sources, which are geodynamically active faults (GDAFs), is insignificant, which also cannot result in attenuation of the electromagnetic wave causing an anomalous value in the measurement zone [11]. From the viewpoint of the fundamental theory, the waveform of the electromagnetic pulse generated by GDAFs, its frequency, period, and other parameters typical for the electromagnetic wave, are studied by various scientific schools, whose representatives are Mulev, S.N. [36], Frid, V. [37], Simpson, F. [38], and Sabaka, T. J. [39]. However, it is worth noting that studies are most often carried out for a specific area with its own geological structure and electromagnetic properties of soils [40,41].
Drainage systems, the most vulnerable from the viewpoint of accidents, represent linear objects of urban communications, which are usually shallow buried (15-90 m [42], in quaternary deposits). The sewers are buried deeper-in the Upper Kotlin clays (in the southern part of the city, this depth corresponds to the Lower Cambrian clays).
From a tectonic point of view, St. Petersburg is situated at the junction of two large structures: the Baltic Shield and the northwestern part of the Russian Plate. In this area, movements of the foundation blocks were recorded even in the modern, quaternary period. In addition, it was found that basement faults are prolonged in sedimentary rocks-clays and sandstones-disintegrating them [34]. The Verkhnekotlinsky clays are also fractured, as evidenced by the data on subway tunnels. The city is characterized by fault tectonics of northeastern, northwestern and sublatitudinal strike. In addition, there are paleovalleys in the city.
Previously, scientists from the Mining University conducted an analysis of sections constructed during the engineering of a subway line, which resulted in a map of tectonic faults, shown in Figure 3 [43,44]. This map is the geological basis for research, as it contains the most recent and complete information among open sources. The measurement pitch and the length of the profile in this study were theoretically determined in advance and will be addressed in greater detail in Section 2.3. This work is mostly focused on the method of data processing, as it is the most relevant issue. At the same time, the task of the adequate spatial visualization of data comes first, including the selection of interpolation tools and the way of presenting the results of the analysis. In addition, determining a method of spatial interpolation for constructing such a map will allow for the interpolation of unknown profiles in the future.

Substantiation of the Processes of Electromagnetic Field Change in Zones of Tectonically Active Faults
According to the research, a pulsed electromagnetic field, which is a superposition of fields of different natures and frequencies, is recorded on the Earth's surface. As is known, the Earth's natural electromagnetic field has a frequency of about 10 Hz and is a low-frequency field, while fields of different natures will have frequencies that fall not only within the low-and mid-frequency ranges [45].
The main, non-interconnected, reasons for these fields are: 1. Mechanoelectric phenomena in rocks, which occur during elastic or plastic deformation and destruction (piezoelectric effect, electrification of cracks) [46]. Deforming forces, resulting in mechanoelectric phenomena, arise due to changes in mechanical For large cities it is not possible to cover the entire area with a dense line of measurements, especially at the initial stages of the study when it is fraught with errors in the choice of the measurement pitch, the length of profiles, as well as in the method of data processing. The above-mentioned reasons behind the inaccuracy of measurements are associated with the fact that, to date, not a single measurement technique exists for taking measurements by the method of electromagnetic emission.
The measurement pitch and the length of the profile in this study were theoretically determined in advance and will be addressed in greater detail in Section 2.3. This work is mostly focused on the method of data processing, as it is the most relevant issue. At the same time, the task of the adequate spatial visualization of data comes first, including the selection of interpolation tools and the way of presenting the results of the analysis. In addition, determining a method of spatial interpolation for constructing such a map will allow for the interpolation of unknown profiles in the future.

Substantiation of the Processes of Electromagnetic Field Change in Zones of Tectonically Active Faults
According to the research, a pulsed electromagnetic field, which is a superposition of fields of different natures and frequencies, is recorded on the Earth's surface. As is known, the Earth's natural electromagnetic field has a frequency of about 10 Hz and is a low-frequency field, while fields of different natures will have frequencies that fall not only within the low-and mid-frequency ranges [45].
The main, non-interconnected, reasons for these fields are: 1.
Mechanoelectric phenomena in rocks, which occur during elastic or plastic deformation and destruction (piezoelectric effect, electrification of cracks) [46]. Deforming forces, resulting in mechanoelectric phenomena, arise due to changes in mechanical and hydrogeological properties, the stress state of the rock mass, and soils, etc. It was found that GDAFs are sources of such mechanoelectric phenomena [47].

2.
Thermal processes that result in exoelectronic emission from minerals. This phenomenon is caused by the release of weakly bound water, giving rise to an increase in electrical conductivity and a decrease in the activation energy of current carriers, followed by an increase in the number of emitted electromagnetic pulses [48].

4.
Man-made load-cable and overhead transmission lines, transformer substations, subway lines, and electric transport. However, the operating frequency, typical for man-made load, as a rule, is 50 Hz [51]. It is worth noting that there are no subway lines or industrial buildings in the studied part of the city. Due to the fact that electric transport is widely used in the city, measurements were stopped when a trolleybus or tram passed, so the impact of man-made load was minimal.

5.
Pulse signals from remote sources [52][53][54]. The main reasons behind their occurrence are atmospheric phenomena, namely lightning discharges. As a rule, among the remote sources, atmospherics and whistling atmospherics are distinguished, which are characterized by an extremely low frequency starting from a fraction of Hertz [55]. 6.
Ionospheric and induced currents in the Earth's deep interior at a depth of 100-600 km [56].
Since the principle of superposition is valid for the electromagnetic field, it is logical to assume that the occurrence of a new process in rocks or a change in another (existing) process will result in a change in the resulting field. Therefore, in the place of measurement, to which any of the above-mentioned phenomena is confined, the parameters of the electromagnetic field will differ upward or downward. This paper, as already mentioned, discusses geodynamically active faults. The physics of the formation of changes in the parameters of the electromagnetic field in the zone of tectonically active faults is explained by the fundamental laws of electromagnetism, as well as by the theory of static and dynamic magnetic fields.
Within the framework of this study, the method of electromagnetic emission (EME), which was developed in the middle of the 20th century, was applied. The essence of the method is to measure the parameters of the electromagnetic field, including the intensity of the electromagnetic field (the pulse repetition rate of the field), the dominant frequency, signal amplitude, energy, and spectral density [57]. In this work, the derivative of only one parameter, namely the signal amplitude, expressed in EMF, was fixed. This approach has proved its efficiency in previous long-term studies [11]. Over the course of fieldwork, direct measurements of the EMF induced in the receiving antenna (magnetic rod with a ferrite core), described in detail in [11], were made when it entered a pulsed electromagnetic field. In order to minimize the man-made load (electric transport, road works, construction works, etc.), measurements were made early in the morning.
GDAFs are sources of pulsed electromagnetic waves (as sources of mechanoelectric phenomena), which result in a change in the resulting electromagnetic field on the surface of the measurement zone. The specific features of the passage of an electromagnetic pulse in an elastic medium, which rocks belong to, are disclosed in [58,59]. For most media, which include sands, clays and loams, granite, concrete, and asphalt, with varying degrees of humidity (these media represent the rock mass between GDAF and the measurement profile), the propagation of electromagnetic pulses depends mainly on the dielectric constant of the medium, and, to the greatest extent, on its imaginary component and humidity [60]. With an increase in the dielectric constant of the medium, the speed of signal propagation decreases. The penetration depth of the signal depends inversely on the specific attenuation and humidity. To minimize the impact of additional water cut, the studies were conducted in dry weather. Since the soil in which the electromagnetic wave propagates has a nonzero conductivity, the wave will attenuate when passing through it. However, due to the facts that the frequencies of the pulses under consideration lie in the range of 1-18 kHz and the wavelengths are more than 100 m with attenuation taken into account, after passing through the surface layer, an electromagnetic pulse can be registered with a receiving antenna [11].

Monitoring Network
The studied profile was made earlier within the framework of the study on correlating the locations of active tectonic faults and accidents on pipelines [11], and the data showed high correspondence. In total, five closed profiles with a total length of more than 25 km were made (observing the average length of every single profile of about 5 km). Table 1 shows the measurement data of the ENPEMF magnetic characteristic, which were obtained in 2014. The statistical description of the spatial distribution [61,62] is used in this work due to the fact that the data obtained during the experiment have random distribution and do not allow obtaining a formulaic distribution law based on the physical processes that cause the phenomenon of EMF in the receiving antenna (see Section 2.1). In this study, the spatial variable is the magnitude of voltage U, recorded by the multimeter in the voltmeter mode (APPA-109N). The variable is measured in µV. This value is: (1) The data array comprises a set of 210 points; (2) It is a continuous variable; (3) It is spatially referenced; (4) It is spatially correlated. Statistical parameters are presented in Table 2.  Min. 45 5 Per. 25 322. 5 6 Per. 50 630 7 Per. 75 1100 8 Max. 6250 All profiles crossed tectonic faults with different strike azimuths. Thus, when processing every single street, no problems were encountered with regard to identifying anomalies and correlating them with tectonic faults or man-made sources. However, when trying to perform large-scale processing and modeling data of the ENPEMF measurements, there was a difficulty in interpolating data for different directions. The selected profile, due to the absence of intersections of tectonic faults of different strikes directly within the profile, allows convenient dividing of a common database into sub-bases without additional manipulations.
At the same time, as mentioned earlier, the task was no longer set to prove a correlation between the ENPEMF anomalies and the location of tectonic faults, since the existence of this correlation was confirmed over the previous years of research [11,60]. To further substantiate this relationship, it is also necessary to combine the already existing separate databases, increase the volume of measurements and correctly select interpolation methods that are characterized by the greatest efficiency in identifying tectonic faults.
Closed profiles allow the intersection of tectonic faults at two locations. Thus, there was no need to make additional profiles through the same tectonic faults that had already been crossed. The streets selected for profiling look optimal in terms of the ratio of labor costs and the amount of data received ( Figure 4).
The study used the technique of discrete observations on the surface with a pitch of 20 m. With such a pitch, there will be at least 2 measurements in each zone of geodynamically active faults (GDAFs), the widths of which vary from 50 to 200 m. The total length of the profile was about 4.5 km.
The orientation of the magnetic moment of the antenna was chosen to be vertical. It is this direction of the antenna that has shown itself as the most efficient way in previous studies, and is also the only "geological" one, taking into account the slight slope feature in this part of the city.
The profile crosses several tectonic faults, which are grouped into two types: the upper part of the polygon is characterized by the development of tectonic faults of the northeastern strike, the lower part, by the northwestern one. Furthermore, in the lowest part of the polygon, at its southern end, the profile intersects and goes through 14 measurement points along a large tectonic fault of the northeastern strike. Thus, based on previous studies, two directions of anisotropy, typical for these measurements, should be expected.
In addition, this area is characterized by the following important geological features: • The profile is not complicated by the presence of rock contacts, the intersection of which also causes changes in the structure of the ENPEMF [34]; • The profile partially crosses the ancient riverbeds [34], which makes it possible to analyze the impact of this factor on measurements; • The thickness of quaternary deposits does not exceed 70 m [34].
Thus, taking into account the above-mentioned points and the general focus of the study, this profile makes it possible to analyze the possibilities of interpolation of the ENPEMF data by area and assess the prospects of visualization for the data type under consideration. substantiate this relationship, it is also necessary to combine the already existing separate databases, increase the volume of measurements and correctly select interpolation methods that are characterized by the greatest efficiency in identifying tectonic faults.
Closed profiles allow the intersection of tectonic faults at two locations. Thus, there was no need to make additional profiles through the same tectonic faults that had already been crossed. The streets selected for profiling look optimal in terms of the ratio of labor costs and the amount of data received (Figure 4).

Applied Methods
In this work, methods of statistical and geostatistical analysis were used, including methods of spatial modeling. First of all, a statistical analysis of the measured data array was carried out. The following parameters were determined: average, maximum, and minimum values, and frequency distribution histograms were constructed. The abovementioned parameters were determined according to the classical theory of statistical analysis [63]. Construction of frequency distribution histograms is necessary to determine the anomalous values corresponding to active tectonic zones. The purpose of the geostatistical analysis was to obtain the parameters necessary for building a spatial model-the main and secondary directions of data correlation, and determining the maximum distances at which the correlation is maintained. The spatial model firstly should not overload the map with data, but, at the same time, it is necessary to provide all its users with sufficient information to make it feasible to reveal possible sources of anomalies. The kriging method was used to build the spatial model.
Kriging is one of the methods used to create grids in the Surfer 13 program [64], which is designed specifically for spatial modeling. When choosing a method, a comparative analysis of their capabilities and limitations was carried out. One of the main requirements was the possibility of modeling anisotropy; the following methods were all suitable for this condition:
Triangulation with Linear Interpolation.
Among these methods, kriging was chosen because of its ability to work effectively with unevenly distributed data, data that vary significantly over short distances. The method is also not inclined to make bull's-eye patterns and does not smooth out heterogeneous data. Figure 5 shows a simplified block diagram of the research methodology.
Sustainability 2022, 14, x FOR PEER REVIEW 11 of 23 method is also not inclined to make bull's-eye patterns and does not smooth out heterogeneous data. Figure 5 shows a simplified block diagram of the research methodology. In previous studies [11], data were represented in graphs with characteristic features of the profile (electromagnetic (EM) sources, faults, georeferencing, etc.), plotted along the abscissa axis. However, this method does not make it possible to fully assess the nature of the distribution of values over the area, the presence of patterns, local anomalies that have no geological explanation, etc.
The most promising modern method of visualization is spatial modeling. In this case, one should not overload the map with data, but, at the same time, it is necessary to provide all its users with sufficient information to make it feasible to reveal possible sources of anomalies.
In accordance with all the earlier-mentioned features of the method used, it was decided that the model should comprise the following: tectonic faults, palaeovalleys, and, in general, the depth of bedrock, water bodies, and subway lines. Water bodies and subway lines are marked on the city map, therefore a web mapping service, Yandex.Maps, was used in the study, but any other open source can also be used.
Thus, a fragment of the city map was selected as a base layer, to which the rest of the data were subsequently linked layer-wise. The lower left edge of the map is linked to the zero point (0,0), while the third coordinate is provided only for the ENPEMF data set and contains information on the level of the measured signal. A rectangular coordinate system, wherein the position of a point in space is determined by the distance to the axes, was used. At this stage of the study, there was no need to link measurements to generally accepted coordinate systems, since general access to database filling has not yet been provided.
The following structure was chosen for the resulting map:  [18]. After digitizing and assigning the corresponding Z coordinates to the isolines, the data were presented in the form of a contour map ( Figure 6). In previous studies [11], data were represented in graphs with characteristic features of the profile (electromagnetic (EM) sources, faults, georeferencing, etc.), plotted along the abscissa axis. However, this method does not make it possible to fully assess the nature of the distribution of values over the area, the presence of patterns, local anomalies that have no geological explanation, etc.
The most promising modern method of visualization is spatial modeling. In this case, one should not overload the map with data, but, at the same time, it is necessary to provide all its users with sufficient information to make it feasible to reveal possible sources of anomalies.
In accordance with all the earlier-mentioned features of the method used, it was decided that the model should comprise the following: tectonic faults, palaeovalleys, and, in general, the depth of bedrock, water bodies, and subway lines. Water bodies and subway lines are marked on the city map, therefore a web mapping service, Yandex.Maps, was used in the study, but any other open source can also be used.
Thus, a fragment of the city map was selected as a base layer, to which the rest of the data were subsequently linked layer-wise. The lower left edge of the map is linked to the zero point (0,0), while the third coordinate is provided only for the ENPEMF data set and contains information on the level of the measured signal. A rectangular coordinate system, wherein the position of a point in space is determined by the distance to the axes, was used. At this stage of the study, there was no need to link measurements to generally accepted coordinate systems, since general access to database filling has not yet been provided.
The following structure was chosen for the resulting map: 1. City map; 2.
Processed data of EMF measurements. reintegrated into the model.
The basis, demonstrating the depth of bedrock, was a map of the paleovalleys of St. Petersburg [18]. After digitizing and assigning the corresponding Z coordinates to the isolines, the data were presented in the form of a contour map ( Figure 6). Furthermore, the lines of tectonic faults and measurement points were integrated into the model (Figure 7).
The measurement points at this stage are visualized as five classes by means of a color scale, corresponding to intervals with an equal number of values. It is worth noting that this format does not allow one to fully assess the results of measurements to determine whether tectonics are active and whether there are anomalies of non-geological origin.  Furthermore, the lines of tectonic faults and measurement points were integrated into the model (Figure 7).
The measurement points at this stage are visualized as five classes by means of a color scale, corresponding to intervals with an equal number of values. It is worth noting that this format does not allow one to fully assess the results of measurements to determine whether tectonics are active and whether there are anomalies of non-geological origin.   The measurement points at this stage are visualized as five classes by means of a color scale, corresponding to intervals with an equal number of values. It is worth noting that this format does not allow one to fully assess the results of measurements to determine whether tectonics are active and whether there are anomalies of non-geological origin.

Statistical Description of Data and Construction of Variograms
A frequency analysis of the data distribution was applied to determine the regular and abnormal values of the voltages (µV). To construct histograms, the results obtained were divided into classes. A total of 14 classes were defined, as shown in Table 3. After collecting the data, Exploratory Data Analysis was carried out. The analysis included the determination of absolute (F A ), relative (F R ), and cumulative (F C ) frequencies, followed by the construction of the corresponding histograms ( Figure 8).

Statistical Description of Data and Construction of Variograms
A frequency analysis of the data distribution was applied to determine the regular and abnormal values of the voltages (µV). To construct histograms, the results obtained were divided into classes. A total of 14 classes were defined, as shown in Table 3.
After collecting the data, Exploratory Data Analysis was carried out. The analysis included the determination of absolute (FA), relative (FR), and cumulative (FC) frequencies, followed by the construction of the corresponding histograms ( Figure 8).   , and the values from the third class can be referred to as anomalous U values due to local changes in the ENPEMF parameters. Thus, it was found that the anomalous profile values start from 931 µV.
Due to the fact that a significant deviation in the EMF value was revealed in different parts of the profile, it was decided to divide the profile into two parts-northeastern (NE) and northwestern (NW), the analysis results of which are presented in Table 4. Dividing was carried out on the basis of the main directions of tectonic faults in such a manner that each of the sub-bases contain information on measurements made through the tectonic faults of one direction (Figure 9). This profile allows doing this without additional complex manipulations. Due to the fact that a significant deviation in the EMF value was revealed in different parts of the profile, it was decided to divide the profile into two parts-northeastern (NE) and northwestern (NW), the analysis results of which are presented in Table 4. Dividing was carried out on the basis of the main directions of tectonic faults in such a manner that each of the sub-bases contain information on measurements made through the tectonic faults of one direction (Figure 9). This profile allows doing this without additional complex manipulations. The methodology for dividing the database into several sub-bases was not the ultimate goal of this study, since initially the goal was only to determine the efficiency of such an approach.  The methodology for dividing the database into several sub-bases was not the ultimate goal of this study, since initially the goal was only to determine the efficiency of such an approach.
As it was found during the study, at frequency variations typical for the main profile and its northwestern part (NW), anomalous values start from the third frequency-from 931 µV and 886 µV and above, respectively. (The values for the NW part of the profile are inverted, since the first frequency does not correspond to the minimum voltage value). It is worth noting that the values belonging to the NW part of the profile are characterized by a high, most-common value, higher than in the entire profile and its NE part. The NE part of the profile is characterized by a different distribution of values-anomalous values start from the second frequency, at voltage values from 620 µV and above. Moreover, for the NE part of the profile, there is a significant increase in the voltage value on one of the linear parts. The operators also noticed this feature during the measurements. However, the attempts made to repeat the measurements gave similar results. Among the most likely sources of such a signal level, the government buildings located in this place, i.e., the Prosecutor's Office and the Ministry of Justice, should be listed. Nonetheless, the nearby building of the Constitutional Court is not geographically connected with the area of high values. On the other hand, this place is characterized as difficult from a geological point of view, which is known in connection with St. Isaac's Cathedral, which is experiencing a slow immersion and is a constant object of study for both geologists and surveyors [65,66].

Geostatistical Analysis of Spatial Continuity: Variograms Estimation and Modeling
In this step, a geostatistical study of the data was carried out to characterize the spatial distribution of variable U for the entire study area. For the purpose of analysis and modeling spatial continuity, the variograms were calculated and modeled, for subsequent application in the kriging interpolation models [67]. The variograms summarize the main characteristics of the spatial continuity pattern of a physical phenomenon, the spatial dispersion of variable U, in particular:

•
The large-scale spatial trend; • The local spatial variogram anisotropy; • The ranges of the spatial correlation of variable U.
It is necessary to determine the anisotropy ratio and main directions of spatial continuity as accurately as possible to ensure the quality of the spatial model. As is known, small ranges, which mean high variability at short distances, will result in the appearance of smooth interpolated images, and too large a range will result in artificial patterns, mainly in areas with a lack of data. The optimal choice for the anisotropy range's value is based on the analysis of the variogram conditioned by expert knowledge of physical phenomena [68].
Based on the analysis, three variograms were constructed (Figures 10-12): a variogram for a full data array covering the entire area, and variograms for two distinct areas, the north and south parts.
Sustainability 2022, 14, x FOR PEER REVIEW 15 of 23 931 µV and 886 µV and above, respectively. (The values for the NW part of the profile are inverted, since the first frequency does not correspond to the minimum voltage value.) It is worth noting that the values belonging to the NW part of the profile are characterized by a high, most-common value, higher than in the entire profile and its NE part. The NE part of the profile is characterized by a different distribution of values-anomalous values start from the second frequency, at voltage values from 620 µV and above. Moreover, for the NE part of the profile, there is a significant increase in the voltage value on one of the linear parts. The operators also noticed this feature during the measurements. However, the attempts made to repeat the measurements gave similar results. Among the most likely sources of such a signal level, the government buildings located in this place, i.e., the Prosecutor's Office and the Ministry of Justice, should be listed. Nonetheless, the nearby building of the Constitutional Court is not geographically connected with the area of high values. On the other hand, this place is characterized as difficult from a geological point of view, which is known in connection with St. Isaac's Cathedral, which is experiencing a slow immersion and is a constant object of study for both geologists and surveyors [65,66].

Geostatistical Analysis of Spatial Continuity: Variograms Estimation and Modeling
In this step, a geostatistical study of the data was carried out to characterize the spatial distribution of variable U for the entire study area. For the purpose of analysis and modeling spatial continuity, the variograms were calculated and modeled, for subsequent application in the kriging interpolation models [67]. The variograms summarize the main characteristics of the spatial continuity pattern of a physical phenomenon, the spatial dispersion of variable U, in particular: • The large-scale spatial trend; • The local spatial variogram anisotropy; • The ranges of the spatial correlation of variable U.
It is necessary to determine the anisotropy ratio and main directions of spatial continuity as accurately as possible to ensure the quality of the spatial model. As is known, small ranges, which mean high variability at short distances, will result in the appearance of smooth interpolated images, and too large a range will result in artificial patterns, mainly in areas with a lack of data. The optimal choice for the anisotropy range's value is based on the analysis of the variogram conditioned by expert knowledge of physical phenomena [68].
Based on the analysis, three variograms were constructed (Figures 10-12): a variogram for a full data array covering the entire area, and variograms for two distinct areas, the north and south parts.    According to the results obtained, it can be concluded that the anisotropy directions are displayed partially or completely in the direction that is adjacent to the direction of GDAFs.
An exponential model was chosen to fit all the variograms, since spatial autocorrelation decreases exponentially with increasing distance. The variogram anisotropies were modeled by a common geometric anisotropy, defined, in 2D, by the ranges of the main direction and the perpendicular direction to the main one. Table 5 presents the results of the analysis of the obtained variograms and the spatial direction of GDAFs (positive directions were taken into account counterclockwise according to the trigonometric circle).  According to the results obtained, it can be concluded that the anisotropy directions are displayed partially or completely in the direction that is adjacent to the direction of GDAFs.
An exponential model was chosen to fit all the variograms, since spatial autocorrelation decreases exponentially with increasing distance. The variogram anisotropies were modeled by a common geometric anisotropy, defined, in 2D, by the ranges of the main direction and the perpendicular direction to the main one. Table 5 presents the results of the analysis of the obtained variograms and the spatial direction of GDAFs (positive directions were taken into account counterclockwise according to the trigonometric circle). According to the variograms, the direction of anisotropy is displayed for the following directions: 1.
For the entire area, the main direction of anisotropy corresponds to the azimuth of 31 degrees (counterclockwise rotation), Figure 10; this direction corresponds to three directions (27,25,26 degrees) of GDAFs. Other directions of GDAFs are not revealed by the variograms.

2.
After dividing the profiles into the two areas, the main directions of anisotropy are: 2.1. NE-65 degrees; northeastern part, Figure 11; 2.2.
According to the results obtained, it can be concluded that the anisotropy directions are displayed partially or completely in the direction that is adjacent to the direction of GDAFs.
An exponential model was chosen to fit all the variograms, since spatial autocorrelation decreases exponentially with increasing distance. The variogram anisotropies were modeled by a common geometric anisotropy, defined, in 2D, by the ranges of the main direction and the perpendicular direction to the main one. Table 5 presents the results of the analysis of the obtained variograms and the spatial direction of GDAFs (positive directions were taken into account counterclockwise according to the trigonometric circle).

Kriging Maps
Based on the variogram models, the following kriging maps were built (Figures 13  and 14) by using the variogram models of Section 3.2 and the U data. The Surfer [64] and GEOMS2 software packages were used for the estimation of kriging maps. GEOMS2 is a free software product (GNU General Public License version 3.0 (GPLv3)).

Kriging Maps
Based on the variogram models, the following kriging maps were built (Figures 13  and 14) by using the variogram models of Section 3.2 and the U data. The Surfer [64] and GEOMS2 software packages were used for the estimation of kriging maps. GEOMS2 is a free software product (GNU General Public License version 3.0 (GPLv3)). As one can see from the kriged maps (Figure 13), the main spatial patterns follow the main variogram spatial structures, which do not correspond to the directions of faults. According to the results obtained, the following conclusions can be drawn: 1. The construction of a spatial model of profiles with an irregular grid with a pronounced geometric shape and high concentration of data along the profile perimeter is difficult to implement. There is only one spatial trend that coincides with only one global direction of GDAFs (northeast 26°). This can be explained by the difference in anomalous values. In addition, in the northeast direction, there is a dense accumulation of anomalous values in the area of St. Isaac's Cathedral (see the description above). between values of 600 µV and 900 µV) may be difficult to find without additional analysis, and potentially the subsequent selection of the correction coefficient allows us to make the transition from the actual measurement data to more convenient values associated with the classes. Moreover, anomalous values of class three and higher (depending on the profile, classes may vary) should have a 100% correlation in terms of displaying GDAFs, which also requires the introduction of additional correction coefficients for sample data.

Discussion
Despite the fact that dividing the database into sub-bases by the direction of the strike of tectonic faults has proved to be a promising method, the proposed approach still requires the addressing of some related problems. In particular, it is necessary to take into account those EMF values, the overestimated level of which is associated with the technogenic source. The solution to this problem can be the replacement of the amplitude data in µV with a special dimensionless correction coefficient. Furthermore, the introduction of such a correction coefficient will help to solve the problems of the interpolation clarity As one can see from the kriged maps (Figure 13), the main spatial patterns follow the main variogram spatial structures, which do not correspond to the directions of faults. According to the results obtained, the following conclusions can be drawn:

1.
The construction of a spatial model of profiles with an irregular grid with a pronounced geometric shape and high concentration of data along the profile perimeter is difficult to implement. There is only one spatial trend that coincides with only one global direction of GDAFs (northeast 26 • ). This can be explained by the difference in anomalous values. In addition, in the northeast direction, there is a dense accumulation of anomalous values in the area of St. Isaac's Cathedral (see the description above).

2.
Dividing the profile into northeastern and northwestern sectors did not result in the initially expected results ( Figure 14). The reason for the resulting trends is a non-standard, open, irregular grid of data with a big, empty, no-data space within the study area. It makes it very hard to assume the sample data values are spatially representative of the study area. It is necessary to conduct additional studies to collect more data in the "empty" sections of the profile measured by the EMP method, for which it is possible to build a more accurate spatial model without changing the approach to the method implementation. 3.
Significant spatial variation in anomalous values gives rise to problems during modeling. According to the conducted studies, the U values that correspond to GDAFs may belong to different classes, and a correlation between these classes (for example, between values of 600 µV and 900 µV) may be difficult to find without additional analysis, and potentially the subsequent selection of the correction coefficient allows us to make the transition from the actual measurement data to more convenient values associated with the classes. Moreover, anomalous values of class three and higher (depending on the profile, classes may vary) should have a 100% correlation in terms of displaying GDAFs, which also requires the introduction of additional correction coefficients for sample data.

Discussion
Despite the fact that dividing the database into sub-bases by the direction of the strike of tectonic faults has proved to be a promising method, the proposed approach still requires the addressing of some related problems. In particular, it is necessary to take into account those EMF values, the overestimated level of which is associated with the technogenic source. The solution to this problem can be the replacement of the amplitude data in µV with a special dimensionless correction coefficient. Furthermore, the introduction of such a correction coefficient will help to solve the problems of the interpolation clarity and displaying GDAFs in the model in conditions of a significant difference in measurements of the EMF value associated with tectonic faults having different coordinates. The features of filtering a man-made signal and anomalies of different levels must be considered separately.
Within the framework of this study, such a scope of work cannot be presented. Thus, the development of such a correction coefficient may become the topic of further research.
In addition, for a widespread application of the method, it is necessary to adapt it to work with profiles in tectonically highly disturbed areas, as well as in areas where it is difficult to group tectonics into sub-bases.
When analyzing variograms for which constant coefficients corresponding to only two values, regular and anomalous, were introduced, the main direction of anisotropy did not correspond to the direction of GDAFs as well, but was co-directional with the direction of the spatial trend for all classes of values. If the value belonged to the range of 0 to the anomalous value, then a coefficient of 0 was assigned; if the value was equal to the anomalous value and exceeded it, then the coefficient value was equal to 1. Regular and anomalous values were determined based on frequency diagrams and correspond to 931 µV, 620 µV and 886 µV for the main profile and sub-bases.
The discrepancy between the main correlation direction when dividing the database into sub-bases can be explained by the following reasons: 1.
the lack of data on one of the linear parts of the profile. The solution to this problem can only be the introduction of additional profiles; 2.
the lack of correlation between values belonging to different classes. The solution to this problem, as mentioned earlier, can be the introduction of additional correction coefficients reflecting the belonging of the value to normal or anomalous levels or classes.
Moreover, in the future, it is necessary to take into account man-made sources (subway, substations, etc.) in the model, as well as welding sites outside GDAF zones, because repeated pipe breaks often occur in such places.
Thus, the main directions of further research are the following: the development of a correction coefficient that allows the separating of anomalous values of all levels from regular ones; work on filling the model; development of a measurement technique taking into account newly identified data features; further statistical processing of data.

Conclusions
In this research, an attempt was made to find an approach to the spatial modelling of data characterized by location along a perimeter with a complex geometric shape and two directions of data anisotropy, as well as by step transition from high values to low.
In the course of the study, the tentative attempt to create a spatial model of the EN-PEMF measurement data did not allow for determining an unambiguous methodological approach to such data visualization, but made it possible to reveal the main challenges of the existing databases. The solution to which will allow us to approach the solution to the main issue.
First of all, the main problem of the processed data is the absence of measurement points inside the profile with their high density around the perimeter. It was established that it is necessary to introduce amendments to the technique of detecting the location of GDAFs by the EMP method. There are two potentially possible options: changing the monitoring network for this method; introducing limitations for the method's applicability, which are based on defining the maximum permissible distance between measurement points. The second option requires additional research to find the maximum permissible distances. It is worth noting that these amendments should be developed specifically for the urban environment and large areas with a complex building system. At the same time, options must meet both the requirements for identifying GDAFs (equally with the adopted profiling technique using a dense line of measurements across tectonic faults), and the requirements for auxiliary points when constructing a kriging network. Both possible options require additional research.
Secondly, the splitting of the profile into sub-bases did not provide the expected result. The expected azimuths of anisotropy did not coincide with the calculated ones. The most probable explanation is the high heterogeneity of the data. Therefore, establishing the method of dividing the profile into sectors also requires additional research aimed at finding the optimal division criteria and limitations of the applicability of the dividing method. Furthermore, to solve this problem, it is necessary to conduct research on determining the coefficients, which will allow the interpolation of data between anomalous values of different levels (classes).
Further development of the proposed methods will facilitate the further modelling of profiles such as those obtained in the linear profiling mode, characterized by step differences and described by multiple directions of anisotropy. Similar data are typical for Earth sciences.
The results of this research will contribute to improving the efficiency of the designing of underground utilities and reduce the accident rate involving underground utilities by means of the early repair of pipeline sections located in zones of geodynamically active faults. The obtained spatial models can be used as the basis for planning the repairs in underground utilities, which in turn will reduce the probability of accidents, the consequences of which may not only result in changes in traffic patterns, transport collapses, and road accidents, but also human and economic losses. The circumstances listed above correspond to one of the 17 Sustainable Development Goals, namely Goal 9, "Build resilient infrastructure, promote sustainable industrialization and foster innovation", Target 9.1: "Develop quality, reliable, sustainable and resilient infra-structure, including regional and transborder infrastructure, to support economic development and human well-being, with a focus on affordable and equitable access for all". It should also be noted that, in addition to the positive impact on infrastructure, monitoring of pipelines will reduce the negative impact on the environment, which is especially relevant for St. Petersburg due to the proximity of the Gulf of Finland. It also corresponds to Goal 11 "Sustainable cities and communication", Target 11.6 "By 2030, reduce the adverse per capita environmental impact of cities, including by paying special attention to air quality and municipal and other waste management". Funding: This study was carried out at the expense of the subsidy for the state assignment in the field of scientific activity for 2021, No. FSRW-2020-0014. This research work by Prof Anton Rassõlkin has been supported by the Estonian Research Council under grant PSG453 "Digital twin for propulsion drive of autonomous electric vehicle".