Calculation of Radar Probability of Detection in K-Distributed Sea Clutter and Noise Stephen Bocquet Joint Operations Division Defence Science and Technology Organisation DSTO-TN-1000 ABSTRACT The detection performance of maritime radars is usually limited by sea clutter. The K distribution is a well established statistical model of sea clutter which is widely used in performance calculations. There is no closed form solution for the probability of detection in K-distributed clutter, so numerical methods are required. The K distribution is a compound model which consists of Gaussian speckle modulated by a slowly varying mean level, this local mean being gamma distributed. A series solution for the probability of detection in Gaussian noise is integrated over the gamma distribution for the local clutter power. Gauss-Laguerre quadrature is used for the integration, with the nodes and weights calculated using matrix methods, so that a general purpose numerical integration routine is not required. The method is implemented in Matlab and compared with an approximate solution based on lookup tables. The solution described here is slower, but more accurate and more flexible in that it allows for a wider range of target fluctuation models. RELEASE LIMITATION Approved for public release Published by Joint Operations Division DSTO Defence Science and Technology Organisation Fairbairn Business Park Department of Defence Canberra ACT 2600 Australia Telephone: (02) 6265 9111 Fax: (02) 6265 2741 © Commonwealth of Australia 2011 AR 014-977 April 2011 APPROVED FOR PUBLIC RELEASE Calculation of Radar Probability of Detection in K-Distributed Sea Clutter and Noise Executive Summary Modelling of the radar returns from the sea is required for operations analysis of maritime patrol, in order to calculate probabilities of detection for targets of interest. The detection performance of maritime radars is usually limited by sea clutter. A statistical model of the clutter is normally used, and the K distribution has become standard for this purpose. There is no closed form solution for the probability of detection in K-distributed clutter, so numerical methods are required. A method for calculation of the probability of detection in K-distributed clutter and noise is described, and compared with other approaches. It is faster than direct numerical integration, and more flexible and more accurate, but slower, than an approximate method based on interpolation. Faster evaluation is highly desirable for simulation models used in operations research. The method described here could be used to create tables for interpolation, which should produce the appropriate mix of speed and accuracy for operations research simulation models. UNCLASSIFIED Contents 1. INTRODUCTION...............................................................................................................1 2. PROBABILITY OF DETECTION.....................................................................................1 2.1 The K Distribution...................................................................................................1 2.2 Probability of Detection Calculation....................................................................2 2.3 Integration over the Clutter Power........................................................................4 2.4 Normalisation............................................................................................................5 2.5 Probability of False Alarm......................................................................................6 3. SAMPLE CALCULATIONS..............................................................................................7 3.1 Probability of False Alarm......................................................................................7 3.2 Probability of Detection..........................................................................................7 3.3 Target Fluctuation Models......................................................................................9 3.4 Convergence of the Integral over the Clutter Power..........................................9 3.5 Rayleigh Limit.........................................................................................................12 3.6 Comparison with the Method of Watts and Wicks..........................................12 3.7 Interpolation Scheme.............................................................................................13 4. POSSIBLE EXTENSIONS...............................................................................................14 4.1 Constant False Alarm Rate Processing...............................................................14 4.2 KK Distribution......................................................................................................14 4.3 Clutter Correlations................................................................................................14 4.4 Target Models..........................................................................................................14 5. CONCLUSION..................................................................................................................15 6. ACKNOWLEDGEMENTS..............................................................................................15 7. REFERENCES....................................................................................................................16 APPENDIX A: THE SERIES EXPANSION FOR THE PROBABILITY OF DETECTION IN NOISE.........................................................................18 APPENDIX B: CHERNOFF BOUNDS............................................................................19 APPENDIX C: INTERPOLATION SCHEME................................................................22 APPENDIX D: EFFECTIVE SHAPE PARAMETER......................................................26 UNCLASSIFIED UNCLASSIFIED DSTO-TN-1000 1. Introduction Modelling of the radar returns from the sea is required for operations analysis of maritime patrol, in order to calculate probabilities of detection for targets of interest. In this context the backscatter from the sea is referred to as sea clutter. A statistical model of the clutter is normally used, and the K distribution has become standard for this purpose [1]. Sea clutter consists of a rapidly varying speckle component and an underlying mean amplitude which varies more slowly. The K distribution provides a compound representation which includes both components. Thermal noise, which has Gaussian statistics, must also be considered in calculating probabilities of detection. Fluctuations in the target return also need to be taken into account, and this is normally done with statistical models based on the gamma distribution. The calculation of the probability of detection in K-distributed clutter and noise is quite difficult. An approximate method was developed by Watts and Wicks [2,3]. This method is fast, because the probability of detection is obtained by linear interpolation in a table, following some preliminary calculations which need only be done once for each set of parameters. However, the result is only approximate, and the necessary tables of coefficients [3] are only available for the Swerling 1 and 2 target fluctuation models, as well as a non- fluctuating target (Swerling 0). Alternatively, the probability of detection can be calculated by numerical integration over the K distribution. This will give accurate results, but it is slow, and therefore impractical for use in simulations where the probability of detection must be updated frequently. Another approach is described in [1], and further explored in this report, in which the K distribution is separated into its components and only the integration over the local clutter power is carried out numerically. 2. Probability of Detection 2.1 The K Distribution The K distribution for the clutter intensity can be expressed as ∞ P(z)=∫ P(z|x)P(x)dx (1) c 0 where 1 P(z|x)= exp(−z x) (2) x and bνxν−1 P(x)= exp(−bx) (3) c Γ(ν) Equation (3) is the Gamma distribution for the local clutter power x. The mean power is UNCLASSIFIED 1 UNCLASSIFIED DSTO-TN-1000 p = x =ν b (4) c Equation (1) evaluates to 2b ( )ν−1 ( ) P(z)= bz K 2 bz (5) Γ(ν) ν−1 The modified Bessel function K in equation (5) gives the distribution its name. However, for calculating the probability of detection it is advantageous to keep the speckle and modulation components separate, assuming that the local clutter power x remains constant over the beam dwell time. 2.2 Probability of Detection Calculation The probability of detection in Gaussian noise and speckle is ∞ P (Y |x,N)=∫ P (µ|s,N)dµ (6) d R Y The threshold Y is normalised by the noise and local clutter power x. The sum of N radar returns, normalised by the noise and local clutter power is 1 N µ= ∑z (7) x+ p i n i=1 The sum of the target powers from the N pulses is 1 N s = ∑A2 (8) x+ p i n i=1 The probability density function for the sum µ of N radar returns from a fixed target in Gaussian noise is the multilook Rice distribution, (N−1)/2 ⎛µ⎞ ( ) P (µ|s,N)= e−(µ+s)I 2 µs ⎜ ⎟ R ⎝ s ⎠ N−1 (9) ∞ e−µµN+i−1 e−ssi =∑ (N +i−1)! i! i=0 The series form of equation (9) is obtained using the series expansion for the Bessel function [4, §8.445]: ∞ 1 ⎛ z⎞n+2i I (z)=∑ ⎜ ⎟ (10) n i!Γ(n+i+1)⎝2⎠ i=0 UNCLASSIFIED 2 UNCLASSIFIED DSTO-TN-1000 The integral in equation (6) can be evaluated to give the double sum formula for the probability of detection for a fixed target [5,6,7]: ∞ ∞ e−ssi 1 ∞ ∫ P (µ|s,N)dµ=∑ ∫ e−µµN+i−1dµ Y R i! (N +i−1)! Y i=0 ∞ e−ssi N+i−1Y j =∑ e−Y ∑ i! j! i=0 j=0 ∞ e−ssi ⎛N−1Y j N+i−1Y j ⎞ =∑ e−Y ⎜∑ + ∑ ⎟ (11) i! j! j! i=0 ⎝ j=0 j=N ⎠ N−1 Y j ∞ e−ssi N+i−1 Y j =∑e−Y + ∑ ∑ e−Y j! i! j! j=0 i=0 j=N N−1 Y j ∞ Y j ⎛ j−N e−ssi ⎞ =∑e−Y + ∑e−Y ⎜1− ∑ ⎟ j! j!⎝ i! ⎠ j=0 j=N i=0 The manipulation of the series in equation (11) makes use of the result ∞ e−ssi ∑ =1 (12) i! i=0 Target fluctuations are modelled using the Gamma distribution for the sum of target returns s: sk−1 ⎛ k ⎞k P (s|S,k)= ⎜ ⎟ e−ks/S (13) T Γ(k)⎝S ⎠ where N A2 S = (14) x+ p n The Gamma distribution encompasses all the Marcum-Swerling and Weinstock target models [8],[9, Section B] as shown in the table below: Weinstock 0 < k < 1 Swerling 1 k = 1 Swerling 2 k = N Swerling 3 k = 2 Swerling 4 k = 2N Swerling 0 k → ∞ A non-fluctuating target (Swerling 0) corresponds to k → ∞ or s = S. UNCLASSIFIED 3 UNCLASSIFIED DSTO-TN-1000 The probability of detection for a fluctuating target is [9, Section B] ∞ ∞ P (Y |x,S,k,N)= ∫ P (s|S,k)∫ P (µ|s,N)dµds d T R 0 Y =∑N−1e−Y Y j + ∑∞ e−Y Y j ⎛⎜1− ∑j−N 1 ⎛⎜ k ⎞⎟k ∫∞e−s(1+k/S)si+k−1ds⎞⎟(15) j! j!⎜ Γ(k)i!⎝S ⎠ 0 ⎟ j=0 j=N ⎝ i=0 ⎠ N−1 Y j ∞ Y j ⎛ j−N Γ(k+i)⎛ k ⎞k ⎛ S ⎞i⎞ =∑e−Y + ∑e−Y ⎜1− ∑ ⎜ ⎟ ⎜ ⎟ ⎟ j! j!⎜ Γ(k)i! ⎝S+k ⎠ ⎝S+k ⎠ ⎟ j=0 j=N ⎝ i=0 ⎠ Numerical calculation of the probability of detection from equations (11) and (15) is not straightforward, due to underflow and loss of significance. Appendix A gives some insight into the nature of the series in these equations, and shows why it is not always possible to begin the summation with the first term, even with high precision arithmetic. Shnidman [5,6,7,9,10] has developed methods to carry out these calculations efficiently, and coded the resulting algorithms in Matlab 1. Chernoff bounds (Appendix B) are used to avoid unnecessary calculations. 2.3 Integration over the Clutter Power The probability of detection in K-distributed clutter with a single pulse threshold y and local n clutter power x is ∞ P (y ,N)=∫ P (Y |x,N)P(x)dx (16) d n d c 0 where P (Y|x) is calculated as for detection in noise, equation (6), and P(x) is the Gamma d c distribution (3). The threshold y is normalised by the sum of the noise power p and the mean n n clutter power p. At first sight, numerical integration over the local clutter power seems c unlikely to provide a fast and efficient method for calculating the probability of detection. However, if we make the substitution t = bx in the Gamma distribution (3) this becomes tν−1e−t P(x)dx= P(t)dt = dt (17) c c Γ(ν) In terms of the mean clutter power p, t = νx/p. Equation (16) becomes c c 1 ∞ P (y ,N)= ∫ P (Y |t,N)tν−1e−tdt (18) d n Γ(ν) 0 d This integral is readily evaluated using Gauss-Laguerre quadrature [11, Section 4.5]. The Laguerre polynomials are generated from a recurrence relation, and the nodes and weights are calculated from the eigenvalues and eigenvectors of a symmetric tridiagonal matrix. Gautschi [12] has written Matlab code to do this. In most cases sufficient accuracy in P is d 1 Private communication. UNCLASSIFIED 4 UNCLASSIFIED DSTO-TN-1000 obtained with only 10 quadrature points. The exception is low thresholds and ν < 1, where many more points are needed to cope with the singularity at the origin. However, this situation is not of much practical importance, because a high threshold is needed to achieve a satisfactory probability of false alarm when ν is small. The maximum value of ν for which this method of integration can be used in Matlab is 171, because it requires Г(ν) to be less than the largest floating point number which can be represented, i.e. Г(ν) < 10308. As ν → ∞ the clutter becomes noise like, so integration over the Gamma distribution is not required. There is a very small, but perceptible, difference between the probability of detection calculated for ν = 170 and the Rayleigh limit where the clutter is treated as noise. 2.4 Normalisation The threshold Y is normalised by the noise and local clutter power. In terms of a single pulse, unnormalised threshold y Ny Y = (19) x+ p n We would like to express this in terms of the clutter to noise ratio CNR and the integration variable t. y N(p + p ) Y = c n p + p tp c n c + p ν n (20) N(1+CNR) = y n t 1+ CNR ν where y is the single pulse threshold normalised by the mean clutter plus noise power. n In the same way we would like to express the signal S in terms of the signal to noise ratio SNR, the clutter to noise ratio CNR and the integration variable t. From equation (14), N A2 S = x+ p n N A2 = (21) tp c + p ν n N SNR = t 1+ CNR ν UNCLASSIFIED 5 UNCLASSIFIED DSTO-TN-1000 For the clutter only case we have p = 0 or CNR → ∞ so n Nν Nν Y = y and S = SCR (22) n t t where SCR is the signal to clutter ratio. 2.5 Probability of False Alarm The probability of false alarm is obtained from equation (11) with s = 0 (note that 00 ≡ 1) N−1 Y j Γ(N,Y) P (Y |x,N)=∑e−Y = (23) fa j! Γ(N) j=0 In terms of the Matlab incomplete gamma function, P = 1 – gammainc(Y,N). fa In order to calculate the probability of detection for a given probability of false alarm, the threshold must be determined from the probability of false alarm. This can be done using the Matlab root finding function fzero. This function requires either a starting estimate for the root, or an interval which brackets the root. In the first case fzero searches for the starting interval itself, which tends to cause it to fail, because the probability of false alarm is undefined for negative thresholds, and underflows to zero for large positive thresholds. A starting interval of 0.1 to 1100 for y was found to be satisfactory. Integration over the local n clutter power does not impede the calculation much, because the nodes and weights depend only on ν, so they don’t need to be recalculated during the root finding process. In the case of a single pulse return from clutter, the integral over the clutter power can be evaluated to give 1 ∞ P (y )= ∫ P (Y |t,1)tν−1e−tdt fa n Γ(ν) 0 fa (24) 2(νy )ν/2 ( ) = n K 2 νy Γ(ν) ν n This solution is used in Section 3.4 to assess the accuracy of the numerical integration with different numbers of quadrature points. UNCLASSIFIED 6