Optical remote sensing method for detecting surface traces as indicators of underground structures

Abstract

The invention relates to a computer-implemented-method for detecting underground structures through analysis of time-series optical remote-sensing data. Unlike approaches that rely on single images, vegetation indices, or pairwise change detection, the method analyzes temporal spectral behavior of the ground surface as an information-bearing signal indicative of subsurface anomalies. The method comprises obtaining a plurality of georeferenced optical images acquired over an extended time-period, spatially aligning the images, extracting spectral values at defined spatial sampling loci, and converting the extracted spectral values into standardized single-channel values. A multi-stage normalization process is applied to mitigate natural soil heterogeneity and transient environmental effects while amplifying persistent temporal patterns associated with underlying subsurface structures. The method generates an output spatial dataset, optionally visualized as a synthetic map, to highlight locations indicative of underground anomalies. The approach provides a deterministic, training-independent, and computationally-efficient-solution suitable for applications including archaeology, underground infrastructure monitoring, and related fields.

Description

[0001] T I T L E

[0002] Optical Remote Sensing Method for Detecting Surface Traces as Indicators of Underground Structures

[0003] D ESC RI PT I O N

[0004] The present invention addresses a critical challenge in remote sensing: the reliable and efficient detection of subsurface structures through indirect surface indicators. Many existing advanced techniques require specialized expertise in remote sensing and complex computational methods, limiting their accessibility to non-specialized users. It leverages innovative image processing techniques applied to time-series data from satellites, drones, and other aerial platforms. The invention is particularly applicable for identifying traces of subsurface constructions — whether archaeological remains or modern infrastructures — by analyzing temporal variations and changes in ground surface properties.

[0005] The primary objective of the invention is to detect subsurface structures through the indirect identification of subtle surface changes, such as variations in soil color and other indicators, which are typically associated with underground anomalies. By employing a robust timeseries analysis combined with a novel soil color normalization process, the method compensates for natural variability in soil characteristics and seasonal changes. This enables the reliable extraction of persistent signals that may indicate the presence of buried structures. Consequently, the invention provides an effective and computationally efficient tool for applications in archaeology, underground infrastructure monitoring, and related fields, facilitating the automated detection of subsurface features using remotely sensed data.

[0006] Level of Prior Art and Its Evaluation

[0007] Remote sensing techniques for investigating subsurface anomalies have been reported across multiple application domains, including archaeology, geotechnical engineering, and infrastructure monitoring. These techniques include conventional image processing, spectral analysis, and machine-learning-based approaches applied to aerial or satellite imagery. However, the prior art does not disclose a systematic, deterministic methodology that treats temporal variation of soil-surface spectral behavior as an information-bearing signal and quantitatively exploits that variation to infer the presence of subsurface structures, which is the core subject matter of the present invention.Paz-Pujalt et al. disclose a method for identifying surface failures in infrastructure by comparing images acquired at different times to detect visible morphological changes, such as cracks or surface disruptions. The method relies on detecting differences between images representing distinct states of the surface. This approach is directed to the identification of overt surface damage and does not address subsurface structures whose influence manifests through indirect spectral changes rather than visible morphological alterations. Moreover, the method operates through first-order image comparison, evaluating differences between two images without analyzing the temporal behavior of pixel values across an extended time series. It therefore does not disclose or suggest computing temporal baseline statistics, local neighborhood statistics, or higher-order metrics derived from the temporal evolution of spectral values, as required by the present invention.

[0008] WO 2013 / 070945 A1 discloses image analysis systems and methods for analyzing timesequence images in order to detect and quantify changes occurring between images acquired at different times. The disclosed techniques employ registration, alignment, and comparison of reference and test images to identify surface-level changes resulting from deformation, movement, or the appearance or evolution of visible features. The methodology is primarily directed to detecting explicit changes between images or image pairs and to highlighting such changes once they occur, rather than synthesizing information from temporal variability itself to infer a stable underlying subsurface structure.

[0009] In terms of subsurface infrastructure and construction detection, recent years have seen many methodological reports in the literature that attempt to detect subsurface structures using images obtained from aerial platforms, particularly satellite images in the visible and hyperspectral ranges (Melilos, 2019). These methodologies have significant results, especially in archaeological applications (El-Behaedi, 2022). Generally, the underlying principle of these methodologies is that subsurface structures indirectly affect their environment (e.g., by altering local soil moisture content, drainage rates, soil composition and geometry, and variability in various indices that assess the qualitative characteristics of vegetation). These changes, when observed on the surface, may indicate the presence of subsurface structures. For instance, Pascucci et al. (2010) assess the suitability of remote sensing data from CASI (Compact Airborne Spectrographic Imager) and ATM (Airborne Thematic Mapper) for detecting subsurface archaeological structures byanalyzing the spectral footprint resulting from different plant cover. Kaimaris and Patias (2014) and Kaimaris (2024) acknowledge that soil and vegetation marks can be indicative of subsurface structures, recognizing that subtle variations in soil color or reflectance may serve as markers. However, the above methods remain primarily observational and qualitative. Crucially, they do not incorporate a systematic, quantitative analysis of timeseries data to normalize soil color differences or enhance the detectability of temporary soil anomalies. Their recognition of soil markers is mainly descriptive, lacking a defined methodology for extracting and processing temporal soil discoloration data. They fail to capture the dynamic "fingerprint" of soil discoloration changes over time, which is the central innovation of the present invention. Their approaches are therefore limited to identifying potentially visible markers, and do not provide a robust, automated, and temporally-sensitive method for detecting underground structures.

[0010] Menze and Ur propose a multitemporal fusion methodology aimed at detecting human-altered soils by classifying individual images and combining results to reduce temporal variability. The method is designed to emphasize static soil properties by suppressing temporal fluctuations.

[0011] This approach treats temporal variability as an artifact to be removed and therefore does not suggest, nor motivate, the use of temporal variability as an information-bearing signal. By design, it suppresses the transient spectral behavior that the present invention identifies as indicative of subsurface structures. Accordingly, this methodology teaches away from the present invention, which explicitly exploits temporal variation rather than eliminating it.

[0012] Kharate et al. describe the use of convolutional neural networks for archaeological site detection using multi-date or multi-band imagery. Such approaches rely on supervised or semi-supervised learning paradigms and require training datasets containing labeled examples of known archaeological features.

[0013] These methods are probabilistic and data-driven, with performance dependent on the availability, quality, and representativeness of training data. They do not disclose a deterministic algorithmic transformation that operates directly on spectral time-series data to derive a quantitative metric without prior training. In contrast, the present invention is unsupervised and deterministic, requiring no labeled data and no prior knowledge ofsubsurface structures. It derives information solely from the temporal behavior of spectral values at and around each spatial sampling locus, enabling deployment in previously undocumented regions.

[0014] Other methodologies for determining subsurface structures, are based on semi-automatic systems that can process large volumes of geospatial data to identify potential archaeological sites (Orengo, 2021, Winton, 2021). However, these methodologies also face challenges mainly related to the quality of input data and the representativeness of the machine learning training data. The effectiveness of machine learning models depends heavily on the analysis and quality of the satellite images used. Low-resolution or poorquality images can lead to misclassifications. Additionally, environmental conditions, such as vegetation, seasonal changes, and human activities, can affect the results, adding noise and reducing model accuracy. Finally, the complexity and lack of transparency of machine learning algorithms can make it difficult for interested parties to understand and trust the results.

[0015] Furthermore, the computational demands of deep learning methods like those used by Kharate et al. render them less accessible and scalable compared to the algorithmic approach of the present invention, which is specifically designed for efficient and robust detection using readily available resources. Moreover, the complexity of implementing and interpreting deep learning models often requires specialized expertise, making these methods less accessible to non-expert users in fields like archaeology or local infrastructure management. In contrast, the present invention provides a more user-friendly and readily deployable solution.

[0016] CN100492050C describes a method for reflectance-based inversion interpretation of geophysical formations using ground reflectance distribution images. The method is directed to static reflectance analysis and does not disclose a time-series workflow for evaluating temporal soil-surface spectral behavior across multiple acquisition dates. It therefore does not teach the normalization and metric-derivation steps central to the present invention.In contrast to the limitations of the prior art, the present invention introduces a comprehensive time series approach that specifically delivers the following technical breakthroughs:

[0017] A. systematically acquires and aligns time-series optical images and treats temporal spectral variation as an information-bearing signal rather than noise;

[0018] B. standardizes pixel-level spectral information into single-channel numeric values suitable for temporal analysis;

[0019] C. applies a spatial sampling framework that supports multiple spatial representations while enabling quantitative tracking of temporal behavior;

[0020] D. computes temporal baseline statistics and local neighborhood statistics to suppress background variability and transient environmental effects; and

[0021] E. derives a quantitative change-emphasis metric from normalized temporal behavior that correlates with the likelihood of underlying subsurface structures.

[0022] None of the cited prior art discloses or suggests a method that is simultaneously deterministic, unsupervised, independent of training data, and configured to derive a quantitative metric from the temporal spectral behavior of the soil surface. The present invention therefore represents a distinct and non-obvious technical contribution to subsurface detection using remote sensing imagery.

[0023] Objective Technical Problem and Inventive Step Starting from prior-art methods that analyze static soil properties, vegetation indices, or suppress temporal variability to identify persistent surface patterns, the objective technical problem addressed by the present invention is to provide a training-independent and rulebased remote-sensing method capable of identifying subsurface structures by deriving a quantitative indicator from temporal spectral behavior of the soil surface, without reliance on training data or prior site-specific knowledge. None of the cited prior art suggests treating temporal soil-discoloration variability as an information-bearing signal or deriving a quantitative change-emphasis metric from standardized spectral time-series data using temporal baseline and local neighborhood statistics. On the contrary, approaches such as multitemporal fusion explicitly suppress temporal variation, while machine-learning-based methods rely on supervised training and probabilistic inference. Accordingly, the claimedsolution is not obvious to the skilled person starting from the cited prior art and involves an inventive step within the meaning of Article 56 EPC.

[0024] Disclosure of the Invention

[0025] The method generates a georeferenced output dataset that aggregates information from different acquisition dates and enables identification of surface traces that may not be detectable from analysis of a single image or sporadically acquired images. The processing reveals patterns and differences on the ground surface that may be related to subsurface constructions and that would otherwise go unnoticed. In some implementations, the output dataset is rendered as a synthetic map to support interpretation and prioritization of locations for further on-site investigation. This approach increases the probability of detecting subsurface structures and provides a computationally efficient tool for archaeology, geology, critical infrastructure management, and related fields.

[0026] Below is a detailed presentation of the stages of the methodology and the procedures applied in each one:

[0027] A. Image Collection:

[0028] Initially, a plurality of optical images of the area of interest corresponding to a selected time period are obtained. The time period may extend over multiple acquisition dates. The images meet the following conditions:

[0029] • The images are georeferenced in a common coordinate system and cover substantially the same area. Spatial misalignment due to platform motion, camera tilt, acquisition angle, and / or topographic effects is corrected, or the images are obtained in a form in which such alignment has been performed, so that pixel-level comparisons across dates are supported.

[0030] • Where applicable, the images are subjected to atmospheric correction, or the images are obtained in a form that has already been atmospherically corrected, to reduce radiometric distortions due to varying atmospheric conditions.

[0031] Each image corresponds to a distinct acquisition date.

[0032] B. Determination of Basic Coordinates:Next, the coordinates of the specific area for which the analysis is requested (basic coordinates) are determined.

[0033] C. Creation of a Grid of Points:

[0034] Based on the previously defined basic coordinates, the invention defines a plurality of spatial sampling loci over the area of interest. Each spatial sampling locus corresponds to a spatial location or spatial unit from which pixel-level color information is extracted from the georeferenced optical images for subsequent algorithmic processing.

[0035] The spatial sampling loci may be defined using any suitable spatial sampling or discretization scheme, and the invention is not limited to a particular geometric representation or data structure, provided that each sampling locus can be associated with pixel values from the time-series images. Depending on the implementation, the spatial sampling loci may correspond, for example, to individual pixels in a full-resolution raster, to regularly or irregularly spaced sampling points, to randomly or pseudo-randomly selected locations, or to region-based spatial units derived from image segmentation or similar techniques.

[0036] In one preferred implementation, the spatial sampling loci are realized as a structured grid of sampling points, which provides a systematic and computationally efficient framework for extracting and organizing pixel-level information. This preferred implementation is described in detail below. Definition of basic coordinates: A set of basic coordinates is defined, which consists of reference points from which the additional grid points will be generated.

[0037] Preferred Vector-Grid Implementation

[0038] In this preferred implementation, a grid of multiple points is generated, which forms the sampling points from which pixel information will be extracted. The data from these pixels undergo algorithmic processing in a subsequent step. The process of creating the grid is described as follows:

[0039] 1. Definition of basic coordinates:

[0040] A set of basic coordinates is defined, which consists of reference points from which the additional sampling points will be generated.

[0041] 2. Preparation for grid creation:For each basic coordinate, a direction vector is calculated to define the path along which new points will be created. This vector is generated after considering the spatial relationship between successive reference points, ensuring that the new points are produced along a direction consistent with the arrangement of the basic coordinates.

[0042] 3. Generation of sampling points and determination of point density:

[0043] New sampling points are created along the predefined paths at predetermined spatial intervals. The step size of these intervals controls the distance between successive sampling points and thus the spatial resolution of the analysis.

[0044] 4. Determination of density of additional points around the sampling points (enhancement points):

[0045] Additional sampling points are generated within a predetermined maximum distance from each sampling point (Figure 1.A), thereby increasing local sampling density and sensitivity to spatially confined subsurface anomalies.

[0046] 5. Organization, classification, and storage of sampling points:

[0047] As new sampling points are generated, they are stored in a structured format, for example in arrays or lists classified according to one or more coordinate values (such as the y-coordinate in ascending or descending order), to facilitate efficient data handling and subsequent processing. 6. Structuring the output format of the sampling points:

[0048] The set of sampling points generated in this manner is compiled into a complete dataset, forming the sampling-point output file (Figure 1.B and Figure 1.C), which is used in the subsequent pixel-sampling and temporalanalysis stages of the method.

[0049] D. Sampling Pixel Values on Georeferenced Images and Spectral Dimensionality Reduction

[0050] Based on the spatial sampling loci defined in Step C, the pixel sampling and spectral processing procedure begins, as described below. The objective of this step is to extract pixel-level spectral information from the time-series images and to transform that information into one or more standardized numeric values suitable for subsequent temporal analysis.

[0051] 1. Accessing sampling loci:The spatial sampling loci defined in Step C are accessed. Each sampling locus is checked for validity, and, in cases of incorrect formatting or out-of-bounds coordinates, the sampling locus is rejected and an error is recorded.

[0052] 2. Extraction of spectral values:

[0053] For each sampling locus, spectral values are extracted from each georeferenced optical image obtained in Step A. The extracted spectral values may correspond to one or more spectral bands available in the image data, including, but not limited to, visible bands (e.g., red, green, blue), near-infrared (NIR), short-wave infrared (SWIR), or other spectral channels provided by the imaging sensor. The pixel position within each image is determined based on the image geometry and the spatial location of the sampling locus.

[0054] 3. Recording spectral values in a color information database:

[0055] The extracted spectral values are stored in a database (color information database) in association with the corresponding sampling locus and image acquisition date. The stored spectral values are then transformed into one or more standardized single-channel numeric values by applying a spectral dimensionality-reduction operation, as described below.

[0056] 3.1. Spectral dimensionality reduction using weighted band combinations: In this step, the multi-band spectral information associated with each sampling locus is reduced to a single scalar value by applying a weighted combination of spectral components. In a general form, the standardized numeric value V is computed as:

[0057] V= Z (i=1 - n) wi Ci

[0058] where (^represents the value of the i-th spectral band at the sampling locus, wTepresents a corresponding weighting coefficient, and nrepresents the number of spectral bands used in the combination. The weighting coefficients may be selected to emphasize soil-related reflectance characteristics, suppress spectral noise, or optimize sensitivity to surface variations associated with subsurface anomalies.The selection of spectral bands and weighting coefficients may depend on the imaging sensor, environmental conditions, and the application context, and may be predetermined or user-defined.

[0059] 3.2. Exemplary grayscale embodiment:

[0060] In one preferred embodiment, the standardized numeric value V is a grayscale value derived from red (R), green (G), and blue (B) spectral components using predetermined weighting coefficients. An exemplary, non-limiting formula is:

[0061] V=0.299R+0.587G+0.114B (Eq. 1)

[0062] While the examples herein utilize grayscale values to demonstrate the inventive principle of statistical temporal normalization, it is understood that other single-channel spectral metrics exhibiting similar statistical distributions may be employed. In this embodiment, the resulting grayscale numeric value is a standardized single-channel scalar derived from the pixel color components and used as a numerical basis for subsequent temporal analysis. Values approaching the upper end of the grayscale range correspond to brighter surface traces, while lower values correspond to darker surface areas. Other weighting schemes or spectral band combinations may be used without departing from the scope of the invention.

[0063] 3.3. Rounding and standardization:

[0064] The resulting standardized numeric values may be rounded to a predetermined number of decimal places, such as three decimal places, to standardize data format and facilitate subsequent analysis.

[0065] 4. Recording spectral values in a color information database:

[0066] The standardized numeric values are organized and stored in the color information database for each processed image, such that values are indexed by sampling locus and acquisition date. The data corresponding to each sampling locus are aggregated to form a time-series dataset containing the standardized numeric values for the entire set of images acquired over the selected time period, enabling efficient access during the subsequent normalization and temporal-analysis stages.E. Generation of an Output Spatial Dataset and Synthetic Map:

[0067] The standardized single-channel numeric values produced in Step D are processed to compute normalization metrics and to generate an output spatial dataset associated with the spatial sampling loci defined in Step C. In some implementations, the output spatial dataset is visualized as a synthetic map to facilitate interpretation and prioritization of candidate locations for further investigation.

[0068] 1. Data normalization:

[0069] The normalization process is configured to mitigate natural spatial heterogeneity and transient environmental effects while amplifying persistent variations in the standardized single-channel values that may be correlated with subsurface anomalies. The normalization includes the following steps.

[0070] • Calculation of a temporal baseline statistic:

[0071] For each spatial sampling locus, a temporal baseline statistic is computed across the entire set of images corresponding to the selected time period. In one preferred embodiment, the temporal baseline statistic is a cumulative mean of the standardized single-channel values forthat sampling locus.

[0072] • Calculation of a local neighborhood statistic:

[0073] For each spatial sampling locus and for each image acquisition date, a local neighborhood statistic is computed using additional sampling loci located within a predetermined spatial neighborhood around the sampling locus. In one preferred embodiment, the local neighborhood statistic is a mean value computed from enhancement points surrounding the sampling locus.

[0074] • Generation of first-order normalized data:

[0075] Normalized data are generated by comparing the local neighborhood statistic to the temporal baseline statistic. In one preferred embodiment, this comparison is performed by calculating a ratio between the local neighborhood mean and the cumulative temporal mean, thereby producing first-order normalized values.

[0076] • Storage of normalized data:The normalized values are stored for each spatial sampling locus and acquisition date, forming a normalized time-series dataset suitable for subsequent analysis and visualization

[0077] 2. Derivation of a Change-Emphasis Metric

[0078] To emphasize persistent spatial variations in the normalized data, a change-emphasis metric is derived from the normalized values.

[0079] • Determination of a maximum normalized value:

[0080] For each spatial sampling locus, the maximum value among the first-order normalized values within the associated local neighborhood is identified and denoted as Hmax,i

[0081] • Calculation of deviation measures:

[0082] For each normalized value H within the local neighborhood, a first deviation measure is calculated as: (AHmax,i= H - Hmax,i). A second deviation measure is calculated using the minimum normalized value / / min lwithin the same neighborhood as: (AHmax,2= Hmax,i -Hmjn,l).

[0083] • Calculation of second-order normalized difference:

[0084] A second-order normalized difference metric A / / "is computed as:

[0085] AH" = AH max,1 / AH max, 2 / Eq. (2)

[0086] The value of H "expresses the relative deviation of a normalized value with respect to the local maximum and is bounded between 0 and 1.

[0087] 3. Generation of Visualization Values

[0088] In implementations where visualization is desired, the change-emphasis metric may be transformed into visualization parameters.

[0089] Exemplary gray-level mapping:

[0090] In one preferred, non-limiting embodiment, a visualization intensity value T is computed using an exponential mapping:

[0091] Tgray = lnt(255*exp(-k*AH”) Eq. (3)Where:

[0092] • Int: denotes extraction of the integer part,

[0093] • 255: represents the maximum grayscale value,

[0094] • e: is the base of the natural logarithm, and

[0095] • k: is a positive scaling coefficient controlling contrast enhancement.

[0096] The parameter values may be adjusted or optimized depending on application requirements, sensor characteristics, or visualization preferences.

[0097] Adjustment of marker attributes: In visualization embodiments, one or more graphical attributes associated with each spatial sampling locus — such as brightness, color intensity, or marker size — may be adjusted according to the derived visualization value to enhance the visual prominence of locations exhibiting higher change-emphasis metrics.

[0098] 4. Interpretation and Evaluation

[0099] The resulting output spatial dataset and / or synthetic map may be interpreted by an operator or processed automatically to identify and prioritize candidate locations indicative of subsurface anomalies for further on-site investigation.

[0100] Computer Implementation:

[0101] The method steps described herein are computer-implemented and may be executed by a computing system configured to obtain the georeferenced optical images, store intermediate values (including the color information database and normalized metrics), and generate the output spatial dataset and / or synthetic map for display, printing, storage, and / or transmission. In some implementations, the computing system executes program instructions stored in memory, and such instructions may be provided on a non-transitory computer-readable medium.DESCRIPTION OF OPERATION EXAMPLES

[0102] The operation of the invention is described below for the following case examples:

[0103] A. Case of Known Infrastructure Location and Targeted Detection

[0104] This case involves the use of the invention for analyzing the impact of a known underground infrastructure, such as a water pipeline, on the surface of the ground.

[0105] B. Case of Hypothetical Presence of Underground Structure and Targeted Detection

[0106] In this case, there is a hypothesis regarding the presence of an underground structure, such as the western walls of ancient Pella. The invention is used in a targeted manner to confirm the existence of the structure and to determine its location and extent.

[0107] C. Case of Unknown Target Search and Scanning of an Entire Area

[0108] This case pertains to the detection of a wide area without specific information about what is being sought, as in the area south of ancient Pella. The invention is used to scan the area and identify any surface traces that may indicate the possible existence of archaeological interest or other underground infrastructures.

[0109] A. Example of Detecting Ground Traces for a Known Infrastructure Location

[0110] In this case, the methodology was applied to record the traces observed on the ground surface under which a water supply pipeline passes, providing drinking water to a large urban center. As shown in Figure 2, after collecting satellite images (Figure 2.A) in the near-infrared spectrum, a grid of points was created based on which the sampling was conducted (Figure 2.B). Specifically, the process includes the following stages:

[0111] 1. Collection of Satellite Images: Images in the near-infrared spectrum were collected at regular intervals, ensuring coverage of the same area with georeferencing in the same coordinate system.

[0112] 2. Creation of Grid of Points: Based on the collected images, a grid of points covering the entire area of interest was created.3. Sampling of Points: The RGB values of the pixels from the satellite images were converted to grayscale values, allowing for data analysis.

[0113] 4. Data Analysis: The data were analyzed using the method proposed in the present invention, producing a synthetic map that depicted the traces of the pipeline (Figure 2.C).

[0114] On-site inspections have confirmed that the points the method displays as brighter on the synthetic map indeed show significant color deviations. These observations support the accuracy of the method and its ability to detect color anomalies on the ground surface through the analysis of remote sensing data. These deviations are likely related to changes in surface moisture or other environmental factors influenced by the interaction of the underground pipeline with the subsoil. Identifying these points where deviations are detected suggests further on-site investigation with test trenches at the locations indicated by this method.

[0115] B. Example of Detecting Hypothetical Presence of Underground Structure and C. Example of Scanning an Entire Area

[0116] (B) The method was also tested in areas of archaeological interest, prompted by the passage of the water transfer pipeline from the Aravissos springs to the city of Thessaloniki through the archaeological site of Ancient Pella. Specifically, traces were detected on the ground surface north of the pipeline that refer to the western walls of Ancient Pella, and south of the pipeline that refer to the ancient islet of Phakos. As shown in Figure 3 and 4, after collecting satellite images in the near-infrared spectrum (Figure 3. A and 4. A), a grid of points was created (Figure 3.B and 4.B) based on which the sampling was conducted. The process includes the following stages:

[0117] 1. Collection of Satellite Images: Images in the near-infrared spectrum were collected at regular intervals, ensuring coverage of the same area with georeferencing in the same coordinate system.

[0118] 2. Creation of Grid of Points: Based on the collected images, a grid of points covering the entire area of interest was created.3. Sampling of Points: The RGB values of the pixels from the satellite images were converted to grayscale values, allowing for data analysis.

[0119] 4. Data Analysis: The data were analyzed using the method proposed in the present invention, producing a synthetic map that depicted the traces of the pipeline (Figure 3.0 and 4.C).

[0120] Regarding the area of the western walls of Ancient Pella, the detected traces are attributed to these walls because the trace under examination follows the direction of a section of the walls that has already been revealed in the past by the competent archaeological service.

[0121] (C) Of particular interest is the case of the Phakos islet (Figure 4), where traces were detected that correspond to bibliographic sources hypothesizing the existence of an islet with geometric characteristics similar to the found traces detected by the proposed methodology, as well as a large circular construction that refers to the fort described by ancient authors of the islet. The results from applying the method at this specific location show that the invention is capable of enhancing the depiction of the traces on the ground surface as they are rendered on the synthetic map, to such an extent that traces are revealed which do not appear through simple observation of any available satellite images. The results of the analysis revealed traces that agree with historical and archaeological sources, enhancing the reliability and accuracy of the method. These detections confirm the method's reliability and its ability to contribute to a better understanding of the topography and ancient underground structures in areas of archaeological interest, offering valuable data for further research and on-site study.

[0122] Independent archaeological observations conducted in the examined areas were found to be consistent with known archaeological and historical information, supporting the technical plausibility of the detected surface traces.Y1

[0123] D RAWI N G S CA TA L O G U E

[0124] Figure 1: Generation of Sampling Point Grid

[0125] Figure 1.A: Visualization of the sampling area around each grid point (enhancement points), with a specified maximum distance.

[0126] Figure 1.B: Example of an output file with the sampling points organized and sorted. Figure 1.C: Visualization of the complete grid of points to be used for pixel sampling. Figure 2: Trace Detection for Known Underground Infrastructure (Water Supply Pipeline) Figure 2.A: Area of interest for detecting the water supply pipeline.

[0127] Figure 2.B: Grid of sampling points for the pipeline area.

[0128] Figure 2.C: Synthetic map depicting the traces of the pipeline on the ground surface. Figure 3: Detection of a Hypothetical Underground Structure (Western Walls of Ancient Pella)

[0129] Figure 3.A: Area of interest for detecting the western walls of Ancient Pella.

[0130] Figure 3.B: Sampling point grid for the area of the western walls.

[0131] Figure 3.C: Synthetic map depicting the traces of the western walls.

[0132] Figure 4: Scanning of the Entire Area (Fakos Islet)

[0133] Figure 4.A: Area of interest for scanning the Fakos Islet.

[0134] Figure 4.B: Sampling point grid for the Fakos Islet area.

[0135] Figure 4.C: Synthetic map depicting traces of archaeological interest, including the circular structure.

Claims

ClaimsIndependent Claim1. A computer-implemented method for detecting patterns on a ground surface indicative of subsurface anomalies using time-series optical remote sensing image data, the method being configured to analyze temporal variations in soil discoloration as an indicator of the presence of said subsurface anomalies, and comprising the steps of:(a) Acquiring Multi-Temporal Images: obtaining a plurality of georeferenced optical images of a predefined area of interest over a selected time period, each image corresponding to a distinct acquisition date and having a spectral sensitivity and spatial resolution sufficient for capturing subtle variations in soil color;(b) Aligning and Correcting Images: processing the optical images, or using optical images that have been processed, to ensure pixel-level spatial congruence across the time series and to reduce radiometric and geometric distortions, wherein said processing comprises registering the optical images and applying, or using optical images that have been subjected to, atmospheric and geometric corrections, thereby providing a spatially aligned and radiometrically consistent image dataset suitable for comparison of soil color across multiple dates;(c) Defining Spatial Sampling Loci: defining a plurality of spatial sampling loci distributed across the area of interest according to a spatial sampling or discretization scheme, and further generating additional local reference sampling points (enhancement points) located within a predetermined proximity to each sampling locus, the spatial sampling loci being configured to facilitate systematic pixel-level extraction of soil spectral data;(d) Extracting and Standardizing Pixel Spectral Values: for each spatial sampling locus, automatically extracting pixel spectral values from one or more spectral bands including visible and / or infrared bands of each aligned and corrected image in the time series at the spatial sampling locus and at the enhancement points, and converting the extracted pixel spectral values into one or more standardized single-channel numeric values (spectral metrics) configured to reduce the influence of natural spectral variability and enhance the comparability of subtle soil color changes over time;(e) Temporal Trend Analysis of Soil Discoloration: analyzing the sequence of standardized single-channel values at each spatial sampling locus across the entire time series to identify and quantify persistent value alterations, said analysis distinguishing short-term environmental effects from enduring soil discoloration changes indicative of potential subsurface anomalies;(f) Multi-Stage Normalization: performing a multi-stage computational normalization process to mitigate natural soil heterogeneity and transient environmental effects while amplifying persistent value changes, the process comprising:(i) computing a temporal baseline statistic across the time series for each spatial sampling locus;(ii) computing a local neighborhood statistic for the enhancement points on a per-image basis;(iii) generating normalized data by comparing the local neighborhood statistic to the temporal baseline statistic; and(iv) deriving a change-emphasis metric from the normalized data to emphasize soil discoloration changes correlated with possible subsurface anomalies; and(g) Generating an Output Spatial Dataset: generating an output spatial dataset based on the change-emphasis metric that associates the spatial sampling loci with the emphasized soil discoloration changes so as to identify candidate locations indicative of subsurface anomalies for guiding further on-site investigation.Dependent Claims2. The method of claim 1, wherein obtaining the plurality of georeferenced optical images comprises acquiring, via a remote sensing system, the plurality of georeferenced optical images.

3. The method of claim 2, wherein the remote sensing system for acquiring optical images comprises a source selected from the group consisting of satellites, unmanned aerial vehicles, and other aerial platforms, wherein each platform provides multi-temporal optical data at intervals sufficient to capture transient soil color variations.

4. The method of claim 1, wherein, in one implementation, defining the spatial sampling loci comprises:(a) defining a set of basic coordinates that serve as reference points within the area of interest;(b) algorithmically creating additional grid points at predetermined spatial intervals along directions determined from the basic coordinates; and(c) generating said enhancement points around each grid point within a predetermined distance, thereby increasing local sampling density and detection sensitivity for spatially confined subsurface anomalies.

5. The method of claim 1, wherein applying, or using optical images that have been subjected to, atmospheric and geometric corrections comprises:(a) applying geometric corrections to compensate for movements of an image acquisition means, acquisition angle, and topographic variations; and(b) applying atmospheric corrections to reduce radiometric distortions due to varying atmospheric conditions, thus ensuring reliable pixel-level color comparison throughout the time series.

6. The method of claim 1, wherein the standardized single-channel numeric values comprise grayscale values obtained by converting the extracted pixel spectral values through a mathematical transformation of the spectral components.

7. The method of claim 6, wherein converting pixel spectral values to standardized grayscale values comprises applying a pre-selected weighting function to each spectral channel, said weighting function being chosen to emphasize soil-related reflectance properties and to provide a consistent numerical basis for subsequent temporal analysis.

8. The method of claim 6, wherein, in one embodiment, the multi-stage computational normalization process comprises:(a) calculating cumulative mean grayscale values across the time series for each spatial sampling locus as the temporal baseline statistic;(b) determining local mean grayscale values for the enhancement points on a perimage basis as the local neighborhood statistic;(c) generating first-order normalized data by computing ratios of the local mean grayscale values to the cumulative mean grayscale values; and(d) deriving a second-order difference metric from the first-order normalized data as the change-emphasis metric.

9. The method of claim 1 , wherein the multi-stage normalization process further includes applying a data-filtering or smoothing technique configured to suppress noise arising from sporadic image artifacts or environmental disturbances, and to preserve meaningful temporal patterns linked to enduring subsurface anomalies.

10. The method of claim 1, wherein the output spatial dataset is optionally rendered as a synthetic map that visually represents the spatial sampling loci in proportion to the changeemphasis metric, such that locations exhibiting prominent or enduring value variations are visually highlighted.

11. The method of claim 1, wherein the method is unsupervised and does not require training data, labeled examples, or prior knowledge of subsurface structures and wherein the change-emphasis metric is derived exclusively from the statistical properties of the obtained plurality of images without reference to external training datasets.

12. The method of claim 1, wherein the method is deterministic and produces identical outputs for identical input datasets under identical processing parameters.

13. The method of claim 10, wherein the synthetic map visually adjusts one or more of the following attributes for each spatial sampling locus in proportion to the change-emphasis metric: brightness, color intensity, and marker size.

14. The method of claim 10, wherein the synthetic map can be displayed on a computer monitor or printed output and is configured to overlay additional geographic or archaeological information and facilitate manual or semi-automatic interpretation of hotspots with elevated likelihood of subsurface anomalies.

15. The method of claim 1, wherein defining the spatial sampling loci is performed according to claim 4, converting the pixel spectral values comprises generating grayscale values according to claim 6 and applying a pre-selected weighting function according to claim 7, the multi-stage normalization comprises the steps of claim 8, and generating the output spatial dataset comprises rendering the output spatial dataset as a synthetic map that visually represents the spatial sampling loci in proportion to the change-emphasis metric.

16. A computer-implemented system comprising one or more processors and a memory storing instructions which, when executed by the one or more processors, cause the system to perform the method of any one of claims 1 to 15.

17. A non-transitory computer-readable medium storing instructions which, when executed by one or more processors, cause performance of the method of any one of claims 1 to 15.

Images