Zooming Into the Darknet: Characterizing Internet Background Radiation and its Structural Changes

Table 1 summarizes the features used to represent Darknet activity. As mentioned earlier, Darknet events in this study are aggregated on a daily basis. For example, a scanner identified $n$ times within a given day—possibly scanning multiple ports and protocols—will be represented in the clustering input as a single entry annotated with the features of Table 1. The chosen features fall into the following categories: 1) intensity of the activity (e.g., total number of packets captured in the Darknet for a given scanner within the day, total duration of scanning, number of ports scanned, etc.), 2) scanning intentions (i.e., the set of ports and protocols involved), 3) scanning strategy (i.e., information about the unique destinations reached in the Darknet, both in terms of IPs and /24 prefixes), and 4) external threat intelligence for the scanner itself (i.e., information gathered from Censys.io). We also annotate each event with BGP routing information, organizational and geolocation features, DNS-based features and others to help the analyst interpret the clustering results.

Notice here that several features used are extremely high dimensional; e.g., the number of unique TCP/UDP ports is $2^{16}$. Therefore, we chose an one-hot encoding scheme for these high-dimensional features where only the top $u$ values (ranked according to the number of distinct source IPs involved) of each feature are encoded for the task at hand.

We emphasize here that the feature set could be easily expanded (e.g., to include time-series data for each Darknet scanner or encoding DNS information, when available) and the proposed framework could easily accommodate the added features. We leave this feature expansion as part of future work.

Motivated by the recent success of deep representation learning, we employ the idea of autoencoders (pmlr-v70-yang17b, ; 10.1145/3278532.3278545, ; 2013arXiv1312.6114K, ) to learn low-dimensional numerical embeddings of the input data. The resulting heterogeneity of the input data features, their high dimensionality, and the need to cope with potentially non-linear interactions between features motivated us to employ deep representation learning to address these challenges. Figure 2 illustrates how the proposed approach can “learn” meaningful features and map them into a low-dimension latent space, while keeping representations of similar points close together in the latent space. These low-dimension embeddings are then passed as input to a traditional K-means clustering algorithm to get the sought clusters. The workflow is shown in Figure 3.

One can employ fully-connected MLP neural networks for the implementation of both mapping functions $f(\cdot;\theta),g(\cdot;\gamma)$. In order to learn representations that preserve the information of input data, we consider minimizing the reconstruction loss:

where the $\ell(\cdot):\mathbb{R}^{P}\rightarrow\mathbb{R}$ is a loss function that quantifies the reconstruction error. For simplicity, we choose the Euclidean distance $\ell(\bm{x},\bm{y})=\left\lVert\bm{x}-\bm{y}\right\rVert_{2}^{2}$. $R(\cdot)$ is a regularization term for the model parameters to help us avoid “overfitting” the data. In this work, we adopt the squared $\ell_{2}$ norm, such that $R(\theta)=\left\lVert\theta\right\rVert_{2}^{2}$. $\lambda\geq 0$ is the regularization coefficient. All model parameters, i.e., $\{\theta,\gamma\}$, can be jointly learned using standard stochastic gradient-based optimization methods, such as Adam (kingma2014adam, ). The dimension $Q$ of the latent space, the number and size of inner layers of the MLP architecture and other hyper-parameters (such as the regularization coefficient) are tuned using the procedure outlined in the sequel (see Section 5.2).

We perform standard K-means clustering directly on the low-dimensional representation to partition the data. Formally, in this step, we aim to minimize the following clustering loss:

where $M$ is the clustering assignment matrix, the entries of which are all binary. $C$ is the matrix of clustering centers that lie in the representation space. $1_{K}$ is a $K$-dimensional column vector of ones. We utilized broadly available libraries for the implementation of the representation learning and clustering steps, such as PyTorch and scikit-learn.

To evaluate the quality of the clustering results, we utilized the silhouette score, which is the standard measure for clustering quality, along with two additional evaluation metrics to obtain a more comprehensive assessment: the Jaccard score (halkidi2001survey, ) that is based on a set of external labels (created for evaluation purpose, and not used by clustering) based on domain knowledge (e.g., known fingerprints of Mirai malware, known port vulnerabilities such as port 445 associated with the attacks on the SMB protocol, known heavy scanners such as shodan.io, Censys.io, etc.) to determine the similarity of the clustering outcomes with the partition formed by the labels; and a stability score that evaluates the stability of clustering results by comparing the similarity of clusters generated by different sub-samples.

where $M_{11}$ is the total number of pair of points that belong to the same group in $C$ as well as the same group in $P$, $M_{01}$ is the total number of pair of points that belong to the different groups in $C$ but to same group in $P$ and $M_{10}$ is the total number of pair of points that belong to the same group in $C$ but to different groups in $P$.

This cluster evaluation metric incorporates domain knowledge (such as Mirai, Zmap and Masscan scanners, that can identified by their representative packet header signatures (MiraiUSENIX2017, ; Durumeric:2014:IVI:2671225.2671230, ), and other partitions as outlined earlier) and measures how compliant the clustering results are with the known partitions. Jaccard score decreases as we increase the number of clusters used for clustering. This decrease is drastic at the beginning and slows down eventually forming a “knee” (see Figure 5). The “knee” where the significant local change in the metric occurs reveals the underlying number of groupings in the data (halkidi2001survey, ).

In order to analyze the stability of the clusters, we generate multiple subsampling versions of the data by using bootstrap resampling. These samples are clustered individually using the same clustering algorithm. The cluster stability score is, then, the average of the pairwise distances between the clustering outcomes of two different subsamples. For each cluster from one bootstrap sample, we identify its most similar cluster among clusters from another bootstrap sample using Jaccard index as the pairwise distance metric. In this case, the Jaccard index is simply the ratio of the intersection and union between the clusters. The average of these Jaccard scores across all pairs of samples provides a measure of how stable the clustering results are.

Systematic tuning of the MLP hyper-parameters is critical to ensure that the network learns an information-preserving low dimension representation of the high dimensional input data. Following best practices for hyper-parameter tuning in machine learning, the set of hyper-parameters is adjusted to find parameter values that result in low autoencoder loss based on evaluation data not used in the training phase of the autoencoder. Hyper-parameters tuned include the size $Q$ of the embeddings (i.e., latent space), the regularization coefficient, the “learning rate” used in stochastic gradient descent, the “dropout” probability employed to prevent over-fitting the training data, the number of “epochs” and the “batch size”.

We employed the Ray[Tune] framework (raytune, ) to tackle the important task of hyper-parameter tuning. Based on the tuning results, we choose the dimension of the latent space to be 50. The “best model” for both MLP and thermometer encoding MLP have a hidden layer of dimension 1000. Using 100 epochs, a batch size of 2000 data points and a value of 0.001 for learning rate, dropout probability and regularization weight provided optimal performance. All results presented in the following sections were generated using these tuned hyper-parameters settings.

To evaluate the clustering outcomes of the two autoencoder architectures we utilize the 3 metrics outlined earlier. We employ a Darknet dataset compiled for the day of January 9th, 2021, which includes about 2 million scanners. Figure  5 shows the clustering performance for different number of clusters ($K$). Following the suggestions in (halkidi2001survey, ), we thus select the number of clusters for clustering Darknet data to be 200. A random sample of 500K scanners were used to perform 50 iterations of training autoencoders and k-means clustering, using 50K scanners in each iteration. The mean and standard deviation of the three clustering evaluation metrics, as well as the mean and standard deviation of the loss function (L2 for MLP, Hamming distance for thermometer-encoding-based MLP (TMLP)), are shown in Table 2.

The results indicated that the TMLP autoencoder led to better clustering results based on the silhouette and stability scores. However, a smaller Jaccard score was reported when compared to the MLP autoencoder. By inspecting the clusters generated, we noticed that this is probably due to the fact that TMLP tended to group scanners into smaller clusters that are similar in size. I.e., it generated multiple fine-grained clusters that correspond to a common large external label used for external validity measure (i.e., the Jaccard score). Because our current Jaccard score computation does not take into account the hierarchical structure of external label, fine-grained partition of external labels are penalized, even though they can provide valuable characteristics of subgroups in a malware family (e.g., Mirai). Henceforth, though, we present results using the MLP architecture that scored very well on all metrics and provided more interpretable results.

We propose the use of decision trees (hastie01statisticallearning, ) for clustering interpretation. Decision trees are conceptually simple, yet powerful, for supervised learning tasks (i.e., when labels are available) and their simplicity makes them easily understandable by human analysts. Specifically, we are interested in classification trees.

In a classification tree setting, one is given $N$ observations that consist of $p$ inputs, that is $x_{i}=(x_{i_{1}},x_{i_{2}},\ldots,x_{i_{p}})$, and a target variable $y_{i}$. The objective is to recursively partition the input space and assign the $N$ observations into a classification outcome taking values $1,2,\ldots,K$ such that the classification error is minimized. For our application, the $N$ observations correspond to the $N$ Darknet events we had clustered and the $K$ labels correspond to the labels assigned by the clustering step. The $p$ input features are closely associated with the $P$ features used in the representation learning step (see Section 4). Specifically, we still employ all the numerical features but we also introduce the new binary variables / tags shown in Table 3. These “groupings”, based on domain knowledge, succinctly summarize some notable Darknet activities we are aware of (e.g., Mirai scanning, backscatter activities, etc.) and, we believe can help the analyst easily interpret the decision tree outcome.

Traditionally, classification trees are constructed using heuristics to split the input space (Breiman1983ClassificationAR, ; bertsimas17optimal, ; hastie01statisticallearning, ). These greedy heuristics though lead to trees that are “brittle”, i.e., trees that can drastically change even with the slightest modification in the input space and therefore do not generalize well. One can overcome this by using tree ensembles or “random forests” but this option is not suitable for the interpretation task at hand since one then needs to deal with multiple trees to interpret a clustering outcome. Hence, we decided to work with optimal classification trees (bertsimas17optimal, ) which are nowadays feasible to construct due to recent algorithmic advances in mixed-integer optimization and hardware improvements that speed-up computations.

We used the software provided by interpretable.ai (interpretableai, ) to generate all optimal decision trees shown in the sequel. Figure 6 shows an example generated for 467,293 Darknet events for Sept. 14th, 2020. The structure of the tree, albeit minimal, is quire revealing. First, the leaves correspond to the largest 4 clusters (with sizes 14953, 11013, 10643 and 9422, respectively) found for Sept. 14th which means that the clusters with the most impact are captured. Another important observation is that the type of decision rules used to split the input space (namely, scanning, censys:mgtm and orion:remote) are indicative of the main Darknet activities during that day. Comparing with the non-optimal, heuristic-based decision tree of Figure 7, we recognize some important differences: 1) two new clusters have emerged (with labels 100 and 191) that do not rank within the top-4 clusters (they rank 8th and 10th, respectively, with 6977 and 6404 members); and 2) there is some “redundancy” in the decision rules used for splitting when both the tags UDP and “scanning” are present. This is because UDP and scanning (i.e., TCP SYN requests and ICMP Echo Requests) are usually complementary to each other.

One of the important challenges in clustering is identifying characteristics of a cluster that distinguish it from other clusters. While the center of a cluster is one useful way to represent a cluster, it can not clearly reveal the features and values that define the cluster. This is even more challenging for characterizing clusters of high-dimensional data, such as the scanner profiles in the Darknet. One can address this challenge by defining “internal structures” based on the decision trees learned.

Given a set of clusters ${C_{1},C_{2},\ldots,C_{K}}$ that form a partition of a dataset $D$, a disjunctive normal forms (DNF) $S_{i}$ is said to be an internal structure of cluster $C_{i}$ if any data items in D satisfying $Si$ are more likely to be in $C_{i}$ than in any other clusters. Hence, an internal structure of a cluster captures characteristics of the cluster that distinguishes it from all other clusters. More specifically, the conjunctive conditions of a path in the decision tree to a leaf node that predicts cluster $C_{i}$ forms the conjunctive (AND) component of the internal structure of $C_{i}$. Conjunctive path description from multiple paths in the decision tree that predict the same cluster (say $C_{i}$) are combined into a disjunctive normal form that characterizes the cluster $C_{i}$. Hence, the DNF forms revealed by decision tree learning on a set of clusters expose the internal structures of these clusters. We plan to pursue detection of structural (di)similarities in the Darknet based on DNF forms applied on decision tree outcomes as part of future work.

Figure 8 (bottom row, left plot) shows a significant increase in the EMD distance between the clustering outcome of Sept. 5th and Sept. 6th. Figure 8 (top row, left plot) shows the aggregate number of unique scanners for the total month of September and a detailed view of the top-3 countries for that month (regarding unique scanners). With a bare eye we observe an increase in scanners originating from Egypt, but looking at a mere scanning volume cannot let us diagnose what really happened. We therefore look at the top clusters extracted for both days for an answer.

Tables 4 and 5 tabulate the top-10 clusters for each day. Notice that the column “Mirai Coverage” (which denotes the portion of Darknet scanners within that cluster that include the well-known Mirai fingerprint (MiraiUSENIX2017, ) in the packet header of the packets they emitted) is zero for all top clusters of Sept. 5th but a relatively large Mirai-related cluster emerges on Sept. 6th. Clearly, we now have a good indication about what might have caused the shift in the clustering structure between Sept. 5th and 6th.

Looking at our data more closely we found that there were 15970 total Mirai victims on Sept. 5 and 28261 on Sept. 6th. Table 6 shows the top Mirai-related clusters for Sept. 5th and Table 7 the top ones for the next day. As observed, the Mirai clusters have increased significantly. Notice also that Censys can provide some valuable information for interpreting each cluster; the largest Mirai cluster on Sept. 6th seems to be heavily associated with devices having the FTP port open.

Examining the distribution of countries affected by the Mirai malware (see Table 8, recall that we geo-annotate all scanners for result interpretation, but we don’t use this information as a clustering feature), we see Egypt emerging as the number one affected country. A few days later, India will be facing a similar Mirai outbreak (the large increase in unique scanners seen in Figure 8 for Sept. 10th is India’s Mirai outbreak with scanners having a similar profile like the ones observed in Egypt).

The second case study we examine is attributed to the large EMD spike detected on Nov. 27th (i.e., denoting a large difference in the clustering outcomes between that day and the day before). Figures 9 and 10 illustrate the distribution of the cluster sizes for these 2 days. There are two extremely large clusters formulated on Nov. 27th with sizes 26662 (cluster-1) and 11262 (cluster-19). The largest cluster on Nov. 26th had only about 6000 members (cluster-2). We here start to speculate that these 2 large clusters might be the culprits for the dramatic shift in the EMD distances, but we need to look deeper into the clustering results.

Figure 8 does not indicate that the change is associated with any of the top-5 ports for that month. Similarly, Table 9 and 10 that include the top-5 ports zoomed into both days do not provide any further insights. This illustrates that “group by” operations (performed frequently by forensic analysts when trying to troubleshoot an incident) might sometimes provide limited information. Further, when there are multiple categorical features (as is the case here, with categorical features being the ports scanned, protocols, Censys tags, Censys ports, countries, etc.) and multiple numerical features, forensics analysis and result interpretation might be challenging.

However, by reflecting on the decision tree paths depicted in Figure 11 we start to gain some understanding about what might be responsible for the structural clustering changes between the 2 days. The rule on the root for Nov. 26th immediately tells us that the top 2 clusters are events that are not associated with “Censys:Web” (meaning that ports 80 or 443 were not found open when Censys scanned the scanners belonging in these clusters; see Table 3). At the same time, the decision rule “Censys:Web” is active for the top-2 clusters for Nov. 27th. Hence, the top-2 very large clusters on Nov. 27th and the top-2 clusters on Nov. 26th are characterized by a distinctive feature (i.e., the “Censys:Web” grouping) and this very likely explains the clustering differences observed. Notably, the top cluster for Nov. 27th and the second-top cluster only differ on their last rule of their decision paths (i.e., the rule about “Censys: embedded”) which means the clusters are very similar except when it comes with the “embedded tag” assigned to them by Censys.

Examining Table 11, which ranks the scanners by their Censys-reported “open ports”, seems to validate that the culprits for the Nov. 27th events are scanners with ports 80 and 443 open that became compromised between the 2 days. Closer investigation of these scanners showed that the behavior was UDP-related (shown also in the decision paths as “!scanning”, i.e., not TCP or ICMP scanning; see Table 3) and the scanners were targeting a variety of randomly looking, high-value ephemeral ports. Further investigation is required to fully understand the intentions behind this event that was originating from several countries and BGP routing prefixes. We note, though, that only 606 source IPs from the top-cluster of Nov. 27th match with any of the scanners from the previous day suggesting the event of Nov. 27th involved newly infected hosts.

The last case study involves heavy SSH scanning associated primarily with devices/hosts attributed to a US tier-1 ISP. Consulting Figure 8, we see a change-detection in the EMD chart between Jan. 6th and Jan. 7th. The SSH scanning activity peaks on Jan. 9th, as seen in the top-right panel of Figure 8, but the change in unique scanners between Jan. 6th and 7th is not that obvious when visualizing that plot. Indeed, looking at the top-5 ports of Table 12, we see that port 22 (SSH) is absent (port 22 ranked top-44th on Jan. 6th and top-13th on Jan. 7th).

Nevertheless, the clustering approach is able to capture this onset of this important event. Looking at Table 13, we observe that there was only a single cluster with port 22 as the top port within the cluster. However, the clustering structure changes significantly on Jan. 6th as evidenced by Table 14. We see several SSH-related ports emerging and many of the clusters are associated with the CWMP tag. This suggests that the infected devices might be ISP-managed modems that were susceptible and got compromised. We note here that, although not captured by the EMD detector, a similar SSH incident occurred back on Nov. 6th, 2020 that involved the same tier-1 ISP. In both cases we notified their “network abuse” team about the incidents.

Figure 12 illustrates a pruned version of the optimal decision tree obtained for Jan. 7th (the full tree is of depth 7). It again demonstrates the utility of decision trees for clustering interpretation. Five out of six clusters associated with SSH activity are extensively explained by the decision tree path highlighted with the red nodes that leads to Node 90 in the tree. That set of rules explain the following fraction of scanners (see Table 14): for cluster-20: 2468/2647, for cluster-35: 5621/5739, for cluster-44 192/958, for cluster-66: 1367/1441, for cluster-74: 2/2023 and for cluster-105: 4233/4325.