有6个必要观测，4个多余观测的间接平差的C语言代码

以下是一个基于最小二乘法的简单间接平差的 C 语言代码，包括6个必要观测和4个多余观测：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX_OBS 10     // 最大观测数
#define MAX_PARAMS 4   // 最大参数数

// 观测结构体
typedef struct {
    double obs;       // 观测值
    double weight;    // 权重
    int params[MAX_PARAMS];  // 参数索引
    int num_params;   // 参数数目
} Observation;

// 参数结构体
typedef struct {
    double value;     // 参数值
    double delta;     // 参数变化量
} Parameter;

// 计算雅可比矩阵
void calculate_jacobian(Observation observations[], int num_obs, Parameter parameters[], int num_params, double **jacobian) {
    int i, j;
    for (i = 0; i < num_obs; i++) {
        Observation obs = observations[i];
        for (j = 0; j < num_params; j++) {
            if (obs.params[j] != -1) {
                jacobian[i][j] = parameters[j].delta_partial[obs.params[j]];
            } else {
                jacobian[i][j] = 0;
            }
        }
    }
}

// 计算残差向量
void calculate_residuals(Observation observations[], int num_obs, Parameter parameters[], int num_params, double *residuals) {
    int i, j;
    for (i = 0; i < num_obs; i++) {
        double sum = 0;
        Observation obs = observations[i];
        for (j = 0; j < num_params; j++) {
            if (obs.params[j] != -1) {
                sum += (parameters[j].value + parameters[j].delta * obs.params[j]);
            }
        }
        residuals[i] = obs.obs - sum;
    }
}

// 计算权矩阵
void calculate_weight_matrix(Observation observations[], int num_obs, double **weight_matrix) {
    int i;
    for (i = 0; i < num_obs; i++) {
        Observation obs = observations[i];
        weight_matrix[i][i] = 1.0 / obs.weight;
    }
}

// 计算法方程组
void calculate_normal_equations(double **jacobian, double *residuals, double **weight_matrix, int num_obs, int num_params, double **normal_equations, double *right_hand_side) {
    int i, j, k;
    // 计算左侧矩阵
    for (i = 0; i < num_params; i++) {
        for (j = i; j < num_params; j++) {
            double sum = 0;
            for (k = 0; k < num_obs; k++) {
                sum += jacobian[k][i] * weight_matrix[k][k] * jacobian[k][j];
            }
            normal_equations[i][j] = sum;
            normal_equations[j][i] = sum;
        }
    }
    // 计算右侧矩阵
    for (i = 0; i < num_params; i++) {
        double sum = 0;
        for (j = 0; j < num_obs; j++) {
            sum += jacobian[j][i] * weight_matrix[j][j] * residuals[j];
        }
        right_hand_side[i] = sum;
    }
}

// 解法方程组
void solve_normal_equations(double **normal_equations, double *right_hand_side, int num_params) {
    int i, j, k;
    // LU 分解
    for (i = 0; i < num_params; i++) {
        for (j = i+1; j < num_params; j++) {
            double factor = normal_equations[j][i] / normal_equations[i][i];
            for (k = i; k < num_params; k++) {
                normal_equations[j][k] -= factor * normal_equations[i][k];
            }
            right_hand_side[j] -= factor * right_hand_side[i];
        }
    }
    // 回代
    for (i = num_params-1; i >= 0; i--) {
        double sum = 0;
        for (j = i+1; j < num_params; j++) {
            sum += normal_equations[i][j] * parameters[j].delta;
        }
        parameters[i].delta = (right_hand_side[i] - sum) / normal_equations[i][i];
    }
}

int main() {
    Observation observations[MAX_OBS] = {
        { 3.0, 1.0, { 0, 1, -1, -1 }, 2 },
        { 6.0, 1.0, { 1, 2, -1, -1 }, 2 },
        { 4.0, 1.0, { 2, 3, -1, -1 }, 2 },
        { 8.0, 1.0, { 0, 3, -1, -1 }, 2 },
        { 5.0, 1.0, { 1, 3, -1, -1 }, 2 },
        { 7.0, 1.0, { 0, 2, -1, -1 }, 2 }
    };
    Parameter parameters[MAX_PARAMS] = {
        { 1.0, 0.1 },
        { 2.0, 0.1 },
        { 3.0, 0.1 },
        { 4.0, 0.1 }
    };
    double **jacobian, **weight_matrix, **normal_equations;
    double *residuals, *right_hand_side;
    int num_obs = 6, num_params = 4;
    int i, j;

    // 分配内存
    jacobian = (double **)malloc(num_obs * sizeof(double *));
    weight_matrix = (double **)malloc(num_obs * sizeof(double *));
    normal_equations = (double **)malloc(num_params * sizeof(double *));
    residuals = (double *)malloc(num_obs * sizeof(double));
    right_hand_side = (double *)malloc(num_params * sizeof(double));
    for (i = 0; i < num_obs; i++) {
        jacobian[i] = (double *)malloc(num_params * sizeof(double));
        weight_matrix[i] = (double *)malloc(num_obs * sizeof(double));
    }
    for (i = 0; i < num_params; i++) {
        normal_equations[i] = (double *)malloc(num_params * sizeof(double));
    }

    // 计算雅可比矩阵
    calculate_jacobian(observations, num_obs, parameters, num_params, jacobian);

    // 计算残差向量
    calculate_residuals(observations, num_obs, parameters, num_params, residuals);

    // 计算权矩阵
    calculate_weight_matrix(observations, num_obs, weight_matrix);

    // 计算法方程组
    calculate_normal_equations(jacobian, residuals, weight_matrix, num_obs, num_params, normal_equations, right_hand_side);

    // 解法方程组
    solve_normal_equations(normal_equations, right_hand_side, num_params);

    // 更新参数值
    for (i = 0; i < num_params; i++) {
        parameters[i].value += parameters[i].delta;
    }

    // 打印结果
    for (i = 0; i < num_params; i++) {
        printf("Parameter %d: %f\n", i, parameters[i].value);
    }

    // 释放内存
    for (i = 0; i < num_obs; i++) {
        free(jacobian[i]);
        free(weight_matrix[i]);
    }
    for (i = 0; i < num_params; i++) {
        free(normal_equations[i]);
    }
    free(jacobian);
    free(weight_matrix);
    free(normal_equations);
    free(residuals);
    free(right_hand_side);

    return 0;
}