{'fraction_non_alphanumeric': 0.029609644479355213, 'fraction_numerical': 0.00495201666609747, 'mean_word_length': 4.634404464113912, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust Oil-spill Forensics and Petroleum Source Differentiation using Quantized Peak Topography Maps

Identification and classification of environmental forensics, with the petroleum forensics as the main application, requires an effective technology or method to distinguish between the closely located forensics as they share many main biomarkers. Two-dimensional gas chromatography is one of these technologies with which a petroleum forensic is separated into its chemical compounds, resulting in a three-dimensional image, GCXGC image. Therefore, distinguishing between two petroleum forensics is equivalent to the comparison between their corresponding GCXGC images. In this paper, we present a technique, called Quantized Peak Topography Map (QPTM), which results in a better separation between the GCXGC images. We validate our proposed method on a model dataset, consisting of thirtyfour GCXGC images, extracted from the different parts of the world.


I. INTRODUCTION
In this work, we introduce a robust signal processing approach for petroleum fingerprinting, in particular, source forensics in the aftermath of a major oil spill such as the Deepwater Horizon disaster, Gulf of Mexico, April 2010. Petroleum fingerprinting has played a crucial role in apportioning the environmental impact of major oil spills and industrial leaks. Furthermore, robust source differentiation between closely spaced oil reservoirs enable determination of underground connectivity of petroleum reservoirs. Such sourcespecific knowledge can significantly reduce the environmental, financial, and safety risks of multiple drilling points in offshore rigs. Therefore, robust source fingerprinting for closely spaced petroleum sources is of paramount importance to environmental health, public safety and resource management for the petroleum industry. Petroleum identification based upon a pixel-driven compression technique was first introduced in [1] where a representation of the GC×GC Chromatography image called ρ τ -map was applied for the classification of different GC×GC images using their compressed versions. In [1] the compressed version of a GC×GC image of a petroleum source from one specific geographical region is achieved by identifying the portion of the image which is common among a bunch of sources from that region. In this scenario, in case there is a noise in one the sources, it could affect the whole representative compressed version of the region. In this paper a compressed version of each GC×GC image is achieved using Symbolic Aggregate ApproXimation (SAX) method and just based upon its own image and hence a robust technique is presented.
We present here a signal processing perspective to this challenging environmental application that combines recent research in peak topography mapping with well-known techniques in high-volume data quantization. The key innovation developed in this article is (Q)uantized (P)eak (T)opography (M)ap: a scaled quantized topography mapping technique that achieves two important signal processing objectives: (i) Peak Topography Mapping: We significantly extend the seed idea of peak topography maps, introduced in [2], [3] to include a wider range of peaks, algorithmically isolated from large-scale chromatographic data. Novel extensions include disentangling low signal-to-noise (SNR) ratio peaks otherwise lost against baseline noise and other perturbations typical in chromatographic instruments. (ii) Quantization of peak information based on prioritization of peaks based on their local SNR and relevance to source signature within the joint peak profile of petroleum hydrocarbons. This step is critical in robust pattern recognition of petroleum source signatures, which may span hundreds of well-known (target) and unknown (nontarget) biomarkers.
In synopsis, the QPTM technique proposed here enables efficient source differentiation of petroleum sources, with broader applicability to a wide range of peak profile data, typically generated in gas chromatographic and mass spectrometric instruments. Thus, from a signal processing perspective, the proposed technique is applicable far beyond the application focus, i.e., petroleum forensics, to most chemical interpretation of instrumental raw signals. Typical applications that may benefit from this signal processing method include, but are not limited to, air pollutant studies (.g. fingerprinting key pollutants sampled in urban air), water sustainability investigations (e.g. isolating key toxins in ground and surface water and fingerprinting them to local industrial agents such as a factory or common household cleaning products), among many others.
We have chosen to focus on petroleum forensics as the target application in this work due to its well-known relevance to arXiv:1807.07484v1 [eess.SP] 19 Jul 2018 environmental studies, public health and national economy. The rest of this paper is organized as follows. Section II provides background and related work in signal processing and high-volume data informatics. Sections III and IV describe the mathematical formulation of the problem and the proposed approach. Sections VI and VI-B present the experimental results and related discussions. Conclusions are presented in Section VII.

II. BACKGROUND AND RELATED WORK
A common analytical approach to petroleum forensics is to separate an oil sample into an information-rich distribution of constituent compounds, which span across an intricate network of hundreds, sometimes thousands, of hydrocarbons. In particular, from a pattern recognition perspective, most petroleum forensic techniques focus on the hydrocarbon "biomarker" groups, hopanes and steranes [4], which are well-known to withstand environmental degradation. In practical terms, this means that oil sampled from grass blades in Louisiana coastline share very similar biomarker distribution with oil sampled from the Macondo well, source of the Deepwater Horizon spill. Therefore, robust source fingerprinting of an oil sampled hundreds of miles away from a spill site reduces to accurate and efficient matching of the biomarker peak distribution between sampled crude oil and oil taken from the spill site.
A. Two-dimensional Gas Chromatography (GC × GC) Two-dimensional Gas Chromatography (GC × GC) is a technology used to separate any petroleum forensic into its constituent components. A petroleum source is injected into a GC×GC system which consists of two stages (columns) where each of the components interact differently with these columns and will pass these columns in different times (retention time). Note that, we have two stages for a better resolution as those components which have not been separated in the first stage can get separated in the second stage. Hence, for each component we will have two retention times corresponding to the time they passed the first and the second columns and also its concentration. The resulting pattern would be a threedimensional image where the first and second dimensions are the first and second retention times and the third dimension is the concentration of the corresponding component or its peak point ( Figure 1).

B. Mathematical model of biomarker distribution and related challenges
Section II-A provides a basic description of twodimensional gas chromatography (GC×GC), the analytical technology that performs these separation of complex mixtures such as crude oil into constituent compounds. The outcome of GC×GC separations are high-resolution images of hydrocarbon peaks. Mathematically, the peak distribution of petroleum hydrocarbons in a crude oil sample can be modeled as: where I(x, y) denotes the pixel intensity at the pixel location (x, y) of the GC×GC image, (x i , y j ) denotes the peak location for the j th peak along the i th column of the GC×GC image, and p ij denotes the peak height and ψ ij (x − x i , y − y j ) denotes the peak concentration distribution function for this peak centered at (x i , y j ). From a modeling perspective, the peak concentration distribution ψ ij (x − x i , y − y j ) roughly follows a Gaussian or bell-shaped distribution, with peak summit of p ij located at (x i , y j ). However, the exact model of each peak profile and the joint peak distribution in Equation (1) is challenging, sometimes impossible, to determine effectively for three related and compelling reasons: (i) Experimental variability introduces baseline uncertainty, which changes non-linearly across time, and thus, renders an effective model for ψ(·) difficult to frame across the entire time period of the experimental run; (ii) Co-eluting peaks, i.e., peaks with overlapping distribution functions are challenging to deconvolve; (iii) Precise modeling of peak profiles is impractical, if not impossible, across hundreds of hydrocarbons, where each intricate joint peak distribution can uniquely fingerprint thousands of closely related oil reservoirs in a petroleumrich area (Gulf of Mexico alone has ∼ 3000 active oil sources within geographic proximity of the Deepwater Horizon spill. Therefore, adopting a model-based approach to petroleum forensics is rife with statistical assumptions that are impossible to assert across a broad variety of test cases. This motivates to adopt a method that is not explicitly dependent on an underlying peak distribution model, but rather, exploit the rich diversity of the peak profile. We aim to prioritize on major peaks that dominate the peak topography, while accounting for the minor peaks that are individually less significant, but cumulatively contribute to the overall distribution. A recent attempt to achieve this using peak manifolds defined across clusters of GC×GC pixels is documented in [3] along with data compression implications in [5], [6]. However, these pixel-driven techniques did not achieve robust source differentiation due to disambiguation and dimensional challenges discussed below.

C. Computational and signal processing challenges to robust pattern recognition
Despite phenomenal advances in separation technology [7]- [11] the fundamental challenges in petroleum forensics are well-known source disambiguation issues in the signal processing and pattern recognition community. Specifically, the challenges may be iterated as follows.
(i) Closely spaced oil reservoirs share bulk of their biomarker peak distributions due to the common regional fingerprint, thus effectively creating "near-far" challenge for disambiguation of weaker source-specific fingerprint against the stronger regional fingerprint. (ii) The high dimensionality of the data limits the performance and scalability of data mining algorithms due to the "curse of dimensionality" [12]. The answer to this question has been sought in two complementary directions, target compound analysis [13]- [15] and statistical multivariate approaches [16]- [18]. While both approaches provide insight into what the source fingerprint of a petroleum reservoir might look like, neither offer a clear path to robust fingerprinting that is cognizant of the network of biomarker compounds that constitute crude oil. This is a critical and compelling gap to fill as closely spaced petroleum sources in an oil-rich locale, e.g. the Gulf of Mexico, often share regional commonalities, which translates to significant similarities in the biomarker distribution. From a pattern recognition perspective, This poses a significant problem where the regional fingerprint masks the source-specific fingerprint.

III. PROBLEM STATEMENT
Consider a source library D, with K GC × GC template images, D = {I 1 , I 2 , . . . , I K } where each image represents the unique GC × GC biomarker fingerprint of an oil source. The geographical regions of these GC × GC images are saved in a set R = {r 1 , r 2 , . . . , r K } where the i th element in R indicates the geographical region of i th image of D.
We assume K to be large enough that D has the requisite data diversity to robustly differentiate source fingerprints, despite significant overlap between neighboring sources. The petroleum fingerprinting problem then reduces to matching the GC×GC peak distribution Φ test (x, y) of a newly-extracted unknown sample, based on the observed GC×GC image I test , to the closest source template image I k ∈ D. Mathematically, the closest match I opt ∈ D for I test is determined using the similarity criterion S such that: The key contribution of this work is designing this similarity criterion such that it is robust to statistical uncertainties in the peak profile distribution. We define the similarity criterion S(·) in Section (IV-A) based on the relative peak distribution of the reference and test samples.

IV. TECHNICAL APPROACH
We present the solution to the problem stated in Section III with the framework shown in Figure 2. It is noteworthy that the computational framework in Figure 2, though similar in architecture to an artificial neural network (ANN) [19], it is not implemented as such, but rather employed to identify peak profile features for source fingerprinting in an unsupervised setup. The proposed framework consists of the three layers: (i) Input, (ii) Comparison and (iii) Decision.
We now provide some details on the different parts of the framework in Figure 2. The box on the left, has two parts, first the GC×GC image of the injection under the test (I test ) and a block for the index of the reference injections from the library D to be tested against I test , which this index is initialized to one (RefIndex=1). The test image is an M × N matrix, so in the input layer we have M × N inputs, shown as X 1 to X M ×N where each of these inputs indicates one pixel of the test image. These inputs are transmitted by the weights W τ to the Comparison layer. In Section IV-B we will discuss on the procedure to construct W τ .
In the Comparison layer the G function is implemented between I test and the reference image. Finally, in the Decision layer, the function F has the role to determine the geographical region of I test . F has K memories. First, the Reference Index is one, so I test is compared against the first element of D and the correspondingd Ref Index is saved in one memory of F . Note thatd Ref Index is a criterion of deviation of the test image from the reference image chosen from D. Once this comparison is done completely andd Ref Index is saved, F sends an "Increment Reference Index" command to the leftmost box and in the next time the test image should be compared against the second element of D and so on. Finally once Reference Index exceeds the size of the dictionary, K, F outputs the index having the lowestd Ref Index which is the index of the geographical region of I test from the set R.

A. Similarity Criterion
As discussed comprehensively in [3], [20] comparing two GC ×GC images is the comparison between their corresponding peaks at the same location of their image.
We sweep the image by going along each of its N columns and save the local maxima, detected along each of the columns as peaks. Suppose we are to compare two images, I ref and I test . At one location, like (x 0 , y 0 ) the amplitude of the peaks in I ref and I test are p ref and p test , respectively. In [3], [20] the similarity between these two peaks is defined as: As known, the function max(α, α −1 ) has a value greater or equal to one for all α ∈ R, α = 0. Therefore, if the function Sim(·) has a value of one, then α = α −1 or in our case p ref = p test . Any difference between the values of p ref and p test results in a value greater than one for Sim(·). We can call two peaks as similar once Sim(·) for their corresponding location is unity. However, to account for experimental variability during the process of producing the GC × GC image, we relax the definition of similarity to within an uncertainty bound . Hence, we define a peak-ratio parameter τ as: We claim two peaks as similar, once the Sim(·) function for their corresponding location is less than or equal to τ (Sim(·) ≤ τ ). The parameter indicates the amount of deviation that is acceptable for us to consider two peaks as similar; the more the value of is, the less strict we are in the definition of similarity between the two peaks. In an extreme case, if we set = ∞, all of the peaks will be considered as similar because Sim(·) is always less than ∞.
B. How to set W τ There are M × N inputs, X 1 to X M ×N , but as discussed in the previous section, just the local maxima of the image are saved and used as the comparison. Therefor, the weights for the non-peak pixels or inputs are set to zero. For calculating the weight of one peak like p test at the location (x 0 , y 0 ), once the peak in the same location in the reference image is p ref , while the peak-ratio is τ , we perform the following two-stage process: 1. Calculate the Sim(·) for the location of the peaks as in Equation 3. 2. a. If Sim(·) ≤ τ then replace both p ref and p test by a common peak, p common : so the weight of p test will be: b. Otherwise, keep the exact values of p ref and p test which means w test (x 0 , y 0 ) = 1. The reason which we replace p ref and p test by p common is that, we have assumed once Sim(·) is less than equal to τ for that location, these two peaks are similar, and therefore we manually assign an equal value to them. We use the minimum value as the common peak for both of the reference and test images, because the minimum peak amplitude is common between the two peaks.

C. How to set
Suppose we have the library R = {I k1 1 , I k2 2 , . . . , I k N K }, where the i th element of the library, I ki i , or i th family, means we have k i number of GC ×GC images for the region indexed by i. For each of these image families, we also learn an (i) with some training members of the family. For each of these training set injections, let's say I train , we first construct its maximum-along-Interval representation as following: 1. We pass each column of I train to the proposed network and apply its corresponding W τ with respect to a reference image of the family. 2. For each column of I train , named as I c train , which is a time series, we do the following: In other words, for each column we have n w intervals, and we replace each interval with its maximum value. Therefore, we have constructed the maximum-along-Interval representation of I train . Once, we have done the above-mentioned two stages for all of the images in the training set, we compute the average L 2 norm distance of the all of the images in the training set. This average L 2 norm distance is equivalent to the deviation of the training images after applying the proposed two stage method. Then we sweep over an increasing range of 's and compute its corresponding deviation. It is noteworthy that the amount of the deviation is disproportional to the value of , i.e, as we increase , the deviation among the images decrease. This is because, as denoted in equation 4, the more the value of , the more similar two arbitrary peaks are to each other and the less strict we are in our definition of similarity.
This can be seen in Figure 3, where the deviation is decreasing by increasing , as we are to set the for the family of images from Macondo well. In this figure, we have set the first seven Macondo injections as the training set for calculating the . But we wish to set the value of not too far from zero, since, at least in an extreme case, setting = ∞ would indicate any two peaks with different values as being similar. So we wish to choose a good choice for where the deviation is low and it is in its minimum possible value. As can be seen in Figure 4, the deviation drops until = 0.5 but remains the same for higher values. So the best choice for would be half.

V. FUNCTION G
As discussed earlier, in our proposed framework, the Comparison layer which is implemented by function G, computes the distance between the images under the test and the reference image using their SAX representations. In the following, we briefly touch upon the SAX method.

A. SAX algorithm
The comprehensive study on Symbolic Aggregate ApproXimation (SAX) algorithm can be found in [21]- [23]. Here, we touch upon it briefly, as we will use it in our proposed framework in Figure 2 as function G in the Comparison layer. Formally, SAX is a method to represent time series of size n using a string of arbitrary size w(w < n). Applying SAX, one can control the dimensionality reduction imperative when dealing with large time series. As formally defined in [21], a time series C = c 1 c 2 c 3 . . . c n of length n can be represented in a w-dimensional space by a vectorC =c 1c2c3 . . .c w . The i th element of C is calculated by the following equation: In other words, the time series is divided into w intervals of the same size, then the data in each of these intervals is replaced by the mean of the data. This representation is called the Piecewise Aggregate Approximation (PAA) representation of the time series. Now each of these intervals is mapped into one symbol. These samples are chosen out of some set with the cardinality of α. After applying the SAX method, the final symbol representation of the times series C will be: Note that C is the original times series,C is its PAA representation andĈ is its SAX representation.
B. Calculating the distance between two time series using their PAA and SAX representations Suppose we have two time series C 1 and C 2 . The Euclidean distance between the PAA representation of the two time series is :d Wherec i 1 andc i 2 are the i th element of theC 1 andC 2 , respectively. After constructing the symbol representations of C 1 and C 2 , we wish to calculate the symbol distance of their representations. Therefore, we need a look-up table in order to have the distance between the symbols. Such a table is given in Table I. Then, The Euclidean distance between the two SAX representation of C 1 and C 2 is given as [21]: Whereĉ i 1 andĉ i 2 are the i th element of theĈ 1 andĈ 2 , respectively.
L 2 norm distance (Euclidean distance) between C 1 and C 2 is: The two-dimensional correlation between the two images I ref and I test images are computed as following: where m and n represent the first and the second dimension of the image, respectively.Ī ref andĪ test also represent the mean value of the reference and test images, respectively. We applied PCA using two principle components and the similarity between I ref and I test is calculated as the correlation between the scores in their corresponding principle components. PTM score is also completely explained in [2], [3].
Basically, we quantized the image treated as a time series using SAX symbols and then apply the similarity criterion in PTM, therefore we call our proposed network or method, Quantized Peak Topography Map or briefly QPTM.
In order to evaluate QPTM, we test it against the other established comparison criteria, L 2 norm distance, Correlation, PCA, PTM score and SAX. Note that in these methods, we don't use the proposed network, we just apply these methods directly and use them as the difference criteria. Specifically, we have used SAX algorithm both individually as the comparison criterion, which is named "SAX" method in the figures and also as the Comparison layer in function G in QPTM. Also in the method named "PAA" in the figures, we have applied the dataset into the proposed network but in the Comparison layer, we just computed the Euclidean distance between the PAA representations of the reference and test images.

A. Data description
We verify the proposed method in this work with a dataset of thirty-four injections with GC ×GC patterns from different parts of the world. Of particular interest, there are seventeen injections from the Gulf of Mexico from which fourteen of them are from Macondo well, one injection from Eugene Island, one injection from Southern Louisiana and one injection is a natural seep in the Gulf of Mexico. There also exist three standard NIST injections unrelated to the Gulf of Mexico in our model dataset. The rest fourteen injections are from the other parts of the world.
As stated in III for a given test image, we should find the most similar image from the library in which the similarity is also above a threshold. But it requires to have at least one image from each of the geographical region, in which we don't have in our model dataset. Therefore, we just consider the case where we have twenty-seven test samples and test them against the first element of D, which is an injection from Macondo well. Note that our model dataset has thirty-four injections, but because we have set aside seven of the injections from Macondo well as the training set for computing , we end up having a dataset of size twenty-seven. Once, we compare our test images against this Macondo injection, the question will turn into the detection of the Macondo samples from the test images. In other words, we are to detect either each of the test images is from Macondo well or not. For better illustrating the comparison between QPTM and the other methods, we have shown the results of comparison of our twenty-seven injection with respect to the reference injection from the Macondo well (the first element of D) in four consecutive figures, Figure 5, 6, 7 and 8, in each of them comparing QPTM against some of the other methods.

B. Discussion over the results
In each of these figures, we are testing the other injections against the first injection in our dataset which is from Macondo well. We have used the reference injection along with six other injections as the training set for learning the parameter for the family of the images from Macondo well, which has been explained comprehensively in Section IV-C.
As can be seen, in all of the Figures 5 to 8, the first thirteen injections have a high percentage of match with respect to the reference injection. The first six injections are the Macondo injections used for learning of Macondo well and the rest seven injections are some other injections from the Macondo well. Therefore, these figures show that the proposed method have performed well in detecting the Macondo injections. In these figures, we have also shown the standard NIST injections. These injections have almost exactly similar percentage of match to the reference injection. In Figure 9 we have also shown the cross QPTM score for the thirty-four injections using = 0.5.
The problem statement is detecting the injections from different geographical regions and classifying them accordingly. Therefore, our proposed network serves as a classifier. The commonly-used metric to observe the accuracy of a classifier is its confusion matrix. In table II we have shown the confusion matrix for detecting the injections from Macondo well, using different methods in matrix C. C i,j element of C represents the number of the injections being from region indexed by i and being detected as j. In our case the confusion matrix is a two by two matrix, as our classifier has two classes, being from Macondo or not. The first and second rows of C correspond to the injections being from Macondo well and some other source rather than Macondo, respectively. Likewise, the first and second column of C represents the injections being detected as Macondo well and some other source, respectively.
As can be seen in the figures 5 to 8, we have chosen θ = 96% as the classification threshold. In Table III we have shown the values for Accuracy, precision, Sensitivity, Specificity and F 1 scores for the different methods.
As can be seen QPTM has the highest scores in Accuracy, Sensitivity and F 1 among the others. In Table IV we

VII. CONCLUSION
Source differentiation in petroleum fingerprinting using hydrocarbon biomarkers is a classic "near-far" disambiguation challenge from the signal processing perspective, where "near" and "far" represent the stronger regional fingerprint and weaker source-specific fingerprint respectively. In this paper, we propose a signal representation method called quantized Peak Topography Mapping (QPTM) along with related differentiation techniques aimed at better separation between the petroleum injections. The proposed method provides a layered interpretation of two-dimensional gas chromatography images. The main contributions of this work is interpreting intricate peak profile distributions as a multi-dimensional time series, which is quantized resulting in its symbolic aggregate representation. Thus, the distance between two time series is reduced to the difference between their symbolic representations, which is performed using a look-up table (Table I). QPTM out-performed the currently established methods for better classification of the petroleum injections in terms of accuracy, sensitivity and F 1 score. Fig. 9. Cross-QPTM score for the thirty-four injections using QPTM for = 0.5. The white box on upper left are injections from the Macondo well which have a very high percentage of similarity. Note that, in each row, we set the first injection as the reference sample, the rest six injections as the training injections to evaluate the and then the rest would be the testing set. NIST injections are also shown at the middle of the image.