01542475

Transcript

4472 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 Detection of Particle Sources With Directional Detector Arrays and a Mean-Difference Test Zhi Liu, Student Member, IEEE, and Arye Nehorai, Fellow, IEEE Abstract—In this paper, the problem of detecting far-field particle sources, such as nuclear, radioactive, optical, or cosmic, is considered. This problem arises in applications including security, surveillance, visual systems, and astronomy. The authors propose a mean-difference test (MDT) with cubic and spherical detector arrays, assuming Poisson distributed measurement models. Through performance analysis, including computing the probability of detection for a given probability of false alarm, the authors show that the MDT has a number of advantages over the generalized likelihood-ratio test (GLRT), such as computational efficiency, higher probability of detection, asymptotic constant false-alarm rate (CFAR), and applicability to low signal-to-noise ratio (SNR). For each array, the authors also present an estimator to find the source direction. Finally, Monte Carlo numerical examples are conducted that confirm the analytical results. Index Terms—Directional detector array, generalized likelihood-ratio test, mean-difference test, nuclear detection, particle source, security. I. INTRODUCTION W E consider the detection of particle sources using directional detector arrays originally proposed in [1]. The sources of interest may be nuclear, radioactive, optical, or cosmic, emitting particles such as rays, neutrons, photons, particles, or particles. An example of potential application is the detection of nuclear sources, which has increasingly become an important issue in security and defense. Another example is the detection of photon illuminants in artificial vision systems. (For details on radiation particle detectors, see [2]–[6], and on photon detectors see [7]–[10].) In [1], we have introduced several directional detector arrays and used them to localize particle sources. We presented projection-based sensor models and assumed Poisson distributed data to derive lower bounds on the mean-square angular error [(MSAE), see also [11]] of the source direction estimation. In this paper (see also [12]), we consider two detector arrays similar to [1], namely, cubic (wherein each side of the cube is a detector) and spherical (with a large number of small-size detectors completely covering a sphere surface). These detectors are directional since their surfaces have different orientations with Manuscript received April 19, 2004; revised January 31, 2005. This work was supported by the Air Force Office of Scientific Research Grant F49620-02-10339 and the National Science Foundation Grant CCR-0330342. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Feng Zhao. The authors are with the Department of Electrical and Computer Engineering, University of Illinois at Chicago, Chicago, IL 60607 USA (e-mail: [email protected]; [email protected]). Digital Object Identifier 10.1109/TSP.2005.859334 respect to the source direction. We employ the measurements of these directional arrays to detect the existence and estimate the direction of the source. In contrast, without directional arrays, many existing techniques can detect particle sources only indirectly, such as through image analysis [13], but may not be capable of estimating the source direction. Estimating the source direction is important in many applications and should also improve the detection performance. We use Poisson distributions to model the detectors’ measurements of the source and noise particles. We first consider the generalized likelihood-ratio test [(GLRT), see [14, Sec. 6.4.2]] for the cubic array. We find that the GLRT has a complex structure in this case (see Section III-A), which makes it difficult to analyze its performance, i.e., find analytically the probability of detection ( ) for a given probability of false alarm ( ) [14]. Moreover, our numerical examples show that the performance of the GLRT may not be acceptable at low signal-to-noise ratio (SNR) values. We then present a new test, the mean-difference test (MDT), based on the difference between the average of the temporal sample means over two groups of detectors, namely, source-exposed and noise-only detectors. We derive the asymptotic distribution of the MDT for the cubic array and analyze its performance. Compared with the numerical results of the GLRT, the MDT has higher probability of detection, especially at low SNR values. In addition, the MDT requires estimating only one parameter, the source direction, which is of interest in its own right. Next, we address the detection using spherical arrays. We define a spherical array as multisurface, wherein each surface corresponds to a detector, the area of each detector approaches zero, and the number of detectors approaches infinity [1]. The detectors fully cover the sphere surface. For spherical arrays, the GLRT may not be analytically computable due to the complexity of the joint likelihood function. In contrast, we derive and analyze the performance of the MDT and show that it outperforms the cubic array with equal surface area. It is computationally simple, applicable at low SNR values, isotropic to the source direction (detection performance remains the same for arbitrary source direction), and has asymptotic constant falsealarm rate (CFAR) [14]. The remainder of the paper is organized as follows. In Section II, we present the basic mathematical model. In Sections III and IV, we derive and analyze the aforementioned tests for both arrays. In Section V, we present numerical examples that confirm our analysis. In Section VI, we discuss possible extensions of our work. 1053-587X/$20.00 © 2005 IEEE LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT 4473 II. BASIC MODELING To simplify the presentation, we consider a single source assumed to be time invariant in the far field. Our model is based on the projection-based arrays originally proposed in [1], which employ the fact that the number of particles impinging on each detector surface in the array is proportional to its cross section. Denote by the number of detectors in the array. Each detector counts the particles impinging on its surface. Assume the count of source particles is a random variable with Poisson distribution; in addition, assume the count of noise particles of each detector is also Poisson distributed and independent of the source particle count. Further assume each detector in the array has a congruent shape, i.e., it is a face of an -surface convex regular polyhedron, and let its surface area be . We also follow the assumptions A1–A4 in [1], namely, narrow-band detectors, far-field source, and isotropic efficiency. Let be the expected flux of the source particles at the detector array, and be the expected rate of the noise counts on the detector surface. More specifically, is the expected count of source particles per unit time crossing a unit area perpendicular to the particles’ direction of arrival at the detector, whereas is the expected count of the noise particles per unit time per unit area on the detector surface. We assume that and are constant for all the detectors in the same array, due to the time-invariance . and far-field assumptions. Define the SNR as Consider the detector array measurements of particles over nonoverlapping time-unit intervals. The th temporal unit measurement of the th detector is denoted as . Under the above can be expressed as modeling and assumptions, (1) where and denote, respectively, the source and noise counts of the detector. Both are random variables, mutually independent, and have Poisson distributions with expected rates and , respectively, where denotes the cross section of the th detector with respect to the source direction. Due and , the measurement is also to the independence of Poisson distributed with expected rate [1] (2) Let be the unit vector at the origin pointing toward the source in some reference frame. Denote , where and are, respectively, the azimuth and elevation angles of . Thus, and . The unknown parameter vector in this case is defined as Fig. 1. Cubic array. III. DETECTION USING CUBIC ARRAYS Consider a projection-based array with a cubical surface wherein each face is an individual detector’s surface. In this , and the area of each face of the cube is . In case general, if a source is present, three detectors are exposed to the source particles. We label the indexes of these detectors as 1, 2, 3, and their opposite detectors as 4, 5, 6, respectively. Without loss of generality, we choose the cube’s axes such that they coincide with the axes of the reference coordinating system, as shown in Fig. 1. To be more specific, denoting the ), unit vector normal to the cube’s th face as ( we choose the cube’s axes and the reference coordinate system are the positive , , and axes, and such that , , and , , and are the negative , , and axes of this system, respectively. Recall the assumption that there are time-unit mea( ; surements of each detector, and here ) denotes one of the temporal measurements. Note that these measurements are statistically independent and have Poisson distributions, and the measurements of an individual detector are independent and identically , where distributed (i.i.d.). In this case, , and “ ” is the inner product. From (2), it follows that obeys the Poisson distribution with parameter . In the unknown source direction case, we do not know which detectors face the source, which causes a problem when labeling the faces of the cube as described above. To solve this problem, we first arbitrarily label the six faces with opposite pairs 1 and 4, 2 and 5, and 3 and 6. Then, we employ the following algorithm to estimate the source direction , based on the relationship between the temporal sample means of the opposite detector pairs [1]: (3) The source detection problem of the detector array can now be simply defined as hypothesis test where the null hypothesis represents the case where there is represents the case no source, and the alternative hypothesis where a source is present. (4) where , ( ) denotes the temporal sample mean of the th deis the Euclidean norm. Finally, using the estitector, and mated source direction, we relabel the faces such that the source is in the first octant. 4474 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 From Fig. 1, we have , , . Thus, , and the two hypotheses can be redefined as where where We thus construct a mean-difference (MD) statistic MD , which is defined as the difference between the temporal sample means averaged over both groups of detectors (7) is unknown where is the temporal sample mean averaged over source-exposed detectors. Therefore, the MD test for the cubic array is constructed as is unknown MD A. GLRT for Cubic Arrays Consider the GLRT for the cubic array GLR (5) where , and and denote the joint likelihood functions of under and , respectively. Through the derivations in Appendix A, we show that the GLRT for the cubic array is equivalent to (8) We analyze the performance of the MDT in (8) as follows, by under deriving first the asymptotic distributions of the MD both hypotheses. : Note that both and Under have the same Poisson distribution with ; by the central limit theorem (CLT) the parameter and should have the same Gaussian sample means , asymptotically (as approaches distribution infinity). Moreover, due to the mutual independence of all the measurements, and should be statistically independent. Then, the asymptotic distribution of MD under is MD GLR (6) where is the temporal sample mean averaged over noise-only detectors, and is the temporal sample mean averaged over all detectors. Observe that it is difficult to analyze the performance of the GLRT, due to its complex structure and non-Gaussian nature of the measurements. As an alternative, in Section V, we present numerical examples to investigate its empirical performance. The results show that the GLRT can detect a source at relatively high SNR values; however, the computations are intensive. In contrast to the GLRT, a novel test with the following properties is desirable: simplicity in computation, tractability in finding the distribution of the test statistic, and having high probability of detection for a given probability of false alarm. In the remainder of this paper, we will present a new test that has these properties and analyze its performance. (9) Under : In this case, has the same asymptotic as under . However, Gaussian distribution ( ,2,3; ) is Poisson distributed with expected rate . Therefore, has Poisson distribution with parameter . has asymptotic Gaussian distribution Thus, by the CLT, . Therefore, the asymptotic distribution of MD under is MD (10) Using the distributions in (9), we are able to compute the threshold in (7) that guarantees the preset probability of false . Since alarm MD MD (11) B. Mean-Difference Test for Cubic Arrays The proposed test is motivated by the fact that since the first three detectors measure signal and noise particles, whereas the last three detectors measure only noise particles, the source can be detected by comparing the averages of the sample means over these two groups of detectors. If the two averages differ by more than a prespecified threshold , we decide that a source is present, and vice versa. where is the probability of an event ; then and (12) where is the inverse function of . LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT Fig. 2. Receiver operating characteristics of the MDT for a cubic array and different source directions u ; SNR = 0:75. Using the distribution in (10), we derive the asymptotic probability of detection for the MDT in (7) as MD 4475 detectors discussed in Section III, if a source is present, would detect both the source and noise particles, and the other would detect only noise particles. For convenience, denote the former detectors as source exposed and the latter as noise only. Then, the detection procedure consists of two steps: first identify the source-exposed detectors, and then construct the MDT accordingly. In the first step, we define a half-array as a half of the array, conand its measurement as the total particle counts of the tiguous detectors in it. In principle, the source-exposed halfarray (wherein all detectors are source exposed) can be identified by choosing the half-array with the highest temporal sample is sufficiently large. Due to the finite length of mean when the temporal measurements, the identified source-exposed halfarray may have a small “bias” from its true version. However, the probability of having this bias will vanish as approaches infinity, since the temporal sample mean is a consistent estimate of the expected rate. In the second step, consider the difference between the temporal sample means of the source-exposed and noise-only halfand arrays, which is expected to be relatively large under small under . Consequently, we construct the MDT as MD MD (13) is proportional to the cube’s cross where section with respect to the source direction. is a function of , we conclude from (13) that the Since performance of the MDT varies with the source direction, i.e., the performance of the cubic array is not isotropic with respect is monoto the source direction. In fact, as shown in (13), tonically increasing with respect to . Subject to the constraint , reaches its peak when and is minimized when only one entry of is 1 and the other two are 0 (in this case, only one face of the cube is exposed to the source particles). Fig. 2 shows the analytical receiver operating characteristics (ROCs) [14] for different values of . As can be seen, the MDT for the cubic array has the optimal detection sensitivity when the source is at the direction . IV. DETECTION USING SPHERICAL ARRAYS We examine the detection of particle sources using spherical arrays. First, we briefly present multisurface arrays as an introduction. A. Multisurface Array Detection We extend the analysis of the cubic array to a multisurface array with congruent faces, wherein each face is a detector. For symmetry, assume to be an even positive integer. Due to the increased complexity of deriving the GLRT, we consider only the MDT for this array. Similar to the cubic array (14) When the detectors in a multisurface array become extremely approaches infinity and the size densely distributed, i.e., of each detector approaches zero, the array asymptotically becomes spherical. We analyze the detection using a spherical array in detail in the following section. B. Spherical Array Detection We consider a spherical array with all its detectors located on the surface of a sphere, wherein a large number of small-size detectors completely cover the sphere surface [1]. The center of the sphere is chosen as the origin of the reference coordinate system. As shown in Fig. 3, the sphere surface can be divided into two halves in the presence of a source: the source-exposed and its complement, the noise-only hemisphere surface . Similar to the multisurface array hemisphere surface discussed previously, the detection using spherical arrays consists of two steps: searching the source-exposed hemisphere surface and then applying the mean-difference testing. 1) Searching the Source-Exposed Hemisphere Surface: Under the previously mentioned assumptions, the measurements for the spherical array are spatial densities of particle counts at all points on the sphere surface, since the area of each detector approaches zero. To clarify further, we redefine the model of the spherical array as follows. and Let be the radius of the sphere, and be the azimuth and elevation angles of a point on the sphere surface, respectively. Therefore, for a given , the doublet defines an arbitrary point on the sphere surface. at a neighborConsider a small area on the sphere surface, whose area hood of a point . Let be the is (approximately) temporal measurement of in the th time unit, and 4476 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 rate of should be the maximum among all hemispheres surfaces on the sphere surface (under ), i.e., Fig. 3. Spherical array. be the temporal sample time-unit intervals. We mean of measurements of during then define the spatial density of the temporal sample mean at as the point (15) and define the corresponding density of the expected rate at as (16) where is the mathematical expectation. on the Further consider an arbitrary hemisphere surface sphere surface. Observe that it is uniquely identified by its center (on its surface). Depoint, i.e., the “symmetric center” of the hemisphere surface with center point at note by . Define the measurement of as the sum of measurements and of all the detectors in . For clarification, let be the temporal sample mean (stochastic, based on ) and the expected rate (deterministic, based on ) of particle counts of , respectively. It is shown in Appendix B can be expressed as1 that (18) Note that the center point of lies on the source direction, . i.e., Furthermore, since the temporal sample mean is a consistent estimate of the expected rate , the source-exposed hemisphere surface should also maxwith high probability as approaches inimize finity, according to (18). In other words, the identified sourcecan be exposed hemisphere surface by maximizing arbitrarily close to the actual one, as approaches infinity. Also provides an estimate note that the procedure of finding of the source direction, which is of interest in its own right. is an integration of the stochastic Observing that variable , in Appendix C we examine the asymptotic distribution of and show that is unbiased and consistent with when is sufficiently large. Below, we present the steps of a gradient-based iterative al(hence also the source direction). Other gorithm to find details of this algorithm are given in Appendix D. Step 1) Let . Denote the point with the maximum as . Then, choose the initial value and . Set a small as positive quantity as the threshold of convergence. . If so, stop the iteraCheck if tion; if not, go to Step 2). Step 2) In the th iteration, we do alternate searching on and , respectively. First, fixing , find a scalar such that Second, fixing , find a scalar such that Here we employ the Fibonacci method2 in deciding and . Then, let and (17) . Step 3) In the and where . is the only hemisphere surface wherein all Obviously, the detectors measure source particles. In principle, the expected 1Without loss of generality, here we assume < 0. A similar expression can be obtained when 0.  and Step 4) If claim otherwise, let th iteration, let . , stop the iteration and is the desired point ; and go back to Step 2). 2http://math.fullerton.edu/mathews/n2003/FibonacciSearchMod.html LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT The definition of is given in (D1), and gradiand are expressed in ents (D2) and (D5). We denote the searching result of the source, an estiexposed hemisphere surface as mate of , where is the corresponding estimate of . Let be the complement of the source direction . 2) MDT for Spherical Arrays: We construct the MDT for the spherical array as MD 4477 then (25) Thus, the asymptotic probability of detection is MD (19) where and are the temporal sample means of particle counts of and , respectively. In the following part of this section, we analyze the performance of the MDT for the spherical array in two cases, as follows. test for known a) Performance analysis of the MD source direction: The performance of the MDT in this ideal case can be used as an upper bound of its performance in and be the temporal sample general cases. Let and . Since in this means of particle counts of case and , the MDT in (19) becomes MD (20) Under , since only noise particles are measured, within a time-unit interval, the measurements of and obey the same Poisson distribution with parameter . Then asymptotically, by the CLT it follows (21) and are statistically independent. Hence, Note that in this case is under , the asymptotic distribution of MD (22) MD , since measures only noise particles, the Under asymptotic distribution of remains the same as in (21). with respect to On the other hand, the cross section of . Taking into account the noise parthe source direction is ticles, within a time-unit interval, the measurement of obeys a Poisson distribution with parameter . has the asymptotic Gaussian distriThus, by the CLT, . In addition, bution and are statistically independent. Therefore, under , the asymptotic distribution of MD in this case is MD (23) (26) Observe first from (22) that the asymptotic distribution of the has as the only unknown parameter, MD statistic under which makes it difficult to compute in (25). Hence, if we replace in (22) with its maximum likelihood estimator (MLE) , this distribution will not depend on any unknown parameters. Moreover, we can expect that the performance of the MDT after this replacement would not differ much is large, since is a consistent estifrom that in (24) when mate of . Then, the threshold required to maintain a constant can be found from (25) with replaced by , i.e., the MDT for the spherical array is asymptotically CFAR, which illustrates the capability of the MDT to detect the source in unknown environments. Second, observe that as expected, the detection with does not depend on the the spherical array is isotropic since source direction. test for unknown b) Performance analysis of the MD source direction: In this case, we take into account the pos. Denote by the error sible small error in searching for angle between the actual source direction and its estimate . We assume that is stochastic and has probability . From the isotropic nature of the density function (PDF) spherical array, we further assume that does not depend on or . Now consider the conditional distribution of MD in (19) given . , since no source is present, for arbitrarily chosen Under hemisphere surface and its complement , the MD statistic has the same distribution as in (22). Note that this distribution does not depend on . In this case, the estimation error has no effect on the distribution of MD . Under , the cross section of , the estimated , with respect to the source direction is . Therefore, the expected count of source particles that impinge is , whereas that of the noise particles is on . According to the Poisson data model, the PDF given is From (22), it follows that MD MD Since (24) particle count of the spherical array 4478 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 the conditional PDF of given is shown in the equation at the bottom of the page. Observe that the conditional distriMD is also Poisson bution of with parameter . By the CLT, as approaches infinity, the asymptotic conditional distribution of given is MD MD (27) and . where In the following, we analyze the performance of the MDT in this case, according to the distributions in (22) and (27). Since depends only on the distribution of MD under , which is the same in both known and unknown source direction cases, the threshold in (19) and (20) should also be the same, which , we calculate as (28), is given in (25). Then, for a given shown at the bottom of the page. MD Fig. 4. Numerical histograms of the log-GLR for a cubic array for various SNR values: (a) SNR = 0:5; (b) SNR = 0:75; (c) SNR = 1; and (d) SNR = 2. MD (28) In practice, however, usually it is difficult to know thoroughly , which makes it the distributions of the estimation errors in (28) in this case. difficult to evaluate Despite the difficulty in analytically computing , we observe that since the distribution of the MD statistic remains the same as the known source direction case, the MDT for the spherical array is also CFAR when the source direction is unknown. is independent of and , we Moreover, due to the fact that does not depend on the source direction due to the claim that nonappearance of and in (28), i.e., the MDT for the spherical array also has isotropic detection performance in this case. V. NUMERICAL EXAMPLES In this section, we first employ Monte Carlo (MC) simulations to examine the experimental performance of the GLRT for the cubic array. Then, we compute the analytical performance of the MDT for both arrays and compare them with the experimental performance. In these MC simulations, we conduct a total of 10 000 experiments; in each single experiment, 100 time-unit measurements for each detector are generated 100). In the following numerical examples, we choose ( and . We also let the Fig. 5. Numerical receiver operating characteristics of the GLRT for the cubic array for various SNR values. cubic and spherical arrays have equal surface area to make their performance comparable. A. Numerical Examples for Cubic Arrays Consider first the numerical performance of the GLRT for the cubic array. We select the threshold in (5) using MC simulations for a given . Figs. 4 and 5 show the experimental histograms and ROCs of the GLRT, respectively, for various SNR values: 0.5, 0.75, 1, and 2. As shown in Fig. 4, the GLRT statistic has non-Gaussian numerical distributions under both hypotheses, as expected. We further observe that the distributions of the GLRT statistic under both hypotheses are too similar . From Fig. 5, we to enable good detection when SNR LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT 4479 Fig. 8. Probability of miss P of the GLRT and the mean-difference test as functions of the SNR for the cubic array; P = 0:05. Fig. 6. Numerical and analytical distributions of the mean-difference statistic for the cubic array for various SNR values: (a) SNR = 0:5; (b) SNR = 0:75; (c) SNR = 1; and (d) SNR = 2. Fig. 7. Receiver operating characteristics of the mean-difference test for the cubic array for various SNR values: analytical and numerical results. observe that the GLRT has a relatively poor performance for an SNR smaller than 1. These experimental results illustrate that the performance of the GLRT for the cubic array may not be satisfactory, especially at low SNR values. The analytical and numerical performance results of the MDT for the cubic array are shown in Figs. 6 and 7. Fig. 6 shows a close match between the numerical and theoretical distributions of MD , which verifies the asymptotic Gaussian distributions given in (9) and (10). Through the ROC curves at different SNR values shown in Fig. 7, we can see that the MDT for the cubic array performs quite well even at low SNR values such as 0.5 ) and 0.75. By depicting the probability of miss ( , Fig. 8 compares the versus SNR curves for a given performances of the MDT and the GLRT, wherein the difference of the detection performances between both tests can be observed. Fig. 9. Numerical and analytical distributions of the mean-difference statistic for the spherical array and known source direction, for various SNR values: (a) SNR = 0:5; (b) SNR = 0:75; (c) SNR = 1; and (d) SNR = 2. We then conclude from the above results that the MDT obviously outperforms the GLRT, especially at low SNR values. For and , of the GLRT is instance, as SNR about 0.6, whereas that of the MDT is 0.45. (The measurement data used for the above analysis are the same.) B. Numerical Examples for Spherical Arrays 1) Known Source Direction Case: Consider the performance of the MDT for the spherical array in the known source direction case. We first generate temporal measurements , and , . Then, the temporal sample means are and . Construct the MDT according to (20). Fig. 9 shows the analytical and experimental distribution of MD , in which we can see that the numerical histograms 4480 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 Fig. 10. Receiver operating characteristics of the MDT for the spherical array and known source direction, for various SNR values: analytical and numerical results. match the analytical distributions (given in (22) and (23)) closely, under both hypotheses. It can be further observed that at relatively low SNR values such as 0.75, the distribution of and are separate enough to the MDT statistic under allow a relatively high . As shown in Fig. 10, the analytical and numerical ROCs of the MDT for the spherical array match closely to each other.3 We also observe that increases quickly with and becomes . The high of the MDT at low close to one when SNR values illustrates its applicability to detection in a noisy environment. 2) Unknown Source Direction Case: We present the numerical performance results of the MDT for the spherical array with unknown source direction. Here, we approximate a spherical 60 faces [15]. array with a multisurface array with Fig. 11 presents the numerical histogram of the MD statistic for the spherical array in this case, whereas Fig. 12 shows the ROCs of the MDT at different SNR values. In Fig. 13, we compare the probability of miss of the MDT between known and . As we unknown source direction cases, for a fixed can see, compared with the MDT in the known source direction case, the MDT with unknown source direction performs simi. Also observe that when SNR , larly well when SNR it does not perform as well due to the increase in the error of the source direction estimation. VI. CONCLUDING REMARKS We presented mean-difference tests for detecting particle sources using cubic and spherical arrays and analyzed their performance. For the cubic array, we showed that the MDT outperforms the GLRT by its higher probability of detection for at low SNR values. For the spherical array, whereas a fixed the GLRT may be too complex to derive, we have shown that the MDT has desired properties such as computational simplicity, 3The deviation between the numerical and analytical ROC curves as P is less than 0.1 is due to the “quantization error” in the numerical analysis. This deviation can be reduced by increasing N . Fig. 11. Numerical histograms of the mean-difference statistic for the spherical array in the unknown source direction case, for various SNR values: (a) SNR = 0:5; (b) SNR = 0:75; (c) SNR = 1; and (d) SNR = 2. Fig. 12. Numerical receiver operating statistics of the mean-difference test for the spherical array in the unknown source direction case, for various SNR values. being asymptotically CFAR, applicability to low SNR cases, and isotropic detection performance. As part of the detection, we also estimated the source direction, which is of interest in many cases and should improve the detection performance. In the performance analysis of the MDT for spherical arrays, we compared the analytical with the experimental results and confirmed that they match very well. Further research will include analyzing the performance of the proposed source direction estimators, investigating arrays with lenses, and optimally designing arrays by maximizing their detection performance. Some practical aspects, such as including the effect of a specific shape of a mounting part, will LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT 4481 (A4) (A5) where . From (A5), we have (A6) Substituting (A6) into (A2), we get the maximized likelihood as function under Fig. 13. Probability of miss P as a function of SNR, of the mean-difference test for the spherical array with both known and unknown source direction; = 0:05. P be considered in the future. The results of this paper can serve as a benchmark of the best performance when this part does not affect the performance. APPENDIX A DERIVATION OF THE GLRT FOR THE CUBIC ARRAY (A7) Then, the GLRT follows according to (5), (A1), and (A7): GLR (A8) Further simplification of (A8) leads to the GLRT-equivalent test for the cubic array given in (6). To derive the generalized likelihood-ratio test (GLRT) for cubic arrays, we maximize the likelihood function under both hypotheses. As assumed, each unit-time measurement of the th detector is Poisson distributed with parameter , i.e., . A. Under All the measurements are i.i.d. The likelihood function . Letis , we get , where ting . Then, the maximized likelihood becomes function under (A1) B. Under APPENDIX B DERIVATION OF Here, we derive the expression of in (17). As shown in Fig. 14, under the assumption that , there is a point on the surface of which has the largest coordinate, . It is whose azimuth angle is . Denote the point as . easy to show that For simplification, we project the boundary of onto the plane, and let denote the resulting ellipse, as shown in Fig. 15. In a two-dimensional (2-D) polar coordinate system, let and an points and be the projections of the point arbitrary point on onto the plane, respectively. one of the points of intersection between ’s Denote by plane. Obviously, and are boundary circle and the the semimajor and semiminor axes of , respectively. In addition, , and . From the property , , of the orthogonal projection, we have . and be the coordinate of point in a 2-D polar system, Let where is the distance from the origin and is the azimuth angle. It then follows that (A2) (B1) Taking the derivative of , , we get with respect to , , , where (A3) . The two solutions of are 4482 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 or and (B3) Integrating over the region subject to the constraints in (B3), we get the time-averaged measurement of hemisphere in (17). surface APPENDIX C ASYMPTOTIC DISTRIBUTION OF Fig. 14. Hemisphere surface H with center point at ( ; We derive the asymptotic distribution of to investigate the conditions under which it is consistent with . For a on the source-exsmall area in a time-unit posed hemisphere surface, its measurement interval obeys a Poisson distribution with parameter , where is the angle between the source diand the radius vector pointing to . Note that rection is a function of and for fixed and . By the CLT, the temporal sample mean of measurement of during time-unit intervals has the following asymptotic distribution: ). From (15), the asymptotic distribution of is (C1) Fig. 15. Projection of the boundary of H onto the XY plane. is And the expected particle counting density at ; here, belongs to the source-ex. posed hemisphere surface and thus is asymptotically unbiased to As can be seen in (C1), , while the variance of depends on the product of and . Assume , where is a positive real quantity, then we find the following. and Then, the two solutions of are a) (B2) stant. b) From (B2), for any point that belongs to the hemisphere surface , we have and . Also observe from Fig. 14 that when . In summary, a point belongs to the if and only if it satisfies the folhemisphere surface lowing constraints: and : c) , where is a conhas a finite variance : variance, which is not desired. : . . has an infinite . The variance of is 0, asymptotically. This is the case that we desire. is large From the above discussion, we conclude that if enough such that , the density of the temporal sample mean is unbiased and consistent with the density of the expected rate . LIU AND NEHORAI: DETECTION OF PARTICLE SOURCES WITH DIRECTIONAL DETECTOR ARRAYS AND A MDT APPENDIX D OTHER DETAILS OF THE ITERATIVE ALGORITHM TO FIND where 4483 , , , and . Hence, it follows from (D3) that In Section IV-B-1), we present the steps of an iterative alby maximizing .4 Here, we gorithm to find discuss other details of the algorithm, including initializing the , and how to preparameters, deriving the gradient of vent the algorithm from converging to local extrema. is a single peak function and First note that . From the consistency reaches its maximum when is likely to have a single peak as property, becomes large. Therefore, we can maximize it by searching for . the zero of the gradients of , as , we Since choose the initial value of as the point with the maximum value, since is a consistent estimate of when is sufficiently large. is a suffiAssume the temporal sample mean and , such that its first partial ciently “smooth” function of derivatives exist almost everywhere on the sphere surface. This assumption is valid when is large enough. Define (D1) Obviously (D2) To derive the partial derivative of recall that with respect to , we (D3) and are differentiable functions of . Then, we where rewrite (17) into the following form: (D4) 4This algorithm depends on the assumption of < 0, to be consistent with the assumption of < 0 in (17). This will not cause loss of generality since we can first assume  < 0 and apply the following algorithm, and then assume  0 and apply a similar algorithm. Combination of both results provides a thorough search of the source direction.  (D5) 4484 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 12, DECEMBER 2005 The notations and indicate that they are and . functions of To avoid converging to local extrema (which is the case where has multiple peaks), we repeat Steps 1)–4) with multiple initial values of , and then compare the rewith the maximum . sults and find the point ACKNOWLEDGMENT The authors would like to thank Dr. E. Paldi for his help in simplifying (6) as well as helpful comments on presentation. REFERENCES [1] A. Nehorai and E. Paldi, “Localization of particle sources with detector arrays,” in Proc. 27th Asilomar Conf. Signals, Systems, Computers, Pacific Grove, CA, Nov. 1993, pp. 981–985. [2] J. B. Bricks, The Theory and Practice of Scintillation Counting. Oxford, U.K.: Pergamon, 1964. [3] E. Fenyves and O. Haiman, The Physical Pronciples of Nuclear Radiation Measurements. New York: Academic, 1969. (German). [4] Semiconduct Detectors, G. Bertolini and A. Coche, Eds., Elsevier North Holland, Amsterdam, The Netherlands, 1968. [5] K. Kleinknecht, Detectors for Particle Radiation. Cambridge, U.K.: University Press, 1986. [6] G. F. Knoll, Radiation Detection and Measurement, 3rd ed. New York: Wiley, 2000. [7] R. H. Kingston, Detection of Optical and Infrared Radiation. Berlin, Germany: Springer-Verlag, 1978. [8] R. W. Boyd, Radiometry and the Detection of Optical Radiation, 2nd ed. New York: Wiley, 1983. [9] B. E. Saleh and M. C. Teich, Fundamental of Phontonics. New York: Wiley, 1991. [10] A. Yariv, Optical Electronics, 4th ed. New York: Holt, Rinehart and Winston, 1991. [11] A. Nehorai and E. Paldi, “Vector-sensor array processing for electromagnetic source localization,” IEEE Trans. Signal Process., vol. 42, no. 2, pp. 376–398, Feb. 1994. [12] Z. Liu and A. Nehorai, “Detection of particle sources with directional detector arrays,” in Proc. 3rd IEEE Sensor Array Multichannel Signal Processing Workshop, Sitges, Spain, Jul. 2004, p. 5. [13] P. C. Schaich et al., “Automatic image analysis for detecting and qualifying gamma-ray sources in coded-aperture images,” IEEE Trans. Nucl. Sci., vol. 43, no. 4, pp. 2419–2426, Aug. 1996. [14] S. M. Kay, Fundamentals of Statistical Signal Processing, Vol. II: Detection Theory. Englewood Cliffs, NJ: Prentice-Hall, 1998. [15] Matlab 6.5 Help File: The Bucky Ball. [Online]. Available: http://www. mathworks.com/access/helpdesk/help/techdoc/math/sparse12.html Zhi Liu (S’04) received the B.Sc. degree in electrical engineering from the North China University of Technology (NCUT), Beijing, China, in 1996 and the M.Sc. degree in information and control engineering from Tongji University, Shanghai, China, in 2001. He is currently working toward the Ph.D. degree in the Electrical Engineering Department at the University of Illinois at Chicago (UIC). His research interests are in statistical signal processing and its application to sensor-array processing. He is also interested in image processing, medical imaging (e.g., PET and SPECT), and biological vision. Mr. Liu has been a member of the Phi Kappa Phi honor society by election of the UIC Chapter since April 2004. He received the Siemens Prize for outstanding graduate study in 2000. Arye Nehorai (S’80–M’83–SM’90–F’94) received the B.Sc. and M.Sc. degrees in electrical engineering from The Technion—Israeli Institute of Technology, Haifa, and the Ph.D. degree in electrical engineering from Stanford University, Stanford, CA. From 1985 to 1995, he was a Faculty Member with the Department of Electrical Engineering at Yale University. In 1995, he joined the Department of Electrical Engineering and Computer Science at the University of Illinois at Chicago (UIC), as a Full Professor. From 2000 to 2001, he was Chair of the department’s Electrical and Computer Engineering (ECE) Division, which is now a new department. In 2001, he was named University Scholar of the University of Illinois. Dr. Nehorai has been a Fellow of the Royal Statistical Society since 1996. He is Vice President (Publications) of the IEEE Signal Processing Society (SPS). He is also Chair of the Publications Board and a Member of the Board of Governors and the Executive Committee of this Society. He was Editor-in-Chief of the IEEE TRANSACTIONS ON SIGNAL PROCESSING from 2000 to 2002 and is the Founding Editor of the special columns on Leadership Reflections in the IEEE SIGNAL PROCESSING MAGAZINE. He was co-recipient of the IEEE SPS 1989 Senior Award for Best Paper with P. Stoica, as well as coauthor of the 2003 Young Author Best Paper Award and of the 2004 Magazine Paper Award with A. Dogandzic. He was elected Distinguished Lecturer of the IEEE SPS for the term 2004 to 2005.