Signal Processing Techniques for Maritime Surveillance with Skywave Radar
Mike D.E.Turley
Defence Science and Technology Organisation
Edinburgh,SA,Australia
mike.turley@v.au
Abstract—Detection and tracking of maritime targets using skywave radar is influenced by the propagation medium,in-terference environment and target scenario.Acquired data dis-play distortion,fading,non-stationarity,and heterogeneity.Brief examples of data are given,then signal processing techniques are developed to provide robust adaptive Doppler processing, rejection of impulsive noi,improved CFAR using the Weibull distribution with robust two-parameter estimation,and a simple track-before-detect scheme for enhancing small SNR target detection performance.
I.I NTRODUCTION
Skywave radar prents a means to survey large areas of a for maritime targets[1].Targets of interest include defen platforms,commercial activities,and illicit incursions and trafficking.Radar operations will typically also monitor associated air movements.In this paper signal processing techniques are prented for improving the radar detection capability of surface targets.The radar cross ction of the targets is quite varied and depending on their radial velocity component the targets are to be detected against clutter and/or noi.
The performance of a skywave radar ship detection system is highly dependent on the received clutter.Thefirst step is thus to lect the optimal operating waveform and scan strategy for a given target scenario[2].Once the waveform is lected a number of signal processing objectives may be considered.Ideally the receiver operating characteristic for the complete signal processing detection and tracking system should be optimized.This is however somewhat intractable given the number of available algorithms and parameters,so instead,the focus is on building parable signal processing components that adapt to the target scenario and radar en-vironment.Rudimentarily the signal processing components are range,Doppler,beamforming,peak lection and track-ing.Each of the components are adaptively procesd to match signals,counter distortion,and mitigate clutter and interference.In some situations performance may be improved by joint processing.In this paper we do not address spatial processing(e for example[3]).
Section II briefly describes the data characteristics,with the following ctions describing appropriate temporal,spec-tral,constant fal alarm rate(CFAR)and track-before-detect (TBD)signal processing techniques.
II.D ATA C HARACTERISTICS
Figure1shows two stacked example range-Doppler(RD) maps collected by the Jindalee radar[4]using conventional ta-pered Fourier transform Doppler spectral analysis with0.02Hz analysis bandwidth.The dominant feature is a backscatter. The a backscatter appears as positive and negative Doppler shifted lines,known as Bragg lines caud by resonant scat-tering from advancing and receding a waves of wavelength equal to half the radar carrier wavelength[5].The Bragg lines are smeared due to ionospheric distortion,and have range dependent Doppler shift impod by ionospheric effects and ocean currents.Second order a backscatter is prent decreasing the potential to detect targets near to and between the Bragg lines.Two target-like respons are obrved outside the advance Bragg line.The impulsive nature of a meteoric trail echo produces a Doppler steak at near range.The prence of a cond ionospheric propagation path produces a cond t of weaker highly range dependent Bragg lines.Detection of targets removed in velocity from the surface backscatter is obviously more
tractable.
Fig.1.Example RD maps collected using the Jindalee radar. Figure2shows the spectrogram of the range containing the candidate target.The spectrogram spans eleven conc-utive radar scans.Temporal gaps occur between scans.The spectrogram shows the dominant Bragg lines and target signal have associated time varying Doppler shift.The Bragg lines are poorly resolved indicating perhaps additional ionospheric distortion,multipath propagation,or insufficient spatial reso-lution.Parent describes spectral broadening mechanisms[6]. Hence the goal is to detect targets against a clutter and noi
background that is non-stationary and
heterogeneous.
Fig.2.Spectrogram of the range containing the candidate target.Eleven concutive scans are plotted.
III.A DAPTIVE D OPPLER P ROCESSING Adaptive Doppler processing is an area of rearch aimed at coherent processing,mitigating spread Doppler clutter,re-jecting interference,and high resolution processing.Adaptive Doppler processing in HF surface wave radar is described by Fabrizi
o et al.[7].Doppler processing is gregated here as a temporal pre-conditioning process and as spectral processing.
A.Temporal Pre-processing
The main objectives of temporal processing are to reject im-pulsive interference,and do ionospheric distortion correction (IDC).Examples of impulsive noi rejection techniques for skywave radar are given in[8],[9],[10],[11],[12].Other techniques are also suitable,with for example the Sacchi hnique[13]developed further here in Section III-B.IDC rves two purpos,namely a pre-processing step to coherent Doppler processing by removing the intra-scan frequency variation[11],[12],[14],[15],and to remove the mean ionospheric Doppler shift thereby improving scan-to-scan Doppler registration prior to subquent tracking.
To register the Doppler shift the full Doppler spectrum was analyd for peak frequency and land or a classification. Smoothing of estimates over nearby range and azimuth cells was employed byfitting a low order2D polynomial using a least squares estimate with sample weights bad on outlier likelihood.Boundary constraints were included together with a bias towards lower orderfitting.A parameter to model ocean current is also necessary when the range-azimuth map spans a land-a interface.
B.Spectral Processing
High resolution spectral analysis methods have the potential to either improve detection of poorly resolved signals,or alter-natively achieve a desired performance with reduced coherent integration time(CIT).The latter incentive has huge impact on radar operations by freeing the precious radar resource for additional scanning.High resolution methods such as auto-regressive(AR)models[16]have been considered,but are nsitive to parameter lection and suffer performance degradation at low SNR.The following subctions describe robust high resolution processing bad on parametric and non-parametric data extrapolation techniques.Strictly the are also temporal pre-processing steps to conventional Fourier spectral analysis.
1)Parametric Data Extrapolation:Data extrapolation can be ud to improve Doppler resolution by extrapolation of the radar sweeps into the past and future.Doppler spectra are characterized by band-limitedness and high dynamic range. The spectra are well suited to robust resolution enhance-ment.Figure3provides a schematic indicating the profit of the techniques.Consider acquisition of N sweeps that with conventional Fourier processing are tapered to reduce sidelobe leakage.The red taper is an N=64sample80dB Taylor example taper.If the sweeps are extrapolated into the past and future,followed by conventional Fourier processing then better u of the acquired
data are made.The blue/yellow taper shows the ca for an extrapolation factor of2.The blue portion shows incread gain from the acquired data,and the yellow portion shows the extrapolation contributes little energy,but is however critical in controlling the sidelobe leakage from strong
signals.
Fig.3.Temporal tapers.Red-taper if only half data was ud.Blue-the taper applied to the actual data.Yellow-continuation of the blue taper into the extrapolated sweeps.
Data extrapolation methods such as DATEX[17],[18]were introduced to skywave radar application by Turley and V oigt in1992[19].Gadwal and Krolik[20]make comparison of DATEX with other techniques such as Hoppler[21].
The DATEX method extrapolates the data ends using linear prediction coefficients.The coefficients are obtained from the full data quence by assuming a P th order AR model:
x n=
P
p=1
a p x n−p+e n N
where x n are the data samples,a p are stable AR coefficients and e n is a white noi process.Then the linear prediction forward estimates are:
ˆx F n=
P
p=1
a p x n−p n>N
and the backward estimates are:
ˆx B n=
P
p=1
a∗p x n+p n<1
Figure4(i)shows a conventionally procesd RD map for 128acquired sweeps.This map is considered as‘truth’data. The DATEX method with extrapolation factor2was applied using only the central64sweeps of the original data.The RD map produced using DATEX,Fig.4(ii)is almost identical to the original,indicating a3dB processing gain.Even the broadband signals(impulsive meteoric clutter),for which the AR model is mismatched,are also replicated.This leads to our suggestion to halve the radar scan CIT,and apply the DATEX technique.This frees half the radar resource and thus allows doubling of the scan rate.The technique is robust to model order lection since only the strong clutter signals require modelling and extrapolation.
To improve robustness the following modifications to the DATEX method were made:
a)Impulsive noi outlier corruption of the AR coeffi-
cients was mitigated using the Nuttall[22]technique.
b)To account for non-stationary clutter(a potential
problem for long CITs)the forward and backward
predictions were bad on AR coefficients calculated
from respective halves of the acquired data.
c)Extrapolation from bad data samples was avoided
using the minimumfilter error technique
[19].
Fig.4.RD map Doppler processing(i)original‘truth’generated from128 sweeps,(ii)generated from central64sweeps with DATEX extrapolation to128sweeps,(iii)Sacchi extrapolation,and(iv)Sacchi extrapolation& interpolation.
2)Non-parametric Data Extrapolation:Non-parametric techniques may also be employed to extrapolate/interpolate the N dimensional time ries data x,to form a high resolution K dimensional spectrum s.Our objective is to produce an extended time quence that is like the original data samples, subject to maximizing the spectral sub-clutter visibility(SCV). SCV is defined as the ratio of total energy to average noi energy.This objective is basically satisfied by the Sacchi et al.[13]method where the following cost function is minimized: J(s)=(x E−F s)H P t(x E−F s)+λ2
K
k=1
ln
1+
s∗
k
s k
2σ2s
where F is the inver discrete Fourier transform,H denotes Hermitian transpo,x E is the zero padded extended version of x,λ2is a weighting term with the cond term as a regu-larizer impod by modelling the spar spectral components as Cauchy distributed.The P t term is introduced here as a temporal lection matrix.The solution is:
s=(λ2Q(s)+F H P t F+γ2I)−1F H P t x E(1) whereγ2is introduced to add a small amount of diagonal loading.The diagonal matrix Q is basically a bandlimiting projection that is biad towards the noi spectra:
Q kk=
1+
s∗
k
s k
2σ2s
−1
with the sparness controlled by theσs parameter.Note the projection matrix Q is a function of the spectrum s,so the solution is found by calculating(1)iteratively.This method is referred to here as the‘Sacchi’method.The RD map produced using the Sacchi method shown in Fig.4(iii)is also almost identical to the original‘truth’data.The formulation given in (1)allows interpolation or replacement of missing or corrupt sweep data,identified in P t.Impulsive noi was detected using a clutter whitening AR technique[8];the corrupt sweeps were then replaced using(1)with results shown in Fig.4(iv). Good rejection of the impulsive meteoric clutter was achieved.
A similar algorithm bad on the l p-norm method[23]offers comparable performance.
IV.T WO-P ARAMETER CFAR
In skywave radar the receiver operating characteristic de-pends critically on the CFAR algorithm.The RD map data is heterogeneous in range,Doppler and azimuth,and also varies with propagation medium,scattering sources and system gains. The objective of the CFAR algorithm is transform each test cell such that a global threshold will produce a constant fal alarm rate.It is desirable that the CFAR transformation has the properties of localization,edge performance,and is robust to target outliers.Due to the heterogeneous nature the scale estimates must be calculated from a small local data sample. The sample support size is a function of range and Doppler, and in particular the support size may be larger in the noi zone than the surface clutter zone.The requirements lead to using order statistic methods[24]with sample kernels matched to typical skywave radar clutter and interference profiles[25]. Figure5plots the cumulative distribution function of 146,000samples from veral RD maps,one of which is shown in Fig.6(i),on Weibull probability paper.The data values were normalized to a noi energy estimate taken as the 25th percentile.The Weibull probability distribution is given as:
p(x)=
C
B
x
B
C−1
exp
−
x
B
C
x 0
where x is the magnitude of the complex data sample,B is a distribution scale parameter,and C is a shape parameter.For the example data of Fig.6(i)the low data values(background noi)follow a shape C=2distribution(Rayleigh),whilst the higher clutter values follow a shape C=0.33distribution. T
his plot reprents the global statistics of the RD maps, and that the local statistics for any test cell may contain contributions from both clutter and noi.
To form a CFAR output each test cell x was transformed to a unit Rayleigh distribution using a local shape and scale
parameter:
x Rayleigh=
x
B(x)
C(x)2
(2)
Fig.5.Weibull paper plot of CDF of ARD data.
Prior to CFAR processing the data is assumed to be Doppler shift corrected to align the dominant Bragg lines across the range and azimuth dimensions.
A.Scale Estimation
Local scale was estimated using order statistics.For range spread clutter,the kernel is declared narrow in Doppler with exclusion of the cell under test.The median of the kernel samples is calculated.The median statistic x med= L}of the L local samples has good properties near clutter edges,is robust to target outliers and also range dependent spread Doppler clutter caud by residual meteor echoes.The median values were then smoothed to form a smoothed estimate x smooth.Thefinal scale estimate was taken as a linear combination of the median and smooth estimates:
ˆB=w x
med
+(1−w)x smooth(3) where the sigmoid like function(0 w 1):
w=erf(4log10x med−2)+1
2
provides a convenient methodology to adaptively increa sample size in the more homogeneous noi zones(thereby reducing unnecessary CFAR loss),whilst maintaining edge performance in the high powered clutter zones.
As an example,the scale parameter(3)was estimated for the raw data of Fig.6(i),and is shown in Fig.6(ii).Application of the scaling is displayed in Fig.6(iii),showing suppression of clutter and maintenance of a potential target. This process would then be followed by a cond stage with a processing kernel aimed at reducing fal alarms caud by spread Doppler clutter.It is noted that peak detection may be performed on either the CFAR Fig.6(iii))or the original raw Fig.
6(i)).
Fig.6.Range scale processing,(i)example RD map post IDC processing, (ii)estimated range scale with kernel,and(iii)scaled RD map.
B.Shape Estimation
A number of methods have been reported for estimating the shape parameter of a Weibull distribution.Levanon et al.
[26]suggest using Dubey’s two-sample order statistics shape estimate[27].For optimal parameter lection the standard deviation in the shape estimate is:
ˆσC=
0.916
L
C
So for example the standard deviation is0.43for C=2 and L=20.Unfortunately the estimation variance makes this technique unsuitable for the prent small sample sizes,and would cau additional CFAR loss.
Alternatively,examination of Fig.5reveals that for this data the Weibull shape parameter is a function of the data scale (at least in a global n).Let us lect a smooth transition function to model the Weibull shape distribution as a function of the estimate scale parameter:
C=
C hi−C lo
2
erf
4
B−B lo
B hi−B lo
−2
+1
+C lo(4)
Such a transform limits the Weibull shape parameter to reasonable values.For example tting C hi=2and B lo= 20dB assumes the data is Rayleigh distributed for scales B<20dB.The transformation curve(4)was estimated by lecting the shape that produced a Rayleigh output,averaged over124rad
ar scans.
The CDF of the Fig.5data after Weibull CFAR processing (2)is shown in Fig.7.The data now follows a Rayleigh
distribution,apart from a few outliers that are attributed to either target candidates or residual clutter.
Figure8(i)shows pre-procesd(no signal conditioning) peak history for124radar scans as a function of Doppler. Each pixel reprents the maximum over all procesd ranges and azimuths for that particular Doppler and scan.The peak spectra are dominated by Bragg clutter.The Bragg lines show the Doppler spectra have a time-varying Doppler shift. Occasional jumps are due to either irregular time steps or radar carrier frequency changes.In the background noi zone targets are masked by broad spectra caud by impulsive noi.Figure8(ii)shows the peak history after processing that included IDC registration,impul rejection(bad on the linear interpolation technique[8])and Weibull CFAR.
A number of potential target candidates are obrved.Some clutter residuals also exist indicating the CFAR algorithm does not quite produce the desired uniform noi distribution in the Doppler
dimension.
Fig.7.Weibull paper plot of CDF of data after CFAR processing.
(i)
(ii)
Fig.8.Peak history as a function of Doppler for(i)pre-processing,and(ii) post IDC,impul rejection,and CFAR.
V.T RACK-B EFORE-D ETECT
For HF radar waveform parameters target structure is unre-solved and candidate detections are typically made by locating waveform ambiguity spreading function peaks.At low SNRs fal alarm rates increa and noi introduces a location measurement error.Thus prior to tracking the candidate peak detections are amplitude thresholded,and low SNR targets are lost.A method to improve detection and tracking of the low SNR targets is TBD processing.A host of TBD methods and applications exist.In skywave radar Perlovsky et al.[28] introduced the Maximum Likelihood Adaptive Neural System (MLANS)method for the purpo of aircraft tracking.This iterative method incorporated a function that modeled clutter, targets and noi.As introduced the MLANS method clutter models are not suitable for the ship detection mission.Wallace [29]described a TBD scheme for pul-Doppler radar bad on the Viterbi algorithm with claims of up to10dB SNR gain.Like SIFTER[30]the Wallace method us the whole volume of radar data,rather than propagating target models. The Wallace method is applied here to skywave radar ship detection.The method process RD map data to produce a TBD score s bad on the previous TBD score,the logarithm of current RD map D,and the supported target models: s(d,r,t)=α·max{s(d m,r m,t−1),D(d,r,t)}
+β·D(d,r,t)(5) where r,d,and t are range,Doppler,and scan-time indices respectively,the combination r m,d m denote the M list of RD cells from which a model target is calculated to arrive from,and
βandαcontrol growth and normalization. This scheme has the advantage that the output score dimen-sions are constant,only a single scan history is maintained, and processing load is negligent.
Figure9shows conventional thresholded peak detections on a range-time display for a24minute period of actual skywave radar data.The blue dots indicate estimated peak position for each radar scan,the red vectors indicate range prediction bad on the measured Doppler,and the vertical green lines indicate the depth of a single range cell.Other examples(not shown)have been obrved where the range-rate does not match the Doppler;for the cas it is desired that the TBD scheme defaults to conventional processing. Without a priori knowledge of ionospheric effects the control parameters of(5)were lected to ensure target and fal alarm persistence was limited.Figure10shows the peak detections after TBD processing.The peak detection threshold level was t to provide a similar fal alarm rate to the conventional output of Fig.9.Initialization and maintenance of potential targets1,2,3&4has improved over conventional processing, and potential fal alarms5&6expire as desired.
VI.C ONCLUSIONS
Detection and tracking of maritime targets using skywave radar is influenced by the propagation medium,interference environment and target scenario.The data characteristics were
Fig.9.Conventionally procesd range-time peak
detections.
Fig.10.TBD procesd range-time peak detections.
briefly described,then signal processing techniques were de-veloped to provide robust adaptive Doppler processing,rejec-tion of impulsive noi,improved CFAR using the Weibull distribution with robust two-parameter estimation bad only on scale,and a practical track-before-detect scheme for en-hancing small SNR target detection performance.
A CKNOWLEDGMENT
The author would like to thank the US ROTHR Program Office for the collaborative provision of veral of the data ts analyzed and displayed.
R EFERENCES
[1]J.Headrick and M.Skolnik,“Over-the-horizon radar in the HF band,”
Proc.IEEE,vol.62,no.6,pp.664–673,1974.
[2]R.Barnes,“Automated propagation advice for OTHR ship detection,”
Radar,Sonar and Navigation,IEE Proceedings,vol.143,no.1,pp.
53–63,Feb.1996.
[3]G.Fabrizio,A.Gershman,and M.Turley,“Robust adaptive beamform-
ing for HF surface wave over-the-horizon radar,”IEEE Transactions on Aerospace and Electronic Systems,vol.40,no.2,pp.510–525,Apr.
2004.
[4] D.Sinnott,“Jindalee–DSTO’s over-the-horizon radar project,”Digest
IREECON87,21st International Electronics Convention and Exhibition, Sydney,pp.661–664,Sept.1987.
[5] D.Crombie,“Doppler spectrum of a echo at13.56Mc./s.”Nature,
vol.175,pp.681–682,1955.
[6]J.Parent,“Statistical study of the spectral broadening of skywave
signals backscattered by the a surface:Application to RMS wave height measurement with skywav
e radar,”IEEE Trans.Anten.And Prop., vol.37,no.9,Sept.1989.
[7]G.Fabrizio,L.Scharf,A.Farina,and M.Turley,“Robust adaptive
doppler processing for HF surface wave OTH radar,”Proc.Defen Applications of Signal Processing2004.
[8]M.Turley,“Impulsive noi rejection in HF radar using a linear
prediction technique,” RADAR2003,Sept.2003. [9]M.Turley and D.Netherway,“OTHR signal reconstruction for data
corrupted by impulsive noi,”Proc.Radarcon’90,Adelaide,Apr.1990.
[10]R.Jarrott,D.Netherway,and S.Anderson,“Signal processing for ocean
surveillance by HF skywave radar,”IEEE Australian Symposium on Signal Processing and Applications,Adelaide,Apr.1989.
[11] D.Netherway,G.Ewing,and S.Anderson,“Reduction of some envi-
ronmental effects that degrade the performance of HF skywave radars,”
IEEE Australian Symposium on Signal Processing and Applications, Adelaide,Apr.1989.
[12]Y.Abramovich,S.Anderson,G.Fabrizio,G.Frazer,I.Solomon,and
M.Ringer,“Adaptive signal processing techniques and architectures for HF skywave radar,”Proc.Defence Applications of Signal Processing, 2001.
[13]M.Sacchi,T.Ulrych,and C.Walker,“Interpolation and extrapolation
using a high-resolution discrete Fourier transform,”IEEE Trans.On Signal Processing,vol.46,no.1,Jan.1998.
[14]J.Parent and A.Bourdillon,“A method to correct sky-wave backscat-
tered signals for ionospheric frequency modulation,”IEEE Trans.Anten.
And Prop.,vol.36,no.1,1988.
[15]S.Anderson,“A unified approach to detection,classification,and cor-
rection of ionospheric distortion in HF skywave radar systems,”Radio Science,vol.33,no.4,1998.
[16]J.Barnum,“Ship detection with high-resolution HF skywave radar,”
IEEE Journal of Oceanic Engineering,vol.OE-11,no.2,Apr.1986.
[17] D.Swingler and R.Walker,“Linear-predictive extrapolation for narrow-
band spectral estimation,”Proc.IEEE,vol.76,no.9,Sept.1988. [18]——,“Line-array beamforming using linear prediction for aperture
interpolation and extrapolation,”IEEE Trans.ASSP,vol.37,pp.16–30,1989.
[19]M.Turley and S.V oigt,“The u of a hybrid AR/classical spectral
analysis technique with application to HF radar,”Proc.ISSPA92, International Symposium on Signal Processing and its Applications, Gold Coast,Australia,Aug.1992.
[20]V.Gadwal and J.Krolik,“A performance evaluation of autoregressive
clutter mitigation methods for over-the-horizon radar,”Thirty-Seventh Asilomar Conference on Signal
s,Systems&Computers,vol.1,Nov.
2003.
[21]G.Frazer,“High–resolution Doppler(Hoppler)processing for skywave
radar,”2001,private communication.
[22] A.Nuttall,“Spectral analysis of a univariate process with bad data
points,via maximum entropy and linear predictive techniques,”Naval Underwater Systems Center,Tech.Doc.5303,New London,CT,1976.
[23]M.Cetin,D.Malioutov,and A.Willsky,“A variational technique for
source localization bad on a spar signal reconstruction perspective,”
Proc.IEEE ICASSP2002,Orlando,May2002.
[24]H.Rohling,“Radar CFAR thresholding in clutter and multiple target
situations,”IEEE Trans.Aerosp.Electron.Syst.,vol.AES-19,pp.608–621,July1983.
[25]M.Turley,“Hybrid CFAR techniques for HF radar,”IEE RADAR97,
Oct.1997.
[26]N.Levanon and M.Shor,“Order statistics CFAR for Weibull back-
ground,”IEE Proceedings,vol.137,Pt.F,no.3,June1990.
[27]S.Dubey,“Some percentile estimators for Weibull parameters,”Tech-
nometrics,vol.9,pp.119–129,1967.
[28]L.Perlovsky,V.Webb,S.Bradley,and C.Hann,“Improved ROTHR
detection and tracking using the MLANS algorithm,”Studies in proba-bilistic mult-hypothesis tracking and related topics,ed.R.L.Streit,vol.
SES-98-01,Feb.1998.
[29]W.Wallace,“The u of track-before-detect in pul-Doppler radar,”
Radar2002,pp.315-319,Oct.2002.
[30]L.Nickisch,S.Fridman,and M.Hausman,“SIFTER:Signal inversion
for target extraction and registration–coherent processing of IQ data,”
Naval Surface Warfare Center Final Technical Report MCR/MRY-R-111, July2003.
