The centered calibration approach involves converting a nonlinear calibration problem into a linear one using a centering approximation. This technique allows for the use of linear least squares to solve for magnetometer biases, scale factors, and non-orthogonality corrections. Here is a summary of the process and an example of how it can be implemented in MATLAB:

Summary of the Centered Calibration Approach

1. Measurement Model: The magnetometer measurements can be modeled as:

   $B_k = (I_{3x3} + D)^{-1} (O^T A_k H_k + b + \epsilon_k)$

   where:

   • B_k is the magnetometer measurement at time t_k.
   • H_k is the geomagnetic field in Earth-fixed coordinates.
   • A_k is the unknown attitude matrix of the magnetometer.
   • D is the matrix of scale factors and non-orthogonality corrections.
   • O is an orthogonal matrix.
   • b is the bias vector.
   • \epsilon_k is the measurement noise.

2. Attitude-Independent Observation: Using the norms of the body-measurement and geomagnetic-reference vectors, we obtain an attitude-independent observation:

   $V_k = \|B_k\|^2 - \|H_k\|^2$

   This observation is used to form a quadratic loss function that can be minimized.

3. Centering Approximation: To linearize the problem, a centering approximation is applied. Define weighted averages and centered variables:

   $L = \frac{1}{N} \sum_{k=1}^N L_k, \quad z = \frac{1}{N} \sum_{k=1}^N z_k$

   $\tilde{L}_k = L_k - L, \quad \tilde{z}_k = z_k - z$

   The optimal estimate \theta^* is then given by:

   $\theta^* = (P_{\theta\theta})^{-1} \left( \sum_{k=1}^N \tilde{L}_k \tilde{z}_k \right)$

4. Conversion to D and b: Once \theta^* is determined, convert to D and b using the equations: $D = UWU^T, \quad b = (I_{3x3} + D)^{-1}c$

   where $E = U V U^T$ and $c$ is derived from \theta.

MATLAB Code Implementation

Here is an example of how you might implement the centered calibration in MATLAB:

matlabCopy code
function [D, b] = centered_calibration(B, H)
    % B: Magnetometer measurements (Nx3)
    % H: Geomagnetic reference (Nx3)

    N = size(B, 1);

    % Compute norms
    B_norms = sqrt(sum(B.^2, 2));
    H_norms = sqrt(sum(H.^2, 2));

    % Attitude-independent observation
    V = B_norms.^2 - H_norms.^2;

    % Define L and z matrices
    L = zeros(N, 9); % Adjust size based on the number of terms in L_k
    for k = 1:N
        B_k = B(k, :);
        S_k = [B_k(1)^2, B_k(2)^2, B_k(3)^2, 2*B_k(1)*B_k(2), 2*B_k(1)*B_k(3), 2*B_k(2)*B_k(3)];
        L(k, :) = [2*B_k, -S_k];
    end
    z = V;

    % Compute weighted averages
    L_avg = mean(L, 1);
    z_avg = mean(z);

    % Centered variables
    L_tilde = L - L_avg;
    z_tilde = z - z_avg;

    % Linear least squares solution
    P_theta_theta = L_tilde' * L_tilde;
    theta_star = P_theta_theta \\ (L_tilde' * z_tilde);

    % Convert to D and b
    E = [theta_star(1), theta_star(4), theta_star(5);
         theta_star(4), theta_star(2), theta_star(6);
         theta_star(5), theta_star(6), theta_star(3)];
    [U, V] = eig(E);
    W = diag(sqrt(diag(V)));
    D = U * W * U';

    c = mean((eye(3) + D) * B', 2) - mean(H', 2);
    b = (eye(3) + D) \\ c;
end