Using an MM-principle to achieve fast image data estimation from large image data sets

Abstract

A majorize-minimize (MM) mathematical principle is applied to least squares regularization estimation problems to effect efficient processing of image data sets to provide good quality images. In a ground penetrating radar application, these approaches can reduce processing time and memory use by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse.

Claims

1. A method of creating an image from an initial data set, the method comprising: receiving the initial data set, the receiving comprising receiving data representing transmission site locations of radar pulses applied into ground in a ground penetrating radar application, reception site locations of reception of reflections from the radar pulses, radar-return profiles for pairings of the transmission site locations and the reception site locations, and data samples associated with individual radar-return profiles; processing the initial data set with a processing device by creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem associated with a mathematical model of image data from the initial data set; displaying the image using the estimated image value of individual voxels of the image.

2. The method of claim 1 further comprising iteratively deriving an estimated image value for a given voxel using
for l=1,2, . . . ,L $x_l^{\left(m+1\right)}=\frac{H_l{x}_l^{\left(m\right)}+G_l^{\left(m\right)}}{H_l+\frac{}{2.Math.x_l^{\left(m\right)}.Math.}}$ $\mathrm{where}$ $H_l=\underset{k=0}{\overset{K-1}{.Math.}}{r}_k{A}_{\mathrm{kl}}^2$ $G_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}{A}_{\mathrm{kl}}\left(y_k-{\left[{\mathrm{Ax}}^{\left[m\right]}\right]}_k\right)$ and where x.sub.1.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at an l.sup.th voxel for iterate number (m), A is a KL system matrix associated with the mathematical model of the image data, K=IJL where I is a number of transmission site locations, J is a number of reception sites for each transmit site location, L is a number of pixels in the image, r.sub.k is a number of nonzero elements in a k.sup.th row of the matrix A, and y.sub.k is a k.sup.th sample of the initial data set.

3. The method of claim 1 further comprising calculating terms used to obtain the estimated image value, wherein the calculating reduces processing time and required memory by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse.

4. The method of claim 3 wherein the calculating the terms comprises: computing G.sub.l.sup.(m) by applying a hash-table-based computation to $q_{\mathrm{ij}}\left[k\right]=\underset{l{\mathrm{��}}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}$ where d.sub.ijl.sup.(m)=.sub.ijlx.sub.l.sup.(m) and S.sub.k={l=1, 2, . . . , L: n.sub.ijl=k} and .sub.ijl represents attenuation of the given radar pulse during travel from an i.sup.th transmit location to an l.sup.th voxel and back to a j.sup.th receiver.

5. The method of claim 4 further comprising computing G.sub.l.sup.(m) via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $G_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}{G}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $G_{\mathrm{ijl}}^{\left(m\right)}=\left\{\left(w*\left(y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}\right)\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}{s}_{\mathrm{ij}}^{\left(m\right)}\left[n\right]=\left(q_{\mathrm{ij}}*w\right)\left[n\right]$ and where y.sub.ij[k] is a k.sup.th sample of a radar-return profile associated with an i.sup.th transmit location and j.sup.th receiver, s.sub.ij.sup.(m)[k] is the m.sup.th estimate of a noise-free component of y.sub.ij[k], w is a discretized version of the given radar pulse, and n.sub.ijl is a discrete time-delay corresponding to rounding the quotient of the time taken by the given radar pulse to travel from a transmitter at the i.sup.th transmit location to the l.sup.th voxel and back to the j.sup.th receiver and a sampling interval, and * denotes discrete-time convolution.

6. The method of claim 3 wherein the calculating of the terms comprises: computing H.sub.l by applying a hash-table-based computation to
|S.sub.k| where |S.sub.k| denotes a number of elements in the set S.sub.k.

7. The method of claim 6 wherein the computing H.sub.l further comprises computing H.sub.l via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $H_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}^2{H}_{\mathrm{ijl}}$ $\mathrm{where}$ $H_{\mathrm{ijl}}=\left({}_{\mathrm{ij}}*h\right)\left[k\right]{.Math.}_{k=n_{\mathrm{ijl}}},{}_{\mathrm{ij}}\left[n\right]=v\left(\min \left(n+M,N\right)-v\left(\max \left(0,n-M\right)\right)\right),v\left(m\right)=\underset{k=0}{\overset{m}{.Math.}}.Math.S_k.Math.$ and where h is a discretized version of a squared radar pulse and 2M+1 is a number of non-zero samples in the given radar pulse.

8. The method of claim 1 further comprising: emitting a radar pulse applied into the ground at specified intervals into a scene of interest; detecting magnitude of signal reflections from the scene of interest from the radar pulse; recording position data corresponding to individual radar pulse emissions and individual receptions of the signal reflections; and creating the initial data set from the position data and detected magnitudes of the signal reflections.

9. A method of creating an image from an initial data set, the method comprising: receiving the initial data set, the receiving comprising receiving data representing transmission site locations of radar pulses applied into ground in a ground penetrating radar application, reception site locations of reception of reflections from the radar pulses, radar-return profiles for pairings of the transmission site locations and the reception site locations, and data samples associated with individual radar-return profiles; processing the initial data set with a processing device by: creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve the l.sub.1-regularized least-absolute deviation estimation problem associated with a mathematical model of image data from the initial data set; displaying the image using the estimated image value of individual voxels of the image.

10. The method of claim 9 further comprising iteratively deriving an estimated image value for a given voxel using $\mathrm{for}l=1,2,\mathrm{.Math.},L$ $x_l^{\left(m+1\right)}=\frac{D_l^{\left(m\right)}{x}_l^{\left(m\right)}+N_l^{\left(m\right)}}{D_l^{\left(m\right)}+\frac{}{.Math.x_l^{\left(m\right)}.Math.}}\mathrm{for}l=1,2,\mathrm{.Math.},L$ $\mathrm{where}$ $D_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}\frac{r_k{A}_{\mathrm{kl}}^2}{.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}$ $N_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}{A}_{\mathrm{kl}}.Math.\mathrm{sign}\left(y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)$ and where x.sub.l.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at an l.sup.th voxel for iterate number (m), A is a KL system matrix associated with the mathematical model of the image data, r.sub.k is a number of nonzero elements in a k.sup.th row of the matrix A, and y.sub.k is a k.sup.th sample of the initial data set.

11. The method of claim 9 further comprising calculating terms used to obtain the estimated image value, wherein the calculating the terms reduces processing time and required memory by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse.

12. The method of claim 11 wherein the calculating the terms comprises: computing N.sub.l.sup.(m) by applying a hash-table-based computation to $q_{\mathrm{ij}}\left[k\right]=\underset{l{\mathrm{��}}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}$

13. The method of claim 12 wherein the computing the term N.sub.l.sup.(m) further comprises computing N.sub.l.sup.(m) via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L; $N_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}{N}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $N_{\mathrm{ijl}}^{\left(m\right)}=\left\{w*\left(\mathrm{sign}\left(y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}\right)\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}},s_{\mathrm{ij}}^{\left(m\right)}\left[n\right]=\left(q_{\mathrm{ij}}*w\right)\left[n\right],$ and y.sub.ij is a k.sup.th sample of a radar-return profile associated with an i.sup.th transmit location and a j.sup.th receiver, s.sub.ij.sup.(m)[k] is an m.sup.th estimate of a noise-free component of y.sub.ij[k], w is a discretized version of the given radar pulse, .sub.ijl represents attenuation of the given radar pulse during travel from an i.sup.th transmit location to an l.sup.th voxel and back to a j.sup.th receiver, and n.sub.ijl is a discrete time-delay corresponding to rounding a quotient of time for the given radar pulse to travel from a transmitter at the i.sup.th transmit location to the l.sup.th voxel and back to the j.sup.th receiver and a sampling interval, and * denotes discrete-time convolution.

14. The method of claim 11 wherein the calculating the terms comprises: computing D.sub.l.sup.(m) by applying a hash-table-based computation to |S.sub.k| where |S.sub.k| denotes a number of elements in the set S.sub.k.

15. The method of claim 14 wherein the computing the term D.sub.l.sup.(m) further comprises computing D.sub.l.sup.(m) via
for i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $D_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}^2{D}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $D_{\mathrm{ijl}}^{\left(m\right)}=\left\{\left(\left(\frac{{}_{\mathrm{ij}}}{.Math.y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}.Math.}\right)*h\right)\right\}\left[k\right]{.Math.}_{k=n_{\mathrm{ijl}}},{}_{\mathrm{ij}}\left[n\right]=v\left(\min \left(n+M,N\right)-v\left(\max \left(0,n-M\right)\right)\right),v\left(m\right)=\underset{k=0}{\overset{m}{.Math.}}.Math.S_k.Math.$ and where h is a discretized version of a squared radar pulse and 2M+1 is a number of non-zero samples in the given radar pulse.

16. The method of claim 9 further comprising: emitting a radar pulse applied into ground in a ground penetrating radar application at specified intervals into a scene of interest; detecting magnitude of signal reflections from the scene of interest from the radar pulse; and recording position data corresponding to individual radar pulse emissions and individual receptions of the signal reflections; creating the initial data set from the position data and detected magnitudes of the signal reflections.

17. An apparatus for detecting objects in a scene of interest, the apparatus comprising: a vehicle; a plurality of radar transmission devices mounted on the vehicle and configured to transmit radar pulses into ground in a scene of interest; a plurality of radar reception devices mounted on the vehicle configured to detect magnitude of signal reflections from the scene of interest from the radar pulses; a location determination device configured to detect location of the vehicle at times of transmission of the radar pulse from the plurality of radar transmission devices and reception of the signal reflections by the radar reception devices; a processing device configured to process an initial data set representing transmission site locations of individual ones of the radar pulses, reception site locations of reception of individual ones of the signal reflections, and number of data samples per reception profile by: creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem associated with a mathematical model of image data from the initial data set.

18. The apparatus of claim 17 wherein the processing device is further configured to iteratively derive an estimated image value for a given voxel using
for l=1,2, . . . ,L $x_l^{\left(m+1\right)}=\frac{H_l{x}_l^{\left(m\right)}+G_l^{\left(m\right)}}{H_l+\frac{}{.Math.2x_l^{\left(m\right)}.Math.}}$ $\mathrm{where}$ $H_l=\underset{k=0}{\overset{K-1}{.Math.}}{r}_k{A}_{\mathrm{kl}}^2$ $G_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}{A}_{\mathrm{kl}}\left(y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)$ where x.sub.l.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at the l.sup.th voxel for iterate number (m), A is a KL system matrix associated with the mathematical model of the image data, K=IJL where I is a number of transmission site locations, J is a number of reception sites, and L is a number of pixels in the image, r.sub.k is a number of nonzero elements in a k.sup.th row of the matrix A, and y.sub.k is a k.sup.th sample of the initial data set.

19. The apparatus of claim 17 wherein the processing device is further configured to calculate terms used to obtain the estimated image value, wherein the calculating reduces processing time and required memory by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse.

20. The apparatus of claim 19 wherein the processing device is further configured to calculate terms by computing G.sub.l.sup.(m) by applying a hash-table-based computation to $q_{\mathrm{ij}}\left[k\right]=\underset{l{\mathrm{��}}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}$ where d.sub.ijl.sup.(m)=.sub.ijlx.sub.l.sup.(m) and S.sub.k={l=1, 2, . . . , L: n.sub.ijl=k} and .sub.ijl represents attenuation of the given radar pulse during travel from an i.sup.th transmit location to an l.sup.th voxel and back to a j.sup.th receiver.

21. The apparatus of claim 20 wherein the processing device is further configured to compute G.sub.l.sup.(m) via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $G_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}{G}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $G_{\mathrm{ijl}}^{\left(m\right)}=\left\{\left(w*\left(y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}\right)\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}},s_{\mathrm{ij}}^{\left(m\right)}\left[n\right]=\left(q_{\mathrm{ij}}*w\right)\left[n\right],$ and where y.sub.ij[k] is a k.sup.th sample of a radar-return profile associated with an i.sup.th transmit location and j.sup.th receiver, s.sub.ij.sup.(m)[k] is the m.sup.th estimate of a noise-free component of y.sub.ij[k], w is a discretized version of the given radar pulse, and n.sub.ijl is a discrete time-delay corresponding to rounding a quotient of time for the given radar pulse to travel from a transmitter at the i.sup.th transmit location to the l.sup.th voxel and back to the j.sup.th receiver and a sampling interval, and * denotes discrete-time convolution.

22. The apparatus of claim 19 wherein the processing device is further configured to calculate terms by computing H.sub.l by applying a hash-table-based computation to
|S.sub.k| where |S.sub.k| denotes a number of elements in the set S.sub.k.

23. The apparatus of claim 22 wherein the processing device is further configured to compute H.sub.l by computing H.sub.l via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $H_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}^2{H}_{\mathrm{ijl}}$ $\mathrm{where}$ $H_{\mathrm{ijl}}=\left\{\left({}_{\mathrm{ij}}*h\right)\right\}\left[k\right]{.Math.}_{k=n_{\mathrm{ijl}}},{}_{\mathrm{ij}}\left[n\right]=v\left(\min \left(n+M,N\right)-v\left(\max \left(0,n-M\right)\right)\right),v\left(m\right)=\underset{k=0}{\overset{m}{.Math.}}.Math.S_k.Math.$ where h is a discretized version of a squared radar pulse and 2M+1 is a number of non-zero samples in the given radar pulse.

24. An apparatus for detecting objects in a scene of interest, the apparatus comprising: a vehicle; a plurality of radar transmission devices mounted on the vehicle and configured to transmit radar pulses into ground in a scene of interest; a plurality of radar reception devices mounted on the vehicle configured to detect magnitude of signal reflections from the scene of interest from the radar pulses; a location determination device configured to detect location of the vehicle at times of transmission of the radar pulse from the plurality of radar transmission devices and reception of the signal reflections by the radar reception devices; a processing device configured to process an initial data set representing transmission site locations of individual ones of the radar pulses, reception site locations of reception of individual ones of the signal reflections, and number of data samples per reception profile by: creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve the l.sub.1-regularized least-absolute deviation estimation problem associated with a mathematical model of image data from the initial data set.

25. The apparatus of claim 24 wherein the processing device is further configured to iteratively derive an estimated image value for a given voxel using
for l=1,2, . . . ,L $x_l^{\left(m+1\right)}=\frac{D_l^{\left(m\right)}{x}_l^{\left(m\right)}+N_l^{\left(m\right)}}{D_l^{\left(m\right)}+\frac{}{.Math.x_l^{\left(m\right)}.Math.}}\mathrm{for}l=1,2,\mathrm{.Math.},L$ $\mathrm{where}$ $D_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}\frac{r_k{A}_{\mathrm{kl}}^2}{.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}$ $N_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}{A}_{\mathrm{kl}}.Math.\mathrm{sign}\left(y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)$ where x.sub.l.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at an l.sup.th voxel for iterate number (m), A is a KL system matrix associated with the mathematical model of the image data, r.sub.k is a number of nonzero elements in a k.sup.th row of the matrix A, and y.sub.k is a k.sup.th sample of the initial data set.

26. The apparatus of claim 24 wherein the processing device is further configured to calculate terms used to obtain the estimated image value, wherein the calculating the terms reduces processing time and required memory by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse.

27. The apparatus of claim 26 wherein the processing device is further configured to calculate terms by computing N.sub.l.sup.(m) by applying a hash-table-based computation to $q_{\mathrm{ij}}\left[k\right]=\underset{l{\mathrm{��}}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}.$

28. The apparatus of claim 27 wherein the processing device is further configured to compute N.sub.l.sup.(m) via
i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L; $N_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}{N}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $N_{\mathrm{ijl}}^{\left(m\right)}=\left\{w*\left(\mathrm{sign}\left(y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}\right)\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}},s_{\mathrm{ij}}^{\left(m\right)}\left[n\right]=\left(q_{\mathrm{ij}}*w\right)\left[n\right],$ and y.sub.ij is a k.sup.th sample of a radar-return profile associated with an i.sup.th transmit location and a j.sup.th receiver, s.sub.ij.sup.(m)[k] is an m.sup.th estimate of a noise-free component of y.sub.ij[k], w is a discretized version of the given radar pulse, .sub.ijl represents attenuation of the given radar pulse during travel from an i.sup.th transmit location to an l.sup.th voxel and back to a j.sup.th receiver, and n.sub.ijl is a discrete time-delay corresponding to rounding a quotient of time for the given radar pulse to travel from a transmitter at the i.sup.th transmit location to the l.sup.th voxel and back to the j.sup.th receiver and a sampling interval, and * denotes discrete-time convolution.

29. The apparatus of claim 26 wherein the processing device is further configured to calculate terms by computing D.sub.l.sup.(m) by applying a hash-table-based computation to |S.sub.k| where |S.sub.k| denotes a number of elements in the set S.sub.k.

30. The apparatus of claim 29 wherein the processing device is further configured to calculate terms by computing D.sub.l.sup.(m) via
for i=1,2, . . . ,I; j=1,2, . . . ,J; l=1,2, . . . ,L $D_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}^2{D}_{\mathrm{ijl}}^{\left(m\right)}$ $\mathrm{where}$ $D_{\mathrm{ijl}}^{\left(m\right)}=\left\{\left(\left(\frac{{}_{\mathrm{ij}}}{.Math.y_{\mathrm{ij}}-s_{\mathrm{ij}}^{\left(m\right)}.Math.}\right)*h\right)\right\}\left[k\right]{.Math.}_{k=n_{\mathrm{ijl}}},{}_{\mathrm{ij}}\left[n\right]=v\left(\min \left(n+M,N\right)-v\left(\max \left(0,n-M\right)\right)\right),v\left(m\right)=\underset{k=0}{\overset{m}{.Math.}}.Math.S_k.Math.$ where h is a discretized version of a squared radar pulse and 2M+1 is a number of non-zero samples in the given radar pulse.

31. A method of creating an image from a data set, the method comprising: receiving a DAS image data set created by applying a delay-and-sum (DAS) algorithm to create an initial data set, the receiving comprising receiving data representing transmission site locations of radar pulses applied into ground in a ground penetrating radar application, reception site locations of reception of reflections from the radar pulses, radar-return profiles for pairings of the transmission site locations and the reception site locations, and data samples associated with individual radar-return profiles; processing the DAS image data set with a processing device by creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem that selects a sparse image derived from the DAS image data set; displaying the image using the estimated image value of individual voxels of the image.

32. The method of claim 31 further comprising iteratively deriving an estimated image value for a given voxel using
for l=1,2, . . . ,L $x_l^{\left(m+1\right)}=\frac{x_{\mathrm{DAS},l}}{1+\frac{}{2.Math.x_l^{\left(m\right)}.Math.}}$ and where x.sub.l.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at an l.sup.th voxel for iterate number (m), is a penalty parameter, and X.sub.DAS,l is a value of the l.sup.th voxel of the DAS image data set.

33. The method of claim 31 further comprising: emitting a radar pulse applied into the ground at specified intervals into a scene of interest; detecting magnitude of signal reflections from the scene of interest from the radar pulse; recording position data corresponding to individual radar pulse emissions and individual receptions of the signal reflections; creating the initial data set from the position data and detected magnitudes of the signal reflections; and applying the delay-and-sum (DAS) algorithm to the initial data set to create the DAS image data set.

34. An apparatus for detecting objects in a scene of interest, the apparatus comprising: a vehicle; a plurality of radar transmission devices mounted on the vehicle and configured to transmit radar pulses into ground in a scene of interest; a plurality of radar reception devices mounted on the vehicle configured to detect magnitude of signal reflections from the scene of interest from the radar pulses; a location determination device configured to detect location of the vehicle at times of transmission of the radar pulse from the plurality of radar transmission devices and reception of the signal reflections by the radar reception devices; a processing device configured to process an initial data set representing transmission site locations of individual ones of the radar pulses, reception site locations of reception of individual ones of the signal reflections, and number of data samples per reception profile by: applying a delay-and-sum (DAS) algorithm to the initial data set to create a DAS image data set, creating an estimated image value for each voxel in the image to create an estimated image value data set having most of the estimated image value data set's values at or near zero by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem that selects a sparse image derived from the DAS image data set.

35. The apparatus of claim 34 wherein the processing device is further configured to iteratively derive an estimated image value for a given voxel using
for l=1,2, . . . ,L $x_l^{\left(m+1\right)}=\frac{x_{\mathrm{DAS},l}}{1+\frac{}{2.Math.x_l^{\left(m\right)}.Math.}}$ and where x.sub.l.sup.(m+1) is the estimated image value at an l.sup.th voxel for iterate number (m+1), x.sub.l.sup.(m) is the estimated image value at an l.sup.th voxel for iterate number (m), is a penalty parameter, and X.sub.DAS,l is a value of the l.sup.th voxel of the DAS image data set.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The above needs are at least partially met through provision of the methods and apparatuses for receiving and processing image data as described in the following detailed description, particularly when studied in conjunction with the drawings, wherein:

(2) FIG. 1 comprises a schematic of operation a prior art GPR system;

(3) FIG. 2 comprises a schematic of an example radar system as configured in accordance with various embodiments of the invention;

(4) FIG. 3 comprises a perspective view of an example implementation of a prior art ultra wide-band (UWB) synchronous reconstruction (SIRE) radar system;

(5) FIG. 4 comprises a schematic demonstrating a prior art approach to obtaining initial data from a SOI;

(6) FIG. 5 comprises a graph demonstrating application of the mathematical majorize-minimize principle;

(7) FIG. 6 comprises a flow diagram of an example algorithm applying the M-M principle to an L1-LS estimation as configured in accordance with various embodiments of the invention;

(8) FIG. 7 comprises a graph displaying accuracies of various algorithms as applied to processing a given data set;

(9) FIG. 8 comprises a displayed image of objects in a SOI using image data processed according to an M-M application to an L1-least squares method as configured in accordance with various embodiments of the invention;

(10) FIG. 9 comprises a displayed image of the objects in the SOI of FIG. 8, in this case using image data processed according to a prior art DAS method;

(11) FIG. 10 comprises a displayed image of the objects in the SOI of FIG. 8, in this case using image data processed according to a prior art RSM method;

(12) FIG. 11 comprises a flow diagram of an example algorithm applying the M-M principle to an L1-LAD estimation as configured in accordance with various embodiments of the invention;

(13) FIG. 12 comprises a graph displaying accuracies of various algorithms as applied to processing a given data set having an outlier;

(14) FIG. 13 comprises a graph displaying accuracy of an example algorithm applying the M-M principle to an L1-LAD estimation as applied to processing the given data set having an outlier of FIG. 12;

(15) FIG. 14 comprises a graph displaying accuracies of various algorithms as applied to processing a given data set without outliers;

(16) FIG. 15 comprises a graph displaying a cost function for an example algorithm applying the M-M principle to an L1-LAD estimation as configured in accordance with various embodiments of the invention;

(17) FIG. 16 comprises a displayed image of objects in a SOI using image data processed according to a prior art DAS algorithm;

(18) FIG. 17 comprises a displayed image of the objects in the SOI of FIG. 16, in this case using image data processed according to an M-M application to an L1-LAD method as configured in accordance with various embodiments of the invention;

(19) FIG. 18 comprises a flow diagram of an example algorithm applying the M-M principle to an L1-least squares estimation applied to an image data set from a DAS algorithm (L1-SIR) as configured in accordance with various embodiments of the invention;

(20) FIG. 19 comprises a displayed image of objects in a SOI using image data processed according to a prior art DAS algorithm;

(21) FIG. 20 comprises a displayed image of the objects in the SOI of FIG. 19, in this case using image data processed according to the L1-SIR algorithm as configured in accordance with various embodiments of the invention;

(22) FIG. 21 comprises a displayed image of the objects in the SOI of FIG. 19, in this case using image data processed according to example algorithm applying the M-M principle to an L1-LS estimation as configured in accordance with various embodiments of the invention.

(23) Skilled artisans will appreciate that elements in the figures are illustrated for simplicity and clarity and have not necessarily been drawn to scale. For example, the dimensions and/or relative positioning of some of the elements in the figures may be exaggerated relative to other elements to help to improve understanding of various embodiments of the present invention. Also, common but well-understood elements that are useful or necessary in a commercially feasible embodiment are often not depicted in order to facilitate a less obstructed view of these various embodiments. It will further be appreciated that certain actions and/or steps may be described or depicted in a particular order of occurrence while those skilled in the art will understand that such specificity with respect to sequence is not actually required. It will also be understood that the terms and expressions used herein have the ordinary technical meaning as is accorded to such terms and expressions by persons skilled in the technical field as set forth above except where different specific meanings have otherwise been set forth herein.

DETAILED DESCRIPTION

(24) Referring now to the drawings, and in particular to FIG. 2, an illustrative apparatus that is compatible with many of these teachings will now be presented. In a GPR application, a vehicle 202 includes a plurality of radar transmission devices 204 and 206 mounted on the vehicle 202. The radar transmission devices 204 and 206 are configured to transmit radar pulses 208 and 210 into a scene of interest 212. A plurality of J radar reception devices 214 are mounted on the vehicle 202 and configured to detect magnitude of signal reflections from the scene of interest 212 from the radar pulses 208 and 210. The vehicle 202 can be any structure able to carry the radar transmitters and receivers to investigate a scene of interest.

(25) FIG. 3 illustrates one example implementation: a truck 302 mounted ultra wide-band (UWB) synchronous reconstruction (SIRE) radar system developed by the US Army Research Laboratory (ARL) in Adelphi, Md. This system includes a left transmit antenna 304, a right transmit antenna 306, and 16 receive antennas 314. Other systems may have different numbers of and arrangement of transmit and receive antennas. The transmit and receive antennas are mounted to a support structure 321, which supports these elements on the truck 302. The truck 302 drives through a scene of interest while the transmit antennas 304 and 306 alternately transmit radar pulses and the receiving antennas 314 receive reflections of the transmitted radar pulses from backscatter objects in the scene of interest.

(26) Referring again to the example of FIG. 2, the positions along the vehicle's 202 path at which a radar pulse is transmitted are referred to as the transmit locations, I. As the vehicle 202 moves, the two transmit antennas 204 and 206 alternately send respective probing pulses 208 and 210 toward the SOI 212, and the radar-return profiles reflected from the SOI 212 are captured by multiple receive antennas 214 at each transmit location.

(27) A location determination device 220 detects the location of the vehicle 202 at times of transmission of the radar pulse from the plurality of radar transmission devices 204 and 206 and reception of the signal reflections by the radar reception devices 214. In one example, the location determination device is a global positioning system (GPS) device as commonly known and used, although other position determination devices can be used. Accordingly, the positioning coordinates of the active transmit antenna 204 or 206 and all the receive antennas 214 are also logged. When using the UWB SIRE system of FIG. 3, there is typically a minimum range of detection of objects in the SOI from the vehicle 2020 of about 8 meters, a maximum range of about 34 meters, and a cross-range of about 25 meters. The UWB SIRE system uses a FORD EXPLORER as the vehicle 202 such that the transmit antennas 204 and 206 have about a two meter separation between 16 receivers.

(28) In one approach, the vehicle 202 includes a processing device 242 in communication with the location determining device 220, the transmit antennas 204 and 206, and the receivers 214 to coordinate their various operations and to store information related to their operations in a memory device 244. Optionally, a display 246 is included with the vehicle 202 to display an image related to the data received from the scanning of the scene of interest 212.

(29) Due to the large size of the scene-of-interest, an initial data set is not generated by processing all voxels at once. Such an image would have cross-range resolution that varies from the near-range to the far-range voxels. The voxels in the near-range would have much larger resolution than those for the far-range ones. To create GPR images with consistent resolution across the scene-of-interest, we use the mosaicking approach discussed in L. Nguyen, Signal and Image Processing Algorithms for the U.S. Army Research Laboratory Ultra-wideband (UWB) Synchronous Impulse Reconstruction (SIRE) Radar, ARL Technical Report, ARL-TR-4784, April 2009, which is incorporated by reference and described with reference to FIG. 4. The steps taken to produce a complete image of the scene-of-interest (using the mosaicking approach) are described as follows. The image space associated with the scene-of-interest is divided into 32 sub-images of size 252 m.sup.2. Each sub-image has 250 voxels in the cross-range direction and 100 voxels in the down-range direction. Thus, each sub-image has L=25000 voxels. The aperture (meaning the distance over which the vehicle travels while accumulating data for a SOI) is divided into 32 sub-apertures corresponding to separate, overlapping distances traveled by the vehicle, where adjacent sub-apertures (or vehicle travel windows) overlap by approximately 2 meters. Each sub-aperture has 43 transmit locations and is approximately of size 122 m.sup.2. The radar-return and location positioning measurements associated with a sub-aperture are used to estimate the reflectance coefficients for the corresponding sub-image. The reconstructed sub-images are merged together to obtain the complete image of the scene-of-interest.

(30) In another approach, referring again to FIG. 2, a separate computing device 260 may receive an initial data set from the vehicle based system to further process to create and optionally display images related to the SOI. The computing device 260 will typically include a processing device 262 in communication with a memory 264 to allow for processing the data according to any of the methods described herein. A display 266 may be included with or separate from the computing device 260 and controlled to display images resulting from the processing of the data received from the vehicle 202. Those skilled in the art will recognize and appreciate that such processing devices 242 and 262 can comprise a fixed-purpose hard-wired platform, including, for example, parallel processing devices, or can comprise a partially or wholly programmable platform. All of these architectural options are well known and understood in the art and require no further description here. Moreover, the memory devices 244 and 264 may be separate from or combined with the respective processing devices 242 and 262. Any memory device or arrangement capable of facilitating the processing described herein may be used.

(31) With respect to the collection of data, consider a single scatterer, with spatial position p.sub.s, located at the center of a voxel (i.e., volume element) within the SOI. The spatial positions of the active transmit antenna and a receive antenna are denoted by p.sub.t and p.sub.r, respectively. If the contributions of measurement noise are momentarily ignored, the relationship between the transmitted signal p(t) and the received signal g.sub.s(t) can be modeled as
g.sub.s(t)=.sub.s.Math.p(t(p.sub.s,p.sub.t,p.sub.r)).Math.x.sub.s,(2)
where x.sub.s is the reflection coefficient of the voxel, (p.sub.s,p.sub.t,p.sub.r) is the time it takes for the pulse to travel from the transmit antenna to the scatterer and back to the receive antenna, and .sub.s is the attenuation the pulse undergoes along the round-trip path.

(32) The single scatterer model in (2) can be generalized to describe all the measurements captured by the SIRE GPR system. The SOI is subdivided into a rectangular grid of L voxels and the unknown reflection coefficient at the l.sup.th voxel is denoted by x.sub.l. Extending the model in (2) to the SIRE system, the output of the j.sup.th receive antenna at the i.sup.th vehicle-stop is given by

(33) $\begin{array}{cc}{s}_{\mathrm{ij}}\left(t\right)=\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.p\left(t-{}_{\mathrm{ijl}}\right).Math.x_l+w_{\mathrm{ij}}\left(t\right),i=1,2,\mathrm{.Math.},I,j=1,2,\mathrm{.Math.},J.& \left(3\right)\end{array}$

(34) In this equation, .sub.ijl is the time it takes for the transmitted pulse to propagate from the active transmit antenna at the i.sup.th transmit location to the l.sup.th voxel and for the backscattered signal to return to the j.sup.th receive antenna. The parameter .sub.ijl is given by

(35) $\begin{array}{cc}{}_{\mathrm{ijl}}=\frac{d_{\mathrm{il}}+d_{\mathrm{ijl}}}{c},& \left(4\right)\end{array}$
where d.sub.il denotes the distance from the active transmit antenna at the i.sup.th transmit location to the l.sup.th voxel, d.sub.ijl denotes the return distance from the l.sup.th voxel to the j.sup.th receive antenna when the truck is at the i.sup.th transmit location, and c is the speed of light.

(36) The notation .sub.ijl is the propagation loss that the transmitted pulse undergoes as it travels from the active transmit antenna at the i.sup.th transmit location to the l.sup.th voxel and back to the j.sup.th receive antenna. The parameter .sub.ijl is given by

(37) $\begin{array}{cc}{}_{\mathrm{ijl}}=\frac{1}{d_{\mathrm{il}}.Math.d_{\mathrm{ijl}}}.& \left(5\right)\end{array}$
The notation w.sub.ij(l) represents the noise contribution.

(38) The above mathematical model defined in (3) is continuous whereas, in practice, the SIRE GPR system only stores discrete and separate sampled versions of the return signals. Thus, we introduce the following discrete-time signals to adpat the above model to the real world application: for i=1, 2, . . . , I, j=1, 2, . . . , J and n=0, 1, . . . , N1,
y.sub.ij[n]s.sub.ij(nT.sub.s)(6)
e.sub.ij[n]w.sub.ij(nT.sub.s)(7)
where T.sub.s is the sampling interval and N is the number of samples per radar return. From (6) and (7) we can write

(39) $\begin{array}{cc}{y}_{\mathrm{ij}}\left[n\right]=\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.p\left({\mathrm{nT}}_s-{}_{\mathrm{ijl}}\right).Math.x_l+e_{\mathrm{ij}}\left[n\right].& \left(8\right)\end{array}$
The corresponding system of N equations can be written in matrix form as
y.sub.ij=D.sub.ijP.sub.ijx+e.sub.ij(9)
where

(40) $\begin{array}{cc}{y}_{\mathrm{ij}}\stackrel{}{=}\left[\begin{array}{c}{y}_{\mathrm{ij}}\left[0\right]\\ y_{\mathrm{ij}}\left[1\right]\\ \mathrm{.Math.}\\ y_{\mathrm{ij}}\left[N-1\right]\end{array}\right],x\stackrel{}{=}\left[\begin{array}{c}{x}_1\\ x_2\\ \mathrm{.Math.}\\ x_L\end{array}\right],e_{\mathrm{ij}}\stackrel{}{=}\left[\begin{array}{c}{e}_{\mathrm{ij}}\left[0\right]\\ e_{\mathrm{ij}}\left[1\right]\\ \mathrm{.Math.}\\ e_{\mathrm{ij}}\left[N-1\right]\end{array}\right].& \left(10\right)\end{array}$
The matrix D.sub.ij is an LL diagonal matrix containing the attenuation coefficients that is given by

(41) $\begin{array}{cc}{D}_{\mathrm{ij}}\left[\begin{array}{cccc}{}_{\mathrm{ij}1}& 0& \mathrm{.Math.}& 0\\ 0& {}_{\mathrm{ij}2}& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& 0& & \mathrm{.Math.}\\ 0& \mathrm{.Math.}& 0& {}_{\mathrm{ijL}}\end{array}\right].& \left(11\right)\end{array}$
The matrix P.sub.ij is an NL matrix containing shifted versions of the transmitted pulse that is defined to be

(42) $\begin{array}{cc}\left(12\right)P_{\mathrm{ij}}\stackrel{}{=}[\begin{array}{cccc}p\left(0.Math.T_s-{}_{\mathrm{ij}1}\right)& p\left(0.Math.T_s-{}_{\mathrm{ij}2}\right)& \mathrm{.Math.}& p\left(0.Math.T_s-{}_{\mathrm{ijL}}\right)\\ p\left(1.Math.T_s-{}_{\mathrm{ij}1}\right)& p\left(1.Math.T_s-{}_{\mathrm{ij}2}\right)& \mathrm{.Math.}& p\left(1.Math.T_s-{}_{\mathrm{ijL}}\right)\\ \mathrm{.Math.}& \mathrm{.Math.}& & \mathrm{.Math.}\\ p\left(\left(N-1\right).Math.T_s-{}_{\mathrm{ij}1}\right)& p\left(\left(N-1\right).Math.T_s-{}_{\mathrm{ij}2}\right)& \mathrm{.Math.}& p\left(\left(N-1\right).Math.T_s-{}_{\mathrm{ijL}}\right)\end{array}].& \end{array}$

(43) In other words, the sampled data vectors (i.e., values for position and signal for transmission and reception of radar pulses) for all transmitters and receivers pairs {y.sub.ij} are concatenated to obtain a K1(K=IJN) data vector y. Extending the model in (9) to account for all I.Math.J transmitter and receiver pairs yields the desired model
y=Ax+e,(13)
where the K1 data vector y, KL system matrix A, and K1 Gaussian noise vector e are given by

(44) $\begin{array}{cc}y\stackrel{}{=}\left[\begin{array}{c}{y}_{11}\\ y_{12}\\ \mathrm{.Math.}\\ y_{1J}\\ y_{21}\\ \mathrm{.Math.}\\ y_{2J}\\ \mathrm{.Math.}\\ y_{I1}\\ \mathrm{.Math.}\\ y_{\mathrm{IJ}}\end{array}\right],A\stackrel{}{=}\left[\begin{array}{c}{P}_{11}{D}_{11}\\ P_{12}{D}_{12}\\ \mathrm{.Math.}\\ P_{1J}{D}_{1J}\\ P_{21}{D}_{21}\\ \mathrm{.Math.}\\ P_{2J}{D}_{2J}\\ \mathrm{.Math.}\\ P_{I1}{D}_{I1}\\ \mathrm{.Math.}\\ P_{\mathrm{IJ}}{D}_{\mathrm{IJ}}\end{array}\right],e\stackrel{}{=}\left[\begin{array}{c}{e}_{11}\\ e_{12}\\ \mathrm{.Math.}\\ e_{1J}\\ e_{21}\\ \mathrm{.Math.}\\ e_{2J}\\ \mathrm{.Math.}\\ e_{I1}\\ \mathrm{.Math.}\\ e_{\mathrm{IJ}}\end{array}\right].& \left(14\right)\end{array}$
Because the SIRE GPR system uses an UWB radar, the duration of the transmitted pulse p(t) is relatively short so that the system matrix A is sparse.

(45) Given the pulse p(t), location (e.g., GPS) data, and observation-vector y, the objective is to estimate the unknown reflection coefficient vector x, which represents the material reflecting radar pulses in the SOI Displaying this reflection coefficient data will correspond to displaying the objects in the SOI.

(46) The Majorize-Minimize Principle

(47) The MM (which stands for majorize-minimize in minimization problems, and minimize-majorize in maximization problems) principle is a prescription for constructing solutions to optimization problems. An MM algorithm minimizes an objective function by successively minimizing, at each iteration, a judiciously chosen objective function that is known as a majorizing function. Whenever a majorizing function is optimized, in principle, a step is taken toward reaching the minimizer of the original objective function. A brief summary of the MM principle is now given with reference to FIG. 5.

(48) Let be a function to be minimized over some domain D.sup.L, i.e., the function's minimum value is to be found within the given domain. A real value function g with domain DD is said to majorize if
g(x,y)(x) for all x,yD(15)
g(x,x)=(x) for all xD.(16)

(49) Suppose the majorizing function g is easier to minimize than the original objective function . Then, the MM algorithm for minimizing is given by

(50) 0$\begin{array}{cc}{x}^{\left(m+1\right)}=\underset{xD}{\arg \min }g\left(x,x^{\left(m\right)}\right),& \left(17\right)\end{array}$
where x.sup.(m) is the current estimate for the minimizer of . The algorithm defined by (17), which is illustrated in FIG. 5, is guaranteed to monotonically decrease the objective function with increasing iteration. In other words, a further minimal or smaller value for the function is obtained with each iteration of the algorithm, stepping between successive functions g. To see this result, first observe that, by definition,
g(x.sup.(m+1),x.sup.(m))g(x.sup.(m),x.sup.(m)).(18)
Now from (15) and (16), it follows that
(x.sup.(m+1)g(x.sup.(m+1),x.sup.(m))g(x.sup.(m),x.sup.(m))=(x.sup.(m)).(19)
In other words and as illustrated in FIG. 5, the function g(x,x.sup.(n)) intersects with the function (x) at x.sup.(n) and also has a further minimum at point x.sup.(n+1). That further minimum point is used in the next iteration as a new g(x.sup.(n)) from which a new minimum at a new x.sup.(n+1) may be determined.

(51) MM-Based Image Reconstruction Using L1-Regularization

(52) Using the above parameters, in a first approach to providing a fast and accurate image by exploiting the known sparsity of the scatterers, the object data represented by the reflection coefficient vector is estimated using the well-established l.sub.1-LS estimation method

(53) $\begin{array}{cc}\hat{x}=\underset{x}{\arg \min }{.Math.y-\mathrm{Ax}.Math.}_2^2+{.Math.x.Math.}_1,& \left(20\right)\end{array}$
where is the regularization parameter or penalty parameter. In contrast to previous approaches, the optimization problem in (20) is solved using the above described MM framework, which leads to an iterative algorithm that is efficient, straightforward to implement and amenable to parallelization. Additionally, the algorithm is guaranteed to monotonically decrease the objective function in (20) to guarantee coming to a final result through the iteration.

(54) Recall that the objective function to be minimized is of the form
(x)=.sub.1(x)+.sub.2(x)(21)
where .sub.1(x)yAx.sub.2.sup.2 and .sub.2(x)x.sub.1, and the regularization or penalty parameter is strictly positive. To find a majorizer for the function .sub.1, we use a result from DePierro (A. R. De Pierro, A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography, Medical Imaging, IEEE Transactions on, vol. 14, no. 1, pp. 132 to 137, 1995, which is incorporated by reference herein) outlined as follows. First, .sub.1(x) is expressed as

(55) $\begin{array}{cc}{}_1\left(x\right)=\underset{k=0}{\overset{k-1}{.Math.}}{y}_k^2-2\underset{k=0}{\overset{K-1}{.Math.}}{y_k\left[\mathrm{Ax}\right]}_k+\underset{k=0}{\overset{K-1}{.Math.}}{\left({\left[\mathrm{Ax}\right]}_k\right)}^2\mathrm{where}& \left(22\right)\\ {\left[\mathrm{Ax}\right]}_k\underset{l=1}{\overset{L}{.Math.}}{A}_{\mathrm{kl}}{x}_l& \left(23\right)\end{array}$
is the k.sup.th component of the vector Ax. They then exploit the convexity of the square function and construct a majorizing function for ([Ax].sub.k).sup.2. By denoting r.sub.k as the number of nonzero elements in the k.sup.th row of A and defining

(56) $\begin{array}{cc}{c}_{\mathrm{kl}}\{\begin{array}{cc}{r}_k^{-1},& A_{\mathrm{kl}}0\\ 0,& A_{\mathrm{kl}}=0,\end{array}\mathrm{they}\mathrm{have}& \left(24\right)\\ {\left({\left[\mathrm{Ax}\right]}_k\right)}^2={\left(\underset{l=1}{\overset{L}{.Math.}}{c}_{\mathrm{kl}}\left(r_k{A}_{\mathrm{kl}}{x}_l-r_k{A}_{\mathrm{kl}}{x}_l^{\left(m\right)}+{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)\right)}^2& \left(25\right)\end{array}$
for any vector x.sup.(m) in .sup.L. Because

(57) $\begin{array}{cc}\underset{l=1}{\overset{L}{.Math.}}{c}_{\mathrm{kl}}=1& \left(26\right)\end{array}$
for all k, it follows from the convexity of the square function that
([Ax].sub.k).sup.2q(x,x.sup.(m)),(27)
where

(58) $\begin{array}{cc}q\left(x,x^{\left(m\right)}\right)\underset{l=1}{\overset{L}{.Math.}}{c_{\mathrm{kl}}\left(r_k{A}_{\mathrm{kl}}{x}_l-r_k{A}_{\mathrm{kl}}{x}_l^{\left(m\right)}+{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)}^2.& \left(28\right)\end{array}$
From Equation (27) and the fact that q(x.sup.(m),x.sup.(m))=([Ax].sub.k).sup.2, it follows that q is a majorizing function for ([Ax].sub.k).sup.2. Thus replacing ([Ax].sub.k).sup.2 by q(x,x.sup.(m)) in (22) produces

(59) $\begin{array}{cc}{Q}_1\left(x,x^{\left(m\right)}\right)=\underset{k=0}{\overset{K-1}{.Math.}}{y}_k^2-\underset{k=0}{\overset{K-1}{.Math.}}{y_k\left[\mathrm{Ax}\right]}_k+\underset{k=0}{\overset{K-1}{.Math.}}q\left(x,x^{\left(m\right)}\right).& \left(29\right)\end{array}$
the desired majorizing function for .sub.1.

(60) A quadratic majorizer for the absolute value function |x| was derived by de Leeuw and Lange (J. de Leeuw and K. Lange, Sharp quadratic majorization in one dimension, Computational statistics and data analysis, vol. 53, no. 7, pp. 2471 to 2484, 2009, which is incorporated by reference herein) where

(61) $\begin{array}{cc}z\left(x,y\right)=\underset{l=1}{\overset{L}{.Math.}}\frac{x^2}{2.Math.y.Math.}+\frac{1}{2}.Math.y.Math..& \left(30\right)\end{array}$
It follows readily from this result that a majorizing function for the function .sub.2 is

(62) $\begin{array}{cc}{Q}_2\left(x,x^{\left(m\right)}\right)=\underset{l=1}{\overset{L}{.Math.}}z\left(x_l,x_l^{\left(m\right)}\right).& \left(31\right)\end{array}$
Because is positive, a majorizing function for the l.sub.1-LS objective function is
Q(x,x.sup.(m))=Q.sub.1(x,x.sup.(m))+Q.sub.2(x,x.sup.(m))(32)
From the general expression in provided by G. Davis, S. Mallat, and M. Avellaneda, Adaptive greedy approximations, Constructive approximation, vol. 13, no. 1, pp. 57 to 98, 199, which is incorporated by reference herein, it follows that the next iterate is given by

(63) $\begin{array}{cc}{x}^{\left(m+1\right)}=Q\left(x,x^{\left(m\right)}\right).& \left(33\right)\end{array}$

(64) We obtain the desired iterative algorithm by setting to zero the derivative of Q(x,x.sup.(m)) with respect to the components of x. Straightforward calculations show that the partial derivative of Q(x,x.sup.(m)) with respect to x.sub.1 is given by

(65) 0$\begin{array}{cc}\frac{Q\left(x,x^{\left(m\right)}\right)}{x_l}=-2\underset{k=0}{\overset{K-1}{.Math.}}{y}_k{A}_{\mathrm{kl}}+2\underset{k=0}{\overset{K-1}{.Math.}}\left(r_k{A}_{\mathrm{kl}}^2{x}_l-r_k{A}_{\mathrm{kl}}^2{x}_l^{\left(m\right)}+{A_{\mathrm{kl}}\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)+\frac{x_l}{.Math.x_l^{\left(m\right)}.Math.}.& \left(34\right)\end{array}$
Setting the derivatives in (34) to zero leads to the following MM algorithm for the l.sub.1-LS estimation problem, representing application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem associated with a mathematical model of image data from the initial data set:

(66) $\begin{array}{cc}{x}_l^{\left(m+1\right)}=\frac{H_l.Math.x_l^{\left(m\right)}+G_l^{\left(m\right)}}{H_l+\frac{}{2.Math.x_l^{\left(m\right)}.Math.}}\mathrm{for}l=1,2,\mathrm{.Math.},L\mathrm{where}& \left(35\right)\\ H_l=\underset{k=0}{\overset{K-1}{.Math.}}{r}_k{A}_{\mathrm{kl}}^2,& \left(36\right)\\ G_l^{\left(m\right)}=\underset{k=0}{\overset{K-1}{.Math.}}{A}_{\mathrm{kl}}\left(y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right).& \left(37\right)\end{array}$

(67) In other words, one may iteratively derive an estimated image value for a given voxel using the equation (35). Accordingly, a processing device may be configured to process an initial data set by creating an estimated image value for each voxel in the image by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem associated with a mathematical model of image data from the initial data set.

(68) One readily observes in (35) that the estimation of individual reflectance coefficients is decoupled because, for a given pixel, the computation of the next estimate x.sub.1.sup.(m+1) only depends on the current estimate x.sub.1.sup.(m). The proposed algorithm is thus amenable to parallel, distributed, and/or graphics processing unit (GPU) processing, which can further expedite processing of the image data over other processing approaches that cannot be implemented using such processing techniques. Moreover, this approach can be readily applied to a variety of applications where datasets are collected using synthetic aperture imaging measurement principles.

(69) A further benefit of this approach is the stability of the algorithm because of its convergence properties. To illustrate this benefit, certain theoretical results on the convergence of MM algorithms are described below to analyze the convergence properties of the above described MM-based l.sub.1-LS algorithm.

(70) First, we re-state here the so called Condition C2 and Theorem 3 of F. Vaida, Parameter convergence for EM and MM algorithms, Statistica Sinica, vol. 15, no. 3, p. 831, 2005, which is incorporated herein by reference, in the context of the MM algorithms, where they apply with minor modifications. These modifications include: (1) the regularity condition R4, which concern the missing data distribution, is not necessary; and (2) for the regularity condition R5 and the condition C2, the expected log-likelihood function of the augmented data is now replaced by the majorizing function.

(71) For the Condition C2 aspect as discussed in the Vaida reference, let be the set of stationary points defined as

(72) $\begin{array}{cc}{S}_{}=\left\{x^{*}\text{:}\frac{}{x}\left(x^{*}\right)=0\right\},& \left(38\right)\end{array}$
where is the l.sub.1-LS cost function. For all x, there exists a unique global minimizer of the majorizing function Q.

(73) For the Theorem 3 aspect as discussed in the Vaida reference, consider an MM iteration sequence {x.sup.(m)} that is defined by the starting point x.sup.(0) and iteration x.sup.(m+1)=(x.sup.(m)). If Condition C2 holds, then for any starting point x.sup.(m).fwdarw.x* as m.fwdarw., for some stationary point x* in . Moreover, x*=(x*) and, if x.sup.(m)x* for all m, the sequence of cost function values (x.sup.(m)) is strictly decreasing to (x*).

(74) Theorem 3 gives a simple condition to test the convergence of MM algorithms; that is, if the global minimum of the majorizing function Q(,x*) is unique for all x*, then the sequence the MM iterates {x.sup.(m): m=0, 1, . . . } will converge to a stationary point. We now show that the majorizing function Q for the l.sub.1-LS cost function is strictly convex and, thus has a unique global minimizer for all stationary points.

(75) The second-order partial derivatives of the function Q(,x*) are given by

(76) $\begin{array}{cc}\frac{{}^{2}Q\left(x,x^{*}\right)}{x_l^2}=\{\begin{array}{cc}2\underset{k=0}{\overset{K-1}{.Math.}}{r}_k{A}_{\mathrm{kl}}^2,\frac{}{\left|x_l^{*}\right|},& l=k\\ 0,& lk.\end{array}& \left(39\right)\end{array}$
Thus, the Hessian matrix of Q(,x*) is diagonal. By the principles disclosed by T. T. Wu and K. Lange, The MM alternative to EM, Statistical Science, vol. 25, no. 4, pp. 492 to 505, 2010, which is incorporated by reference herein, x*.sub.i must be non-zero for all l. Therefore, from (39) it follows that

(77) $\begin{array}{cc}0<\frac{{}^{2}Q\left(x,x^{*}\right)}{x_l^2}<,\mathrm{for}\mathrm{all}x,x^{*},& \left(40\right)\end{array}$
which implies that the Hessian matrix of Q(,x*) is strictly positive definite. Consequently, Q(,x*) is a strictly convex function and thus has a unique global minimum for all x*. Finally, by Theorem 3, the MM-based l.sub.1-LS algorithm is guaranteed to converge to a stationary point.

(78) Description of the Fast Implementation

(79) When applied in a typical GPR context, the computation of the term G.sub.l.sup.(m) in (35) requires the KL matrix A where K=IJN. For the above described UWB SIRE radar system, these parameters are I=43 transmit locations, J=16 receive antennas, N=1350 data samples per return-profile, and L=25000 pixels.

(80) These parameter settings require 173 gigabytes (GB) of memory to merely store the system matrix A. Because A has many zero-elements, however, the data could be more efficiently stored as a sparse matrix. Nevertheless, a sparse representation for A would still require approximately 16 GB of memory. With such a large memory size requirement, the construction of the A matrix in the current formulation of the algorithm is not feasible or practical for typical computing platforms, especially in field deployment where on site imaging would be advantageous. Indeed, virtually any other GPR image formation method that requires explicitly constructing the system matrix would have comparable requirements.

(81) In addition to memory size challenges, computational cost would also be an issue for the current format of the MM-based l.sub.1-LS algorithm. At each iteration, the computation of G.sub.l.sup.(m) would require the matrix multiplication Ax.sup.(m), which has (KL) time complexity. This operation is thus not practical for large-scale implementations where the parameters K and L are expected to be relatively large. For example, in our case, we have K=27520 and L=25000. Additional costs include the computation of the term H.sub.l where the number of non-zero elements in each of the K rows of A is needed. To arrive at a fast and memory-efficient implementation of the algorithm in (35), the following acceleration techniques may be implemented.

(82) Fast Implementation of G.sub.l.sup.(m)

(83) In a GPR context, the mathematical expressions at equations (36) and (37) above can be modified to reduce processing time and required memory by accounting for a symmetric nature of a given radar pulse, accounting for similar discrete time delays between transmission of a given radar pulse and reception of reflections from the given radar pulse, and accounting for a short duration of the given radar pulse. Accordingly, the equation for determining estimates for values x representing reflectance coefficients of the objects in the SOI, (35), involves calculation of the terms G.sub.l.sup.(m) and H.sub.l, which calculation can be streamlined according to the above assumptions. In application, a processing device is configured to calculate terms used to obtain the estimated value. Pursuant to these aspects, the expression for G.sub.l.sup.(m) in (37) can be written as

(84) $\begin{array}{cc}{G}_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}\underset{n=0}{\overset{N-1}{.Math.}}{}_{\mathrm{ijl}}.Math.p\left({\mathrm{nT}}_s-{}_{\mathrm{ijl}}\right).Math.g_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)& \left(41\right)\end{array}$
where, for n=0, 1, . . . , N1,

(85) $\begin{array}{cc}{g}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)=y_{\mathrm{ij}}\left[n\right]-\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.p\left({\mathrm{nT}}_s-{}_{\mathrm{ijl}}\right).& \left(42\right)\end{array}$

(86) To facilitate a fast implementation, we approximate the quantity G.sub.l.sup.(m) by

(87) $\begin{array}{cc}{\hat{G}}_l^{\left(m\right)}\stackrel{}{=}\underset{i=i}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}\underset{n=0}{\overset{N-1}{.Math.}}{}_{\mathrm{ijl}}.Math.p\left({\mathrm{nT}}_s-n_{\mathrm{ijl}}{T}_s\right).Math.{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right),\mathrm{where}& \left(43\right)\\ n_{\mathrm{ijl}}\stackrel{}{=}\mathrm{round}\left(\frac{{}_{\mathrm{ijl}}}{T_s}\right)& \left(44\right)\\ {\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)\stackrel{}{=}{y}_{\mathrm{ij}}\left[n\right]-{\hat{s}}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)& \left(45\right)\\ {\hat{s}}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)\stackrel{}{=}\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.p\left({\mathrm{nT}}_s-n_{\mathrm{ijl}}{T}_s\right).& \left(46\right)\end{array}$
We will refer to the set of values {n.sub.ijl} as the discrete-time delays. We can write .sub.l.sup.(m) as

(88) $\begin{array}{cc}{\hat{G}}_l^{\left(m\right)}=\underset{i=i}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}.Math.{\hat{G}}_{\mathrm{ijl}}^{\left(m\right)}\mathrm{where}& \left(47\right)\\ \begin{array}{c}{\hat{G}}_{\mathrm{ijl}}^{\left(m\right)}=\underset{n=0}{\overset{N-1}{.Math.}}p\left({\mathrm{nT}}_s-n_{\mathrm{ijl}}{T}_s\right).Math.{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left({\mathrm{nT}}_s\right)\\ =\underset{n=0}{\overset{N-1}{.Math.}}w\left[n-n_{\mathrm{ijl}}\right].Math.{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left[n\right]\left(49\right)\\ =\left\{\underset{n=0}{\overset{N-1}{.Math.}}w\left[n-k\right]{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left[n\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\left(50\right)\end{array}& \left(48\right)\end{array}$
with w[n]p(nT.sub.s). Because the transmitted pulse p(t) is symmetric, w[nk]=w[kn] holds for all n and k, and thus

(89) $\begin{array}{cc}\begin{array}{c}{\hat{G}}_{\mathrm{ijl}}^{\left(m\right)}=\left\{\underset{n=0}{\overset{N-1}{.Math.}}w\left[k-n\right]{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\left[n\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\\ =\left\{\left(w*{\hat{g}}_{\mathrm{ij}}^{\left(m\right)}\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\left(52\right)\\ =\left\{\left(w*\left(y_{\mathrm{ij}}-{\hat{s}}_{\mathrm{ij}}^{\left(m\right)}\right)\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}.\left(53\right)\end{array}& \left(51\right)\end{array}$

(90) It is readily observed that computing .sub.ijl.sup.(m) requires the convolution of the discrete pulse w[n] with the m.sup.th iteration of the error-term sequence (y.sub.ij[n].sub.ij.sup.(m)[n]). The sequences w[n] and y.sub.ij[n] are given. Hence, to efficiently compute .sub.ijl.sup.(m), a computationally efficient way for generating the sequence .sub.ij.sup.(m)[n] is needed.

(91) First, we note that the collection of discrete-time delays {n.sub.ijl} is expected to have repeated values. Let k.sub.min and k.sub.max denote respectively the minimum and maximum discrete-time delays. The sifting property of the unit impulse function can be used to write

(92) 0$\begin{array}{cc}\begin{array}{c}{\hat{s}}_{\mathrm{ij}}^{\left(m\right)}\left[n\right]=\underset{k=k_{\min }}{\overset{k_{\max }}{.Math.}}\left\{\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.w\left[n-n_{\mathrm{ijl}}\right].Math.\left[k-n_{\mathrm{ijl}}\right]\right\}\\ =\underset{k=k_{\min }}{\overset{k_{\max }}{.Math.}}\left\{\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.w\left[n-k\right].Math.\left[k-n_{\mathrm{ijl}}\right]\right\}\left(55\right)\\ =\underset{k=k_{\min }}{\overset{k_{\max }}{.Math.}}\left\{\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.\left[k-n_{\mathrm{ijl}}\right]\right\}w\left[n-k\right]\left(56\right)\\ =\underset{k=k_{\min }}{\overset{k_{\max }}{.Math.}}{q}_{\mathrm{ij}}\left[k\right].Math.w\left[n-k\right]\left(57\right)\\ =\left(q_{\mathrm{ij}}*w\right)\left[n\right]\left(58\right)\end{array}\mathrm{where}& \left(54\right)\\ q_{\mathrm{ij}}\left[k\right]\stackrel{}{=}\underset{l=1}{\overset{L}{.Math.}}{}_{\mathrm{ijl}}.Math.x_l^{\left(m\right)}.Math.\left[k-n_{\mathrm{ijl}}\right].& \left(59\right)\end{array}$
The term q.sub.ij[k] can then be expressed as

(93) $\begin{array}{cc}{q}_{\mathrm{ij}}\left[k\right]\stackrel{}{=}\underset{l{\mathrm{��}}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}& \left(60\right)\end{array}$
where d.sub.ijl.sup.(m)=.sub.ijl.Math.x.sub.l.sup.(m) and .sub.k={l=1, 2, . . . , L: n.sub.ijl=k}.

(94) In other words, the term q.sub.ij[k] is computed by accumulating all elements of
d.sub.ijl.sup.(m)=[d.sub.ij1.sup.(m),d.sub.ij2.sup.(m), . . . ,d.sub.ijL.sup.(m)](61)
for which associated discrete-time delay indexes n.sub.ijl have the same integer value k. Consequently, q.sub.ij[k] can then be computed in a very efficient manner using the hash table data structure concept. The indexes of a hash table, typically referred to as keys, are the integers between k.sub.min and k.sub.max, and the record associated with the k.sup.th key is the set of values
{d.sub.ijl.sup.(m):l=1,2, . . . ,L; n.sub.ijl=k}.(62)
By one approach, the hash-table-based computation of q.sub.ij[k] is implemented using a processing device configured to use MATLAB using the accumarray function. The variables d, n and q store the following sequences:
dd.sub.ijl.sup.(m)=[d.sub.ij1.sup.(m),d.sub.ij2.sup.(m), . . . ,d.sub.ijL.sup.(m)](63)
nn.sub.ijl=[n.sub.ij1,n.sub.ij2, . . . ,n.sub.ijL](64)
qq.sub.ij=[q.sub.ij[1],q.sub.ij[2], . . . ,q.sub.ij[k.sub.max]](65)
The variable q is computed via the command q=accumarray(n,d) where k.sub.minn.sub.ijlk.sub.max for all l, q.sub.ij[k]=0 for all indexes k<k.sub.min. An example of pseudocode to be run by the processing device for implementation of the proposed algorithm for efficiently computing G.sub.l.sup.(m) is given below.

(95) TABLE-US-00001 .square-solid. Subroutine 1: Pseudocode for computing G.sub.l.sup.(m) for l = 1, 2, . . . , L for i = 1, 2, . . . , I do for j = 1, 2, . . . , J do SET q.sub.ij[k] = 0 for 0 k < k.sub.min for k = k.sub.min, k.sub.min + 1, . . . , k.sub.max do S.sub.k = {l = 1, 2, . . . , L:n.sub.ijl = k} $q_{\mathrm{ij}}\left[k\right]=\underset{l{𝒮}_k}{.Math.}{d}_{\mathrm{ijl}}^{\left(m\right)}\left(\mathrm{hash}\text{-}\mathrm{table}\text{-}\mathrm{implementation}\right)$ end for .sub.ij.sup.(m)[n] = (q.sub.ij * w)[n] for l = 1, 2, . . . , L do .sub.ijl.sup.(m) = {(w * (y.sub.ij .sub.ij.sup.(m)))[k]}|.sub.k=n.sub.ijl end for end for end for ${\hat{G}}_l^{\left(m\right)}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{}_{\mathrm{ijl}}.Math.{\hat{G}}_{\mathrm{ijl}}^{\left(m\right)}$

(96) In addition to being more computationally efficient, the proposed implementation does not require constructing the large system matrix A. A tangible benefit of this fact is the size of data (i.e., the number of transmit locations) that can be used to form an image is no longer limited. It is also readily observed from the pseudocode that the computation of .sub.l.sup.(m) is parallelizable such that faster processing techniques such as parallel or GPU based processing can be used to process the data.

(97) Fast Implementation of H.sub.l

(98) An alternative expression for H.sub.l in (6) is

(99) $\begin{array}{cc}{H}_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}\underset{n=1}{\overset{N}{.Math.}}{\left({}_{\mathrm{ijl}}\right)}^2.Math.r_{\mathrm{ijn}}.Math.\left({\mathrm{nT}}_s-{}_{\mathrm{ijl}}\right),& \left(66\right)\end{array}$
where (t)p.sup.2(t) and r.sub.ijn is the number of non-zero elements in the n.sup.th row of the NL sub-matrix A.sub.ij=P.sub.ijD.sub.ij. To facilitate a fast implementation, we approximate the quantity H.sub.l by

(100) $\begin{array}{cc}{\hat{H}}_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}\underset{n=1}{\overset{N}{.Math.}}{\left({}_{\mathrm{ijl}}\right)}^2.Math.r_{\mathrm{ijn}}.Math.\left({\mathrm{nT}}_s-n_{\mathrm{ijl}}{T}_s\right),& \left(67\right)\end{array}$
We write .sub.l as

(101) $\begin{array}{cc}{\hat{H}}_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{\left({}_{\mathrm{ijl}}\right)}^2.Math.{\hat{H}}_{\mathrm{ijl}}\mathrm{where}& \left(68\right)\\ \begin{array}{c}{\hat{H}}_{\mathrm{ijl}}=\underset{n=0}{\overset{N-1}{.Math.}}{r}_{\mathrm{ijn}}.Math.\left({\mathrm{nT}}_s-n_{\mathrm{ijl}}{T}_s\right)\\ =\underset{n=0}{\overset{N-1}{.Math.}}{}_{\mathrm{ij}}\left[n\right].Math.h\left[n-n_{\mathrm{ijl}}\right]\left(70\right)\\ =\left\{\underset{n=1}{\overset{N}{.Math.}}{}_{\mathrm{ij}}\left[n\right].Math.h\left[n-k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\left(71\right)\\ =\left\{\underset{n=1}{\overset{N}{.Math.}}{}_{\mathrm{ij}}\left[n\right].Math.h\left[k-n\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\left(72\right)\\ =\left\{\left({}_{\mathrm{ij}}*h\right)\left[k\right]\right\}{.Math.}_{k=n_{\mathrm{ijl}}}\left(73\right)\end{array}& \left(69\right)\end{array}$
with h[n](nT.sub.s), and r.sub.ijn is now represented by the n-indexed sequence .sub.ij[n]r.sub.ijn. For the sake of convenience and consistency, we assume here that rows of a matrix are counted starting from a zeroth row. The computation .sub.ijl requires the convolution of the squared and discretized pulse h[n] with the sequence .sub.ij[n]. The computation of .sub.ijl is significantly accelerated with the introduction of a fast procedure for generating .sub.ij[n].

(102) First, we recall that the sample .sub.ij[n] is the number of non-zeros entries in the n.sup.th row of the NL sub-matrix A.sub.ij. Because the radar system has an ultra wide band, the transmitted pulse p(t) is short. Consequently, the samples of the length-N sequence w[n] are zero (or practically zero) for indexes n such that |n|<M and non-zero, otherwise. The parameter M is even and significantly smaller than N. The l.sup.th column of A.sub.ij coincides with the length-N vector
[.sub.ijl.Math.p(0.sub.ijl),.sub.ijl.Math.p((T.sub.s.sub.ijl), . . . ,.sub.ij.Math.p((N1)T.sub.s.sub.ijl)].sup.T.(74)
The (n,l)-entry of A.sub.ij is thus non-zero if
|nT.sub.s.sub.ijl|MT.sub.s.(75)
Using (44), the above rule in (75) can be approximated by
|nn.sub.ijl|M.(76)

(103) A computed delay index n.sub.ijl is such that 0n.sub.ijlN. Consequently, for computational convenience, we write that the (n,l)-entry of A.sub.ij is non-zero if
max(0,nM)n.sub.ijlmin(n+M,N).(77)
The number .sub.ij[n] of non-zeros entries in the n.sup.th row of A.sub.ij is thus equal to the number of elements in the n.sup.th row that satisfy (77). A more convenient definition is
.sub.ij[n]=|.sub.n|(78)
where |.sub.n| denotes the number of elements in the set
.sub.n={l=1,2, . . . ,L|max(0,nM)n.sub.ijlmin(n+M,N)}.(79)

(104) The parameter .sub.ij[n] can be efficiently computed by taking advantage of the hash-table-based fast implementation concept used in (60). First, we write
.sub.ij[n]=|.sub.n.sup.+||.sub.n.sup.|(80)
where
.sub.n.sup.+={l=1,2, . . . ,L|n.sub.ijlmin(n+M,N)}(81)
.sub.n.sup.={l=1,2, . . . ,L|n.sub.ijlmax(0,nM1)}.(82)
The expression in (80) is further expanded as

(105) $\begin{array}{cc}{}_{\mathrm{ij}}\left[n\right]=\underset{k=0}{\overset{\min \left(n+M,N\right)}{.Math.}}.Math.{\mathrm{��}}_k.Math.-\underset{k=0}{\overset{\max \left(0,n-M\right)}{.Math.}}.Math.{\mathrm{��}}_k.Math.& \left(83\right)\end{array}$
Finally, we have
.sub.ij[n]=[min(n+M,N)][max(0,nM)](84)
where

(106) $\begin{array}{cc}v\left[m\right]=\underset{k=0}{\overset{m}{.Math.}}.Math.{\mathrm{��}}_k.Math.& \left(85\right)\end{array}$
with .sub.k={l=1, 2, . . . , L: n.sub.ijl=k}. The inner summation in (85) (and, hence the computation of [m]) is efficiently computed using the hash-table-based fast implementation previously discussed and used in (60). Example pseudocode to be run by the processing device for implementation of the proposed algorithm for efficiently computing H.sub.l is given below.

(107) TABLE-US-00002 .square-solid. Subroutine 2: Pseudocode for computing H.sub.l for l = 1, 2, . . . , L for i = 1, 2, . . . , I do for j = 1, 2, . . . , J do SET q[k] = 0 for 0 k < k.sub.min for k = k.sub.min, k.sub.min + 1, . . . , k.sub.max do S.sub.k = {l = 1, 2, . . . , L:n.sub.ijl = k} q[k] = |S.sub.k| (hash-table-implementation) end for for m = 0, 1, . . . , N do $\left[m\right]=\underset{k=0}{\overset{m}{.Math.}}q\left[k\right]$ end for for n = 0, 1, . . . , N do .sub.ij[n] = w[min(n + M, N)] w[max(0, n M)] end for for l = 1, 2, . . . , L do .sub.ijl = {(.sub.ij * h)[k]}.sub.k=n.sub.ijl end for end for end for 0${\hat{H}}_l=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{\left({}_{\mathrm{ijl}}\right)}^2.Math.{\hat{H}}_{\mathrm{ijl}}$

(108) Putting together the results for calculating terms G.sub.l.sup.(m) and H.sub.l, example pseudocode for the l.sub.1-LS algorithm follows.

(109) TABLE-US-00003 .square-solid. Pseudocode for computing l.sub.1-LS algorithm for m = 1, 2, . . . , num.sub.it Initialization: x.sup.(0) = {x.sub.1.sup.(0), x.sub.2.sup.(0), . . . , x.sub.L.sup.(0)} for l = 1, 2, . . . , L do Compute .sub.l (via Subroutine 2) end for for m = 1, 2, . . . , num.sub.it do for l = 1, 2, . . . , L do Compute .sub.l.sup.(m) (via Subroutine 1) end for $x_l^{\left(m+1\right)}=\frac{{\hat{H}}_l.Math.x_l^{\left(m\right)}+{\hat{G}}_l^{\left(m\right)}}{{\hat{H}}_l+\frac{}{2\left|x_l^{\left(m\right)}\right|}}$ end for

(110) So configured, the described MM-based l.sub.1-LS algorithm is applicable to large-scale, real applications. Although the proposed algorithm effectively estimates reflection coefficients of scenes-of-interest using GPR datasets, the algorithm could be readily applied to a variety of applications where datasets are collected using synthetic aperture imaging measurement principles. When compared to images produced by the DAS or RSM algorithms, the image obtained using the MM-based l.sub.1-LS algorithm is more accurate, is less noisy, and captures the main scatterers in the scene-of-interest while effectively suppressing shadows and side lobes. Although the proposed algorithm is still more computationally expensive than the DAS algorithm, a derived acceleration technique produces a fast-implementation version that is very fast and requires substantially less memory. Moreover, because the algorithm decouples the estimation of individual reflectance coefficients, further computational speed gains are achievable via parallel and GPU processing implementations.

(111) By one approach, the method described above can be implemented as illustrated in FIG. 6, where a processing device receives 605 an initial data set and processes 610 the initial data set by creating an estimated image value for each voxel in the image by iteratively deriving the estimated image value through application of a majorize-minimize principle to solve an l.sub.1-regularized least-squares estimation problem associated with a mathematical model of image data from the initial data set. These basic steps can be applied to achieve fast and computationally efficiently prepared image using the estimated image value of individual voxels of the image that can be displayed 615, wherein the initial data set may be sourced from a variety of applications where datasets are collected using synthetic aperture imaging measurement principles. In the GPR context, the method may further include emitting 620 a radar pulse at specified intervals into a scene-of-interest and detecting 625 magnitude of signal reflections from the scene of interest from the radar pulse. Position data is recorded 630 corresponding to individual radar pulse emissions and individual receptions of the signal reflections. The initial data set in this application is created 635 from the position data and detected magnitudes of the signal reflections. Where the method is carried out remote from the vehicle, it is sufficient where the receipt of the initial data set to be processed includes receiving data representing transmission site locations of radar pulses, reception site locations of reception of reflections from the radar pulses, radar-return profiles for pairings of the transmission site locations and the reception site locations, and data samples associated with individual radar-return profiles.

(112) Results for the MM-Based l.sub.1-LS Algorithm

(113) The performance of the MM-based l.sub.1-LS algorithm can be evaluated using a numerical experiment and using a real dataset as obtained using a UWB SIRE apparatus and provided by the US Army Research Laboratory.

(114) With reference to FIG. 7, the numerical experiment demonstrates that the accuracy of the proposed algorithm is comparable to that of existing standard l.sub.1-LS algorithms despite the described computational efficiencies over these algorithms. In the numerical experiment, a length-N data vector y is generated using the standard additive white Gaussian (AWGN) noise model y=Ax+w, where x is a length-L vector of regression coefficients and w is an AWGN vector with variance .sup.2. The numerical values for the parameters are N=500, L=23 and .sup.2=1. The 50023 system matrix A is randomly generated. FIG. 7 illustrates the estimation accuracy of the proposed MM-based l.sub.1-LS algorithm as compared to that of previously used algorithms including the LASSO, the shooting, and the standard l.sub.1-LS algorithms. As illustrated, the accuracies of these various algorithms are substantially overlapping meaning that the accuracies are approximately the same. As discussed above, while all algorithms give comparable performance results where the data size is relatively small, when the data size is large such as in a GPR application, only the proposed l.sub.1-LS can produce a result using reasonable computing resources, i.e., processing and memory resources. For instance, the LASSO, the shooting, and the l.sub.1-LS algorithms fail due to memory size limitations.

(115) FIGS. 8-10 illustrate the images provided by application of the MM-based l.sub.1-LS algorithm (FIG. 8) to data collected by UWB SIRE system as compared to the images provided by the DAS (FIG. 9) and RSM (FIG. 10) algorithms as applied to the same data set. This test data set corresponds to measurements taken from I=274 consecutive transmit locations using J=16 receive antennas. The scene-of-interest is of size 6525 m.sup.2, and is divided into a grid of 250 voxels in the cross-range direction and 3200 voxels in the down-range direction. The voxels have each 0.1 m in the cross-range direction and 0.02 m in the down-range direction. These parameters and their values are summarized below in Table 1.

(116) TABLE-US-00004 TABLE 1 Parameters for the real data Parameter Value Image dimension 25 m (cross-range) by 65 m (range) Voxel size 0.1 m (cross-range) by 0.02 m (range) Number of transmit locations 274 Transmit locations per sub-aperture 43 Sub-aperture dimension (cross-range) 25 m (cross-range) by 2 m (range)

(117) FIGS. 8-10 show complete images of the scene-of-interest that are generated using the DAS, MM-based l.sub.1-LS, and RSM algorithms. FIG. 9 shows the image obtained using the DAS algorithm. The image has significant side lobes and shadows. Moreover, the presence of background noise is clearly visible. FIG. 10 shows the image obtained using the RSM algorithm. Although, side lobes and shadows are reduced, there is still room for improvement. The image obtained using the proposed l.sub.1-LS algorithm, shown in FIG. 8, is sparser and adequately suppresses both the side lobes and background noise.

(118) Another Approach: MM-Based Least Absolute Deviation (LAD) Algorithm with l.sub.1-Regularization

(119) A second approach to application of the MM principle to processing an image data set includes application of this principle to a different approach to the least squares technique. More specifically, such a method includes creating an estimated image value for each voxel in the image by iteratively deriving the estimated image value through application of an MM principle to solve the l.sub.1-regularized least-absolute deviation (LAD) estimation problem associated with a mathematical model of image data from the initial data set.

(120) The so called l.sub.1-regularized LAD estimation problem is a known approach to the least squares regression analysis. This approach is known to handle outlier data in a better or more robust fashion, but at the cost of increased computational resources. In this approach, the reflectance coefficient vector, which represents objects in the SOI to be detected, is estimated using the l.sub.1-regularized least absolute deviation (l.sub.1-LAD) method:

(121) $\begin{array}{cc}\hat{x}=\underset{x}{\arg \min }{.Math.y-\mathrm{Ax}.Math.}_1+{.Math.x.Math.}_1,& \left(86\right)\end{array}$
where is the regularization parameter. We solve the optimization problem in (86) using the MM principle as shown by D. R. Hunter and K. Lange, A tutorial on mm algorithms, The American Statistician, vol. 58, no. 1, pp. 30 to 37, 2004, which is incorporated herein by reference. The resulting algorithm is straightforward-to-implement, computationally efficient and amenable to parallel (or distributed) implementations.

(122) The MM principle is described above. To solve the GPR image formation problem in (86) using the MM principle, the objective function to be minimized can be written as
(x)=.sub.1(x)+.sub.2(x)(87)
where .sub.1(x)=yAx.sub.1, .sub.2(x)=x.sub.1, and the regularization parameter is positive. A quadratic majorizer for the absolute value function (x)=|x| is given by De Leeuw and Lange as referenced above. From that result, it directly follows that for any real x.sup.(m)0

(123) $\begin{array}{cc}{Q}_2\left(x,x^{\left(m\right)}\right)=\underset{l=1}{\overset{L}{.Math.}}\frac{x_l^2}{2.Math.x_l^{\left(m\right)}.Math.}+\frac{1}{2}.Math.x_l^{\left(m\right)}.Math.& \left(88\right)\end{array}$
is a majorizing function for .sub.2(x)=.sub.l=1.sup.L|x.sub.l| at the point x.sup.(m).

(124) A majorizing function for .sub.1(x) is now constructed by first replacing the absolute value function by De Leeuw and Lange's majorizing function for the absolute value function

(125) $\begin{array}{cc}\begin{array}{c}{}_1\left(x\right)=\underset{k=0}{\overset{K-1}{.Math.}}.Math.y_k-{\left[\mathrm{Ax}\right]}_k.Math.\\ \underset{k=0}{\overset{K-1}{.Math.}}\frac{{\left(y_k-{\left[\mathrm{Ax}\right]}_k.Math.\right)}^2}{2.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}+\frac{.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}{2}.\left(90\right)\end{array}& \left(89\right)\end{array}$
Next, we replace the term (y.sub.k[Ax].sub.k|).sup.2 above by a majorizing function developed by De Pierro as referenced above. More specifically, De Pierro showed that ([Ax].sub.k).sup.2 is majorized by

(126) $\begin{array}{cc}q\left(x,x^{\left(m\right)}\right)\stackrel{}{=}\underset{l=1}{\overset{L}{.Math.}}{C_{\mathrm{kl}}\left(N_k{A}_{\mathrm{kl}}{x}_l-N_k{A}_{\mathrm{kl}}{x}_l^{\left(m\right)}+{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k\right)}^2& \left(91\right)\end{array}$
where N.sub.k is the number of non-zero elements in the k.sup.th row of the system matrix A and C.sub.kl is defined by

(127) $\begin{array}{cc}{C}_{\mathrm{kl}}\stackrel{}{=}\{\begin{array}{cc}{N}_k^{-1},& A_{\mathrm{kl}}0\\ 0,& A_{\mathrm{kl}}=0\end{array}.& \left(92\right)\end{array}$
It then follows that
(y.sub.k[Ax].sub.k).sup.2=y.sub.k.sup.22y.sub.k[Ax].sub.k+([Ax].sub.k).sup.2(93)
y.sub.k.sup.22y.sub.k[Ax].sub.k+q(x,x.sup.(m)).(94)
where q is given by (28). From (90) and (94) it can be seen that a majorizing function for .sub.1 at the point x.sup.(m) is

(128) $\begin{array}{cc}{Q}_1\left(x,x^{\left(m\right)}\right)\stackrel{}{=}\underset{k=0}{\overset{K-1}{.Math.}}\frac{y_k^2-2{y_k\left[\mathrm{Ax}\right]}_k+q\left(x,x^{\left(m\right)}\right)}{2.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}+\underset{k=0}{\overset{K-1}{.Math.}}\frac{1}{2}.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math..& \left(95\right)\end{array}$
Because the penalty factor is positive, a majorizing function for the objective function (x)=.sub.1(x)+.sub.2(x) is
Q(x,x.sup.(m))=Q.sub.1(x,x.sup.(m))+Q.sub.2(x,x.sup.(m))(96)
where Q.sub.2 is given by (31). From the general expression in (17), we obtain the desired iterative algorithm by setting to zero the partial derivatives of Q(x,x.sup.(m)) with respect to the components of x. For l=1, 2, . . . , L,

(129) $\begin{array}{cc}\frac{Q}{x_l}=\frac{}{x_l}\underset{k=0}{\overset{K-1}{.Math.}}\frac{y_k^2-2{y_k\left[\mathrm{Ax}\right]}_k+q\left(x,x^{\left(m\right)}\right)}{2.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.}+\frac{}{x_l}\underset{k=0}{\overset{K-1}{.Math.}}\frac{1}{2}.Math.y_k-{\left[{\mathrm{Ax}}^{\left(m\right)}\right]}_k.Math.+.Math.\frac{}{x_l}\underset{l=1}{\overset{L}{.Math.}}\left[\frac{x_l^2}{2.Math.x_l^{\left(m\right)}.Math.}+\frac{1}{2}.Math.x_l^{\left(m\right)}.Math.\right].& \left(97\right)\end{array}$
Computing the derivatives in (97), simplifying and re-arranging terms gives