Timely poacher detection and localization using sentinel animal movement
Study system and species
This study was performed in Welgevonden Game Reserve (WGR), a privately owned game reserve in the Limpopo province, South Africa (24° 10′ S; 27° 45′ E to 24° 25′ S; 27° 56′ E). The reserve is located in the mountainous Waterberg region. WGR was established on former agricultural lands in the early 1980s and the main occurring vegetation types are Waterberg Mountain Bushveld and Sour Bushveld. The Waterberg region has a temperate climate, with two distinct seasons, characterized by the rainfall regime: a dry season ranging from April to September and a wet season ranging from October to March, with an average annual precipitation in WGR of 634 mm. Our study area is an enclosed breeding camp within WGR, with a size of approximately 1200 ha. Main predator species such as lion, cheetah and spotted hyena were excluded from this study area, as well as elephant and rhino.
WGR equipped 35 impala (Aepyceros melampus), 34 blue wildebeest (Connochaetes taurinus), 35 plains zebra (Equus burchellii) and 34 common eland (Taurotragus oryx) with a GPS and accelerometer sensor equipped custom made collar; an estimated 23% of the individual impalas present in the area, 48% of the eland, 40% of the wildebeest and 40% of the zebra. However, due to malfunctioning and errors made in the sensor development process, only 83 of the sensors yielded data at any point in time, thus lowering the effective density of sentinel animals. During the experimental intrusions (see below), the median number of data-yielding sensors was 47, and minimally 30. The animal movement data were recorded day and night and transmitted wirelessly in near real-time to five long-range low-power LoRa radiocommunication gateways in the study area, from where data packages were routed to an on-line data warehouse via a 3G/4G backhaul. The deployment of these sentinel animals were approved by the board and CEO of WGR as a management action and was performed in accordance with relevant guidelines and regulations (see Supplementary GPS Collaring letter).
Experimental intrusions
Between September 2017 and March 2018, WGR employees performed experimental intrusions (lasting ca. 2 h) on foot and by car through the study area, at varying locations and movement routes through the study area, independent from the locations of the sentinel animals. The movement of the intrusions were tracked by GPS, and the relevant metadata for each intrusion recorded (mode of transport, group size, start time, end time). The intrusions were distributed in a stratified way over the mornings, middays and afternoons (with time slots relative to specific solar positions: sunrise, solar noon and sunset). Furthermore, the intrusions were temporally spread in such a way to avoid a disturbance overflow for the sentinel animals, by performing a maximum of five experiments per week and a maximum of two experiments per day (and then only with one intrusion in the morning and one in the afternoon).
Data gathering
The animal sensors gathered location data via GPS and overall dynamic body accelerations47 (ODBA) via a tri-axial accelerometer (range ± 2 g; sampling frequency 100 Hz, down-sampled to 10 Hz prior to analysis). The GPS was scheduled to record spatial position at irregular intervals depending on the level of activity as gauged by ODBA. All sensors were scheduled to record locations every 15 min in the absence of sufficient activity (given that successive fixes were further than 5 m apart, else a geofence was applied and the new coordinate was omitted to save bandwidth and battery power, thereby assuming that the animal still was at its previous location). The GPS fix rate was increased up to 2- or 10-min intervals (depending on two different sensor settings) when ODBA indicated sufficient activity (after checking for the geofence). ODBA data were sampled continuously and summarized per 15 s window in a mean, maximum and variance value.
The experimentally intruding groups were outfitted with handheld GPS devices that recorded their location every 5 s and these groups logged and timestamped all their pre-defined activities and metadata on a tablet using CyberTracker48 during their intrusion. Most cars traveling through the study area were tracked by GPS as well to filter the animal data for disturbances by cars unrelated to the experimental intrusions.
Weather data (temperature, radiation, precipitation and wind) in the study area were recorded on a 3-min resolution with a weather station in the north of the study area. We assumed the 1200 ha study area to be sufficiently small to assume the weather station data to be representable for the prevailing weather conditions throughout the study area. GIS data of the study area (summarized in Supplementary Table S1) consisted of information on topography, infrastructure (e.g., fences, roads, powerlines, etc.) and vegetation cover (supervised classification of 25 cm resolution aerial imagery into four classes: trees, herbaceous/grass, sand/soil and other/built-up area).
All further data processing and algorithm development was done in the software R3.5.049.
Data pre-processing
To link the animal location data with the intrusion location data, as well as to correct for the substantial level of positional noise present in the animal location data, we modelled the animal location data to regular 1-min resolution trajectories using the following five steps. First, we filtered out large obvious errors (e.g., obvious outliers and irregularities such as locations far outside the study area) from the data. Second, we corrected systematic medium-scale outliers: ‘spikes’ that occurred due to positional outliers. Such spike-like outliers were visible during sensor testing while following known straight-line trajectories along an airstrip, thereby confirming that these spike-like geometries most likely resulted from positional error rather than true animal movement. Points were classified as anomalous spike points when (a) the displacement to and away from this point was high ( > 500 m), (b) when the distance between the locations before and after this point was small, and (c) when the turning angle at this point approached 180 degrees. Therefore, we corrected the locations that were classified as spike-like anomalies by shifting them closer to the straight line between the neighboring points. The extent of this shift was set relative to the degree of spikiness of the points (the spikier the pattern, the larger the shift towards the midpoint of the adjacent coordinates). Third, after filtering and correcting the original locations we smoothed the timeseries of x/y coordinates at each original timepoint with a Kalman smoother using a dynamic linear model. Fourth, we linearly interpolated the locations to a 10 s resolution based on ODBA, where we considered the animal to be stationary between multiple timepoints if the accelerometer signal suggested the animal was not moving. Fifth, we fitted an X-spline through the data, where we gave the linearly ODBA-interpolated locations a smaller weight, and sampled the fitted spline on a regular 1-min resolution. These pre-processing steps resulted in the modelled animal trajectory data, composed of spatial locations every minute, and averaged ODBA statistics per step (i.e., the segments between consecutive coordinates). These data were used as input for the next steps in the analyses. In contrast to the animal data, the raw intrusion data were of a high temporal resolution and spatial accuracy so that we only needed to subset the data in order to acquire 1-min resolution time-synchronized intrusion trajectories.
The first three parts of the data pre-processing were only needed because of firmware issues in our custom-made sensors. Without these issues, a simple denoising technique like a Kalman filter will suffice.
Feature engineering and processing
We computed a plethora of human-engineered features from the animal trajectories, ODBA data, weather data and several GIS layers with environmental data from the study area (summarized in Supplementary Table S1). All features were computed such that they could not directly be linked to specific points in space or time (by computing movement features relative to the environmental variables), so that only behavioral patterns and abnormalities therein could be linked to intrusion presence. After engineering these base features, we transformed certain features (after visual inspection of the histograms) to approximately symmetric distributions using logarithms. Then we truncated the distributions to the lower and upper 0.001 percentile to correct possible outliers. After that, we standardized all computed features to zero mean and unit variance per species. We also computed scaled versions of selected features by subtracting the mean and dividing by the variance of the selected features per reference set to capture deviations from normal behavior: (1) per area (characterized by a 30 by 30 m neighborhood around each grid cell), (2) per time of day (morning, midday, afternoon) in a period of 5 weeks around each intrusion or control, (3) per area per time of day per 5 weeks, and (4) per individual sentinel per time of day per 5 weeks (Supplementary Table S1). Furthermore, after computing and standardizing the features, we computed more features by applying moving window computations (5 min centered, 10 and 20 min lagging, and the difference between these: 5 min centered minus 10 and 20 min lagging) on the standardized features to capture (the change in) the recent history of animal movement descriptors (mean and standard deviation of all features, fitted Mean Squared Displacement exponential function parameters, net-gross distance ratio and variance of log First Passage Times). Finally, we discretized all features to ordinal values to avoid odd-, fat- and heavy-tailed distributions. In total we computed 2117 features describing different aspects of movement geometry of individual trajectories, herd topology and the interactions with landscape variation.
Subsetting and dimensionality reduction
Before analyzing the computed animal movement features, we applied some filtering on the data. We removed all periods with an experimental intrusion during which there were less than 30 active animal sensors in total. We also removed data of both animals and intrusion when they were close to the reserve’s main gate in order to avoid dilution of the data with other known disturbances. This resulted in 57 intrusions that were selected for further analyses. For every intrusion we selected control data of the same period one or 2 days earlier or later during which no intrusion took place, resulting in an approximately balanced intrusion-control dataset. Furthermore, we removed data from animals that were located within 250 m and within 20 min of a vehicle moving through the area that was not part of our experiment.
For each feature, we computed 4 importance metrics based on binary labelled data: records associated to locations within 1 km from the intrusion (subscript 1) versus an equally-sized random selection of data points during control periods (subscript 0): Mahalanobis distance, marginality (computed as (frac{{mu }_{1}-{mu }_{0}}{{sigma }_{0}}), for sample mean (mu) and sample standard deviation (sigma)), specialization (computed as (frac{{sigma }_{1}}{{sigma }_{0}})) and the Mean Decrease Accuracy of a Random Forest classifier (with default hyperparameters). We then ranked the features according to their importance and selected a feature for further analyses if it occurred in the top 125 features for any of the 4 importance measures described above (resulting in a total of 361 selected features). Subsequently, we converted the selected features per main feature class (Supplementary Table S1) to principal components, keeping those principal components that capture the most variation (in total 95%), which resulted in 99 selected components in total. Finally, we transformed these components again via a second principal component analysis, now across all the selected 99 components. In subsequent training of the animal behavior classifier, we optimized the total number of included components as a hyperparameter, which resulted in the first 8 principal components in the best performing classifier.
Labelling
We labelled the sentinel movement data through visual inspection of the animal and intruder trajectories, where we considered the animals’ behavior to be undisturbed when the animal was not near an intrusion, or when the animal was close to an intrusion yet did not visually display a change in behavior. However, when the animal was near the intrusion and displayed a sudden or gradual behavioral change in response to intrusion proximity, we labelled the data as ‘flight’ (changing the movement direction away from the intrusion, possibly with increased speed) or ‘regroup’ (when individuals clustered together). In total, only ca. 1% of the animal data were associated to either flight or regroup behavior (which we will refer to as ‘response’ behavior). A few animals also appeared to exhibit behavior we could label as ‘freeze’, i.e., halting movement in the proximity of the intrusion, yet this class was too underrepresented to be accurately predicted and hence dropped from the final dataset. Furthermore, we assigned a qualitative measure of intensity to each labelled behavioral response (‘low’, ‘medium’, ‘high’) to describe how visually pronounced this response was. Besides the supervised labelling based on visual inspection of behavioral responses via video animations of the trajectories, we also labelled data using an unsupervised k-means nearest neighbor classifier, where we clustered the feature space consisting of the 99 features selected as described above into 25 clusters per species.
Animal behavior classification
We trained an RBF kernel C-classification Support Vector Machine (SVM) with a subsequent moving window over the outputted probabilities to distinguish undisturbed versus response behavior. In the training datasets we only included the data separated by more than 1 km from the intrusion and labelled as ‘undisturbed’, and removed 90% thereof to train algorithms with a more balanced dataset. Furthermore, we only trained and validated on data with intrusions present in the area. We trained another SVM to distinguish the flight response from the regroup response. All computations were done in R 3.5.0 with the e1071 package on the Linux High Performance Cluster of Wageningen University and Research. We optimized the following hyperparameters and model settings during the training phase for the Average Precision via a grid search (with the selected values between brackets):
gamma (undisturbed-response: 10–3.2; flight-regroup: 10–2.0);
cost (undisturbed-response: 10–2.2; flight-regroup: 10–1.5);
number of principal components to include as features (undisturbed-response: 8; flight-regroup: 12);
species-specific models versus one model with species dummy variables included in the features (species-specific models);
specific models for the different times of day versus one model with time of day dummy variables (one model);
response intensities to include in the training data (only medium and high intensities);
weights to assign to the classes (equal weights);
the quantile to be computed of the SVM probabilities by the moving window (100%, i.e., maximum value);
the alignment of the moving window (centered);
the size of the moving window (15 min on both sides).
The best model was selected via a leave-one-intrusion-out cross-validation approach. We summarized the predictive performance by computing the Average Precision of the least occurring class (i.e., ‘response’ for the undisturbed-response model: 46%, Supplementary Fig. S2; and ‘regroup’ for the flight-regroup model: 80%, Supplementary Fig. S3). After having computed these probabilities with an SVM and a temporal window smoother, we tried to improve the predicted performance by including the predicted animal response probabilities of nearby animals. However, this spatial explicit approach hardly improved the predictive performance, indicating that the spatial contextualization of behavioral response was sufficiently captured by the computed features. We therefore did not include this spatial contagion effect of predicted animal response probabilities in the final analysis.
System classification—detection
Based on the predicted SVM response probabilities and feature cluster analysis, we computed summary features per 15 min of each intrusion and control period. These summary features related to the odds ratios of the probability of association of unsupervised clusters with intrusions versus controls, the SVM predicted probabilities of behavioral response, and several features describing the values (and its spatial structure, e.g., clustering or autocorrelation) of these SVM predicted response probabilities. After computing summary features per 15 min, we summarized them even further for the intrusions versus controls using the following eight statistics: mean, standard deviation, minimum, maximum, mean of the lagged differences, standard deviation of the lagged differences, minimum of the lagged differences and maximum of the lagged differences.
After computing the summary features, we build a logistic regression classifier to distinguish intrusions from controls. To create a parsimonious model, we iteratively added features to the model and evaluated its performance after each iteration. We evaluated the performance based on the model accuracy and performed validation through 25 times twofold cross-validation in a stratified way (by 25 times choosing a balanced random sample of intrusions and controls). We determined the sequence of adding features to the model by performing an independent two-sample t-test for each feature between the intrusions and controls. The feature with the largest t-value was then added to the model. After each feature addition, we removed its correlation with the remaining features using linear regressions with the added feature as independent variable and the remaining features as dependent variables, from which we extracted the residuals, standardized them to zero mean and unit variance, and applied the t-tests again. The (original) feature with the largest t-value was then added to the model again. This procedure was repeated until all features were ordered corresponding to their “importance”. We then performed logistic regressions without interactions between the features for an increasing number of features (Supplementary Fig. S4). The model already performed quite accurately with only 7 features (86.1% accuracy ± SD 3.3%, precision 82.6% ± SD 6.9%, recall 89.2% ± SD 5.1%). However, with 20 features and 2-way interactions the model achieved the maximum accuracy (90.9%).
System classification—localization
The data gathered during intrusions that were correctly predicted as such by the detection classifier were used to train the intrusion localization algorithm. The probability surface of the location of the intrusion was fitted relative to that of the sentinel animals using:
$${O}_{i,j} sim frac{{p}_{j}left({f}_{wn}left({theta }_{i,j},{mu }_{j},{rho }_{1}right) {f}_{ln}left({gamma }_{i,j},{mu }_{1},{sigma }_{1}right)right) left(1-{p}_{j}right) left({f}_{wn}left({theta }_{i,j},{mu }_{j},{sigma }_{0}right) {f}_{ln}left({gamma }_{i,j},{mu }_{0},{sigma }_{0}right)right)}{{f}_{wn}left({theta }_{i,j},{mu }_{j},{rho }_{0}right) {f}_{ln}({gamma }_{i,j},{mu }_{0},{sigma }_{0})}$$
where ({O}_{i,j}) is the odds ratio of intrusion presence at location (i) evaluated for individual (j), ({p}_{j}) is the SVM-predicted probability that individual (j) is exhibiting response behavior. The function ({f}_{wn}) is the wrapped normal probability density function, ({theta }_{i,j}) is the direction from location (i) to the location of the focal animal (j), ({mu }_{j}) is the movement direction of individual (j), ({rho }_{1}) and ({rho }_{0}) are the standard deviations of the unwrapped distributions. The function ({f}_{ln}) is the lognormal probability density function, where ({gamma }_{i,j}) is the distance of location (i) to (j), ({mu }_{1}) and ({mu }_{0}) as well as ({sigma }_{1}) and ({sigma }_{0}) are the log-normal distribution parameters (respectively log-mean and log-sd).
The parameters ({mu }_{1}), ({sigma }_{1}) and ({rho }_{1}) capture the geometry of intrusion-animal topology for animals that exhibited a predicted behavioral response to the intrusion. Similarly, ({mu }_{0}), ({sigma }_{0}) and ({rho }_{0}) are the corresponding parameters for animals that were predicted to be undisturbed. The parameters ({mu }_{1}), (mathrm{log}({sigma }_{1})) and (mathrm{log}({rho }_{1})) were fitted to the data assuming a 3rd order polynomial relationship to ({t}_{s}): the time (in minutes) since the start of the predicted behavioral response (using the maximum F1 classification score). Since the behavioral response signature is lost over time, we truncated ({t}_{s}) to 45 min (thus ({t}_{s} >45) min was set to ({t}_{s}=45)). The parameters ({mu }_{0}), ({sigma }_{0}) and ({rho }_{0}) were estimated using the data of the controls and with randomly generated intrusion locations in the study area, in order to correct for the effects of geometry of the study area on the predicted response surfaces. The probability surface ({P}_{i}) was then calculated as:
$${P}_{i}=alpha sum_{j}{O}_{i,j}$$
where (alpha) is a normalization constant so that ({P}_{i}) integrates to 1 over the area covered by the rectangular axis-aligned bounding box around the study area.
To measure the prediction accuracy of each localization surface, we simplified each surface to a point coordinate located at the location of maximum probability, and computed the Euclidian distance to the known true position of the intrusion. We then summarized each experimental intrusion by selecting the 10 prediction surfaces with the most condense highest probability density, i.e., those in which the top 5% probability density is contained in the smallest, most condense, area. The spatial error of the localization prediction associated with these selected predictions was further summarized by taking the average Euclidian distance over the 10 selected predictions. More
