Document Sample

```					               Interferometric
Imaging in CASA
An Introduction for ALMA reduction

Steven T. Myers
Socorro, NM
S. T. Myers                NAASC Charlottesville –08 Sep 2011   1
Imaging
in theory…

S. T. Myers    NAASC Charlottesville –08 Sep 2011   2
• e.g. The EVLA or ALMA
– a bunch (27 to 50) antennas connected together
– independent elements  Earth rotation synthesis

NAASC Charlottesville - 08 Sep 2011   3
Interferometer Baselines
•   Baseline vector B in “aperture plane”
– coherent signal applied to interferometer would produce plane-wave
interference “fringe” on sky with angular period l/B

baseline vector B

q=λ/B
interferometer naturally decomposes sky into plane waves!
NAASC Charlottesville - 08 Sep 2011           4
The Aperture Plane
•   Correlate wavefronts in plane of apertures (Fourier transform of sky)
–   dish optics sum aperture plane at focus
–   visibility is cross-correlation of wavefronts of the 2 apertures

visibility
contains range
of baselines
from
closest to
furthest parts of
apertures

each point on aperture
gets correlated with each                      autocorrelations
point on other aperture                         measure uv
spacings inside
D/l
interferometer cannot measure “zero-spacing” w/o autocorrelations
NAASC Charlottesville - 08 Sep 2011                   5
From sky to uv-plane
The uv-plane is the Fourier Transform of the tangent plane to the sky

aperture xcor
width 2D/λ

F-1
2.5o
baseline u = B/λ
F                       ℓ = 2|u| = 2|B|/λ

Sky Plane x = (x,y)                                       Fourier Plane u = (u,v)
~        ~        2 i vx p
V (u)   d v A(u  v) I ( v) e
2
n
NAASC Charlottesville - 08 Sep 2011                         6
Example: VLA observes Jupiter
• A 6cm VLA observation of Jupiter:

Fourier transform of
nearly symmetric
planetary disk

S. T. Myers              NAASC Charlottesville –08 Sep 2011   7
Reconstruction of the Sky
• Visibilities and the Sky
v = A F-1 s + n
– A known instrumental response, but is not invertible
– instrumental noise n is a random variable
• The issues:
– unknown random noise n
– convolution due to size of A in uv domain
– incomplete sampling of uv-plane by visibilities
• Maximum Likelihood - Optimal Map
mMLE = ( AT N-1 A)-1 AT N-1 v = R-1 d                  d=Hv
No inversion R singular (at best ill-conditioned) for fully sampled s
S. T. Myers                        NAASC Charlottesville –08 Sep 2011             8
The Dirty Map
• Grid onto sampled uv-plane
d = H v = H s + nd
– H should be close to HMLE, e.g.
H = AT N-1                   :        A~A
– AT should sample onto suitable grid in uv-plane
– reminder: need only be approximate for gridding
• Invert onto sky  “dirty image”
d = F d = R s + nd        R = F R F-1
– image is “dirty” as it contains artifacts
• convolution by “point spread function” (columns of R)
• multiplication by response function (diagonal of R)
• noise
S. T. Myers                    NAASC Charlottesville –08 Sep 2011       9
• The VLA is an example of a sparsely filled array
– uv-plane gaps are treated as zeroes, cause “sidelobes” in PSF
– many solutions for sky that fit data, “dirty image” is principal solution
– must use “deconvolution” techniques to “clean” image

snapshot uv-coverage       PSF “dirty beam”                 “dirty” image   “clean” image

Example: VLA 30s snapshot discovery data for gravitational lens CLASS B1608+656
(Myers et al. 1995, ApJL, 447, L5-L8)
S. T. Myers                          NAASC Charlottesville –08 Sep 2011                         10
Image, uv, and Data Spaces
• image plane  uv-plane  visibilities
– operators F , H , A handle these transformations
– not all operators have inverses (H and A do not)
• example: model image m
– first transform sky model to uv-plane
m = F-1 m
– then project onto the visibilities (data space)
vm = A m = A F-1 m
– form residual
dvm = v - vm = A ( s - m ) + n
– finding “best model” will involve minimizing this residual
S. T. Myers                    NAASC Charlottesville –08 Sep 2011        11
Classic Deconvolution
• uv-plane CLEAN algorithm (“major cycles”)
– iterate on residual images removing point models
– initial residual data, and model:                   dv0 = v   m0 = 0
– form dirty image: d0 = F H dv0
– locate peak and residual and put fraction f into model
dm1 = f M d0              mask M : 1 at max, else 0
– increment model: m1 = m0 + dm1
– form cumulative visibilities and residual
v1 = A m1 = A F-1 m1        dv1 = v – v1
– form new dirty residual image: d1 = F H dv1
– and repeat until final residual image df is noise-like
S. T. Myers                    NAASC Charlottesville –08 Sep 2011                  12
image-plane CLEAN
• CLEAN in the image plane (minor cycles)
– initial model:  m0 = 0
– form dirty image: d0 = F H dv0
– locate peak and residual and put fraction f into model
dm1 = f M d0               mask M : 1 at max, else 0
– increment model: m1 = m0 + dm1
– subtract from dirty image to make dirty residual image
d1 = d0 – R dm1 R ~ R (the PSF)
– and repeat until final residual image df is noise-like
– size of kernel R (Hogbom=full-size, Clark=quarter size)
• specified by psfmode in CASA
S. T. Myers                     NAASC Charlottesville –08 Sep 2011      13
CLEAN variations
• Cotton-Schwab (CS) CLEAN
– break into major (uv-plane) & minor (image-plane) cycles
• in CASA initiated by imagermode=‘csclean’
• some options ( ‘mosaic’, ‘mfs’ nterms>1) hardwired cs
– clean in minor cycles to a threshold where max residual
is some fraction of starting max residual
• in CASA will be cyclefactor x max psf sidelobe
– more major cycles = more accurate cleaning but slow
• poor PSF, simple image structure = lower cyclefactor (<1)
• good PSF, complex image structure = raise cyclefactor (>1)
– purpose: correct errors from gridding & minor cycles
• the transform back from model visibilities is as accurate as we
wish to make it

S. T. Myers                       NAASC Charlottesville –08 Sep 2011              14
CLEAN Example
• Jupiter 6cm – interactive cleaning in CASA
as imaging
proceeds, include
regions containing

breaks between
interactive cycles
are major cycles
S. T. Myers              NAASC Charlottesville –08 Sep 2011                        15
Gridding kernels in CLEAN
• This is a choice:  d = Hv
• Default kernel: H0 = N-1
– includes only noise (inverse variance) weighting
• Optimal (mosaic) imaging: H = AT N-1
– uses aperture function (A ~ A)
– will grid mosiacs onto same uv-plane
– can correct for known pointing errors
• Frequency synthesis kernel (MFS)
– add kernels corresponding to Taylor expansion terms
• G0k= 1   G1k= ln(nk/n0) for each channel k
S. T. Myers                       NAASC Charlottesville –08 Sep 2011   16
Mosaicing in the sky and uv planes

phase

S. T. Myers       NAASC Charlottesville –08 Sep 2011   17
Mosaicing kernel
• the usual equation (aperture and offset term)
V (u,v)         dx dy Ax  x , y  y Ix, ye
k              k
2iux x k vy y k 

  dx dy Bx, y,u,v Ix, y e                  2 i ux vy

– note that we assume “phased up” at each pointing xk !
• kernel B and its transform (Fourier shift theorem)
Bx, y,u,v   Ax  x k , y  y k e2 iux k vyk 

Bi u,v  Au  ui,v  vi e2iui xk vi yk 
˜          ˜

          – B = F-1 B F is the “mosaicing” kernel (A-projection)
• if offset xk is unknown, then this is a “pointing error”
• offset for R and L polarizations is the “squint” term

S. T. Myers                               NAASC Charlottesville –08 Sep 2011                                    18
Mosaicing considerations
• In CASA, imagermode=‘mosaic’
– ftmachine=‘mosaic’ uses A-projection kernel
• grids to single uv-plane, optimally weights fields, most efficient
• can be some aliasing, keep mosaic to inner quarter of imsize
– ftmachine=‘ft’ does standard image-plane shift+add
• takes much longer, not recommended
– mosaic uses csclean, watch cyclefactor
– uses approximate single PSF for entire mosaic
• will “correct” in major cycles
– uses POINTING table when present (FIELD for phases)
• Outputs
–   standard: .image, .model, .psf
–   .flux (PB plus weights plus extra PB from A-convolution)
–   .flux.pbcoverage (just the PB effects)
–   e.g. to correct to flux on-sky use .image/.flux

S. T. Myers                            NAASC Charlottesville –08 Sep 2011             19
MEM and CLEAN
• CLEAN
– algorithm: find peak in residual image; add fraction to
model; form new residual data & residual image; iterate
– performance: good on compact emission, difficult for
extended
• Maximum Entropy Method (MEM)
– algorithm: for pixel values p : maximize entropy -S p ln p ;
minimize c2(p) to encourage “smooth” extended emission
– performance: complicated, suppresses spiky emission, but
fast
• CLEAN and MEM use point (pixel) basis
– complete basis – unique representation of image
S. T. Myers                  NAASC Charlottesville –08 Sep 2011         20
Sparse Approximation Imaging
• Problem: find a model to represent the sky as
efficiently as possible, subject to the data
constraints and within the noise uncertainty,
possibly also subject to prior constraints.
– some problems (like ours) cannot be efficiently
reconstructed using orthonormal bases (like pixels or
Fourier modes)
– extensive literature on this!
– use non-orthogonal bases: multiscale (e.g. Gaussians)
– choose dictionary of model elements (atoms)
– efficiency: find a representation that uses the fewest number
of atoms
S. T. Myers                NAASC Charlottesville –08 Sep 2011      21
Example: MEM versus MS-CLEAN

Restored                  Residual              Error

Maximum
Entropy

MS
Clean

S. T. Myers              NAASC Charlottesville –08 Sep 2011           22
The Future of Multiscale Methods
• Algorithms
– mostly iterative, starting from a blank model
– “greedy” methods make locally optimal choices at each step
• MS-CLEAN is a greedy algorithm in this class!
– dictionary of points and Gaussians on different scales
– is essentially a “Matching Pursuit” (MP) algorithm (e.g.
Tropp 2004)
• Key Research Area in CASA
– new arrays are designed for high dynamic range & fidelity
– we need efficient, robust, and accurate multiscale methods
– integrate multiscale (MS) with multispectral (MFS)
S. T. Myers                 NAASC Charlottesville –08 Sep 2011      23
Calibration and Imaging
• Some effects corrupt the visibilities
– most are on a per-antenna basis, other per-baseline
– antenna based effects can be “self-calibrated” out
• The Measurement Equation (ME)
                    
  * *  *  *
Vij  J i si  J j s j  J i  J j si  s j
obs

ideal

– the Jones matrices J contain the corruptions to V

V obs
ij      M ij Bij Gij Dij E ij Pij Tij Fij V                  ideal
ij        Aij
– there are different corruption terms to the J
• gain G, pol leakage D, ionosphere F, parallactic angle P
• can be direction dependent, additive RFI & correlator errors
S. T. Myers                       NAASC Charlottesville –08 Sep 2011                          24
Calibration in Image Plane
• Calibration errors show up as artifacts in
image plane:

Before Correction                         After Correction
– given an approximate model for the image we can solve
for the errors  “self-calibration” through iteration
S. T. Myers                  NAASC Charlottesville –08 Sep 2011            25
Pointing Corrections
• Example of direction-dependent errors:
– VLA antennas have ~10’’ pointing residual
– affects high-dynamic range imaging (>104)
– also “squint” between R and L beams
• Work by Sanjay Bhatnagar (NRAO)
• Simulation of 1.4GHz EVLA observations
• Residual images
– Before correction: Peak 250Jy, RMS 15Jy
– After correction: Peak 5Jy, RMS 1Jy
• Can incorporate into standard self-cal
• Computational cost ok for now
• See EVLA Memos 100 & 84
– Implementing in CASA, testing underway
S. T. Myers                   NAASC Charlottesville –08 Sep 2011   26
Polarization: Vis to Stokes in XY basis
• for antennas with parallactic angles fi and fj
f
 v XX   cos( i  f j )     f
cos( i  f j )  sin(fi  f j ) i sin(fi  f j )  I 
 XY                                                                      
 v    sin(fi  f j ) sin(fi  f j )       f
cos( i  f j )         f
i cos( i  f j )  Q 
 vYX    sin(f  f )   sin(fi  f j ) cos( i  f j )  i cos( i  f j ) U 
f                  f
              i   j
 
 vYY   cos(  f )  cos(  f )  sin(f  f ) i sin(f  f )  V 
fi j           fi j
                                             i     j          i     j    
• and for identical parallactic angles f between antennas:

 v XX           I  Q cos 2f  U sin 2f                      Linear Feeds:
 XY                                                       linear polarization
v               Q sin 2f  U cos 2f  iV                in all hands, circular
 j
 vYX  fi f Q sin 2f  U cos 2f  iV 
f                                             only in cross-hands
                                         
 vYY            I  Q cos 2f  U sin 2f 
                                         

S. T. Myers                          NAASC Charlottesville –08 Sep 2011                            27
Polarization Leakage
• Primary on-axis effect is “leakage” of one
polarization into the measurement of the other (e.g.
X  Y) : X  X + diX Y
– but, direction dependence due to polarization beam!
• Calibrate out on-axis leakage and put direction
dependence in “beam”
– example: expand XY basis with on-axis leakage

ˆ XX  V XX  d X V YX  d *X V XY  d X d *X V YY
Vij     ij     i   ij      j   ij     i    j   ij

ˆ
VijXY  VijXY  diX Vij  d *Y VijXX  diX d Y Vij
YY                         YX
j                j

– remember XX,YY contain IQU and XY,YX contain QUV
– ideally the “d”s are ~1-2% but can be worse, should be stable

S. T. Myers                       NAASC Charlottesville –08 Sep 2011        28
Primary Beam: full field polarization
• VLA primary beams
– Beam squint due to off-axis
system
– Instrumental polarization off-
axis
– AZ-EL telescopes
• Instrumental polarization
patterns rotate on sky
with parallactic angle
– Limits polarization imaging
– Limits Stokes I dynamic range            Green contours: Stokes I 3dB, 6dB, black
contours: fractional polarization 1% and up,
(via second order terms)                 vectors: polarization position angle, raster:
– must implement during imaging            Stokes V

S. T. Myers                     NAASC Charlottesville –08 Sep 2011                              29
Simulations on a complex model
• VLA simulation of ~ 1 Jy
point sources + large           I                                 V
source with complex
polarization (“Hydra A”)
• Long integration with full
range of parallactic
angles
• equivalent to weak
1.4GHz source observed
with EVLA                       Q                                 U
• Antenna primary beam
model by W. Brisken
• See EVLA memo 62

S. T. Myers                   NAASC Charlottesville –08 Sep 2011       30
1-D Symmetric Beam
I                                        V

Dynamic Range
~200 using
symmetric
beam model     Q                                        U

• dynamic range limited by errors from
incorrect approximate primary beam
S. T. Myers             NAASC Charlottesville –08 Sep 2011       31
2-D Polarized Beam
I                                        V

Dynamic Range
~104 using
2-D beam
Q                                        U
model

• need to use accurate polarized beam to
reach high fidelity and dynamic range
S. T. Myers             NAASC Charlottesville –08 Sep 2011       32
CASA CLEAN Imaging Details
• MFS (mode mfs)
–   nterm=2 compute spectral index, 3 for curvature etc.
–   needed for bandwidths ~5% or more (S/N dependent)
–   tt0 average intensity, tt1 alpha*tt0, alpha images output
–   takes at least nterms longer (image size dependent)
• Multiscale (set non-zero multiscale list)
– scales are in units of pixels
– usually set to be multiples of synth. beam size
– e.g. for 3x oversampling of beam:
• multiscale = [0,3,6,12,24]
– takes at least nscales longer (x nterms?)
– can be tricky to get to work right

S. T. Myers                         NAASC Charlottesville –08 Sep 2011    33
CASA CLEAN logistics
• clean can be restarted from current state
– if imagename used before and files exist and same size
– will first recompute residuals from model
• user can input some starting conditions
– previous mask (e.g. from previous clean)
– previous modelimage
• key toggles
–   mode: mfs, channel, velocity, frequency
–   imagermode: ‘’, csclean, mosaic
–   psfmode: clark, hogbom
–   gridmode: advanced stuff, eventually a-projection
S. T. Myers                     NAASC Charlottesville –08 Sep 2011   34
Spectral Cube Considerations
• mode: channel, velocity, frequency
– sets the grid of planes in output image cube
– will apply doppler corrections on the fly
– can set frame info, restfreq, etc.
• data: taken in sky frequency (terrestrial) frame
• velocity: include doppler shifts from:
– Earth rotation: few km/s (diurnal)
– Earth orbit: 30 km/s (annual)
– Earth/Sun motion w.r.t. LSR (e.g. LSRK, LSRD)
– maybe galactic rotation to extragalactic frames
• these shifts are time dependent
– “Doppler Tracking/Setting” adjust sky freq while observing
– can shift and regrid data before imaging : cvel
S. T. Myers                       NAASC Charlottesville –08 Sep 2011             35
Continuum Subtraction
• If you have a viable continuum model
– e.g. via using mode mfs
– can use uvsub to remove from data before imaging
– BUT: must have really good accurate model (not likely)
• In practice, subtract “continuum” in data
– use task uvcontsub or uvcontsub2
• particularly if strong short-baseline emission
– these fit polynomial to each vis spectrum and subtract
• uvcontsub2 can cross spw boundaries
– need to know line-free channel ranges
• if uvcontsub need to fit in each spw
– “continuum” not generally imageable
• e.g. no closure !!!

S. T. Myers                           NAASC Charlottesville –08 Sep 2011   36
Odds and Ends
• Oversampling PSF
– at least 2.5x, I like 3-4x (whichever rounds well), some
use 5x, no strong impact (but can make imaging longer)
• Boxing
– helps guide clean where to pick out real emission,
important when psf sidelobes high/complex and/or when
image has complicated structure
– if initial calibration poor, shallow careful clean with boxing
will get better model for selfcal (don’t burn in artifacts)
– some cases using csclean (w/cyclefactor) you can “turn it
loose”, but be careful before doing this
– we are working on autoclean / autoboxing
S. T. Myers                     NAASC Charlottesville –08 Sep 2011           37
Flotsam and Jetsam
• What, no clean components?
–   CASA clean increments the .model image
–   does not keep track of each iteration component
–   does not report total cleaned flux until end
–   does not keep separate multiscale components
–   might later store cc’s as component lists, but not now
–   you can mask/edit .model image
–   you can supply initial “model” using modelimage
• this can come from anywhere! e.g. single dish
• Stopping clean
– when you get residual that are noiselike :) or a mess :(
– ideally you should get within 2x thermal
• if not, probably calibration errors, bad data, or dyn. range issues
S. T. Myers                         NAASC Charlottesville –08 Sep 2011                 38
Fin
• Final words
– the proof of the quality of the data is in the final imaging,
you generally cannot guess this from metrics on the
calibrators (particularly when selfcal is possible)
• Anything else?

• Documentation