Method and system for analysing distributed acoustic sensing data

Abstract

The invention provides a method of analysing data from a distributed acoustic sensing system. The distributed acoustic sensing system comprises a plurality of spatial channels, where each spatial channel associated with a respective scattering location along an optical fibre of the distributed acoustic sensing system. The method comprises obtaining a set of candidate trajectories corresponding to possible trajectories represented in data from the distributed acoustic sensing system. The method further comprises receiving, from the distributed acoustic sensing system, distributed acoustic sensing data comprising a plurality of measurement signals as a function of time, each of the plurality of measurement signals corresponding to a respective one of the spatial channels. An objective function is computed based on the distributed acoustic sensing data, wherein the objective function is indicative, for each of the candidate trajectories, of a likelihood of that candidate trajectory being represented in the distributed acoustic sensing data. An iterative process is then performed, in which one of the candidate trajectories having a highest likelihood of being represented in the distributed acoustic sensing data is selected. A representation of detected trajectories is output, based on one or more candidate trajectories selected in the iterative process.

Description

[0001]METHOD AND SYSTEM FOR ANALYSING DISTRIBUTED ACOUSTIC SENSING DATA Field of the InventionThe present invention relates to a method of, and system for, analysing data received from adistributed acoustic sensing (DAS) system. In particular, the invention provides a method of, and system for, detecting trajectories that are represented in DAS data. Background Distributed Acoustic Sensing (DAS) is an established technique with several commercialsystems available. In these systems, a pulse or pulses of laser light are launched into a lengthof optical fibre and the light that is scattered within the fibre is analysed in order to derive thenature of the acoustic environment, i.e. any physical vibrations, of the fibre transducer. Thesesystems typically make a measurement of the acoustic strain environment of an optical fibre transducer using an optical time domain reflectometer (OTDR) approach. This gives a differential strain measurement as a function of position along the optical fibre. As an optical fibre is manufactured it is cooled or quenched from a high temperature as it is drawn. This process leads to the presence of small variations in the density of the optical fibre. These tiny variations in density equate to variations in the effective refractive index of the fibre. These discontinuities lead to scattering of laser light passing through the optical fibre, particularly by Rayleigh scattering. The amplitude of the scattering follows a Rayleigh distribution, but the phase angle of the scattering is uniformly distributed around a unit circle, i.e. the phase angle. In the OTDR approach, a coherent light pulse is sent along the optical fibre. As the pulse travels along the optical fibre, it is scattered at scattering locations along the optical fibre, andbackscattered signals are received at a detector stage. Scattered (i.e. backscattered) signalsarising from different scattering locations along the optical fibre will be received at the detectorstage at different times, such that each scattered signal can be assigned to its correspondingscattering location (or “spatial channel”) based on its time of receipt at the detector stage. In thismanner, measuring the scattered signals over time enables a strain, and hence an acoustic field, at the different scattering locations along the optical fibre to be monitored over time.The present invention has been devised in light of the above considerations. Summary of the Invention At its most general, the present invention provides a method of (and system for) analysing distributed acoustic sensing (DAS) data to detect trajectories represented in the data, such as vehicle trajectories. For example, a DAS system can include an optical fibre which extends along (e.g. under) a vehicle path such as a road or track. Vibrations caused by vehicle movement along the path are coupled into the optical fibre, such that data collected from the DAS system can be used to detect vehicle movement along the path as a function of time. The invention enables efficient analysis of DAS data, to allow for rapid detection of vehicle trajectories in the DAS data. According to a first aspect of the invention, there is provided a method of analysing data from a distributed acoustic sensing system, the distributed acoustic sensing system comprising a plurality of spatial channels, each spatial channel associated with a respective scattering location along an optical fibre of the distributed acoustic sensing system, wherein the method comprises: obtaining a set of candidate trajectories corresponding to possible trajectories represented in data from the distributed acoustic sensing system; receiving, from the distributed acoustic sensing system, distributed acoustic sensing data comprising a plurality of measurement signals as a function of time, each of the plurality of measurement signals corresponding to a respective one of the spatial channels; computing an objective function based on the distributed acoustic sensing data, wherein the objective function is indicative, for each of the candidate trajectories, of a likelihood of that candidate trajectory being represented in the distributed acoustic sensing data; performing an iterative process comprising: selecting,based on the objective function, one of the candidate trajectories (from the set of candidatetrajectories) having a highest likelihood of being represented in the distributed acoustic sensingdata; identifying a set of data points in the distributed acoustic sensing data corresponding to the selected candidate trajectory; updating the objective function by selectively removing from the objective function a contribution arising from the identified set of data points; and repeating the iterative process until a termination condition is met; and outputting a representation of detected trajectories, based one or more candidate trajectories selected in the iterative process. The method is a computer-implemented method, which can be implemented using any suitable computing system. For example, the method may be implemented by an analysis system (e.g. as described in the second aspect of the invention, below), which includes a computer memory storing instructions for performing the method, and a processing device (e.g. including one or more processors) configured to execute the instructions to perform the method. The DAS system with which the method is used may comprise any suitable type of DAS system, such as an optical time domain reflectometer (OTDR). The DAS system may be configured to launch a pulsed test signal along the optical fibre, and to receive at a detector stage a plurality of scattered signals, each of which was scattered at a respective location along the optical fibre. Each of the scattered signals thus corresponds to a respective “spatial channel” of the DAS system. Scattered signals corresponding to different scattering locations (spatial channels) can be distinguished on the basis of their time of receipt at the detector stage, e.g. by comparing a time at which a pulse of the test signal was launchedalong the optical fibre and the time of receipt of the scattered signal at the detector stage, takinginto account the speed of light along the optical path. The detector stage may then be configured to output a respective measurement signal for each of the plurality of spatial channels, the measurement signal being a function of the scattered signal for that spatial channel. The DAS data includes a respective measurement signal as a function of time for each of the plurality of spatial channels. Thus, each measurement signal may correspond to a time-series of data points. The measurement signal for each spatial channel may be derived from (i.e. be a function of) the scattered signal received at the detector stage for that spatial channel. The specific nature of the measurement signal may depend on a type of DAS system used and / or a type of measurement performed. For example, each measurement signal may be a function of a phase and / or amplitude of the scattered signal received for the corresponding spatial channel. In some cases, the detector stage may be configured to interfere the scattered signals with a local oscillator signal, in which case the measurement signal may be a function of a phase difference between the local oscillator signal and the scattered signal for the corresponding channel. In some cases, each measurement signal may represent a power level as a function of time in the associated spatial channel in an energy band of interest. This may be obtained, for example, by applying a fast Fourier transform (FFT) to the scattered signal, and averaging the power level over multiple time samples. The optical fibre of the DAS system may extend along a vehicle path, such as a road or a track (e.g. train track). For example, the optical fibre may be buried under the vehicle path. In this manner, vibrations caused by vehicles traveling along the vehicle path will be coupled into theoptical fibre, enabling vehicles moving along the path to be tracked. In particular, as a vehiclemoves along the path (and hence along the optical fibre), it will generate a response in spatial channels of the DAS system at different points in time. Thus, a trajectory of the vehicle may berepresented in DAS data as a response of the spatial channels as a function of time.Herein, a trajectory may refer to a vehicle trajectory. A trajectory may be defined in terms of spatial channel (i.e. position) as a function of time. Prior to analysis of the DAS data, the set of candidate trajectories is obtained. The set of candidate trajectories may herein be referred to as a ‘dictionary’ of possible trajectories for theDAS system. Thus, the set of candidate trajectories may provide a description of possible (e.g.realistic) vehicle trajectories that can be represented in the DAS data, e.g. taking into account parameters such as possible vehicle speeds and directions of travel relative to the optical fibre. For example, the set of candidate trajectories may be computed based on a predetermined range of vehicle speeds (and optionally directions) relative to the optical fibre. Each candidate trajectory may comprise a parametric definition (description) of the candidate trajectory. For example, each candidate trajectory may be a straight line, which is defined by a slope and an offset. This allows for a relatively simple parametric definition of the candidate trajectories, which may enhance computational efficiency of the method.The set of candidate trajectories may be obtained (e.g. accessed, loaded) from a memory inwhich the set of candidate trajectories is stored. Thus, the set of candidate trajectories may bepre-computed and stored in the memory. As the set of candidate trajectories is obtained prior toperforming the iterative process, computing of the set of candidate trajectories does not affect a run time of the iterative process. As a result, a high-density set of candidate trajectories can bepre-computed, which enables high precision trajectory detection without significantly increasingthe run time of the iterative process. The objective function provides, for each candidate trajectory in the dictionary, an indication of likelihood of that candidate trajectory being represented in the DAS data. For example, the objective function may provide a probability or related value for each of the candidate trajectories.The objective function is computed (calculated, determined) using the received DAS data, todetermine the likelihood of each candidate trajectory being represented in the DAS data. For instance, the objective function may be computed by comparing each candidate trajectory to theDAS data, to determine a fit quality (or an error signal) between the candidate trajectory and theDAS data. A candidate trajectory which has a high fit quality (low error signal) with the DAS data may have a relatively high likelihood of being represented in the DAS data, compared toanother candidate trajectory having a lower fit quality (higher error signal).Following computing of the objective function, an iterative process is performed to detect trajectories in the DAS data, by determining which of the candidate trajectories are represented in the DAS data. The iterative process includes a sequence of steps for detecting a trajectory in the DAS data, which is repeated until the termination condition is met. Thus, following termination of the iterative process, one or more trajectories may be detected in the DAS data. In the iterative process, the objective function is used to select the candidate trajectory having the highest likelihood of being represented in the DAS data. For example, this may be achievedby maximising the objective function over the set of candidate trajectories. In other words, thecandidate trajectory for which the objective function has the highest value is selected. Subsequently, the set of data points in the DAS data corresponding to the selected candidate trajectory is identified. In other words, a portion of the DAS data which may represent the selected candidate trajectory is identified. This may be achieved, for example, by mapping the selected candidate trajectory onto the DAS data to determine the set of data pointscorresponding to the selected candidate trajectory. For instance, data points which lie on orwithin a predetermined threshold (distance) of the candidate trajectory may be determined as corresponding to that candidate trajectory. Thus, coordinates (i.e. spatial channel, time) of data points in the DAS data can be checked to see if they lie on or within a predetermined threshold of the selected candidate trajectory. Once the set of data points corresponding to the selected candidate trajectory is identified, the objective function is updated by selectively removing from the objective function a contribution arising from the identified set of data points. This process of updating the objective function is referred to herein as ‘notching’. Thus, the contribution from the set of data points corresponding to the selected candidate trajectory is removed from the objective function, so that it is not takeninto account in subsequent iterations of the iterative process. Accordingly, at each iteration, theobjective function is updated so that previously detected trajectories no longer contribute to the objective function. In this manner, at each iteration, a new candidate trajectory representing a best fit to the DAS data is selected, and then notched from the objective function. Selectively removing from the objective function the contribution arising from the identified set of data points means that only a relatively small portion of the objective function is affected. In thismanner, the entire objective function does not need to be re-computed when updating theobjective function. Rather, by selectively removing the contribution from the identified set of data points, only the likelihoods for candidate trajectories whose representations comprise one or more of the set of data points will be affected. Thus, the objective function need only be re- computed for candidate trajectories whose representations comprise one or more of the set of data points. This process of selectively updating the objective function avoids the computationally intensive task of re-computing the entire objective function. Thus, the selective updating of the objective function contributes to speeding up of the iterative process, by reducing bottlenecks associated with computation of the objective function. The iterative process is repeated until a predetermined termination condition is met. At each iteration, one of the candidate trajectories is selected, providing a potential trajectory detection in the DAS data. Various different types of termination condition may be used for determining when to terminate the iterative process. For example, the iterative process may be repeated until a predetermined number of trajectories is detected. As another example, a quality of detected trajectories may be monitored, and the iterative process may be terminated based on a decrease in quality of detected trajectories. Examples of termination criteria are provided below. Following termination of the iterative process, the method outputs a representation of detected trajectories, based on the output from the iterative process. For example, the representation of detected trajectories may include an indication of the one or more candidate trajectories that were selected during the iterative process. Additionally or alternatively, the representation of detected trajectories may include an indication of one or more sets of data points that were identified as corresponding to candidate trajectories during the iterative process. The representation of detected trajectories may include a graphical representation. For example, the representation may comprise a plot (e.g. spatial channel vs. time) of the detected trajectories. The set of candidate trajectories may comprise a mapping between each candidate trajectory and a corresponding set of data point coordinates. Thus, each candidate trajectory is associated with a set of data points in the DAS data having coordinates lying on or near (e.g.within a threshold distance) a representation of that candidate trajectory. For example, the setof candidate trajectories may comprise, for each candidate trajectory, a mapping between the candidate trajectory and the plurality of spatial channels. The mapping may comprise, for each candidate trajectory, an indication of one or more time coordinates for each spatial channel which correspond to that that candidate trajectory. The mapping may be pre-computed with the set of candidate trajectories.The step of identifying the set of data points may then be performed using the mapping. In thismanner, the pre-computed mapping can be used to facilitate and speed up the process of identifying the data points in the DAS data which correspond to the selected candidate trajectory. In particular, this reduces a run-time of the iterative process by avoiding having to map the selected candidate trajectory onto the DAS data during the iterative process. Thus, when a candidate trajectory is selected, a corresponding set of data points in the DAS data can simply be determined from the mapping. Identifying the set of data points may comprise determining end points and / or a line width of arepresentation corresponding to the selected candidate trajectory in the DAS data. In thismanner, the set of data points may be identified as being data points associated with candidate trajectory which are between the end points and / or within the line width. The end points and / or line width may be determined by analysing a signal intensity associated with data points lyingalong the candidate trajectory, for example using a segmentation process involving edge orboundary detection techniques. This may ensure that only data points which are part of atrajectory that is represented in the DAS data are part of the identified set. Where a mapping for the selected candidate trajectory is used, as mentioned above, the set of data points associated with the selected candidate trajectory in the mapping may first beidentified. Then, only data points located between the determined end points and / or within thedetermined line width are retained and identified as the set of data points corresponding to theselected candidate trajectory. The iterative process further may comprise determining whether the identified set of data pointssatisfies a predetermined condition and, if so, determining the selected candidate trajectory as adetected trajectory. In this manner, a check is performed on the identified set of data points to verify that they actually represent a real trajectory in the DAS data. This may serve to reduce false detections, thereby improving an accuracy of trajectory detection. For example, this may avoid random noise signals in the DAS data being incorrectly detected as a trajectory. The predetermined condition for determining that a trajectory is detected may include any suitable condition that serves to distinguish noise signals from representations of trajectories in the DAS data. For example, the predetermined condition may be based on one or more predetermined characteristics of trajectory representations in DAS data. Thus, if the identified set of data points satisfies the one or more predetermined characteristics, then they can be determined as representing a detected trajectory in the DAS data. Determining whether the identified set of data points satisfies a predetermined condition may comprise: determining if the contribution to the objective function arising from the identified set of data points in the distributed acoustic sensing data satisfies a predetermined uniformity condition; and if the predetermined uniformity condition is satisfied, determining the selected candidate trajectory as a detected trajectory. Thus, the predetermined condition may correspond to a uniformity condition of the set of data points. The inventors have realised that a trajectory is typically represented in DAS data by a set of data points having highly uniform signal intensities. Thus, by checking the set of data points against a predetermined uniformity condition (e.g. threshold), it is possible to verify whether the set of data points represent an actual trajectory or not. Determining if the uniformity condition is satisfied may comprise determining a measure of uniformity of the contribution to the objective function arising from each data point in the identified set (e.g. based on an intensity of each data point). The measure of uniformity can then be compared to a predetermined threshold, to determine if the uniformity condition is satisfied or not. For example, the measure of uniformity may comprise a comparison between an intensityof each data point and an average (e.g. a moving average) across the set of data points. If it isfound that a large proportion (e.g.75% or more) of the set of data points has a low intensity compared to the average, this may be indicative of a non-uniform distribution and the uniformity condition may not be satisfied. Only representations of candidate trajectories which are determined as detected trajectories may be included in the output of the method. The method may further comprise: following an iteration of the iterative process, determining if a re-normalisation condition is met and, if so, re-normalising the distributed acoustic sensing data based on a remaining portion of the distributed acoustic sensing data which does not include one or more sets of data points identified in the iterative process; and re-computing theobjective function based on the re-normalised distributed acoustic sensing data. The inventorshave found that, after multiple iterations of the iterative process, the objective function tends to flatten due to progressive notching of data points from the objective function, which can lead to less accurate detections. This is remedied by re-normalising the DAS data when the re- normalisation condition is met, and re-computing the objective function based on the re- normalised DAS data. Crucially, the re-normalised data does not include contributions from the previously notched data points (i.e. sets of data points which were identified as corresponding to selected candidate trajectories in previous iterations). In this manner, the re-computed objective function will have sharper peaks for the remaining candidate trajectories (i.e. candidatetrajectories which were not selected in previous iterations), which facilitates trajectory detection.The process of re-computing the objective function is done in a similar way to the initial computation of the objective function. As noted above, computing the objective function can be a computationally intensive process. Advantageously, by using a re-normalisation condition to determine when the objective function should be re-computed, the objective function need not be re-computed at each iteration, thus speeding up the iterative process. As an example, the re-normalisation condition may be met if there has been a predetermined number of iterations (e.g. four or more) of the iterative process since the last re-normalisation process. As another example, an efficient detection criterion (EDC) may be defined for the objective function, and the re-normalisation condition may be determined to be met based on the EDC. The step of updating the objective function may comprise: determining one or more affected candidate trajectories which map on to one or more of the identified set of data points; and re- computing components of the objective function corresponding to the affected candidate trajectories to remove a contribution from the identified set of data points. In this manner, only components of the objective function corresponding to candidate trajectories which are affected by the identified set of data points are re-computed, leaving the remainder of the objective function as it is. A candidate trajectory is determined to be affected if it maps on to (e.g. passes through) one or more of the identified set of data points. This can be determined using a pre- computed mapping between each candidate trajectory and a corresponding set of data point coordinates. The components of the objective function can be re-computed, for example, by subtracting from the objective function a signal intensity contributed from each data point in the identified set of data points. The objective function may comprise, for each candidate trajectory, a sum over the plurality of spatial channels of an error signal between the measurement signal and that candidatetrajectory. In other words, for each spatial channel, an error signal between the measurementsignal for that spatial channel and the candidate trajectory is determined (calculated). The error signal provides an indication of how well the candidate trajectory fits the measurement signal of that spatial channel. The error signals across all the spatial channels are then summed, to provide a likelihood of the candidate trajectory being represented in the DAS data. Computing the objective function may comprise performing a parallel computing process, where a contribution to the objective function associated with each spatial channel is computed in parallel, and then summed together. For example, the error signals for each spatial channel mentioned above may be computed in parallel, and then summed together. This takes advantage of a realisation that the contributions from each spatial channel are independent of one another, and so can be computed in parallel. Parallel computing for each of the spatial channels can significantly reduce a computing time for the objective function. The iterative process may further comprise determining if the identified set of data points satisfies a quality metric, and if the quality metric is not satisfied, determining that the termination condition is met if the quality metric was also not satisfied for the preceding Neiterations of the iterative process, where Ne is a predetermined number. In this manner, if thequality metric is not satisfied in the current iteration and the preceding Neiterations, the iterative process is terminated. The quality metric not being satisfied in the current iteration and the preceding Ne iterations may be an indication that all of the trajectories represented in the DASdata have been detected. This also avoids terminating the iterative process too early as a resultof a single ‘bad’ detection, with termination instead being based on a sequence of ‘bad’ detections. When the iterative process is terminated, the selected trajectories from the current iteration and the preceding Neiterations may be omitted from the detected trajectories output by the method. Any suitable quality metric may be used for assessing a quality of the detection. The predetermined number Nemay be any real number equal to or greater than 1. The quality metric for the identified set of data points may be determined based on acomparison with a preceding iteration of the iterative process. For example, the quality metricmay be configured to compare a quality of the set of data points identified in the current iteration with a quality of sets of data points identified in one or more preceding iterations. In this manner,a decrease in quality of detection results over multiple iterations can be determined, which maybe indicative of there being no further trajectories to detect in the DAS data.Computing the objective function may comprise a pre-processing procedure. The pre-processing procedure may include normalising each of the plurality of measurement signals.This may ensure a consistency of signal amplitude across the plurality of spatial channels, which may improve an accuracy of trajectory detection. For example, different spatial channels may have different noise and sensitivity levels, depending on local conditions along the optical fibre. Thus, normalising each spatial channel may serve to mitigate against differences in response of each channel. Any suitable normalisation process may be used. The pre-processing procedure may additionally or alternatively comprise performing noise reduction process on the distributed acoustic sensing data to improve a signal-to-noise ratio ofthe distributed acoustic sensing data. This may improve an accuracy of trajectory detection inthe DAS data. The noise reduction process may comprise: applying a contrast enhancement algorithm to the distributed acoustic sensing data; applying an edge detection algorithm to an output of the contrast enhancement algorithm; and subtracting an output from the edge detection algorithm from the output of the contrast enhancement algorithm. This process serves to sharpen and reduce a width of any trajectories represented in the DAS data, thus improving SNR and facilitating detection of the trajectories. Obtaining the set of candidate trajectories may comprise pre-computing the set of candidate trajectories, the set of candidate trajectories being pre-computed based on dimensions of a data window of the distributed acoustic sensing system, and one or more predefined trajectory parameters. In this manner, the set of candidate trajectories is specifically computed for the data window of the DAS system, which facilitates mapping the candidate trajectories onto data obtained from the DAS system. Thus, each candidate trajectory can be described in terms of coordinates in the data window of the DAS system. Here, the data window of the DAS system corresponds to a data window (sample) having a spatial dimension corresponding to the spatial channels of the DAS system, and a time dimension corresponding to a sample measurement period of the DAS system. The one or more predefined trajectory parameters correspond to parameters which are used to compute the candidate trajectories, so that the candidate trajectories are realistic trajectories that may be represented in DAS data. For example, the one or more predefined trajectory parameters may comprise a range of possible vehicle speeds and / or directions of travel relative to the optical fibre. These can then be used to calculate parameters such as slopes and offsets of the candidate trajectories. The method of the first aspect may further include steps performed with the DAS system to acquire the DAS data. For example, the method may include steps of launching a pulsed test signal along the optical fibre, receiving at a detector stage a plurality of scattered signals, each of which was scattered at a respective location along the optical fibre, and outputting, by the detector stage, DAS data comprising a plurality measurement signals, each measurement signal being a function of a respective one of the scattered signals. According to a second aspect of the invention, there is provided an analysis system for analysing distributed acoustic sensing data, the analysis system comprising a processing device, and a memory storing instruction which, when executed by the processing device, cause the processing device to perform the method of the first aspect of the invention. The analysis system is a computer-implemented system, which can be implemented using any suitable computer system or network of computer systems. The processing device may correspond to one or more computer processors which are coupled to the memory and configured to execute the instructions stored in the memory. The memory may comprise a non- volatile storage medium, such as a local memory drive and / or a cloud-based storage system. According to a third aspect of the invention, there is provided a distributed acoustic sensing system comprising: a pulse generator configured to transmit a pulsed test signal along an optical fibre; a detector stage configured to receive a plurality of scattered signals from the optical fibre, wherein each scattered signal corresponds to a respective spatial channel associated with a respective scattering location along the optical fibre, and wherein the detector stage is further configured to output a respective measurement signal as a function of each of the plurality of scattered signals; and an analysis system configured to receive the respective measurement signals from the detector stage, and to perform the method of the first aspect of the invention. The DAS system of the third aspect may be used to implement the method of the first aspect of the invention. Accordingly, features described above in relation to the first aspect may be shared with the third aspect of the invention (and vice versa). The analysis system in the DAS system of the third aspect may correspond to the analysis system described in relation to the second aspect of the invention above. Accordingly, any features described in relation to the second aspect may be shared with the third aspect (and vice versa). The pulse generator may comprise an optical modulator for generating the pulsed test signal for a received light signal. For example, the DAS system may comprise a coherent light source (e.g. laser), which provides a continuous wave light signal to the pulse generator, which generates the pulsed test signal. The pulse generator may be coupled to an end of the optical fibre to launch the pulsed test signal into the optical fibre. The detector stage may be coupled to the optical fibre to receive the scattered signals. The detector stage may comprise an optical detector for detecting the scattered signals. For example, a square law detector may be used. In some cases, the detector stage may further be configured to receive a local oscillator signal, and to interfere the local oscillator signal with the received scattered signals on the detector. In this manner, a signal output by the detector may be indicative of an interference between the local oscillator signal and the scattered signals. The invention includes the combination of the aspects and preferred features described except where such a combination is clearly impermissible or expressly avoided. Summary of the Figures Embodiments and experiments illustrating the principles of the invention will now be discussed with reference to the accompanying figures in which: Fig.1 shows a schematic diagram of a distributed acoustic sensing (DAS) system and analysis system according to an embodiment of the invention; Fig.2 shows a schematic diagram of an example use of the DAS system for detecting a vehicle trajectory; Fig.3 shows an example graph of a vehicle trajectory represented in DAS data; Fig.4 shows an example waterfall plot of DAS data showing multiple vehicle trajectories; Fig.5 shows a diagram of a method of detecting vehicle trajectories in DAS data according to an embodiment of the invention; Fig.6 shows plots of notch power distribution across spatial channels of the DAS system, for different iterations in the method of Fig.5; Fig.7 shows examples of DAS data in different steps of a pre-processing procedure which can be included in the method of the invention;Fig.8 shows example plots of an objective function across successive iterations in the methodof Fig.5; and Fig.9 shows an example plot of candidate trajectories that may be used in the method of the invention. Detailed Description of the Invention Aspects and embodiments of the present invention will now be discussed with reference to the accompanying figures. Further aspects and embodiments will be apparent to those skilled in the art. All documents mentioned in this text are incorporated herein by reference. Fig.1 shows a schematic diagram of a distributed acoustic sensing (DAS) system 10, according to an embodiment of the invention. The system 10 is arranged to interrogate an optical path, inparticular an optical fibre 1000, which may be of any desirable length for a given purpose. Thesystem 10 comprises a light source which produces coherent light, which is given here as a laser 12, and is used in continuous wave (CW) operation. The light produced by the laser 12 can be directed into an optical isolator to ensure that light is not passed back to the laser 12. Light from the laser 12 is split into two paths by an optical coupler 16 or beam splitter. The firstpath, from which light is directed into the fibre 1000 is known as the launch path. The secondpath, from which light is passed directly to a detector stage 50 (discussed below), is known as the local oscillator path. The light is split between the two paths by the optical coupler 16, for example such that 90% of the incoming light is directed into the launch path, and 10% of the incoming light is directed into the local oscillator path. Of course, the ratio of incoming light directed into each path may be chosen by the operator depending on the nature of the operation for which the OTDR system 10 is used. The laser light which is directed into the launch path then passes through a pulse generator 18, such as an acousto-optic modulator (AOM) 18. The AOM 18 is a device which can simultaneously generate an optical pulse as well as upshift or downshift the frequency of light by an amount equal to the radiofrequency which drives the AOM 18. This frequency shift, F, may be known as the intermediate frequency or the difference frequency. In this way, the AOM 18 is able to generate a pulsed test signal which may be between 5 ns and 100 ns in duration, but not limited to this range. Of course, any preferred method of generating a pulse of light may be used, such as an electro-optic modulator (EOM). The pulsed test signal may also be referred to herein as a launch pulse. The pulse of light can then be amplified using an optical amplifier. The light pulse is thenintroduced to the optical fibre 1000 via an optical circulator 22, which has three ports. Theamplified light pulse enters the circulator 22 through a first port, where it is passed to a secondport in order to enter the optical fibre 1000. As the pulse of light passes through the fibre 1000,a fraction of the light is backscattered from the fibre 1000, e.g. by Rayleigh scattering, and afurther fraction captured and guided back towards the circulator 22. The scattered light, which may be referred to herein as a scattered signal, enters the circulator 22 at the second port, and leaves the circulator 22 to enter a detector stage 50 via a third port. The detector stage 50 has two inputs. The first input is the scattered laser light from the third port of the circulator 22. The second input is the laser light taken directly from the local oscillator(LO) path mentioned above. The scattered light is then mixed with the LO light at an opticalcoupler 28. The light the optical coupler 28 is then allowed to interfere on an optical detector 30 (e.g. square law detector). An output signal from the detector is then taken and measured at an analog-digital-converter (ADC) 32, which outputs a corresponding measurement signal. The measurement signal can then be further analysed, e.g. to determine a strain and / or acoustic field at the corresponding scattering location along the optical fibre 1000. It will be appreciated that the invention is not limited to use of the specific DAS system shown in Fig.1, and that various alternative systems may be used. For example, additional or alternative components may be included in the DAS system. In some cases, a polarisation diverse detector stage 50 is used, e.g. where the scatter signal and the LO signal are each split into horizontal and vertical polarised states, to enable polarisation diverse detection. The system 10 described above can make use of a heterodyne sensing approach, wherein the frequency of the local oscillator and of the launch pulse are shifted relative to one another by the AOM 18. The difference in these two frequencies should be larger than the bandwidth required to represent the scattering without allowing crosstalk between the carrier and the DC terms which are also generated, allowing the phase and amplitude information of the scattering to be recovered using a real carrier. Another method employs a complex carrier detector stage, replicating the polarisation diverse detector stage for two copies of the local oscillator shifted by 90 degrees relative to each other. This allows detection via a complex carrier, allowing either the positive sidelobe or the negative sidelobe of the resulting interference signal to be recovered independently. This allows homodyne operation whereby the local oscillator signal and launch pulse operate at the same optical frequency. As the pulsed test signal travels along the optical fibre 1000, the test signal will be scattered at a plurality of scattering locations distributed along a length of the optical fibre 1000. Thus, for each pulse launched along the optical fibre, a plurality of scattered signals with be received at the detector stage 50, each of the scattered signals corresponding to a respective scattering location in the optical fibre. The scattered signals are received sequentially in time, with the time of receipt for a given scattered signal depending on its scattering location in the optical fibre 1000. Thus, the scattering location corresponding to a received scattered signal can be determined based on its time of receipt at the detector stage 50, e.g. comparing a time at whicha pulse of the test signal was launched along the optical fibre 1000 and the time of receipt,taking into account the speed of light along the optical path. Accordingly, the detector stage 50 can output a respective measurement signal for each of the plurality of scattering locations along the optical fibre 1000, the measurement signal being derived from interference of the scattered signal for that scattering location with the local oscillator signal. Each of the scattering locations in the optical fibre 1000 for which a measurement signal is output by the detector stage 50 may be referred to as a spatial channel of the system 10. Returning to Fig.1, an analysis system 60 is connected to an output of the detector stage 50, to receive DAS data comprising the plurality of measurement signals output by the detector stage 50, i.e. by the ADC 32. The analysis system 60 may be directly coupled to the ADC 32 as shown in Fig.1. Alternatively, the analysis system 60 may be arranged to receive the plurality of measurement systems from the detector stage 50 via a computer network or intermediate storage system which is configured to store data output by the detector stage 50. The analysis system 60 is configured to analyse the DAS data received from the detector stage 50, in order to detect trajectories represented in the DAS data, as discussed in more detail below. The analysis system 60 is a computer-implemented system which includes a memory 62 that stores computer instructions, and a processor 64 configured to execute the instructions to analyse the received DAS data. In practice, the analysis system 60 can be implemented using any suitablearrangement of computer hardware, such as a personal computer (e.g. desktop or laptopcomputer), computer server (e.g. cloud-based server), or a network of computing devices. Fig.2 illustrates a use of the DAS system 10 for detecting a vehicle trajectory along a vehicle path. The optical fibre 1000 extends along the path, e.g. under or next to a surface 204 of the path. As a vehicle 206 travels along the path, it generates an acoustic wave which propagatesin the ground to the optical fibre 1000. In the example shown, the vehicle 206 is a road vehiclemoving along a road or a path. However, the method is equally applicable to other types of vehicle, e.g. a train moving along a track. The vehicle 206 acts as a moving source for the acoustic wave, such that progress of the vehicle 206 along the path can be detected by monitoring arrival of the acoustic wave in the spatial channels of the optical fibre 1000. An example of this is shown in Fig.3, which is an example graph of measurement time as a function of spatial channel in the DAS system 10. The vertical dashed lines in Fig.3 indicate locations of spatial channels that are measured by the DAS system 10. The points 300 in Fig.3 indicate an arrival time of the acoustic wave generated by the vehicle 206 in each spatial channel. As can be seen, an arrival time of the acoustic wave increases for increasing spatial channel number, indicating that the vehicle 206 is moving along the path (and hence along the optical fibre 1000), e.g. in a positive direction. The arrival times of the acoustic wave follow a line 302, which is indicative of a speed and direction of motion of the vehicle 206 relative to the optical fibre 1000. When the acoustic wave reaches a spatial channel, this will result in amodulation of strain in the optical fibre 1000 at that spatial channel (scattering location),resulting in a corresponding response (e.g. modulation or variation) in the measurement signalcorresponding to that spatial channel. Thus, the trajectory of the vehicle (i.e. position as afunction of time) is represented in the DAS data as a sequence of responses in the plurality ofspatial channels.Returning to Fig. 2, an acquisition rate of the system 10 depends on a length L of the sensor, asthe pulse and the related reflections have to reach the detector stage 50 before a new pulse canbe sent without interferences. A typical maximum data rate for a 50 km system is ^^ = 1 kHz,although shorter systems enable a higher acquisition rate. The spatial resolution d (shown inFig.2) between adjacent spatial channels is determined by the duration of the transmitted pulse, which defines the amount of light fired into the fibre at each time step ^. The spatial channels (sensing units) are equally-spaced along the length of the optical fibre 1000, and eachspatial channel is defined by a channel index ^ = ^^⁄ ^ , ^ ∈ {1, 2, … , ^}, where ^^ is the distanceof the spatial channel from the laser 12 and ^ = ^⁄ ^ is the number of spatial channels in thesystem. At each time step ^, the DAS system 10 measures the acoustic energy at each location ^^∈[^, ^]. The data recorded in a time window (^, ^ + ^) of N data points is usually represented in a‘waterfall’ plot as shown in Fig.4, where each column, indexed by ^, is a measurement signal of a respective spatial channel. The measurement signal for each spatial channel may be given as: The entire waterfall dataset Y represents the ^ channel signals in the time window, ^ =[^^, … , ^^], with data point (^, ^) corresponding the value of ^^(^ + ^) and ^ ∈ [1, ^]. In thiscontext, bottom rows in Fig.4 represent older data points, while newer data is displayed at thetop of the waterfall, moving towards the bottom part as new data is collected.The waterfall DAS data may be assumed to contain ^ trajectories which can be approximatedas straight-line trajectories in short time windows. Each line (i.e. trajectory) may be defined by aparameter ^^ = (^, ^) ∈ ^, where ^ is a parameter space, ^ is a slope of the line and ^ is anoffset of the line. The entire set of trajectories represented in the DAS data is described by thevector of parameters ^ = [^^ , … , ^^]^, which the method of the invention aims to estimate. The measurement signal of each spatial channel potentially contributes to representing some of the trajectories depicted in the DAS data. Therefore, a measurement signal for a spatial channel^^ can be represented as a linear combination of ^ + 1 basis functions ^^,^ associated with thespatial channel ^ and defined by ^:(2) (3)where ^^,^ is a ^ × (^ + 1) matrix, ^^ = [^^,^ , ^^,^ , … , ^^,^]^is an unknown vector of amplitudes,and ^^ ∈ ℝ^is a noise component, such as zero-mean Additive White Gaussian Noise(AWGN) with covariance matrix ^^^^×^. The basis function of the parameter ^^ is a columnvector of ones, ^^^,^ = 1^×^, representing the average intensity ^^,^ of the measurement signal^^ over the time window. For each ^^ , ^ ∈ [1, ^], the basis function ^^^ ,^ represents the trajectorydefined by the parameter ^^ at the distance ^^ in the chosen time window. Each basis functioncan be defined as a linear combination of standard basis ^^in an ^-dimensional space of the measurement signals:where ^^ ∈ [1, ^] is a time instance at which the trajectory is detected (represented, depicted) inthe spatial channel ^ and ^ = [0.5^^^^] − 1 defines the duration of the trajectory in terms of datapoints, with ^^the shortest duration of the trajectory signal across all channels. In other words,each basis function ^^^,^ is defined by assigning a value of 1 to the time point at which thetrajectory signal is detected, and a value of 0 for all other time points. Of course, the invention isnot limited to use of the specific basis function shown in (4) and equivalent or alternative basis functions could be used to model trajectory signals. The waterfall DAS data is composed of a collection of measurement signals, each of which may be affected by different noise levels, which depend local conditions at the spatial channels along the optical fibre 1000, such as fibre-soil coupling and the presence of external noise sources. The insertion of the basis function ^^^,^for the measurement signal of each spatial channel acts to normalise the DAS data across the plurality of spatial channels. Fig.5 shows a flow diagram of a method 500 according to an embodiment of the invention. The method 500 is for detecting vehicle trajectories represented in DAS data. For example, the method 500 can be implemented by the analysis system 60 described above, in order to analyse the DAS data received from the detector stage 50. For convenience, the method is described in the context of the examples of Figs.1 and 2, although it will be appreciated that the method 500 can be used with other arrangements of DAS system. The method comprises a step 502 of obtaining a set of candidate trajectories corresponding to possible trajectories represented in data from the DAS system 10. The set of candidatetrajectories is referred to as a dictionary ^ comprising a set of parameters ^^ = (^, ^) definingcandidate trajectories that could potentially be represented in DAS data from the system 10.Accordingly, the dictionary ^ is computed so as to provide a substantially complete definition ofrealistic trajectories. The dictionary ^ is computed based on dimensions of a data window (i.e.of the waterfall dataset) from the DAS system, in other words based on the number of spatialchannels and number of data points in each spatial channel. Additionally, a set of trajectory parameters can be used to compute the sets of parameters in the dictionary ^, to ensure that the candidate trajectories are realistic. In particular, the trajectory parameters include a range of possible vehicle speeds and directions relative to the optical fibre 1000. Then, iterating throughdifferent possible combinations of vehicle speeds, directions and positions, the dictionary ^ iscomputed.As an example, the candidate trajectories in the dictionary ^ can be computed on a canvas,which is a space encompassing all of the data points as the DAS data window (i.e. the waterfalldataset), but with a larger number of spatial channels than the DAS data window. A first set oflines with parameters (^, ^ = 0) can then be determined within a waterfall-sized window of thecanvas, corresponding to candidate trajectories passing through the origin. The waterfall-sizedwindow can then be shifted left and right along the spatial channel axis, such that the first set oflines provides a representation of candidate trajectories (^, ^) where ^ is the shift applied to thewindow.As noted above, the candidate trajectories are defined in a continuous parameter space (^, ^),where ^ and ^ are real numbers. In order to compute the dictionary ^, the parameter space(^, ^) can be discretised according to a resolution ∆^ and ∆^, which can be set by a user, inorder to represent a finite number of candidate trajectories in the parameter space. A density (orresolution) of the dictionary ^ depends on the chosen values of ∆^ and ∆^, e.g. smaller valuesof ∆^ and ∆^ yield a higher density (resolution) dictionary, and vice versa. In practical terms,the density of the dictionary ^ may be set to cover the entire waterfall space with candidatetrajectories. For instance, setting ∆^ to 3 data points (pixels) or less, and ∆^ such that thespeed resolution is 3 spatial channels per data point (pixel) or higher. Note that as ^ isproportional to the inverse of speed, ∆^ is not a constant.Additionally, for each candidate trajectory ^^in the dictionary ^, a mapping between the candidate trajectory and data point coordinates in the waterfall are computed. This can be achieved, for example, by mapping (projecting) the candidate trajectory onto the waterfall, todetermine coordinates (^, ^) of data points lying on the candidate trajectory, or within apredetermined distance of the candidate trajectory. The dictionary ^, together with the mapping of the candidate trajectories is stored in a memory, e.g. the memory 62 or another dedicatedmemory. Then, at the start of the method 500, the dictionary ^ and mapping can be obtained byaccessing the memory. The method 500 further comprises a step 504 of receiving DAS data from the DAS system 10. In line with the discussion above relating to Fig.1, a pulsed test signal is launched along the optical fibre 1000. This results in a plurality of measurement signals being output by the detectorstage 50, each measurement signal corresponding to a respective spatial channel ^ of thesystem 10. Accordingly, in step 504, the method 500 comprises receiving (e.g. by the analysis system 60), DAS data comprising the plurality of measurement signals as a function of time. Inother words, a measurement signal ^^ (equation (1)) is received for each spatial channel. Forexample, the DAS data can be received as a data stream, i.e. as it is being output by the detector stage 50. The received DAS data can be stored (or cached) in a memory of theanalysis system 60, such as the memory 62 or another dedicated memory. Alternatively, themethod 500 may be performed ‘offline’, i.e. on DAS data which was measured at an earlier time and stored in the memory 62 for later analysis. Note that the step 502 of obtaining the candidate trajectories may be performed before, during (e.g. in parallel with), or after the step 504 of receiving the DAS data. At step 506, the method comprises computing an objective function based on the DAS data received in step 504 where the objective function is indicative, for each of the candidate trajectories in the dictionary ^, of a likelihood of that candidate trajectory being represented inthe DAS data. The objective function is computed by comparing each of the candidatetrajectories from the dictionary ^ to the DAS data to assess a fit between the candidatetrajectory and the DAS data, in order to determine a likelihood that the candidate trajectory isrepresented in the DAS data. A detailed example of a derivation of the objective function isprovided below. Objective function derivationAssuming that the number of trajectories ^ in the DAS dataset is known, we consider that theset of unknown parameters is Θ^ = where ^^:^ and ^^^:^are, respectively, sets ofamplitudes and noise variances over the ^ spatial channels. Given an Additive White GaussianNoise (AWGN) assumption, the logarithm of the likelihood can be formulated as (Ref.1): where ^^is a constant term. In order to determine a maximum likelihood estimation (MLE),equation (5) can be maximised with respect to the parameter set Θ^. Fixing ^ and ^^:^,maximisation with respect to ^^^:^can be performed for each spatial channel (Ref.1). Substituting (6) into (5), the following is obtained: where is the sum of the log power of estimated error signals ^^^ (. )For a given ^, (9)-(10) can be minimised with respect to ^^:^ by solving a least log squaresproblem for each spatial channel ^: Substituting (11) in (9), the likelihood can be reformulated in terms of ^ for the optical ^^ ^:^ and^^^^:^as where the log power of the estimated error signal is: (13)Introducing the projection operators onto the span (^^,^) and its orthogonal complement,respectively (Ref.2): we can rewrite the log power of the estimated error signal in (13) as Therefore, the parameter estimation is performed by minimising (12) over the parameter space ^^. Every spatial channel signal ^^makes an equivalent contribution to the estimation of theglobal parameter ^, which represents a set of lines (trajectories) spanning through the ^ spatialchannels in the waterfall dataset ^.Now the assumption is taken that the number of vehicles ^ is unknown, and the parameters ofprevious estimations ^^ = ^^^, … , ^^^^^ ^, ^ − 1 < ^, are known at the beginning of iteration ^. Aniterative minimisation of (12) is performed to infer an unknown parameter ^^^ ∈ ^ resulting in anestimated parameter vector equal to the augmented set ^^ = ^^^ , … , ^^^ ^. ^^ can be referred toas a notch vector, with ^ corresponding to the dictionary of candidate trajectories as describedabove.A maximum likelihood (ML) criterion for a model with ^ parameters can be given by (Ref.3): (17) (18) where is defined as a notch power of the measurement signal ^^with respect to the notch vector ^^,and ^^ ∈ ^. The term ^^∥^^^,^is a projector onto a space defined by the residual vector of g^^,^onto span In order to consider the varying lengths of trajectories, the notch power in (18) is regularisedand the objective function ^^(. ) is defined as: where ^^(^^) is the number of spatial channels in which the trajectory ^^ is supposed to exist.Thus, the objective function provides, for each candidate trajectory ^^in the dictionary ^, a sumover the plurality of spatial channels ^ of the error signal ^^(^^) (equation (17)) between themeasurement signal and that candidate trajectory. Returning to Fig.5, at step 506 the objective function can be computed for the DAS data usingthe formula (21) above. Initially, we set ^^ = ∅, assuming there have not been earlier trajectoryestimations on the dataset. It should be noted that equation (21) above is provided as a specific example of an objective function that can be used in the method of the invention. However, the method is not limited to this specific form of the objective function and other equivalent formulations of the objective function can be used. The objective function can be computed using parallel computing (processing) techniques. In particular, as the contribution from each spatial channel is independent from other spatial channels, the terms in equation (18) corresponding to different spatial channels can be computed in parallel, and then summed together to compute the objective function (21). Such parallel computing may significantly reduce a computing time for the objective function. Following computing of the objective function, the method 500 involves an iterative process for detecting trajectories in the received DAS data. The iterative process first involves a step 508 ofselecting one of the candidate trajectories from the dictionary ^ having the highest likelihood ofbeing represented in the DAS data. This can be done by maximising the objective function (21): (22) Next, in step 510, a set of data points in the DAS data corresponding to the estimatedparameter (i.e. the selected candidate trajectory) ^^^ is identified. This may be referred to as asegmentation process. The parameter ^^^ provides a minimal representation of a vehicletrajectory defined by a set of basis functions {g^^^,^}^^^^. In order to effectively notch the contribution from the selected trajectory in the DAS data, the minimal trajectory representation is refined by estimating boundaries and endpoints of the trajectory in the DAS data. From this, an extended representation of the vehicle trajectory can be obtained, defined by a set of basis functions {g^^^^,^}^^^^.In one example, given that the minimal representation of the estimated trajectory ^^^ is defined interms of slope as ^^^ = (^^^^ , ^^^^ ), ordered sets of offsets are introduced:(23) (24) both representing a subset of lines in ^ parallel to ^^^ in a range of offsets defined by ^^, where^^represents lines with a larger offset than ^, a^^^^ ^nd ^^^^represents lines with a smaller offsetthan ^^^ . Offsets on which boundaries of the trajectory represented in the DAS datalie can then be estimated. For instance, evaluating the objective function ^^(. ) on the orderedsets ^^^^^and ^^^^^, the offsets can be defined according to a first order derivative of^^(. ) with respect to the offset ^ and a multiplicative coefficient ^^: where Δ^ is an offset step size in ^^^ ^^^and ^^^^. Data points corresponding to the selectedcandidate trajectory ^ then can be determined to lie between the lines ^ ) and (^ ^ , ^ ^)Accordingly, the extended representation of the trajectory ^^^ can be expressed in each channelby a new set of basis functions: where ^^^ , ^^^ ∈ [1, ^] are time instants at which the lines and (^^^^ , ^ ^^^^ ) intersectspatial channel ^, respectively. In this manner, data points corresponding to the selectedcandidate trajectory ^^^ are identified as those points lying between the parallel lines (^^^^^, ^^^^)and (^^^^^, ^^^^). A correspondent notch power for the selected candidate trajectory defined on {g^^^^^,^}^^^ is denoted as ^^^,^(^^^ ; ^^).In addition to the above, the segmentation process in step 510 can include determining end points of the representation of the selected candidate trajectory. The representation of theselected candidate trajectory can for example be set to zero for points lying beyond thedetermined end points. So far, it has been assumed that a vehicle trajectory spans for an entire length of the time window of the DAS data and between two channels ^^^ , ^^^ ∈ {1, … , ^}, ^^^ < ^^^. Locations at which the vehicle trajectory appears and disappears in the DAS data can beestimated by analysing ^^ (^^^ ; ^ ) along the trajectory for a set of spatial channels ^^^ = {^^^^ , ^^^ ^ ∶ ≠ 0}. The notch power for ^ ∈ ^^^ can then be compared with apredetermined threshold Γ^ = 0.1^(^^ , ^^ ^^ , ^^ , ^^^), where the second term corresponds to mean notch power observed between the channels ^^and ^^^ , a cardinality of the set ^ ^ . Starting from^^^and ^^^and proceeding towards a centre (middle) of the trajectory, an iterative approach can be employed to determine the first channels ^^ ^^ and ^^ ^ ^^ with ^^ > ^^ ^^ and ^^ ^^ < ^^^, such that thevalue of ^ (^^^ ; is greater than Γ^. The end points [^^ ^^ , ^ ^^ ] of the trajectory can then beby the algorithm. Having identified the set of data points corresponding to the selected candidate trajectory in step 510 by determining the extended representation {g^^^^,^}, the method 500 subsequently involves a step 512 of notching the identified data points from the objective function. In other words, a contribution to the objective function arising from the identified set of data points isremoved from the objective function. To do this, the selected candidate trajectory ^^^ is added tothe notch vector ^^, i.e. ^^ = [^^ , … , ^^^ ]. The objective function then needs to be updated inaccordance with this latest iteration of the notch vector, in accordance with equation (21). In order to avoid re-computing the entire objective function, only portions of the objective function which are affected by the notched data points are re-computed. The set of basisfunctions {g^^^^,^} defines data points in each spatial channel corresponding to the selectedcandidate trajectory. Accordingly, using the mapping determined for the dictionary ^ computedin step 506, it is possible to determine each of the candidate trajectories in the dictionary whichis affected by the notched data points {g^^^^,^}, for example by checking which of the candidatetrajectories map on to data points defined by {g^^^^,^}. Thus, if a candidate trajectory maps onto coordinates of one or more data points defined by {g^^^^,^}, then that candidate trajectory is determined to be affected by the notching. Once the affected candidate trajectories are determined, the objective function is re-computed only for those affected candidate trajectories,using equation (21) above. In this manner, only the likelihood for affected trajectories is re-computed, with the remainder of the objective function remaining the same. The method 500 subsequently includes a step 514 of determining if a predetermined termination condition is met. A detailed example of a termination condition for the iterativeprocess is provided further below. If the predetermined termination condition is met (‘Yes’), thenthe method 500 goes to step 516, and a representation of trajectories detected in the iterative process is output. For example, the representation can include an indication of candidate trajectories that were selected in the iterative process, and / or of data points in the DAS datawhich were identified as corresponding to selected candidate trajectories. For instance, themethod may output a plot of the detected trajectories, and / or highlight the detected trajectories in the original DAS data. If at step 514 it is determined that the termination condition is not met (‘No’), then the method proceeds to step 518. At step 518, it is determined if a re-normalisation condition is met and, if so (‘Yes’), the method proceeds to step 520 where the DAS data is re-normalised and the objective function is re-computed based on the re-normalised data. Following step 520, the method returns to step 508 and the iterative process (i.e. steps 508, 510, 512, 514) is performed again using the re-computed objective function. A detailed example of the re- normalisation condition and re-normalisation process is provided below. If at step 518 it is determined that the re-normalisation condition is not met (‘No’), then the method 500 returns to step 508 to repeat the iterative process again using the objective function as updated (notched) in step 512. As contributions to the objective function arising from trajectories detected in previous iterations of the iterative process, a new trajectory in the DAS data can be detected at each iteration. Uniformity condition The method 500 can optionally check that the candidate trajectory selected in step 508corresponds to an actual vehicle trajectory represented in the DAS data. The objective function(21) sums the contribution of all the spatial channels in the DAS data. In some scenarios, fewhigh-intensity data points along with many low-intensity data points could lead to the generationof false detections. These false detections can be characterised by a highly non-uniform notchpower along the vehicle signal, as shown in Fig. 6. Panel (a) of Fig. 6 shows examples of notchpower distribution across the spatial channels (‘Channel space’) for two true detections for ^ ∈{1,2}, and panel (b) of Fig.6 shows examples notch power distribution across the spatialchannels for a false detection for ^ = 10. As can be seen, the true detections in panel (a) arecharacterised by a much more uniform notch power distribution, whilst the false detection of panel (b) has a highly non-uniform notch power distribution.To minimise false detections, following step 510 (and prior to step 512), the distribution of notchpower along the spatial channels between ^^ ^^ and ^^ ^^ is evaluated. A moving average of thenotch power with a window size of ^^spatial channels is considered: where ^^ is an odd integer. The dashed lines Fig. 6 represent the notch power based on theextended representation, whilst the solid lines depict the moving average of the notch powerwith ^^ = 5. The set of channels not satisfying a uniformity condition is defined as: If the condition holds, more than 75% of the channels show very low notch power compared to the average,and the power distribution along the trajectory is classified as non-uniform. The parameter ^^^ isadded to the notch vector ^^, leading to the notching of the associated trajectory power from thewaterfall according to the set of basis functions {g^^^^,^}. However, the corresponding estimatedtrajectory is omitted and does not become part of the set of returned trajectories. Note thehorizontal lines in Fig.6 express the right-hand side of the inequality in (30). Termination condition A termination condition is defined to check if the iterative process should be terminated in step 514 of the method 500. In general terms, and evolution of a quality of trajectories detected over multiple iterations is monitored. If the quality of the detected trajectories is seen to degrade over multiple iterations, then this may be a sign that there are no more trajectories left to be detected in the DAS data, and the iterative process should be terminated. The sum of the notch power on the minimal and extended representations can be defined as: respectively. Two criteria are used to determine when to terminate the iterative process. Thefirst criterion compares the value of ^(. ) with the previous estimations. If the condition is satisfied, the candidate trajectory ^^^ is considered as a candidate to be the last estimation,is a predetermined parameter. The second criterion considers ^^(. ) normalised by theextension of the trajectory, as expressed by ^^^^^^^. If the absolute value of the differencebetween the normalised ^(. ) of subsequent detections is below a predetermined threshold Γ^^ then the estimated trajectory ^^^ is inserted into a set ^^ of detections meeting the terminationcriteria (34), (35).The termination condition is determined to be met at step 514 if ^^^ of the current iteration isinserted into ^^and if the previous ^^detections are in ^^as well, where ^^is a predetermined number. Once the termination condition is determined to be met, the method 500 proceeds to step 516. At step 516, all of the detected trajectories are output, except for the detections included in ^^. Thus, the output detected trajectories include the trajectories defined in the notch vector ^^, omitting trajectories included in the set ^^and omitting any trajectories which were determined not to meet the uniformity condition (31) mentioned above. Pre-processing Pre-processing and normalisation steps may be performed on the received DAS data, following receiving the DAS data in step 504 and prior to computing the objective function in step 506. The pre-processing step aims to increase the signal-to-noise (SNR) ratio of the received DAS data and partially reduce the width of the vehicle trajectories represented in the DAS data, to facilitate trajectory detection. Different stages of the pre-processing procedure are shown in Fig. 7, where panel (a) shows the raw received DAS data, plotted as a waterfall as in Fig.4. The pre-processing procedure first involves applying a contrast enhancement algorithm to thereceived DAS data. For example, a Contrast-Limited Adaptive Histogram Equalization (CLAHE)method (Ref. 4) is applied to the received DAS data to increase its contrast. Panel (b) of Fig. 7shows the result of applying the CLAHE method to the waterfall of panel (a). The CLAHEmethod enhances the core part of the vehicle signals, as well as their blurred edges and somenoisy components. Next, an edge detection algorithm is applied to the output of the contrastenhancement algorithm, to detect edges of the vehicle trajectories in the data. For example, aLaplacian of Gaussian (LoG) filter can be applied to detect and cancel the blurred edges of thevehicle trajectories (see e.g. Ref.5). The Gaussian filter smooths the waterfall, while theLaplacian operator highlights the regions of rapid intensity change, i.e. the edges of the vehiclesignals. Panel (c) shows the result of applying the LoG filter on the contrast-enhanced data ofpanel (b). Following detection of the edges, the blurred edges are removed by subtracting the output from the edge detection algorithm (i.e. the detected edges) from the output of the contrast enhancement algorithm. Panel (d) of Fig.7 shows the result obtained by subtracting the result in panel (c) from the contrast-enhanced data of panel (b). Following this pre-processing procedure, the DAS data has a higher SNR, with more sharply defined trajectories. This pre-processed DAS data can then be used for the rest of the method500. The pre-processed DAS data may then be normalised, e.g. by insertion of the basisfunction ^^^,^for each measurement signal as discussed above. Re-normalisationFig.8 illustrates an example objective function ^^(^^) evaluated in the first four iterations of theiterative process in method 500. In each panel, the horizontal axis represents an index number ^ for the candidate trajectories ^^in the dictionary ^. Thus, the vertical lines plotted in the panels of Fig.8 represent a likelihood associated with a respective candidate trajectory in the dictionary. The markers 802, 804, 806, 808 in each panel of Fig.8 represent a maximum of^^(^^) in each iteration, and hence indicate the candidate trajectory which is selected in thatiteration (in step 508). As discussed above in relation to steps 510 and 512, following the selection of a candidate trajectory, its contribution to the objective function is notched. Thus, for example, in iteration 2 of the example in Fig.8, the contribution arising from the candidate trajectory at marker 802 is notched and no longer appears in the objective function. The same applies for subsequent iterations, with the contributions of preceding iterations being notched. As a result of this progressive notching of contributions from the objective function, this leads to a gradual reduction in sharpness of the objective function ^^(^^), as illustrated by the progression from iteration 1 to iteration 4 in Fig.8. The progressively reduced sharpness of the objective function may result in less accurate trajectory detection as the iterative process progresses. To allow high detection efficiency in each iteration, step 518 checks if a re-normalisation condition is met. If a re-normalisation condition is met (i.e. because the objective function is no longer sufficiently sharp to yield efficient detection results), then the DAS data is renormalised based on portions of the DAS data which have not yet been notched. Then, the objective function is re-computed based on the re-normalised DAS data, e.g. based on (21).By way of example, an Efficient Detection Criterion, EDC (Ref. 6) can be used to determine ifthe re-normalisation condition is met. Taking ^ as the number of iterations since the last re-normalisation, including ^^^ . Where there have not yet been any re-normalisations, ^ correspondsto the number of iterations since the initial normalisation of the DAS data. A version of the EDC is defined as: (36)where ^ is a coefficient that defines a penalty term. If the criterion (36) after the last detection ^is greater than the criterion after the last re-normalisation (37)then the re-normalisation condition is met. Re-normalisation can then be performed at step 520by redefining the basis function ^^^,^on the diagonal of the projection operator ^^^^^,^ , ^^ = where ^^ is the ^-th unit vector of the standard basis. This operation enables the renormalisationof the waterfall based on the part of the channel signals that has not been notched yet. Then,following re-normalisation of the DAS data, the objective function is re-computed using the re- normalised DAS data. The iterative process then returns to step 508, which is performed using the re-computed objective function. Full re-computing of the objective function may be necessary after each re-normalisation of the DAS data, as the re-normalisation involves re-scaling of the intensity of numerous measurement signals. Efficiency of the full re-computing of the objective function can be enhanced on the basis that several terms in the ML criterion (equation (18)) are common to several candidatetrajectories ^^ in the dictionary ^. As the number of possible basis functions ^^^ ,^ is finite, thechannel notch powers can be pre-computed and updated after each-renormalisation. Inparticular, as each basis function represents a portion of a candidate trajectory in a respective spatial channel, the value of the notch power for each basis function can be assigned to one or more corresponding candidate trajectories, to compute the objective function as in (21). As the number of basis functions is typically significantly smaller than the number of candidate trajectories, re-computing the objective function in this manner results in a lower number ofcomputing operations. Thus, in some implementations, computing (and / or re-computing) theobjective function may comprise determining, for each basis function, one or more candidate trajectories that are described in the corresponding spatial channel by that basis function. Then, the notch power for that basis function is computed, and used to compute the objective function for each of the one or more candidate trajectories. An example of this is illustrated in Fig.9, which shows a plot of three candidate trajectories ^^, ^^, ^^. In line with the above, eachcandidate trajectory can be described by a linear combination of basis functions: ^^ ^= ^^^^ ,^^^^^; As can be seen in Fig.9, the trajectories ^^, ^^, ^^intersect inthe spatial channel ^ = 1, such that the basis function in spatial channel ^ = 1 is the same for allthree trajectories, i.e. ^^^ ,^ = ^^^ ,^ = ^^^ ,^ = ^^,^. Accordingly, the notch power ^^,^(^, ^^)corresponding to basis function ^^,^can be computed, and used for each of the trajectories ^^,^^, ^^ in equation (21) to compute the objective function. In this manner, instead of having tocompute a separate notch power for each of the three trajectories in the spatial channel ^ = 1,only a single notch power need be computed. The features disclosed in the foregoing description, or in the following claims, or in the accompanying drawings, expressed in their specific forms or in terms of a means for performing the disclosed function, or a method or process for obtaining the disclosed results, as appropriate, may, separately, or in any combination of such features, be utilised for realising the invention in diverse forms thereof. While the invention has been described in conjunction with the exemplary embodiments described above, many equivalent modifications and variations will be apparent to those skilled in the art when given this disclosure. Accordingly, the exemplary embodiments of the invention set forth above are considered to be illustrative and not limiting. Various changes to the described embodiments may be made without departing from the spirit and scope of the invention.For the avoidance of any doubt, any theoretical explanations provided herein are provided forthe purposes of improving the understanding of a reader. The inventors do not wish to be bound by any of these theoretical explanations. Any section headings used herein are for organizational purposes only and are not to be construed as limiting the subject matter described. Throughout this specification, including the claims which follow, unless the context requiresotherwise, the word “comprise” and “include”, and variations such as “comprises”, “comprising”,and “including” will be understood to imply the inclusion of a stated integer or step or group ofintegers or steps but not the exclusion of any other integer or step or group of integers or steps. It must be noted that, as used in the specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the context clearly dictates otherwise. Ranges may be expressed herein as from “about” one particular value, and / or to “about” another particular value. When such a range is expressed, another embodiment includes from the one particular value and / or to the other particular value. Similarly, when values are expressed as approximations, by the use of the antecedent “about,” it will be understood that the particular value forms another embodiment. The term “about” in relation to a numerical value is optionaland means for example + / - 10%.References [1] J.-K. Hwang and Y.-C. Chen, “A combined detection-estimation algorithm for the harmonic-retrieval problem”, Signal Processing, vol.30, pp.177–197, Jan. 1993.[2] C. R. R. Christian Heumann and H. T. Shalabh, “Linear Models and Generalizations” Springer Berlin Heidelberg, 2010. [3] J. K. Hwang and Y.-C. Chen, “Superresolution frequency estimation by alternating notch periodogram,” IEEE Transactions on Signal Processing, vol.41, no.2, pp.727–741, 1993.[4] K. Zuiderveld, “VIII.5. - Contrast Limited Adaptive Histogram Equalization,” pp.474–485,Academic Press, Jan.1994.[5] Marr, D. and Hildreth, E., "Theory of edge detection", Proceedings of the Royal Society ofLondon. Series B, Biological Sciences 207, 1167, pp.187—217, 1980.[6] L. C. Zhao, P. R. Krishnaiah, and Z. D. Bai, “On detection of the number of signals when thenoise covariance matrix is arbitrary,” Journal of Multivariate Analysis, vol. 20, pp.26–49, Oct.1986.

Claims

Claims:

1. A method of analysing data from a distributed acoustic sensing system, the distributedacoustic sensing system comprising a plurality of spatial channels, each spatial channel associated with a respective scattering location along an optical fibre of the distributed acoustic sensing system, wherein the method comprises: obtaining a set of candidate trajectories corresponding to possible trajectoriesrepresented in data from the distributed acoustic sensing system; receiving, from the distributed acoustic sensing system, distributed acousticsensing data comprising a plurality of measurement signals as a function of time, eachof the plurality of measurement signals corresponding to a respective one of the spatial channels; computing an objective function based on the distributed acoustic sensing data,wherein the objective function is indicative, for each of the candidate trajectories, of a likelihood of that candidate trajectory being represented in the distributed acoustic sensing data; performing an iterative process comprising: selecting, based on the objective function, one of the candidate trajectories having a highest likelihood of being represented in the distributed acoustic sensing data; identifying a set of data points in the distributed acoustic sensing data corresponding to the selected candidate trajectory; updating the objective function by selectively removing from the objectivefunction a contribution arising from the identified set of data points; andrepeating the iterative process until a termination condition is met; and outputting a representation of detected trajectories, based one or more candidatetrajectories selected in the iterative process.

2. A method according to claim 1, wherein:the set of candidate trajectories comprises a mapping between each candidatetrajectory and a corresponding set of data point coordinates; andidentifying the set of data points is performed using the mapping.

3. A method according to claim 1 or 2, wherein identifying the set of data points comprisesdetermining end points and / or a line width of a representation corresponding to the selectedcandidate trajectory in the DAS data.

4. A method according to any preceding claim, wherein the iterative process furthercomprises determining whether the identified set of data points satisfies a predeterminedcondition and, if so, determining the selected candidate trajectory as a detected trajectory.

5. A method according to claim 4, wherein determining whether the identified set of datapoints satisfies a predetermined condition comprises: determining if the contribution to the objective function arising from the identifiedset of data points in the distributed acoustic sensing data satisfies a predetermined uniformity condition; and if the predetermined uniformity condition is satisfied, determining the selected candidate trajectory as a detected trajectory.

6. A method according to any preceding claim, further comprising:following an iteration of the iterative process, determining if a re-normalisation condition is met and, if so, re-normalising the distributed acoustic sensing data based ona remaining portion of the distributed acoustic sensing data which does not include one or more sets of data points identified in the iterative process; and re-computing the objective function based on the re-normalised distributedacoustic sensing data.

7. A method according to any preceding claim, wherein updating the objective functioncomprises: determining one or more affected candidate trajectories which map on to one or more of the identified set of data points; and re-computing components of the objective function corresponding to the affected candidate trajectories to remove a contribution from the identified set of data points.

8. A method according to any preceding claim, wherein the objective function comprises,for each candidate trajectory, a sum over the plurality of spatial channels of an error signalbetween the measurement signal and that candidate trajectory.

9. A method according to any preceding claim, wherein the iterative process furthercomprises determining if the identified set of data points satisfies a quality metric, and if thequality metric is not satisfied, determining that the termination condition is met if the quality metric was also not satisfied for the preceding Neiterations of the iterative process, where Neis a predetermined number.

10. A method according to claim 9, the quality metric for the identified set of data points isdetermined based on a comparison with a preceding iteration of the iterative process.

11. A method according to any preceding claim, wherein computing the objective functioncomprises a pre-processing procedure, the pre-processing procedure including: normalising each of the plurality of measurement signals; and / or performing noise reduction process on the distributed acoustic sensing data to improve a signal-to-noise ratio of the distributed acoustic sensing data.

12. A method according to claim 11, wherein the noise reduction process comprises:applying a contrast enhancement algorithm to the distributed acoustic sensing data; applying an edge detection algorithm to an output of the contrast enhancement algorithm; and subtracting an output from the edge detection algorithm from the output of the contrast enhancement algorithm.

13. A method according to any preceding claim, wherein obtaining the set of candidatetrajectories comprises pre-computing the set of candidate trajectories, and wherein the set of candidate trajectories is pre-computed based on dimensions of a data window of the distributed acoustic sensing system, and one or more predefined trajectory parameters.

14. An analysis system for analysing distributed acoustic sensing data, the analysis systemcomprising a processing device, and a memory storing instruction which, when executed by the processing device, cause the processing device to perform the method of any preceding claim.

15. A distributed acoustic sensing system comprising:a pulse generator configured to transmit a pulsed test signal along an optical fibre; adetector stage configured to receive a plurality of scattered signals from theoptical fibre, wherein each scattered signal corresponds to a respective spatial channel associated with a respective scattering location along the optical fibre, and wherein the detector stage is further configured to output a respective measurement signal as a function of each of the plurality of scattered signals; and an analysis system configured to receive the respective measurement signals from the detector stage, and to perform the method of one of claims 1 to 14.