Near-Field Resolution in Planar Source Reconstructions

—This paper deals with the classical question of estimating the achievable resolution in terms of the conﬁguration parameters in inverse source problems. In particular, the study focuses on the case of a planar surface magnetic current which is to be reconstructed from near-ﬁeld observed over a bounded rectangular aperture parallel to the source domain. Resolution formulas are well known for far-ﬁeld or Fresnel-zone conﬁgu- rations, and also for near-ﬁeld cases, when data are full-view or the measurement aperture is an unbounded plane. For the case of bounded near-ﬁeld observations the resolution estimations mainly rely on some asymptotic arguments in order to mimic far-ﬁeld reasoning. Here, the plan is to improve these results and to work out a resolution estimation that precisely captures the spatially varying behaviour entailed by the near-ﬁeld and aspect-limited conﬁguration. To this end, the pertinent radiation operator is inverted by an adjoint inversion scheme (a back-propagation-like method) and the corresponding point-spread function is analytically estimated. Numerical examples show that the derived resolution estimation clearly points out the role of the geometrical parameters of the conﬁguration and it is more accurate than other literature results.


I. INTRODUCTION
The reconstruction of a current from its radiated field is a classical problem in electromagnetics [1], which besides being theoretically intriguing, is relevant in a number of applications. Just to quote a few of them, we mention antenna synthesis [2], [3] and/or diagnostics [4] - [6], near-field to far-field transformation [7], [8] or near-zone RCS estimation [9], [10]. Inverse source problems are also strictly linked to linearised inverse scattering problems, where similar integral operators have to be dealt with [11].
In this framework, the achievable resolution plays a key role since it says which is the finest detail that can be reconstructed and it is also linked to other important figures like the number of degrees of freedom of the problem [12] and the information content [13] that can be transmitted from a source region to an observation domain [14]. In particular, because of the ill-posedness of inverse source problems, which is basically due to the filtering introduced by the propagator, a trade-off between accuracy and stability must be established [15]; hence the achievable resolution results limited and in general dependent on the noise level. However, far from the Green function's singularity (i.e., when evanescent waves are Maria Antonia Maisto is with the Department of Engineering, Universitá della Campania, Aversa, (Ce), 81031, Italy, e-mail: mariaantonia.maisto@unicampania.it.
Manuscript received April 19, 2005; revised August 26, 2015. negligible), the kernel of the radiation operator behaves like an entire function of exponential type. As a consequence, the corresponding singular values present a distinctive behaviour which is characterised by an abrupt exponentially quick decay beyond a certain critical index [16]. This entails that, unless a very high (often impractical) signal-to-noise ratio (SNR) is available, resolution is weakly dependent on the noise and mainly related to the configuration parameters. This is indeed the mathematical rationale (often implicitly assumed) behind most of the studies reporting resolution estimation in terms of "only" the parameters of the configuration [17]- [24]. Of course, for very near-zone configurations, the evanescent contribution cannot be neglected and so the SNR turns to play a prominent role on the achievable resolution [25], [26]. However, in this contribution, very near-zone configurations are ruled out. The estimation of the achievable resolution for far-field configurations greatly benefits from the k-space formulation which allows to highlight the portion of the source spatial spectrum that can be reconstructed and therefore to estimate the resolution accordingly. In particular, analytical estimations can be easily obtained for full-view data configurations [1]. For aspect-limited cases, the shape of the spectrum region does not in general permit to perform the Fourier transformation in closed form. In this regard, it is worth mentioning the study reported in [27], where the authors succeed in addressing aspect-limited configurations by finding analytical approximations for the resolution, remarkably even by accounting for the full vector nature of the problems.
For near-zone configurations (in the sense explained above), the k-space approach can still be employed when data are collected over an unbounded line [19] or plane [28]. For fullview data, analytical resolution estimations can be obtained by the multipole expansion if canonical measurement curves or surfaces are considered [11], or again by means of the kspace method under the framework of generalized holography and the Porter-Bojarski integral equation inversion, as shown in [29] and [30]. Indeed, more recently it has been shown that the case of general observation curve can be also successfully addressed by recasting the point-spread function as a Fourier type integral operator [10].
For near-zone aspect-limited configurations previous methods cannot be rigorously applied because truncation affects the computation of the Fourier transform stage. In these cases, the retrievable spectrum region can still be approximately determined by resorting to stationary phase asymptotic arguments [22]. What is interesting here is that the spectrum region is spatially varying, which entails that (differently form farfield cases or full-view or unbounded observation near-zone configurations) the achievable resolution is spatially varying. This is a well known fact observed many times in literature [24]. Nonetheless, often far-field estimations are still employed [21] or the spatially varying behaviour is ignored since the union of the spectrum regions is used to estimate the resolution [24].
In this paper we are concerned with the study of the achievable resolution in the reconstruction of planar sources from near-field field aspect-limited observations. In particular, we consider a magnetic planar surface current which radiates in free-space and with its radiated field observed over a rectangular aperture. The main aim is to provide an analytical estimation of the achievable resolution in terms of the geometrical parameters of the measurement aperture which is more accurate than the previous estimation based on asymptotic arguments [22], [24] and that captures its spatially varying behaviour. To this end, the resolution is estimated by finding an analytical approximation of the pointspread function obtained by a back-propagation like inversion scheme. In particular, we use the same approach developed in [18] for the case of strip currents. Basically, thanks to suitable variable transformations, the point-spread function is recast as a spatially varying band-limited function whose spatially varying band is then enclosed within the smallest spatially varying rectangular domain containing it. The derived resolution estimation is numerically checked and shown to be in excellent agreement with the one yielded by a truncated singular value decomposition scheme, which is here used as benchmark. The paper also includes a comparison with some literature results addressing similar configurations.

II. PROBLEM STATEMENT
Consider a magnetic current J of bounded finite planar support SD whose radiated field is observed over another planar domain, the observation domain OD, located in near-zone. For the sake of simplicity, we assume SD and OD being the located at z = 0 and z = z o , respectively. The source is assumed to be directed in the x − y plane whereas only the tangential components of the radiated field are collected. Under this framework, the problem is split as two scalar problems that can be addressed in the same way. Therefore, here we just consider the current directed along the x-axis , i.e, J = J(x, y)x, and to collect the corresponding tangential y-component of the radiated field (see Fig. 1). Also, OD is assumed larger than SD. The problem is described in the frequency domain by the following radiation operator (unless an unessential factor) with Φ(r o , r) = k|r o − r| and K(r o , r) We are interested in the estimation of the achievable resolution while inverting (1). This formally entails considering a Dirac-delta like source whose reconstruction is the socalled point-spread function ( Fig. 2 illustrates this process). In particular, using the singular system of A , the point-spread function is expressed as where * means conjugation. Note that the summation in (2) is stopped at a certain index N , which is in general dependent on the noise level . This reflects the need to regularize the inversion in order to counteract noise propagation. As a consequence, the achievable resolution results limited. Eq. (2) allows to find the point-spread function (and hence the resolution could be deduced form it) in terms of the singular functions of A. However, for the case at hand, those singular functions are not known in closed form. Of course, one can numerically compute them but this would not give an analytical resolution estimation. Therefore, in the sequel, the regularized inverse is approximated by the adjoint operator A † . Note that this method is basically a back-propagation which is reminiscent of migration, time-reversal and other similar approaches, which are very commonly used in Radar imaging [21], [31], [32]. Accordingly, the point-spread function can be estimated by evaluating the following integral Naturally, it is necessary to establish how (3) relates to (2), which is the actual point-spread function. To this end, it is convenient to rewritep sf in terms of the singular system of A, that isp where σ n are the singular values. Since the σ n tends to zero, the adjoint based inversion is a stable procedure. Moreover, since the singular values exhibit a "step-like" behaviour ( Fig.  3 shows some examples) thep sf is practically coincident with psf when N in (2) is chosen in correspondence of the index for which the singular values start to abruptly decay. Note that, as remarked in the introduction, because of this singular value decay, the resolution is actually weakly dependent on the noise and hence we are entitled to look for an analytical expression which only highlights the role of the configuration parameters. Also, by the adjoint inversion, there is no need to estimate such an index since in (3) truncation is implicitly obtained by the windowing imposed by the singular values themselves.
Eventually, the problem of resolution estimation is cast as the evaluation of (3).

III. POINT-SPREAD FUNCTION EXPRESSION
In order to obtain an analytical approximation of the pointspread function in (3), the main idea it to recast it as a Fourierlike transformation. To this end, we rewrite the phase term as such that p(ν) is a curve whose starting and ending points coincide with r and r, respectively, that is p(ν 0 ) = r and p(ν 1 ) = r. Now, the curve p(ν) can be properly chosen in order to let the phase term resemble a Fourier kernel. This can be achieved, for example, in the following way. Considerr ≡ (x, y ) and then perform integration in (5) along the polygonal line with nodes r ,r and r. Accordingly, we have that where · denotes the scalar product, w ≡ (w x , w y ) and dν (8) It can be easily shown that ∀r, r the transformation w : r o → w(r o , r, r ) is injective and the corresponding Jacobiam matrix full rank. This allows to replace integration in r o with integration in w, which in turn allows to rewrite (3) aŝ psf (r, r ) = Ω(r,r ) H(r, r w)e jw·(r−r ) dw (9) with (10) and with ∂(xo,yo) ∂(wx,wy) being the Jacobian determinant of the introduced transformation. It is instructive to have an idea of what Ω(r, r ) looks like. This is particularly simple when r = r . In fact, in this case and Now, if we look at how lines parallel to the y-axis map in the w plane, as the intercept x o changes, we obtain the following family of ellipses (parametric with respect to x o ) . As a result, the strip −X 0 ≤ x o ≤ X 0 of the r o plane maps into the convex domain of the w plane contained between the two ellipses obtained from (14) for Hence, Ω(r) is given by the intersection of the two domains returned by (14) and (15). Some examples of Ω(r) are shown in Fig. 4. It is interesting to highlight what happens when X 0 and Y 0 approach ∞. In this case, Ω(r) is independent on r and becomes the circle of radius k ∀r ∈ SD. What is more, because of (7) and (8), this holds true also for the band described in (10). This is expected, since for non-aspect limited configuration, the spatially varying behaviour is lost and the retrievable spatial spectrum coincides to the visible circle. In this case, the point-spread function can be easily computed and turns out to be ∝ J 1 (k|r − r |)/|r − r |. Note that this result could have also been obtained by resorting to the plane-wave expansion of the propagator [11] which naturally leads to the k-space approach.
Here, indeed, we are interested in aspect-limited configuration for which things are more involved. In order to simplify a little bit the matter, we note that, because H is a constant sign function, (9) clearly shows that the leading order contribution occurs for r − r = 0 [33]. This entitles us to approximate the amplitude factor by its value assumed for r = r , that is H(r, r , w)) ≈ H(r , r , w) = H(r , w). Therefore, in order to compute H(r , w), we have to consider which is given by (12). The corresponding Jacobian transformation then yields from which finally we have H(r , w) = 1. Eventually, the point-spread function is approximated as followŝ Some comments concerning (18) are now in order. First, since the assumption r = r has been used, we mainly expect that (18) works in approximating the point-spread function around its main beam. However, this is what is needed for resolution estimation. Second, (18) shows the point-spread function as a 2D spatially varying band-limited function [34], which basically entails that resolution will be spatially varying. This of course is a distinctive feature of near-zone aspect-limited configuration and has been already observed in literature by using more approximate (and hence less accurate resolution estimation) arguments [24], [35]. Third, we have already followed a similar procedure for strip currents [18]. In that case we were able to solve the integral expressing the point-spread function by introducing a further non-linear warping transformation [36] of the source spatial variable which allowed to cast the point-spread function as a classical band-limited function with respect to this new variable. In this case, the spatially varying behaviour was embodied within the non-linear transformation.
Here, we plan to apply the same approach. However, now the matter is much more difficult because, unlike the strip current case where spatially varying behaviour manifested in terms of the varying size of the band only, here, both the size and the shape of the band have changed.

RESOLUTION ESTIMATION
According to the previous discussion, the main problem to be faced in order to evaluate (18) is the change in shape that Ω(r, r ) undergoes while r and r range over the source domain SD. To cope with this point, herein, we are content to find only an approximation for (18). To this end, in (18), we consider the smallest rectangular domain Ω R (r, r ) which contains Ω(r, r ), with r, r ∈ SD, and denote bŷ psf R the corresponding point-spread function approximation. Accordingly, say 2∆w x (x, r ) and 2∆w y (r, y ) the sides of Ω R (r, r ) and w mx (x, r ) and w my (r, y ) their middle points, the point-spread function assumes the following more convenient expression with w m = (w mx , w my ).
Eq. (19) is the sought after expression for the point-spread function. However, to be used, w m , ∆w x and ∆w y , that is Ω R (r, r ), still remain to be determined. To achieve this end we have to compute w max and w min y (r, y ) = min ro∈OD {w y (r o , r , y )}. As mentioned above, the Jacobian of the transformation w : r o → w(r o , r, r ) is full rank. Accordingly, both w x and w y cannot have stationary points inside OD. Therefore, their maximum and minimum must be looked for over the boundary of the observation domain. After simple but tedious calculations, it results that from which it readily follows that Now, we can estimate the resolution from the main-beam of the point-spread function when (21) is inserted in (19). In particular, by doing so and by further introducing the transformations then eq. (19) can be rewritten in a more convenient way as where we have considered the magnitude of the point-spread function. If now we denote as ∆x and ∆y the resolution along x and y, from (23) it readily follows that and ξ y (y + ∆y) − ξ y (y) = π As could be expected from (19), (20) and (21), equations (24) and (25) express the achievable resolution in a factorised form with ∆ x being dependent on the extent of OD (i.e., X 0 ) along x and ∆ y being solely related to Y 0 . This is not only due to the considered spectral rectangular domain Ω R in the computation of (18), but also to the polygonal path used in (5) and because, since we have assumed , ∆w x (x, x ) (∆w y (y, y )) depends only on x-component (y-component) of r and r . It is also important to note that the shape of SD does not change the obtained resolution estimation.
Eqs. (24) and (25) are indeed implicit expression of the achievable resolution. Indeed, they are equal to the ones obtained for the case of strip currents because of the occurred point-spread function factorization. As such, the same manipulations adopted in [18] can be employed to explicitly solve them for ∆ x and ∆ y . This yields and

3) 4)
(c) The source point (x , y ) is located at (4, 0λ).  Eqs. (26) and (27) predict that ∆x and ∆y are smaller at the centre of the source region whereas they increase when source point moves toward the edge of SD. Moreover, such a spatially varying behaviour becomes more marked when z o increases or the observation domain size decreases. Finally, as X 0 , Y 0 → ∞ both ∆x and ∆y tend to become constant (i.e., the spatially varying resolution behaviour is lost) and approach to λ/2, λ being the wavelength. This is consistent with the point-spread function pertaining to the unbounded observation domain previously reported in Section III. Indeed, in this case Ω(r, r ) tends to the circle of radius k (independently from r and r ) and (26) and (27) just refer to the square spectrum domain of side 2k containing that circle. Of course in that case there is no need to go through our estimation procedure since closed form results can be easily obtained as (actually) reported above.

A. Numerical check
In order to check the accuracy of the resolution estimations (26) and (27), in this section we compare |psf | and |p sf R | as the source location varies in SD. In particular, the psf is obtained by choosing N in (2) in correspondence to the index for which the singular values start to abruptly decay. The example in Fig. 5 shows the normalized point-spread function amplitudes for X s = Y s = 5λ, X 0 = 6λ, Y 0 = 5λ and z o = 4λ. In order to appreciate the expected spatially varying behavior, four source locations (x , y ) are considered, i.e., (0, 0), (0, 4λ), (4λ, 0) and (4λ, 4λ). More in detail, each sub-figure of Fig. 5 reports |psf | (panel (1) ), |p sf R | (panel (2) ) and the comparison of the corresponding cuts (panels (3) and (4)) along y = y and x = x , respectively. As can be seen, the point-spread functions' main-lobes overlap, meaning that (26) and (27) works very well in predicting the achievable resolution. As expected, the side-lobe structure of the estimated point-spread function differs from that of |psf | because we have considered Ω R (r, r ) instead of Ω(r, r ) in (18). However, as remarked above, the point-spread function main lobe is what matters for resolution estimation. Moreover, resolution degrades while moving away from the point (0, 0). In particular, when the source point is moved along y, as expected, ∆ x remains unchanged whereas ∆ y increases; the opposite occurs when the movement occurs along x axes. Note that the increasing of ∆ x appears less marked since the observation domain has been chosen rectangular with X 0 > Y 0 . This behavior is perfectly, qualitatively and quantitatively, captured by our theoretical estimations.

B. Number of degrees of freedom
The resolution can be conveniently used to estimate the socalled number of degrees of freedom (NDF) of the problem as well. Rigorously speaking the NDF is the dimension of the subspace where the unknown can be reliably reconstructed or equivalently the dimension of the subspace where the field data projects more "significantly". Accordingly, it is related to the singular value behaviour of the radiation operator. In particular, the NDF can be computed as the number of singular values that are above a given noise dependent threshold. As discussed above, the properties of the radiation operator are such that, when evanescent waves are negligible, the singular values exhibit a step-like behaviour and hence the NDF becomes weakly dependent on the noise level and practically coincides with those that precede the abrupt decay. However, in general (as for the case at hand), the singular values cannot be determined analytically. However, borrowing from Optics literature, the NDF can be estimated by counting how many resolvable pointspread functions are required to fill the source domain [14], [37]. By doing so, and by working in the ξ domain, which where ξ x (−X s ) = −ξ x (X s ) and ξ y (−Y s ) = −ξ y (Y s ) have been exploited, we readily obtain In order to check the NDF estimation returned by (28), we turn again to consider the singular value behaviours shown in Fig. 3. In this figure the σ 2 m are shown for different configuration parameters. For these examples, the NDF returned by (28) are 60 and 324, respectively, and they adequately approximate the index beyond which the singular values start to decay very quickly. As seen in section III, when X 0 , Y 0 → ∞, the pointspread function turns to be ∝ J 1 (k|r−r |)/|r−r | whose main lobe half-width is 3,8 k . Note that the latter is a factor 3,8 π = 1, 2 larger than one obtained by our estimation. Accordingly, this entails that (28) leads to an overestimation of the actual NDF. In fact, while the latter is 0, 68 [16XsYs In this section, we briefly compare our resolution estimation with some literature results. Since most of the results that we found refer to radar imaging, we advise the reader that we adapted those results to the case at hand. To this end, we consider radar imaging cross-range resolution pertaining to a monostatic configuration, which involves a similar integral operator. In particular, since even under a multi-frequency configuration, cross-range resolution is estimated in correspondence to a single (the highest or the average) frequency, it is sufficient to consider half those estimations (because in radar imaging k is replaced by 2k in order to account for the twoway propagation path) to obtain the benchmark.
Suppose to simplify (26) and (27) by arresting the Taylor series expansions of ξ x (x + ∆ x ) and ξ y (y + ∆ y ) at the first term. Then ∆ x and ∆ y are approximated as with with Clearly, resolutions are still spatially varying and result dependent on the angular sector under which the source point is "seen" by the observation domain. In particular, the spatially varying bandwidth is basically approximated by dξx dx × dξy dy . It is important to highlight that eqs. (29) and (30) could have been obtained by resorting stationary phase arguments as done in [22], [38]. However, in those papers the spatially varying behavior is lost because resolution was estimated in correspondence to the bandwidth obtained by the union of the bands [22] or the averaged band [38] as the source point ranges within the source domain. Eventually, we can conclude that those results coincide with a first order approximation of our estimation when the spatially varying behavior was not neglected.
In other papers, such as [22], [23] and [31], the following resolution formula is derived and with θ rx (0) and θ ry (0), being the angular sectors under which the centre of the source domain is "viewed". In particular, in [23], θ rx (0) and θ ry (0) are supposed small enough so that tanθ rx (0) ≈ sinθ rx (0) and tan(θ ry (0) ≈ sinθ ry (0). These formulas have shown to provide good resolution estimation for example in Fresnel zone. However, in near-field, they return poor results because they approximate (29) and (30), which, as observed above, are in turn first order approximation of our estimation. Summarizing, the literature resolution formulas considered above are approximated version of (26) and (27) when a first order approximation is applied and the Fresnel zone condition is forced.

VI. CONCLUSION
The aim of this paper has been to provide an analytical estimation of the achievable resolution in the reconstruction of planar sources from near-field data collected over a bounded planar measurement surface. This has been pursued by finding an analytical approximation of the point-spread function returned by a back-propagation inversion scheme. The obtained expressions for the resolution ∆ x and ∆ y are spatially variant. This is a distinctive feature of near-zone aspect-limited configuration that has been already observed in some recent papers [17], [18] addressing the same problem for strip currents. The obtained estimations have been numerically checked and shown to be in excellent agreement with the outcome yielded by a truncated singular value decomposition scheme, used here as benchmark. A comparative discussion with some literature results addressing similar configurations has been addressed as well; it is shown that those resolution formulas are indeed approximated versions of (26) and (27) when a first order approximation is applied.
We end this paper by some further considerations. It is worthwhile remarking that the peculiar singular value behaviour (which as mentioned above justifies deriving resolution in term of configuration parameters) in certain cases can show some little dynamic before the knee (i.e, before the abrupt decay). This happens, for example, for electric currents [38]. In those cases, (3) and (2) show some differences. These, however, can be reduced by considering a suitable weighted back-propagation inversion [18]. Also, as mentioned above, while the shape of source domain does not play any role in the obtained the resolution estimation, the one of measurement surface affects the local band Ω(r, r ). In this regard, while the obtained results are specific to the considered observation domain, the proposed procedure is not and the same method can be used to analyze configurations exploiting different shape of measurement surfaces, including non planar ones.