Spatial resolution enhancement of rotational-radar subsurface datasets using combined processing method

Effective visualisation of railway tunnel subsurface features (e.g. voids, utilities) provides critical insight into structural health and underpins planning of essential targeted predictive maintenance. Subsurface visualisation here utilises a rotating ground penetrating radar antenna system for 360° point cloud data capture. This technology has been constructed by our industry partner Railview Ltd, and requires the development of complimentary signal processing algorithms to improve feature localisation. The main novelty of this work is extension of Shrestha and Arai’s Combined Processing Method (CPM) to 360° Ground Penetrating Radar (360GPR) datasets, for first-time application in the context of railway tunnel structural health inspection. Initial experimental acquisition of a sample rotational transect for CPM enhancement is achieved by scanning a test section of tunnel sidewall - featuring predefined target geometry - with the rotating antenna. Next, frequency data separately undergo Inverse Fast Fourier Transform (IFFT) and Multiple Signal Classification (MUSIC) processing to recover temporal responses. Numerical implementation steps are explicitly provided for both MUSIC and two associated spatial smoothing algorithms, addressing an identified information deficit in the field. Described IFFT amplitude is combined with high spatial resolution of MUSIC via the CPM relation. Finally, temporal responses are compared qualitatively and quantitatively, evidencing the significant enhancement capabilities of CPM.


Introduction
Maintaining the structural health of brick-lined railway tunnels ensures the safety of passengers, crew and rolling stock, alongside the continued efficient operation of the network. Structural issues commonly stem from the formation of voids within the tunnel lining that can show no hallmarks on the interior surface (e.g. hidden shafts, ring separation, water ingress) [1, p.42] [2] [3]. Left unaddressed, void growth within the subsurface can debond adjacent layers of brickwork and channel groundwater, compromising structural rigidity. Planning repairs is also complicated by the presence of concealed pipes and other tunnel utilities that must also be detected, either to fix or avoid damaging during maintenance. Development of an effective means to visualise railway tunnel subsurface is therefore of paramount importance.
Ground Penetrating Radar (GPR) comprised of fixed-directional antenna arrays has gained popularity amongst railway tunnel surveyors for non-destructive subsurface imaging. However, units are typically slow and utilise Uniform Linear Arrays (ULAs) that inherently suffer from blind spots when scanning cylindrical tunnel walls [4][5] [6]. The rotating GPR assembly utilised Table 1. Comparison of algorithms for temporal response extraction.

Algorithm Advantages Disadvantages IFFT
• High-accuracy return amplitudes.
• Fast transforms based on computationally efficient Radix-2 butterfly operations. • An in-built function in many mathematical programming languages, granting straightforward implementation.
• Implicit windowing can create false-artefacts in output signal. • Addition of noise processing filters can be an involved and time consuming process.

MUSIC
• Suppresses noise and returns sharp, high spatial resolution peaks in the temporal response. • Addition of spatial smoothing can decorrelate coherent GPR signals.
• Subspace approximation can result in low amplitude accuracy. • Multi-step algorithm computationally expensive. • Unlikely to be built-in function as very specialist. Thereby may demand highly involved manual programming.
The Combined Processing Method (CPM) -first formally proposed by Shrestha and Arai in 2003 [7][8] -innovatively combines the desirable amplitude of a temporal response recovered by IFFT with the high spatial resolution attained through application of MUSIC processing. The scheme was originally devised to enhance the spatial resolution of conventional ULA GPR, and derived from a weighted average formulation, for which datapoint proximity to signal peaks is considered to asymptomatically derive the scaling coefficients, as detailed in Section 3. Owing to the direct relevance of CPM for spatial resolution enhancement of GPR, and purely numerical implementation (i.e. required no significant hardware modifications), pursuit of its first-time application to a 360GPR dataset was identified as the primary objective for this study. Moreover, the authors identified several self-inconsistencies in the mathematics presented by the original study. To mitigate the risk of errors compounding in future independent research, and aid researcher understanding, the secondary objective of this study was to rectify these inconsistencies. Together, advances made would naturally build upon Shrestha and Arai's initial work through deployment of CPM in a practical context (i.e. tunnel inspection), which was an identified development point in their proof-of-concept investigation.
Interestingly, a notable lack of subsequent works (investigating practical applications of CPM) have been encountered during our consultation of the literature. This serves to reinforce the novelty of this study. The only direct follow up application aligned with GPR for infrastructure inspection concerned utility pipe detection beneath gravel roads in Japan, but was notably conducted by original authors [9]. All other citations of the CPM study almost exclusively did so for consultation of the mathematics behind MUSIC [10][11] [12], in which the inconsistencies identified were present. We considered if confusion induced here was a contributing culprit behind the lack of observed research interest in practical CPM applications. Cross checks with independent GPR-aligned studies centred on use of MUSIC confirmed the existence of innumerable studies discussing the mathematics underpinning MUSIC. However, only [13, p.60] consolidated a definitive sequence of steps for numerical implementation, which ultimately proved too vague for a researcher entering the field to utilise for programming their own version. We also considered accessibility and versatility of pmusic(), an inbuilt numerical computational tool for the MUSIC algorithm present in MATLAB2021(a). As a 'black-box', end-user adaption of the tool was found to be extremely limited. The implementation steps, whilst performed, were not clearly conveyed (considerably limiting accessibility), whilst direct recovery of temporal response profiles proved to be just as ambiguous a procedure (restricting versatility). Ergo, the final objective for this study was set to consolidate and provide clear instructions for the numerical implementation of the MUSIC algorithm (and in turn CPM), in order to fill the identified knowledge gap.

Investigation Outline
In this study, an initial declaration of key terminology and conventions is presented in Section 1.2. We proceed to perform the first numerical implementation of CPM on 360GPR traces constituting the rotational transect of a railway tunnel wall. Section 2 outlines the experimental setup adopted and workflow for data acquisition. Primary post-processing of captured frequency domain data was performed using windowed IFFT and a spatially smoothed variation of the MUSIC algorithm to recover temporal response profiles and facilitate secondary processing via CPM. A procedural overview of the key numerical methods is presented in Section 3. A comparative assessment of achieved temporal amplitude deviation and spatial resolution between IFFT, MUSIC and CPM processed signals follows, for which results are discussed in Section 4. Finally, conclusions and an appraisal of CPM-processed 360GPR efficacy for tunnel subsurface inspection are outlined in Section 5, alongside several future work opportunities.

Conventions Adopted
Terminology conventions adopted in this study include: • The term 'subsurface' concerns any concealed or inaccessible region of space behind a tunnel's innermost surface, and extends through the lining into the surrounding landmass. • Time and depth on a GPR response profile are proportional, therefore used interchangeably depending on context. • Following the summary of numerical implementation in Section 3, reference to 'IFFT' or 'MUSIC' is for brevity and assumes prior application of windowing and MSSP respectively. • 360GPR refers to data recovery achieved using the rotating GPR antenna. A trace is a response profile from a 360GPR dataset recorded with the antenna at consistent azimuthal elevation. A rotational transect is a collection of 360GPR traces acquired at distinct azimuthal elevations.
Mathematical conventions adopted in this study include: • The imaginary unit is denoted i = √ −1. • A subscript i denotes an arbitrary index in a 1D array.
• Estimations of signal return time assume signal propagation is at light speed, c = 299792458m s −1 . In reality, propagation velocity v < c within the coaxial cable and antenna. • Signal spatial dependency is omitted throughout for brevity. For example, where k denotes the associated wave-vector. • Temporal indexing begins at 0.

Rotational Radar Experimental Setup
Undertaking this study required initial acquisition of a sample rotational transect of a railway tunnel wall subsurface, for which multiple 360GPR traces (1D) at distinct azimuthal elevations of antenna pitch would need to be recorded. To achieve this, experimental data collection in a real tunnel was proposed. However, for safety and simplicity, it was finally decided the first deployment of the Railview system should be in a more controlled environment. Therefore, a test section of wall featuring a curved front face and 5-layers of concentric red-brick masonry was constructed to emulate the subsurface geometry of an archetypical UK railway tunnel sidewall. A photograph and schematic of the experimental setup have been provided in Figure 1.

Data Capture
Plane VNA Noteworthy features incorporated behind the brickwork included an air-filled void which spanned the wall length at depth 0.52m relative to the front face, alongside addition of a metallic strip concealed at relative depth 0.63m. The air-filled gap represented a weak-reflecting subsurface void our system would need to detect and relay to surveyors, whilst the metallic strip was expected to induce a stronger reflection akin to a buried utility also requiring detection (e.g. a metal drainage pipe). We note the close-proximity positioning (∆d ≈ 0.1m) and significant depth (d > 0.5m) of these targets was a judicious selection. From this, CPM performance in a worst-case survey scenario could be assessed. Namely, a situation wherein both IFFT and MUSIC could reasonably struggle to resolve both targets whilst maintaining a desirable (i.e. high intensity and accuracy) return-signal amplitude; foremost due to target proximity and secondarily due to the increased potential for scattering and attenuation of signals inherent to backscatter from deep targets.

Rotary Antenna Motor
A portable version of the Railview system was assembled by mounting the developed rotational GPR antenna onto a motorised shaft, which was affixed to a pushcart traction unit for ease of positioning. Air-launched SFCW GPR was implemented using a Vector Network Analyser (VNA). Response capture involved successive transmission of linearly swept frequency sinusoids across a 2.5GHz operating bandwidth with incrementation of 12.5MHz over a 1ms step cycling interval. Backscatter was recorded directly in the frequency domain as Fourier Coefficients, encoding both the phase and amplitude of returning wavefronts. Parsing by azimuth partitioned the rotational transect into individual traces, each containing 201 complex datapoints. These were then exported in a propriety binary format via a GBIP-USB interface cable to data processing software LabVIEW. Subsequent conversion to .txt plain text permitted usage for analysis in MATLAB2021(a). Table 2. 360GPR trace aperture-wall separation radii (r i ) and azimuthal elevation (θ i ). The capture procedure to obtain a 360GPR rotational transect of the test wall subsurface first required recording of a calibration trace. By directing the antenna skyward, backscatter originating internally from hardware reflections at material interfaces within the circuitry (antenna aperture, VNA connection port) could be detected. Subtraction of this response during post-processing of the main traces (see Section 3) would diminish internal reflection prevalence and in turn emphasise reflections associated with critical environmental features. The main rotational transect was captured using the SFCW GPR setup to record traces aligned with the antenna boresight vector as it pitched through 17 discretised angular elevations subject to continuous antenna rotation between azimuths of 1.03rad and 2.11rad. A breakdown of azimuthal elevations tested can be found in Table 2. The pushcart remained static throughout data capture to ensure non-helical (i.e. cylindrical) transect acquisition. We aim to apply CPM to 3D helical datasets in a separate study, so it is not pursued here. Further details can be found in Section 5. Angular elevations tested were not uniformly distributed by constant ∆θ ∈ R >0 in order to ensure the linear spatial sampling density of the scan remained constant across the inspection surface, emulating ULA GPR. Lastly, attention should be drawn to temporal offset. Backscatter recordings begin at twoway travel time t = 0s, which corresponds to the point of signal emission from the VNA output port. For brevity, we shall drop the qualifier when discussing 'time'. The reader should take care not confuse t = 0 with the central axis of antenna rotation (i.e. base of antenna), signal crossing of the antenna aperture or transmission through the front wall surface. Between VNA output port and the wall front surface, propagating signals traversed 0.91m of coaxial cable, followed by 0.35m to reach the antenna aperture, and finally through between 0.15m and 0.33m of free-space subject to pitching of the antenna.

CPM Data Processing
CPM is a secondary data post-processing scheme. For application to our frequency domain dataset X(f ), IFFT was first be applied to recover the temporal response profile x(t) from the Fourier Coefficients output by the VNA. A separate -but equally valid -temporal response y(t) was also be recovered through application of the MUSIC algorithm to X(f ). We direct the reader to Section 3.2 for further details on their numerical implementation. Both returned datasets were also normalised to reside in the amplitude range [0, 1], which avoided skew attributed to dissimilar magnitudes during merger.
In an ideal merger, we would wish to preserve the advantageous characteristics of both temporal response profiles. In reality, combining both fundamentally results in a trade-off. Therefore, we anticipated the CPM temporal response to be bounded to reside between x(t) and y(t) throughout time. This provides a useful check for researchers debugging their own CPM programs.
The CPM relation can be intuitively derived by considering the desired characteristics of a temporal response profile when either close or far away from a peak in the return signal. Close proximity to a peak is most accurately indicated by a sharp increase in the MUSIC response, owing to its higher spatial resolution than the IFFT response. In this region, we desire the CPM response to resemble x(t), as compounded approximations performed during the MUSIC algorithm (see Section 3.2) can limit the amplitude accuracy of y(t). Recalling sharp gradients correspond to high magnitude gradients, we represent this first condition mathematically as: In contrast, when far from a peak, the MUSIC response is near flat compared to the IFFT response, owing to its improved noise suppression. Here, we desire the return signal more closely resemble y(t) to preserve the noise suppression and maintain sharp peaks. Framed mathematically, our second condition takes the form: From this, the CPM relation can be formulated. The initial form resembles a weighted average of both responses, re-standardised by division by magnitude of the weighting coefficients, which are determined by asymptotic fit subject to our two identified conditions. This obtains the CPM temporal response profile, Numerical implementation of the CPM relation following computation of IFFT and MUSIC responses facilitated extraction of the temporal response profiles for the experimentally obtained 360GPR rotational transect, depicted in Figure 2. For further complementary graphical depictions of CPM response extraction, we direct the reader to the figures provided in [8, p.4-6].

Numerical Implementation of IFFT
A built-in IFFT() function does not always feature in high-level mathematical programming languages (e.g. Fortran). Where it does (e.g. MATLAB, Python, Wolfram Language), application of a suitable windowing function for spectral leakage compensation is often an involved and ambiguously conveyed process, discussed solely in the context of FFT when not dismissed as common knowledge. Manually programmed numerical implementation of windowed IFFT is non-trivial, therefore this gap in formal literature presently limits its accessibility, particularly to researchers new to the field. For completeness, we provide a brief summary. Spectral leakage stems from artificial discontinuities introduced to the input signal during Fourier Transform operations, which widen the base of peaks in the temporal response, reducing achievable levels of spatial resolution. This can cause close proximity peaks to merge, alongside false-artefact manifestation. These discontinuities arise because sampled signals rarely contain an integer multiple of periods, resulting in an amplitude offset between the start and end of the waveform profile when the observed signal is assumed to exhibit periodicity. Windowing reduces the prevalence of spectral leakage in the signal response returned by IFFT, improving attainable spatial resolution -albeit not to the same extent as the MUSIC algorithm. The window w(t), similar to CPM, can be regarded as a weighting function. By emphasising or suppressing specific constituent frequencies within a temporal response signal, peak base widening and the prevalence of artefacts originating from spectral leakage can be reduced. However, selecting the most appropriate w[n] is a subjective process, typically informed by experience. In our case, on account of the wide variation in constituent frequency amplitudes in the 360GPR dataset (see Figure 4), this study utilised the Kaiser-Bessel window [14, p.73], Here, each discretised bin of the frequency domain is indexed by the integer n ∈ [1, 2, . . . , N ], where our domain (i.e. trace) contains N datapoints. A key advantage of the Kaiser-Bessel window is that the relative sizes of the main and side lobes can be fine-tuned, here by setting α = 3. We denote the modified zeroth-ordered Bessel Function of the first kind by I 0 .
Windowing must be applied to the dataset before the transform operator (here IFFT) is applied. As a window is a form of digital filter, to 'apply' it, we convolve with the VNA output signal (in the frequency domain) in accordance with the Convolution Theorem. Formulated mathematically, The windowed dataset can now undergo (Discrete) IFFT, where k ∈ [1, 2 . . . , N ] indexes discretised values in the temporal domain. From this, an IFFT temporal response profile with improved spatial resolution can be recovered, as we desired.

Numerical Implementation of MUSIC Algorithm
As a spectral estimation technique widely adopted for signal processing due to its highresolution capability [11][15], the MUSIC algorithm exploits signal-noise orthogonality in received wavefronts to output temporal response profiles exhibiting suppressed noise and sharper peaks. The procedure is considered a generalised signal processing technique [13, p.61], thereby applicable to arbitrary antenna setup topologies, permitting application to our single rotating GPR antenna. Implementing this algorithm numerically reduces to five steps: (i) Approximate the signal covariance matrix, S.
(ii) Decompose S to identify its eigenvalues, λ i and eigenvectors, v i . (iii) Partition v i by magnitude of λ i to approximate the associated Noise Matrix, E W .
(iv) Calculate the associated signal steering matrix, A.
(v) Using E W and A, compute the pseudospectrum to recover the desired temporal response.
We proceed to walk-through the key equations utilised in Steps (i)-(v). During consolidation, dimensional analysis was performed in each step to locate and inform rectification of mathematical ambiguity in the original study. Meaningful elaboration on each step requires formal definition of the problem model. The VNA output signal for each 360GPR trace may be modelled as: where, Here, s is the signal vector, representing the 'true' backscattered signals from K subsurface features. Additive White Gaussian Noise emulating signal clutter and thermal interference is contributed by the noise vector w. The signal steering matrix, A, governs the direction of arrival of received wavefronts. It is comprised from K individual steering vectors which encode the phase delays associated with each wavefront, and by their association with time, we may link them with the discrete time increments of our desired temporal response profiles. The steering vector of wavefront p ∈ [1, 2, . . . , K] takes the form: In Step (i), denoting complex conjugation by * , the ideally required autocorrelation matrix of X is approximated by the N × N sample covariance matrix, As each trace recorded corresponded to one VNA output signal measurement during data acquisition, computation of the expected value is for a sample size of one, permitting its omission here. It should be remarked that this matrix relays the degree of correlation (standardised covariance) between the antenna elements (within a ULA) and the individual backscattered signals which comprise the VNA output [16]. Although 360GPR does not utilise a ULA, we may consider each emitted waveform in the frequency sweep conducted to be representative of a separate constituent 'element' in this approximation. In Step (ii), Eigenvalue Decomposition (EVD) of S abides, For this study, the MATLAB in-built function eig() permitted automatic calculation and sorting of v i into ascending order based on the magnitudes of corresponding λ i (i.e. 0 ≤ λ 1 ≤ λ 2 ≤ · · · ≤ λ N ). We are always permitted to define this ordering as a consequence of S being Hermitian by formulation. For programming environments lacking an in-built function, EVD may be manually performed through computation of the characteristic polynomial, where det() is the determinant operator, I is the N × N Identity Matrix and λ i are recovered by calculation of the roots. Once extracted, back substitution of λ i into Equation 11 return corresponding v i . A low complexity sorting algorithm such as Quicksort() may then be applied to arrange λ i and v i into ascending order of magnitude as described prior. In Step (iii), we proceed to partition the eigenvectors into the constituent Signal Subspace U S and Noise Subspace U W of X in order to exploit the orthogonality of their basis vectors (selective v i ) for noise suppression and scattering centre localisation. Numerical implementation requires horizontal concatenation of v i (as in construction of A) associated with low magnitude λ i to approximate the Noise Matrix E W , which contains the basis vectors of U W . By taking the complement, recovering v i with associated high magnitude λ i , one can also obtain the approximated Signal Matrix E S . For this study, only E W required explicit computation.
We note that what constitutes a low magnitude is inherently ambiguous. In fact, attempts to establish an optimal partition to form E S and E W has lead to development of numerous approaches including the Nyström Method [17], ML Inverse Fitting [18] and more recently, Quantum-Behaved Particle Swarm Optimization [15]. Numerical implementation of these methods is highly involved, and can be over-specification for what is often a visually discernable small proportion of dominating λ i (e.g. [0.01,0.05,10,0.2]). A reasonably reliable and low complexity alternative is to perform an array search over all v i for low magnitude λ i , 10 incorporating a manually defined acceptance threshold parameter β to enable refinement and fine-tuning of subsequent partition attempts. In use of a subjective declaration, we have introduced a small amount of uncertainty into the results, however, the considerable reduction in computation complexity can provide very reasonable justification for this choice.
In accordance with standard practice in statistical confidence tests, this study considered a 10% deviation from the maximal magnitude eigenvalue grounds for inclusion in E W (other common values include 5% and 1% subject to context). This was implemented using the acceptance threshold, where λ b for b ∈ [1, 2, . . . , N ] denotes an arbitrary eigenvalue under consideration. Numerical implementation of Step (iv) centres on definition of discretised frequency values f i and temporal values t i to calculate the Steering Matrix defined in Equation 9. These values are fundamentally linked by the physical relation which is preserved under Fourier Transform operations. Here, T denotes the frame size (i.e. duration) associated with the response signal sampled by the VNA post subsurface excitation. Literature consultation revealed it is customary to plot response profiles to durations up to 10 − 20ns after the final peak, in order to convey the signal has attenuated to the level of background noise (i.e. finished). For example, in this study T ≈ 80ns, but we chose to plot to 40ns in order to focus more clearly on the main peaks. Finally, in Step (v), the pseudospectrum is evaluated to return the desired MUSIC processed temporal response profile: This relation derives from the aforementioned orthogonality of U S and U W , which is shared by E S and E W by design. However, numerical round-off error is introduced during digital computation. Perturbations are also introduced to S, which arise from approximation [13, p.60] (regardless of adopted spatial smoothing -see Section 3.3). Consequently, steering vectors a[t k ] recovered in Step (iv) exhibit near, but not exact, orthogonality with basis vectors in U W , and by extension E W . Definition of a suitable metric returns zero when the orthogonality condition is met between eigenvectors in E W and steering vectors a[t k ] (associated with wavefronts returning from scattering centres). This is equivalent to performing a Fourier Transform, thereby transitioning the returned response to the time domain. Taking the reciprocal results in sharp peaks at scattering centres and near zero response amplitude in-between, thereby returning a noise suppressed high spatial resolution temporal response profile, as desired.

Spatial Smoothing Adjustments
SFCW GPR systems dependant on VNAs are predisposed to generation of coherent returning wavefronts [8, p.3], which can cause issues with MUSIC processing. Specifically, in MUSIC Step (i), it was implicitly assumed that all returning wavefronts represent distinct signal waveforms (incoherence) [13, p.143], this deviation from reality can limit the accuracy of recoverable temporal responses. The same issue pertained to MUSIC processing of 360GPR traces, motivating implementation of spatial smoothing adjusters SSP and MSSP, which were also adopted by Shrestha and Arai to decorrelate coherent return signals. Numerical implementation first required subdivision of X[f ] (a single large snapshot of data) into M overlapping subarrays (i.e. multiple smaller 'snapshots') containing L < N datapoints. Subdivision induces a desirable artificial phase variation between snapshots, which can be exploited for decorrelation. This is achieved through spatial averaging of the covariance matrices associated with each subarray. To formally define these matrices, we must update the VNA output signal model (Equation 7) to accommodate the additional snapshots. For the m th snapshot, this yields where B represents the updated steering matrix, given by Here, D is a diagonal matrix of phase delays, which we may liken to the artificial phase variation we have introduced by considering multiple smaller snapshots (i.e. subdividing the original dataset) to combat coherence. It takes the form: where ∆f is the frequency incrementation of the linear sweep. With this, the now reduced L × L covariance matrix can be approximated: The spatially smoothed covariance matrix, obtained by SSP and MSSP respectively, are in turn and where J is the L × L 'forward-'backward' Exchange Matrix, defined: A downside of spatial smoothing, as with CPM itself, is that implementation involves a native trade-off. By definition of the snapshots, we may relate L, M and N by Therefore, increasing the number of snapshots to improve decorrelation performance must reduce the number of datapoints per snapshot, thereby decreasing the attainable spatial resolution of returned temporal responses (and vice versa). We adopt M = 10 throughout the remainder of our study based on trail and error inspection of output to return peaks resolvable to within than 0.5ns -a scale comparative to the thickness of the air-filled void and metallic strip. MSSP is 12 also exclusively adopted owing to its improved performance capability on data from ultra-wide band GPR systems, as well as facilitate meaningful comparison with results obtained in [8] during design and control testing of a spatially smoothed MUSIC algorithm in MATLAB. Care is also taken to truncate the plotted time frame of temporal responses well in advance of 70ns (in addition to main-peak focusing), as false-artefacts began to form in this particular region when MSSP was utilised. Notably, these artefacts (symmetric reflections of the main peaks where t < 40ns) were not present when only SSP was adopted, or when no spatial smoothing was included. From these observations, we deduced these artefacts likely stemmed from the forward-backward averaging pattern of MSSP. Feasibly, the transient nature (relatively shortduration) of the return waveform placed main peaks closer to the beginning of the time frame than the middle, thereby skewing the average for temporal values symmetrically displaced from the end of the frame.

Numerical Implementation of CPM Relation
Implementation of CPM (Equation 3) required accurate numerical calculation of first-order partial derivative terms. This was achieved using the fourth-order central differencing approximations: where ∆t = T /N is the temporal incrementation value. Ghost-points were defined to enforce physically grounded boundary conditions for recovered temporal response profiles. Considering t = 0, the fixed boundary condition x = y = 0 was enforced for ghost-points t −2 , t −1 . It was reasoned that as no signal backscatter could occur without sufficient elapsed time for the transmission signal to propagate to a subsurface feature, backscatter from it, and return to the receiver antenna. In addition, considering t = T , a further periodic boundary condition was implemented for ghost-points defined at t N and t N +1 : In reality attenuation and signal scattering within the subsurface results in decay of backscatter amplitude intensity to the level of background noise. However, in our finite sample frame, we cannot be completely certain when the return signal has fully attenuated owing to the continual presence of background noise. In other words, we are considering a small sample of a potentially much longer duration signal. This motivated selection of a periodic boundary condition.

Results and Discussion
Post-processing of the acquired 360GPR rotational transect returned the temporal response profiles shown in Figure 2. Initial qualitative inspection highlights Kaiser-Bessel windowing does not completely compensate for spectral leakage. This is evidenced by the considerably wide base of the main temporal response peak -centred on approximately 11ns in Figure 2(a) -which at a glance appears to extend 10ns to either side and merge with neighbouring peaks between 10ns and 13ns. Although distinguishable on the IFFT response with close attention, three more clearly resolvable peaks are exhibited by the MUSIC processed temporal response -in Figure 2(b). This visually demonstrates the significantly improved spatial resolution achievable through exploitation of the orthogonality between the signal and noise subspaces. However, amplitude discrepancy is also apparent, with maximal variation of 20% for the peak at 12ns. We can reasonably assume the amplitude recovered via MUSIC to exhibit lower accuracy than that of the IFFT response, owing to its dependency on compounded (often subjective) estimation of both S and E N , alongside the effective reduction in considered datapoints (N = 201 L = 192) imposed by spatial smoothing.

Qualitative Analysis -CPM Temporal Response
On qualitative inspection of the CPM response profile, we observe in Figures 2(c-d) that the desired amplitude of the windowed IFFT response has been successfully integrated with the enhanced spatial resolution of the MUSIC profile, as initially desired. Consulting the CPM temporal response of the horizontal trace in Figure 2(e), high amplitude peaks near t = 9ns, t = 12ns and t = 13ns indicate the presence of strong reflectors. Direct Arrivals (DA) from initial backscattering at the front air-wall interface are anticipated to mark the first considerable reflection event. Using the fundamental distance-time relation, we may estimate the expected delay time to be: As we have visually estimated location of the peak from the profile and the amplitude of the peak is 40% greater than its nearest neighbours, we may reasonably associate this peak with the direct arrivals, in turn providing a reference point for relating travel time to penetration depth relative to the front wall interface at d = 0m (t ≈ 9ns). We make a side-note of successive small but consistent amplitude peaks observed for t = 5-9ns. It is apparent they are residual artefacts from the IFFT response by visual inspection, most likely attributable to the antenna backlobe interacting with its support mount and the metallic guide rails below the pushcart. The non-zero tails of the Kaiser-Bessel window could permit such artefacts to persist, making this mechanism plausible, but as direct arrivals are readily identifiable, for the present study these artefacts are not of significant concern.
In regard to the close proximity Air Void (AV) and Metallic Strip (MS), at 0.52m and 0.63m respectively, we can again approximate the expected peak delay time: which is in agreement with corresponding peaks present in the CPM temporal response profile, as we desire. This observation evidences successful identification and discrimination of the close proximity targets, for which similar analysis yields unclear interpretations from either IFFT or MUSIC when consulted in isolation. We do note return amplitude of both peaks is low despite their significant associated discontinuities in dielectric constant within the surrounding subsurface media. Signal spreading and attenuation inherent with long duration wave propagation through even slightly inhomogeneous media (bricks, mortar) presents a feasible explanation, for which introduction of a time-variable gain function may permit compensation in a future study. A final point of consideration is the nature of the six peaks recorded between the direct arrivals and the air void. Consulting Figure 1, we can match features by associating successive peaks with consecutive material interfaces, interpreting larger amplitudes with more significant variation in the dielectric constants of the media. Working backwards from the suspected air void, a nearby peak at 12ns of moderate amplitude is indicative of detection of the plastic membrane. Between t = 10.5ns and t = 11.5ns, two significant peaks occur with a third lower amplitude central peak enclosed. This is suggestive of a wide object, for which the two ends would present sharp material discontinuities invoking significant backscatter, whilst multiple partial reflections from within amount to lower amplitude -but still significant -set of central peaks. The vertical backing breeze block matches this description, and is the next substantial material interface, motivating their association. Finally, between the breeze block and wall front face we would expect to observe 4-5 sequential peaks of similar amplitude if individual bricks (or mortar joints) have been resolved. This closely resembles observed amplitude variation between t = 9-11ns. Excluding micro-fractures and fine cracks, resolution of individual bricks could reasonably capture the majority of significant tunnel defects, providing strong evidence that CPM enhanced 360GPR could form the basis of a powerful subsurface inspection method for tunnel structural health surveys.

Quantitative Analysis -Temporal Response Amplitude
To formally assess the amplitude accuracy achieved by CPM in this study, we proceed to make a quantitative comparison with the amplitude accuracy achieved by MUSIC. Taking the IFFT temporal response as a comparative baseline for each trace in the rotational transect, Root Mean Square Error (RMSE) was calculated using amplitudes extracted from both MUSIC and CPM temporal response profiles to compute respective residuals. Results are displayed in Figure 3(a). Subsequent recovery of absolute percentage difference between RMSE output for each trace then permitted calculation of a representative mean average measure of amplitude accuracy, presented in Figure 3(b). It should be remarked that a higher absolute percentage difference equates to a temporal response more closely aligned with the IFFT profile.   Average observed percentage reduction in RMSE across all 17 compared traces was found to be 52.7%, indicating CPM exhibited effectively half the amplitude deviation experienced by MUSIC. We note an apparent rise and fall trend in alignment variation across the rotational transect, with amplitude deviation about the horizontal trace between 20-30% lower (i.e. higher absolute percentage difference) than at extremal azimuthal elevations.
A potential explanation requires consultation of Figure 1, in which we note that the boresight vectors of traces with extremal azimuth, when extended, do not predominately reside within the test wall, but instead inadvertently pass into the ground for θ i≤3 and out of the top of the wall back into free space for θ i≥16 where they encounter few significant reflectors. By Equation 3, deviation between y[t] and z[t] is lowest when the return signal exhibits very few significant peaks, also corresponding to the encounter of few significant reflectors. This consistency complies with our observations. Moreover, the approximate 5% greater percentage difference between traces of low and high azimuth would also be in agreement with their propagation into the ground rather than back into free space. More low amplitude peaks would be returned by ground-entering signals as this medium exhibits greater inhomogeneity than free space, generating more clutterlike backscatter than free space. Despite their low amplitude, such clutter peaks could still cause z[t] to marginally less closely resemble y[t], feasibly resulting in the percentage difference offset we observe against the traces with high azimuthal elevation.
In all however, the rise-fall trend does not undermine that results attained by CPM in this study consistently evidenced at least a 40% reduction in amplitude deviation from IFFT compared to the exclusive application of MUSIC. Excluding extremal traces θ i∈ [4,5,...,15] , this improvement increases to 55.2%. Considering practical application, an operational real tunnel would completely surround the Railview system, preventing traces of high azimuthal elevation re-entering free space unless a void was present. Whilst transmission through the floor would remain possible for low azimuthal elevations, we conclude the origins of this effect are of academic interest, but ramifications of its existence are of negligible consequence for practical deployment.

Quantitative Analysis -Temporal Response Spatial Resolution
Enhancing the spatial resolution of a temporal response results in a new profile containing a higher proportion of sharp peaks, thereby resulting in a corresponding increase to the median frequency of the signal. We denote median signal frequency with an overhead tilde (e.g.x). Performing FFT and comparingx andz across axes standardised by sampling frequency, we can obtain suitable indirect quantitative measure of signal spatial resolution. The reader may have noted thatỹ would significantly exceed both of these responses. However, we do not consider the MUSIC response here owing to its known undesirable amplitude deviations; we have established that CPM represents a more relevant near optimal trade-off in Section 3. Results for extremal and horizontal traces can be found in Figure 4.
Spectra obtained across all 17 responses evidencez is an average 152% greater than thex baseline for the same trace. We also observe an apparent linear growth of approximately 6% per trace in the percentage difference betweenx andz as azimuthal elevation increases.
Similar to trend observation in amplitude deviation, this growth does not undermine the significant practical benefit application of CPM has on spatial resolution of 360GPR. However, this observed trend does instil academic curiosity surrounding its origins. We suspect this behaviour can be jointly attributed to the curvature of the test wall and offset position of the antenna central rotation axis relative to the hypothetical tunnel centreline. Practically, we see in Table 2 that r i more than doubles for the higher elevation traces, increasing two-way travel time for signal propagation. As the transient waveform duration broadens, the IFFT begins to resemble a single large peak, since individual peaks typically do not drop completely to zero amplitude between targets owing to noise prevalence. By contrast, the MUSIC-enhanced spatial resolution of CPM better diminished noise prevalence to resolve the individual peaks,  16 separating them into multiple distinct transient waveforms. As wider waveforms exhibit lower median frequencies, this would explain the growth in percentage difference we have observed as azimuthal elevation increases, becausex is reducing more rapidly thanz in this region.
In practical deployment, improved consistency of r i will be achieved by adaptive repositioning of the Railview system antenna using a robotic armature, enabling it to better follow the tunnel centreline. Here, we chose to position the antenna close to the wall to reduce signal attenuation, as the pushcart setup was too low to viably mount a reflector to boost return signal strength.  Figure 4. Spatial resolution enhancement quantified by consideration of relative difference between median frequency of IFFT and CPM temporal responses respectively. Power spectra are presented above for 360GPR traces of azimuthal elevation (a,b) θ 1 , (c,d) θ 8 and (e,f) θ 17 .

Conclusions
Improved localisation of subsurface features behind railway tunnel walls was successfully demonstrated through practical implementation of CPM on the first 360GPR datasets to be captured by the new Railview system. A comprehensive summary of steps for numerical implementation of IFFT, MUSIC and CPM was also created. In the experimental study conducted, CPM was numerically applied to 17 traces of distinct azimuthal elevation in MATLAB2021(a), which were acquired experimentally from a test section of tunnel wall containing a concealed air-filled void and metallic strip. Analysis of temporal response profiles recovered evidenced an average 52% reduction in amplitude deviation of CPM responses compared to those obtained by MUSIC, whilst spatial resolution was observed to have increased by an average of 152% on the baseline set by IFFT. The numerical implementation summary should provide researchers entering the field of GPR signal analysis with clearer and more comprehensive understanding of essential processing techniques, and has corrected mathematical inconsistencies from previous works. Moreover, the experimental investigation has successfully demonstrated CPM can indeed achieve significant spatial resolution enhancement of 360GPR datasets, and could therefore provide an effective tool in tunnel structural health surveys. Looking ahead, prospects for future work would include: • Generalisation of CPM application to operate on 3D helical datasets obtained by concurrent operation of the rotating antenna and a rail-mounted traction unit. • Replacement of subjective trail-and-error parameter selection with a more consistent selection scheme for numerical implementation of the MUSIC algorithm. Namely, determination of threshold β for estimation of the Noise Subspace U W [18] and refined selection of snapshot parameter M for array subdivision in spatial smoothing. • Processing and analysis of a larger sample of independent 360GPR rotational transects to improve accuracy and reliability of quantitative spatial resolution evaluation. Ideally including authentic datasets obtained from real operating railway tunnels. • Consideration of alternative measures for analysis of spatial resolution enhancement.
For instance, assessment of any observed increase in the associated signal-to-noise ratio, computation of Resolution Matrices [19, p.2] or calculation of metrics derived from the Modulation Transfer Function (e.g. MTF50P [20]). • Comparison between CPM performance and independent spatial resolution enhancement techniques (e.g. Diffraction Stack Migration [21]). • Expansion of relatively computationally lightweight CPM (viable for efficient onsite 360GPR processing) to act as feeder for more intensive deep learning assisted f -k analysis tool for super-resolution GPR processing [22] (viable for detailed 360GPR laboratory processing).