极化相干分解——Cameron分解

MATLAB代码

function target_type = CameronDecomposition(S)
%% 对极化散射矩阵进行Cameron分解，并返回分解类型
% 输入：
%   S               2*2复矩阵
% 输出： target_type       分解目标类型
%
%%
k = reshape(S.',[],1); % 将矩阵直序展开为列向量

% Pauli基
kA = 1/sqrt(2)*[1,0,0,1].';
kB = 1/sqrt(2)*[1,0,0,-1].';
kC = 1/sqrt(2)*[0,1,1,0].';
kD = 1/sqrt(2)*[0,-1,1,0].';


% 1. 互易性判断
Prec = eye(4) - kD*kD.'; % 互易性提取矩阵
krec = Prec*k;
theta_rec = acos(norm(krec)/norm(k)); % 互易性判断参数
if theta_rec > pi/4
%     disp('This target is a Non-reciprocal scatter.');
    target_type = 11;
    return
else
    % 2. 对称性判断
    
    beta = dot(k,kB);  gamma = dot(k,kC); % Pauli基分解参数
    chi = atan2((beta*conj(gamma)+conj(beta)*gamma),(beta*conj(beta) - gamma*conj(gamma)));
    theta = chi/2; % 最大化角度
    kp = cos(theta)*kB + sin(theta)*kC;
    ksym = dot(k,kA)*kA + dot(k,kp)*kp; % 对称部分
    
    tau_sym = acos(norm(dot(krec, ksym))/(norm(krec)*norm(ksym))); % 对称性判断参数
    if tau_sym > pi/8
        
        % 判断是否为螺旋体
        khl = 1/2*[1,1i,1i,-1].'; % 左螺旋体S矩阵
        khr = 1/2*[1,-1i,-1i,-1].'; % 右螺旋体S矩阵
        theta_thl = acos(abs(dot(krec,khl)/norm(krec))); % 目标与左螺旋体的相似性
        theta_thr = acos(abs(dot(krec,khr)/norm(krec))); % 目标与右螺旋体的相似性
        h_th = 0.05; % 判断为螺旋体的门限
        if theta_thl < h_th
%             disp('And this target is a left helix scatter.')
            target_type = 8;
            return
        elseif theta_thr < h_th
%             disp('And this target is a right helix scatter.')
            target_type = 9;
            return
            
        else
%             disp('This target is an asymmetric scatter.')
            target_type = 10;
            return
        end
        
    else
        % 判断与典型对称散射体的对应类型
        
        % 计算旋转角并去取向
        
        
        Ssym = reshape(ksym,[2,2]).'; % 对称散射体矩阵
        [~,Sdiag] = eig(Ssym);
        Sdiag = 1/sqrt(Sdiag(1,1).*conj(Sdiag(1,1)) + Sdiag(2,2).*conj(Sdiag(2,2)))*Sdiag;
        eta = Sdiag(1,1)*conj(Sdiag(2,2)) + Sdiag(2,2)*conj(Sdiag(1,1)) - 1;
        switch eta
            case 0
                psi = 0;
            case -2
                psi = -1/4*chi;
                if psi > pi/4
                    psi = psi-pi/2;
                elseif psi > -pi/4
                    psi = psi;
                else
                    psi = psi+pi/2;
                end
                
            otherwise
                psi = -1/4*chi;
                psi_before_pick = [psi, psi-pi/2, psi+pi/2];
                pic_order = (psi_before_pick > -pi/2) & (psi_before_pick <= pi/2);
                psi = psi_before_pick(find(pic_order==1));
                for l = 1:length(psi)
                    R = [cos(psi(l)),  -sin(psi(l));
                        sin(psi(l)),   cos(psi(l))];
                    Ssym = reshape(ksym,[2,2]).';
                    Sdiag = round(R*Ssym/R,10);
                    if (abs(Sdiag(1,1)) >= abs(Sdiag(2,2)))
                        psi = psi(l);
                        break
                    end
                end
                
        end
        
        R = [cos(psi),  -sin(psi);
            sin(psi),   cos(psi)];
        Sdiag = round(R*Ssym/R,10);
        kdiag = reshape(Sdiag.',[],1);
        
        
        z = [1, -1, 0, 0.5, -0.5, 1i]; % 分别为三面角、二面角、偶极子、圆柱体、窄二面角、四分之一波器件
        typical_symmetric_name = ["1.三面角";"2.二面角";"3.偶极子";"4.圆柱体";"5.窄二面角";"6.四分之一波器件"];
        k_typ = cell(6,1); theta_t = ones(6,1);
        typical_flag = 0;
        s_th = 0.05;
        for j = 1:6
            S_type = typicalSMatrix(z(j));
            k_typ{j} = reshape(S_type.',[],1);
            theta_t(j) = acos(abs(dot(kdiag,k_typ{j})/norm(kdiag)));
            if theta_t(j) < s_th
                typical_flag = 1;
%                 disp(['This target is a typical symmetric scatter and it is ...']);
%                 disp(typical_symmetric_name(j))
                target_type = j;
                break
            end
        end
        
        % 若目标与典型对称散射体都不符合
        if ~typical_flag
%             disp('This target is a symmetric scatter.')
            target_type = 7;
        end
        
    end
end

return



function St = typicalSMatrix(z)
%
kl = [1,0,0,z];
St = 1/norm(kl)*reshape(kl,[2,2]).';

参考资料

[1]Cameron W L , Leung L K . Feature motivated polarization scattering matrix decomposition[C]// IEEE International Radar Conference. IEEE, 2002.
[2]匡纲要, 陈强, 蒋咏梅. 极化合成孔径雷达基础理论及其应用[M]. 国防科技大学出版社, 2011.