Polarimetric calibration (Ainsworth)
This example synthesises fully-polarimetric SAR data with NUFFT projection, applies a known polarimetric distortion (channel imbalance + crosstalk), and recovers it with torchbp.polarimetry.ainsworth [1].
The scene is distributed clutter (many independent scatterers obeying reflection symmetry) plus a single trihedral corner reflector. The corner reflector is used to resolve HH/VV channel imbalance (k), clutter correlation is used to estimate other parameters.
[1] T. L. Ainsworth, L. Ferro-Famil and Jong-Sen Lee, “Orientation angle preserving a posteriori polarimetric SAR calibration,” in IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 4, pp. 994-1003, April 2006.
[1]:
import torch
import numpy as np
import matplotlib.pyplot as plt
from numpy import hamming
import torchbp
import torchbp.polarimetry as pol
[2]:
device = "cuda" if torch.cuda.is_available() else "cpu"
print("Device:", device)
torch.manual_seed(0);
Device: cpu
Radar and imaging geometry
The platform flies along y at constant altitude looking in x.
[3]:
fc = 6e9 # RF center frequency (Hz)
bw = 200e6 # Sweep bandwidth (Hz)
tsweep = 100e-6 # Sweep length (s)
fs = 20e6 # Sampling frequency (Hz)
nsamples = int(fs * tsweep)
nsweeps = 256 # Azimuth measurements
wl = 3e8 / fc
altitude = 100.0
element_spacing = 0.25 * wl
# Projection (scene) grid and imaging (polar) grid
grid_proj = {"x": (60.0, 260.0), "y": (-90.0, 90.0), "nx": 256, "ny": 256}
r0, r1, theta_limit = 90.0, 260.0, 0.35
res = 3e8 / (2 * bw)
nr = int(1.3 * (r1 - r0) / res)
ntheta = int(1 + nsweeps * 1.3 * (element_spacing / wl) * theta_limit / 0.25)
grid_polar = {"r": (r0, r1), "theta": (-theta_limit, theta_limit), "nr": nr, "ntheta": ntheta}
pos = torch.zeros(nsweeps, 3, dtype=torch.float32, device=device)
pos[:, 1] = torch.linspace(-nsweeps / 2, nsweeps / 2, nsweeps) * element_spacing
pos[:, 2] = altitude
print("nr", nr, "ntheta", ntheta, "nsamples", nsamples)
nr 294 ntheta 117 nsamples 2000
True scene
Ainsworth’s only assumption on the scene is scattering reciprocity and reflection symmetry of the distributed clutter, the cross-pol return is uncorrelated with the co-pol returns, \(\langle S_{hh}S_{hv}^*\rangle = \langle S_{vv}S_{hv}^*\rangle = 0\), and \(S_{vh} = S_{hv}\). We build the four channels in [HH, HV, VH, VV] order as independent speckle, with the cross-pol channel drawn independently of the co-pol channels.
A single bright trihedral corner reflector (\(S_{hh}=S_{vv}\), \(S_{hv}=S_{vh}=0\)) is placed at the scene centre to later pin down the HH/VV imbalance.
[4]:
nx, ny = grid_proj['nx'], grid_proj['ny']
def cn(scale):
return scale * (torch.randn(nx, ny, device=device)
+ 1j * torch.randn(nx, ny, device=device)) / np.sqrt(2)
Shh, Svv = cn(1.0), cn(1.0) # co-pol (independent)
Shv = cn(0.4) # cross-pol, independent of co-pol
Svh = Shv.clone() # reciprocity
cx, cy = nx // 2, ny // 2 # trihedral corner reflector
Shh[cx, cy] = Svv[cx, cy] = 200.0
Shv[cx, cy] = Svh[cx, cy] = 0.0
S = torch.stack([Shh, Shv, Svh, Svv]) # [4, nx, ny], order HH,HV,VH,VV
order = ['HH', 'HV', 'VH', 'VV']
Polarimetric distortion
The system distortion mixes the four channels. We use the same parameterisation that ainsworth inverts: channel imbalances alpha (cross-pol), k (HH/VV), and crosstalk u, v, w, z.
Ainsworth is orientation-angle preserving. It recovers only the crosstalk differences \(u-z\) and \(v-w\), because the common-mode part \((u+z)/2\), \((v+w)/2\) is mathematically indistinguishable from a real scene orientation angle. We therefore choose z = -u, w = -v (zero common-mode) so the distortion is fully recoverable.
[5]:
alpha = 1.15 * np.exp(1j * 0.20) # cross-pol channel imbalance
k = 0.85 * np.exp(-1j * 0.35) # co-pol HH/VV imbalance
u = 0.08 + 0.03j
v = -0.05 + 0.06j
w = -v # zero common-mode crosstalk
z = -u
M_true = torchbp.polarimetry.distortion_matrix(alpha, k, u, v, w, z, pol_order=order, device=device)
# Observed (distorted) scene, per channel
O = (M_true @ S.reshape(4, -1)).reshape(4, nx, ny)
Forward project and focus each channel
Each channel is forward projected to raw FMCW data with projection_cart_2d_nufft, range compressed, and backprojected.
[6]:
oversample = 2
data_fmod = -torch.pi * (1 - (oversample - 1) / oversample)
n = nsamples * oversample
r_res = 3e8 / (2 * bw * oversample)
win = torch.tensor(hamming(nsamples)[None, :], dtype=torch.float32, device=device)
fmod_f = torch.exp(1j * data_fmod * torch.arange(n, device=device))[None, :]
def make_image(scene):
data = torchbp.ops.projection_cart_2d_nufft(
scene, pos, grid_proj, fc, fs, bw / tsweep, nsamples,
use_rvp=False, normalization='gamma')[0]
data = torch.fft.ifft(data * win, dim=-1, n=n) * fmod_f
return torchbp.ops.backprojection_polar_2d(
data, grid_polar, fc, r_res, pos, dealias=False, data_fmod=data_fmod)[0]
sar = torch.stack([make_image(O[c]) for c in range(4)]) # [4, nr, ntheta]
print('Image shape:', tuple(sar.shape))
Image shape: (4, 294, 117)
Visualize the measured SAR images
Corner reflector should not reflect cross-polarization, but it is slightly visible also in HV and VH images due to cross-talk.
[7]:
extent = [*grid_polar['r'], *grid_polar['theta']]
fig, axes = plt.subplots(1, 4, figsize=(16, 3.2))
m = None
for ax, c, name in zip(axes, range(4), order):
db = 20 * torch.log10(sar[c].abs() + 1e-9).cpu().numpy()
if m is None:
m = db.max()
ax.imshow(db.T, origin='lower', vmin=m - 60, vmax=m, extent=extent, aspect='auto')
ax.set_title(name); ax.set_xlabel('Range (m)')
axes[0].set_ylabel('Angle (sin radians)')
plt.tight_layout();
Solve k from the corner reflector
The brightest pixel is the corner reflector. Its measured HH/VV ratio fixes the otherwise-ambiguous co-pol imbalance k. We also mask a small window around it with the weight argument so the corner does not bias the distributed-clutter crosstalk statistics.
[8]:
mag = sar[0].abs()
pk = torch.argmax(mag)
pr, pt = int(pk // mag.shape[1]), int(pk % mag.shape[1])
print(f'corner at ({pr},{pt}), peak/median = {(mag[pr,pt]/mag.median()).item():.1f}')
corner_hh_vv = (sar[0, pr, pt] / sar[3, pr, pt]).item()
print(f'measured corner HH/VV = {abs(corner_hh_vv):.3f} <{np.angle(corner_hh_vv):+.3f} rad')
weight = torch.ones(nr, ntheta, device=device)
ww = 6
weight[max(0, pr-ww):pr+ww+1, max(0, pt-ww):pt+ww+1] = 0.0
corner at (121,58), peak/median = 124.0
measured corner HH/VV = 0.964 <-0.290 rad
Calibration
ainsworth returns the normalised correction matrix Minv. A perfect calibration makes Minv @ M_true a scalar multiple of the identity.
Note that there is a possible 180 degree phase ambiguity with corner reflector k solving due to complex square root. A guess for k can also be specified and then the closer solution is used. k can also be specified directly if known and then corner_hh_vv should be None to disable corner reflector k solving.
[9]:
Minv = pol.ainsworth(sar, weight=weight, pol_order=order, corner_hh_vv=corner_hh_vv)
R = Minv @ M_true
R = R / R[0, 0]
eye = torch.eye(4, dtype=R.dtype, device=device)
offdiag = (R * (1 - eye)).abs().max().item()
diagerr = (torch.diagonal(R) - 1).abs().max().item()
print('Minv @ M_true:')
for row in R.cpu().numpy():
print(' ' + ' '.join(f'{x.real:+.3f}{x.imag:+.3f}j' for x in row))
print(f'\nmax off-diagonal (residual crosstalk) = {offdiag:.4f}')
print(f'max diagonal error (residual imbalance) = {diagerr:.4f}')
Minv @ M_true:
+1.000+0.000j -0.000-0.000j +0.000+0.000j +0.000+0.000j
+0.000+0.000j +1.002+0.000j +0.000+0.000j +0.000+0.000j
-0.000-0.000j -0.000+0.000j +1.002+0.000j +0.000-0.000j
+0.000-0.000j -0.000-0.000j +0.000+0.000j +1.004+0.001j
max off-diagonal (residual crosstalk) = 0.0000
max diagonal error (residual imbalance) = 0.0037
Product is almost identity. The small residual error left on diagnal is from corner reflector measurement error.
Applying calibration
torchbp.polarimetry.apply_cal can be used to apply calibration to image. After calibration the corner reflector isn’t anymore visible in HV and VH images.
[10]:
sar_cal = torchbp.polarimetry.apply_cal(sar, Minv)
extent = [*grid_polar['r'], *grid_polar['theta']]
fig, axes = plt.subplots(1, 4, figsize=(16, 3.2))
for ax, c, name in zip(axes, range(4), order):
db = 20 * torch.log10(sar_cal[c].abs() + 1e-9).cpu().numpy()
ax.imshow(db.T, origin='lower', vmin=m - 60, vmax=m, extent=extent, aspect='auto')
ax.set_title(name); ax.set_xlabel('Range (m)')
axes[0].set_ylabel('Angle (sin radians)')
plt.tight_layout();
Orientation-angle and common-mode crosstalk
Ainsworth’s model allows the true scene to carry a non-zero \(\langle S_{co}S_{cross}^*\rangle\), the scene orientation angle, requiring only that it be equal for HV and VH (reciprocity). A common-mode crosstalk \((u+z)/2\), \((v+w)/2\) produces exactly the same observed covariance as such an orientation angle, so the two cannot be separated from distributed targets alone. Ainsworth deliberately leaves the common-mode in place. The check below shows that with a non-zero common-mode crosstalk the residual equals common-mode.
Common-mode crosstalk should also be calibrated for the best calibration quality, but Ainsworth can’t correct it. It can be done with for example corner reflector measurements.
[11]:
# Same distortion but with a non-zero common-mode crosstalk component.
u2, v2 = 0.08 + 0.03j, -0.05 + 0.06j
w2, z2 = 0.04 - 0.07j, 0.06 + 0.02j # u2+z2, v2+w2 != 0
M2 = torch.tensor([
[k*alpha, v2/alpha, w2*alpha, v2*w2/(k*alpha)],
[z2*k*alpha, 1/alpha, w2*z2*alpha, w2/(k*alpha)],
[u2*k*alpha, u2*v2/alpha, alpha, v2/(k*alpha)],
[u2*z2*k*alpha,u2/alpha, z2*alpha, 1/(k*alpha)],
], dtype=torch.complex64, device=device)
O2 = (M2 @ S.reshape(4, -1)).reshape(4, nx, ny)
sar2 = torch.stack([make_image(O2[c]) for c in range(4)])
chv2 = (sar2[0, pr, pt] / sar2[3, pr, pt]).item()
Minv2 = pol.ainsworth(sar2, weight=weight, pol_order=order, corner_hh_vv=chv2)
R2 = Minv2 @ M2; R2 = R2 / R2[0, 0]
resid = (R2 * (1 - eye)).abs().max().item()
print(f'common-mode (u+z)/2 = {abs((u2+z2)/2):.3f}, (v+w)/2 = {abs((v2+w2)/2):.3f}')
print(f'residual off-diagonal = {resid:.3f}')
common-mode (u+z)/2 = 0.074, (v+w)/2 = 0.007
residual off-diagonal = 0.074
Calibration with rotation
Apply rotation to antenna, calibrate, verify that polarimetric orientation angle is still correct after calibration and cross-talk is correctly calibrated. After rotation co- and cross-polarization channels correlate.
[12]:
def copol_xpol_corr(img, weight):
"""Normalized co-pol/cross-pol correlation
~0 for a reflection-symmetric scene"""
d = img[0] - img[3] # HH - VV
x = 0.5 * (img[1] + img[2]) # (HV + VH)/2
num = (weight * d * x.conj()).mean()
den = ((weight * d.abs()**2).mean() * (weight * x.abs()**2).mean()).sqrt()
return (num / den).item()
# Give the clutter a known orientation angle, then distort, image and calibrate.
theta_scene = np.radians(20.0)
S_rot = pol.pol_antenna_rotation(S, theta_scene, pol_order=order)
O_rot = (M_true @ S_rot.reshape(4, -1)).reshape(4, nx, ny)
sar_rot = torch.stack([make_image(O_rot[c]) for c in range(4)])
chv_rot = (sar_rot[0, pr, pt] / sar_rot[3, pr, pt]).item()
Minv_rot = pol.ainsworth(sar_rot, weight=weight, pol_order=order, corner_hh_vv=chv_rot)
cal_rot = pol.apply_cal(sar_rot, Minv_rot)
R_rot = Minv_rot @ M_true
R_rot = R_rot / R_rot[0, 0]
resid = (R_rot * (1 - eye)).abs().max().item()
deg = lambda im: np.degrees(torchbp.polarimetry.orientation_angle(im, weight, order).item())
print(f"scene orientation applied = {np.degrees(theta_scene):+.1f} deg\n")
print(f"{'':30}{'|corr(HH-VV, HV)|':>17}{'orientation':>14}")
for name, im in [("un-oriented, calibrated", sar_cal), ("+20 deg, calibrated", cal_rot), ("+20 deg, uncalibrated", sar_rot)]:
print(f" {name:28}{abs(copol_xpol_corr(im, weight)):>15.3f}{deg(im):>11.1f} deg")
print(f"\nresidual system crosstalk = {resid:.4f}")
scene orientation applied = +20.0 deg
|corr(HH-VV, HV)| orientation
un-oriented, calibrated 0.021 0.5 deg
+20 deg, calibrated 0.504 20.5 deg
+20 deg, uncalibrated 0.511 19.5 deg
residual system crosstalk = 0.0000