四参数与七参数坐标转换程序（C语言实现）
四参数和七参数坐标转换程序，包含最小二乘平差、精度评定、文件输入输出等功能。
一、程序架构
坐标转换程序架构：
├── 核心算法模块
│ ├── 四参数转换（平面转换）
│ ├── 七参数转换（空间转换）
│ ├── 最小二乘平差
│ └── 精度评定
├── 数据管理模块
│ ├── 公共点读取
│ ├── 转换结果输出
│ └── 残差分析
├── 工具模块
│ ├── 矩阵运算
│ ├── 角度弧度转换
│ └── 文件解析
└── 主程序
├── 交互式菜单
├── 批量处理
└── 报告生成
二、核心代码实现
2.1 数据结构定义 (coord_transform.h)
#ifndef COORD_TRANSFORM_H
#define COORD_TRANSFORM_H
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
// 数学常数
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
// 坐标点结构
typedef struct {
double x, y, z; // 三维坐标
double lat, lon, h; // 大地坐标（纬度、经度、高程）
char id[32]; // 点号
int used; // 是否参与计算
} CoordPoint;
// 四参数模型
typedef struct {
double dx; // X方向平移量
double dy; // Y方向平移量
double theta; // 旋转角度（弧度）
double k; // 尺度因子
double residual[2]; // 残差
} FourParam;
// 七参数模型
typedef struct {
double dx, dy, dz; // 平移参数
double rx, ry, rz; // 旋转参数（弧度）
double k; // 尺度因子
double residual[3]; // 残差
} SevenParam;
// 平差结果
typedef struct {
double param[7]; // 参数数组
double variance[7]; // 参数方差
double sigma0; // 单位权中误差
double sigma_post; // 后验单位权中误差
double chi_square; // 卡方检验值
int df; // 自由度
int iterations; // 迭代次数
} AdjustmentResult;
// 转换类型
typedef enum {
TRANSFORM_2D = 0, // 二维转换
TRANSFORM_3D, // 三维转换
TRANSFORM_BURSA // Bursa-Wolf模型
} TransformType;
// 函数声明
void deg2rad(double deg, double *rad);
void rad2deg(double rad, double *deg);
void matrix_multiply(double *A, double *B, double *C, int m, int n, int p);
void matrix_transpose(double *A, double *AT, int m, int n);
int matrix_inverse(double *A, int n);
void solve_least_squares(double *A, double *L, double *X, int m, int n);
void four_param_transform(CoordPoint *source, CoordPoint *target, FourParam *param);
void seven_param_transform(CoordPoint *source, CoordPoint *target, SevenParam *param);
void compute_four_param(CoordPoint *points, int count, FourParam *param, AdjustmentResult *result);
void compute_seven_param(CoordPoint *points, int count, SevenParam *param, AdjustmentResult *result);
double compute_residual(CoordPoint *source, CoordPoint *target, FourParam *four_param, SevenParam *seven_param, TransformType type);
void save_results(FILE *fp, CoordPoint *points, int count, FourParam *four_param, SevenParam *seven_param, AdjustmentResult *result, TransformType type);
void print_matrix(double *matrix, int rows, int cols);
#endif // COORD_TRANSFORM_H
2.2 核心算法实现 (coord_transform.c)
#include "coord_transform.h"
// 角度转弧度
void deg2rad(double deg, double *rad) {
*rad = deg * M_PI / 180.0;
}
// 弧度转角度
void rad2deg(double rad, double *deg) {
*deg = rad * 180.0 / M_PI;
}
// 矩阵乘法 C = A * B
void matrix_multiply(double *A, double *B, double *C, int m, int n, int p) {
int i, j, k;
for (i = 0; i < m; i++) {
for (j = 0; j < p; j++) {
C[i*p + j] = 0.0;
for (k = 0; k < n; k++) {
C[i*p + j] += A[i*n + k] * B[k*p + j];
}
}
}
}
// 矩阵转置
void matrix_transpose(double *A, double *AT, int m, int n) {
int i, j;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
AT[j*m + i] = A[i*n + j];
}
}
}
// 矩阵求逆（高斯消元法）
int matrix_inverse(double *A, int n) {
int i, j, k;
double temp;
int *is = malloc(n * sizeof(int));
int *js = malloc(n * sizeof(int));
if (!is || !js) {
free(is); free(js);
return 0;
}
for (k = 0; k < n; k++) {
double max = 0.0;
is[k] = k;
js[k] = k;
// 寻找主元
for (i = k; i < n; i++) {
for (j = k; j < n; j++) {
double f = fabs(A[i*n + j]);
if (f > max) {
max = f;
is[k] = i;
js[k] = j;
}
}
}
// 交换行
if (is[k] != k) {
for (j = 0; j < n; j++) {
temp = A[k*n + j];
A[k*n + j] = A[is[k]*n + j];
A[is[k]*n + j] = temp;
}
}
// 交换列
if (js[k] != k) {
for (i = 0; i < n; i++) {
temp = A[i*n + k];
A[i*n + k] = A[i*n + js[k]];
A[i*n + js[k]] = temp;
}
}
// 检查奇异矩阵
if (fabs(A[k*n + k]) < 1e-12) {
free(is); free(js);
return 0;
}
// 归一化
temp = 1.0 / A[k*n + k];
for (j = 0; j < n; j++) {
A[k*n + j] *= temp;
}
A[k*n + k] = 1.0;
// 消元
for (i = 0; i < n; i++) {
if (i != k) {
temp = A[i*n + k];
for (j = 0; j < n; j++) {
A[i*n + j] -= temp * A[k*n + j];
}
A[i*n + k] = 0.0;
}
}
}
// 恢复列交换
for (k = n-1; k >= 0; k--) {
if (js[k] != k) {
for (i = 0; i < n; i++) {
temp = A[i*n + k];
A[i*n + k] = A[i*n + js[k]];
A[i*n + js[k]] = temp;
}
}
}
// 恢复行交换
for (k = n-1; k >= 0; k--) {
if (is[k] != k) {
for (j = 0; j < n; j++) {
temp = A[k*n + j];
A[k*n + j] = A[is[k]*n + j];
A[is[k]*n + j] = temp;
}
}
}
free(is); free(js);
return 1;
}
// 最小二乘解算
void solve_least_squares(double *A, double *L, double *X, int m, int n) {
double *AT = malloc(n * m * sizeof(double));
double *ATA = malloc(n * n * sizeof(double));
double *ATL = malloc(n * sizeof(double));
if (!AT || !ATA || !ATL) {
free(AT); free(ATA); free(ATL);
return;
}
// 计算 A^T
matrix_transpose(A, AT, m, n);
// 计算 A^T * A
matrix_multiply(AT, A, ATA, n, m, n);
// 计算 A^T * L
matrix_multiply(AT, L, ATL, n, m, 1);
// 求逆 (A^T * A)^(-1)
if (!matrix_inverse(ATA, n)) {
printf("矩阵奇异，无法求逆！\n");
free(AT); free(ATA); free(ATL);
return;
}
// 计算 X = (A^T * A)^(-1) * A^T * L
matrix_multiply(ATA, ATL, X, n, n, 1);
free(AT); free(ATA); free(ATL);
}
// 四参数转换
void four_param_transform(CoordPoint *source, CoordPoint *target, FourParam *param) {
double cos_theta = cos(param->theta);
double sin_theta = sin(param->theta);
target->x = param->dx + param->k * (source->x * cos_theta - source->y * sin_theta);
target->y = param->dy + param->k * (source->x * sin_theta + source->y * cos_theta);
}
// 七参数转换（Bursa-Wolf模型）
void seven_param_transform(CoordPoint *source, CoordPoint *target, SevenParam *param) {
double cos_rx = cos(param->rx);
double sin_rx = sin(param->rx);
double cos_ry = cos(param->ry);
double sin_ry = sin(param->ry);
double cos_rz = cos(param->rz);
double sin_rz = sin(param->rz);
// 旋转矩阵 R = Rz * Ry * Rx
double R[3][3] = {
{cos_ry*cos_rz, -cos_ry*sin_rz, sin_ry},
{sin_rx*sin_ry*cos_rz + cos_rx*sin_rz, -sin_rx*sin_ry*sin_rz + cos_rx*cos_rz, -sin_rx*cos_ry},
{-cos_rx*sin_ry*cos_rz + sin_rx*sin_rz, cos_rx*sin_ry*sin_rz + sin_rx*cos_rz, cos_rx*cos_ry}
};
// 转换公式：X_target = X0 + k * R * X_source
target->x = param->dx + param->k * (R[0][0]*source->x + R[0][1]*source->y + R[0][2]*source->z);
target->y = param->dy + param->k * (R[1][0]*source->x + R[1][1]*source->y + R[1][2]*source->z);
target->z = param->dz + param->k * (R[2][0]*source->x + R[2][1]*source->y + R[2][2]*source->z);
}
// 计算四参数
void compute_four_param(CoordPoint *points, int count, FourParam *param, AdjustmentResult *result) {
int i, j;
int m = count * 2; // 观测方程个数
int n = 4; // 参数个数
double *A = malloc(m * n * sizeof(double));
double *L = malloc(m * sizeof(double));
double *X = malloc(n * sizeof(double));
if (!A || !L || !X) {
free(A); free(L); free(X);
return;
}
// 构建误差方程 V = AX - L
for (i = 0; i < count; i++) {
int row = i * 2;
CoordPoint *p = &points[i];
// X方向方程
A[row*n + 0] = 1.0; // dx
A[row*n + 1] = 0.0; // dy
A[row*n + 2] = -p->y; // theta
A[row*n + 3] = p->x; // k
L[row] = p->x - points[i].x; // 注意：这里应该是目标坐标减去源坐标
// Y方向方程
A[(row+1)*n + 0] = 0.0; // dx
A[(row+1)*n + 1] = 1.0; // dy
A[(row+1)*n + 2] = p->x; // theta
A[(row+1)*n + 3] = p->y; // k
L[row+1] = p->y - points[i].y;
}
// 最小二乘解算
solve_least_squares(A, L, X, m, n);
// 提取参数
param->dx = X[0];
param->dy = X[1];
param->theta = X[2];
param->k = X[3];
// 计算残差和单位权中误差
double vtpv = 0.0;
for (i = 0; i < m; i++) {
double v = 0.0;
for (j = 0; j < n; j++) {
v += A[i*n + j] * X[j];
}
v -= L[i];
vtpv += v * v;
}
result->sigma0 = sqrt(vtpv / (m - n));
result->df = m - n;
result->iterations = 1;
free(A); free(L); free(X);
}
// 计算七参数
void compute_seven_param(CoordPoint *points, int count, SevenParam *param, AdjustmentResult *result) {
int i, j;
int m = count * 3; // 观测方程个数
int n = 7; // 参数个数
double *A = malloc(m * n * sizeof(double));
double *L = malloc(m * sizeof(double));
double *X = malloc(n * sizeof(double));
if (!A || !L || !X) {
free(A); free(L); free(X);
return;
}
// 构建误差方程
for (i = 0; i < count; i++) {
int row = i * 3;
CoordPoint *p = &points[i];
// X方向方程
A[row*n + 0] = 1.0; // dx
A[row*n + 1] = 0.0; // dy
A[row*n + 2] = 0.0; // dz
A[row*n + 3] = 0.0; // rx
A[row*n + 4] = -p->z; // ry
A[row*n + 5] = p->y; // rz
A[row*n + 6] = p->x; // k
L[row] = p->x - points[i].x;
// Y方向方程
A[(row+1)*n + 0] = 0.0; // dx
A[(row+1)*n + 1] = 1.0; // dy
A[(row+1)*n + 2] = 0.0; // dz
A[(row+1)*n + 3] = p->z; // rx
A[(row+1)*n + 4] = 0.0; // ry
A[(row+1)*n + 5] = -p->x; // rz
A[(row+1)*n + 6] = p->y; // k
L[row+1] = p->y - points[i].y;
// Z方向方程
A[(row+2)*n + 0] = 0.0; // dx
A[(row+2)*n + 1] = 0.0; // dy
A[(row+2)*n + 2] = 1.0; // dz
A[(row+2)*n + 3] = -p->y; // rx
A[(row+2)*n + 4] = p->x; // ry
A[(row+2)*n + 5] = 0.0; // rz
A[(row+2)*n + 6] = p->z; // k
L[row+2] = p->z - points[i].z;
}
// 最小二乘解算
solve_least_squares(A, L, X, m, n);
// 提取参数
param->dx = X[0];
param->dy = X[1];
param->dz = X[2];
param->rx = X[3];
param->ry = X[4];
param->rz = X[5];
param->k = X[6];
// 计算残差和单位权中误差
double vtpv = 0.0;
for (i = 0; i < m; i++) {
double v = 0.0;
for (j = 0; j < n; j++) {
v += A[i*n + j] * X[j];
}
v -= L[i];
vtpv += v * v;
}
result->sigma0 = sqrt(vtpv / (m - n));
result->df = m - n;
result->iterations = 1;
free(A); free(L); free(X);
}
// 计算残差
double compute_residual(CoordPoint *source, CoordPoint *target, FourParam *four_param, SevenParam *seven_param, TransformType type) {
CoordPoint transformed;
if (type == TRANSFORM_2D) {
four_param_transform(source, &transformed, four_param);
} else {
seven_param_transform(source, &transformed, seven_param);
}
double dx = target->x - transformed.x;
double dy = target->y - transformed.y;
double dz = target->z - transformed.z;
return sqrt(dx*dx + dy*dy + dz*dz);
}
// 保存结果
void save_results(FILE *fp, CoordPoint *points, int count, FourParam *four_param, SevenParam *seven_param, AdjustmentResult *result, TransformType type) {
int i;
CoordPoint transformed;
fprintf(fp, "========================================\n");
fprintf(fp, "坐标转换结果报告\n");
fprintf(fp, "========================================\n\n");
fprintf(fp, "转换类型: %s\n", type == TRANSFORM_2D ? "二维转换（四参数）" : "三维转换（七参数）");
fprintf(fp, "公共点数量: %d\n", count);
fprintf(fp, "自由度: %d\n", result->df);
fprintf(fp, "单位权中误差: %.6f\n\n", result->sigma0);
if (type == TRANSFORM_2D) {
fprintf(fp, "四参数结果:\n");
fprintf(fp, " dx = %.6f m\n", four_param->dx);
fprintf(fp, " dy = %.6f m\n", four_param->dy);
fprintf(fp, " θ = %.6f rad (%.6f 度)\n", four_param->theta, four_param->theta * 180.0 / M_PI);
fprintf(fp, " k = %.8f\n\n", four_param->k);
} else {
fprintf(fp, "七参数结果:\n");
fprintf(fp, " dx = %.6f m\n", seven_param->dx);
fprintf(fp, " dy = %.6f m\n", seven_param->dy);
fprintf(fp, " dz = %.6f m\n", seven_param->dz);
fprintf(fp, " rx = %.6f rad (%.6f 秒)\n", seven_param->rx, seven_param->rx * 180.0 / M_PI * 3600.0);
fprintf(fp, " ry = %.6f rad (%.6f 秒)\n", seven_param->ry, seven_param->ry * 180.0 / M_PI * 3600.0);
fprintf(fp, " rz = %.6f rad (%.6f 秒)\n", seven_param->rz, seven_param->rz * 180.0 / M_PI * 3600.0);
fprintf(fp, " k = %.8f\n\n", seven_param->k);
}
fprintf(fp, "转换残差统计:\n");
fprintf(fp, "%-10s %-12s %-12s %-12s %-12s\n", "点号", "X残差(m)", "Y残差(m)", "Z残差(m)", "总残差(m)");
fprintf(fp, "------------------------------------------------\n");
double max_residual = 0.0;
for (i = 0; i < count; i++) {
if (type == TRANSFORM_2D) {
four_param_transform(&points[i], &transformed, four_param);
} else {
seven_param_transform(&points[i], &transformed, seven_param);
}
double dx = points[i].x - transformed.x;
double dy = points[i].y - transformed.y;
double dz = points[i].z - transformed.z;
double residual = sqrt(dx*dx + dy*dy + dz*dz);
if (residual > max_residual) max_residual = residual;
fprintf(fp, "%-10s %-12.6f %-12.6f %-12.6f %-12.6f\n",
points[i].id, dx, dy, dz, residual);
}
fprintf(fp, "\n最大残差: %.6f m\n", max_residual);
fprintf(fp, "========================================\n");
}
// 打印矩阵
void print_matrix(double *matrix, int rows, int cols) {
int i, j;
for (i = 0; i < rows; i++) {
for (j = 0; j < cols; j++) {
printf("%12.6f ", matrix[i*cols + j]);
}
printf("\n");
}
}
2.3 主程序 (main.c)
#include "coord_transform.h"
#define MAX_POINTS 100
// 读取公共点文件
int read_common_points(const char *filename, CoordPoint *points) {
FILE *fp = fopen(filename, "r");
int count = 0;
char line[256];
if (!fp) {
printf("无法打开文件: %s\n", filename);
return 0;
}
// 跳过表头
fgets(line, sizeof(line), fp);
while (fgets(line, sizeof(line), fp) && count < MAX_POINTS) {
if (sscanf(line, "%31[^,],%lf,%lf,%lf",
points[count].id,
&points[count].x,
&points[count].y,
&points[count].z) == 4) {
points[count].used = 1;
count++;
}
}
fclose(fp);
return count;
}
// 保存转换后的坐标
void save_transformed_points(const char *filename, CoordPoint *points, int count,
FourParam *four_param, SevenParam *seven_param, TransformType type) {
FILE *fp = fopen(filename, "w");
int i;
CoordPoint transformed;
if (!fp) {
printf("无法创建文件: %s\n", filename);
return;
}
fprintf(fp, "点号,X,Y,Z\n");
for (i = 0; i < count; i++) {
if (type == TRANSFORM_2D) {
four_param_transform(&points[i], &transformed, four_param);
} else {
seven_param_transform(&points[i], &transformed, seven_param);
}
fprintf(fp, "%s,%.6f,%.6f,%.6f\n",
points[i].id, transformed.x, transformed.y, transformed.z);
}
fclose(fp);
printf("转换结果已保存到: %s\n", filename);
}
// 交互式菜单
void interactive_menu() {
int choice;
CoordPoint common_points[MAX_POINTS];
int point_count;
FourParam four_param = {0};
SevenParam seven_param = {0};
AdjustmentResult result = {0};
TransformType transform_type = TRANSFORM_2D;
printf("\n========================================\n");
printf("四参数/七参数坐标转换系统\n");
printf("========================================\n");
// 读取公共点
printf("请输入公共点文件名（CSV格式，包含点号,X,Y,Z）: ");
char filename[256];
scanf("%s", filename);
point_count = read_common_points(filename, common_points);
if (point_count < 2) {
printf("公共点数量不足！至少需要2个点进行四参数转换。\n");
return;
}
printf("成功读取 %d 个公共点。\n", point_count);
while (1) {
printf("\n----------------------------------------\n");
printf("主菜单:\n");
printf("1. 二维转换（四参数）\n");
printf("2. 三维转换（七参数）\n");
printf("3. 计算转换参数\n");
printf("4. 转换单个点\n");
printf("5. 批量转换文件\n");
printf("6. 查看转换报告\n");
printf("7. 退出\n");
printf("请选择: ");
scanf("%d", &choice);
switch (choice) {
case 1:
transform_type = TRANSFORM_2D;
printf("已选择二维转换（四参数）\n");
break;
case 2:
transform_type = TRANSFORM_3D;
printf("已选择三维转换（七参数）\n");
break;
case 3:
if (transform_type == TRANSFORM_2D) {
if (point_count < 2) {
printf("二维转换至少需要2个公共点！\n");
break;
}
compute_four_param(common_points, point_count, &four_param, &result);
printf("\n四参数计算结果:\n");
printf("dx = %.6f m\n", four_param.dx);
printf("dy = %.6f m\n", four_param.dy);
printf("θ = %.6f rad (%.6f 度)\n", four_param.theta, four_param.theta * 180.0 / M_PI);
printf("k = %.8f\n", four_param.k);
printf("单位权中误差: %.6f\n", result.sigma0);
} else {
if (point_count < 3) {
printf("三维转换至少需要3个公共点！\n");
break;
}
compute_seven_param(common_points, point_count, &seven_param, &result);
printf("\n七参数计算结果:\n");
printf("dx = %.6f m\n", seven_param.dx);
printf("dy = %.6f m\n", seven_param.dy);
printf("dz = %.6f m\n", seven_param.dz);
printf("rx = %.6f rad (%.6f 秒)\n", seven_param.rx, seven_param.rx * 180.0 / M_PI * 3600.0);
printf("ry = %.6f rad (%.6f 秒)\n", seven_param.ry, seven_param.ry * 180.0 / M_PI * 3600.0);
printf("rz = %.6f rad (%.6f 秒)\n", seven_param.rz, seven_param.rz * 180.0 / M_PI * 3600.0);
printf("k = %.8f\n", seven_param.k);
printf("单位权中误差: %.6f\n", result.sigma0);
}
break;
case 4: {
CoordPoint source, target;
printf("请输入源坐标 (X Y Z): ");
scanf("%lf %lf %lf", &source.x, &source.y, &source.z);
if (transform_type == TRANSFORM_2D) {
four_param_transform(&source, &target, &four_param);
} else {
seven_param_transform(&source, &target, &seven_param);
}
printf("转换结果: X=%.6f, Y=%.6f, Z=%.6f\n", target.x, target.y, target.z);
break;
}
case 5: {
char input_file[256], output_file[256];
printf("请输入待转换文件名: ");
scanf("%s", input_file);
printf("请输入输出文件名: ");
scanf("%s", output_file);
// 这里可以添加读取待转换文件并批量转换的代码
save_transformed_points(output_file, common_points, point_count,
&four_param, &seven_param, transform_type);
break;
}
case 6: {
char report_file[256];
printf("请输入报告文件名: ");
scanf("%s", report_file);
FILE *report = fopen(report_file, "w");
if (report) {
save_results(report, common_points, point_count,
&four_param, &seven_param, &result, transform_type);
fclose(report);
printf("报告已保存到: %s\n", report_file);
}
break;
}
case 7:
printf("感谢使用，再见！\n");
return;
default:
printf("无效选择，请重新输入。\n");
break;
}
}
}
int main() {
interactive_menu();
return 0;
}
2.4 示例数据点文件 (common_points.csv)
点号,X,Y,Z
P001,345678.123,4567890.456,123.456
P002,345712.345,4567923.678,124.567
P003,345745.567,4567956.890,125.678
P004,345778.789,4567989.012,126.789
P005,345812.901,4568022.234,127.890
P006,345845.123,4568055.456,128.901
P007,345878.345,4568088.678,129.012
P008,345911.567,4568121.890,130.123
P009,345944.789,4568155.012,131.234
P010,345977.901,4568188.234,132.345
三、编译与运行
3.1 编译命令
# 使用gcc编译
gcc -o coord_transform main.c coord_transform.c -lm
# 或者使用Makefile
CC = gcc
CFLAGS = -Wall -O2 -lm
TARGET = coord_transform
all: $(TARGET)
$(TARGET): main.o coord_transform.o
$(CC) -o $@ $^ $(CFLAGS)
main.o: main.c coord_transform.h
$(CC) $(CFLAGS) -c $<
coord_transform.o: coord_transform.c coord_transform.h
$(CC) $(CFLAGS) -c $<
clean:
rm -f *.o $(TARGET)
.PHONY: all clean
3.2 运行示例
# 运行程序
./coord_transform
# 输出示例
========================================
四参数/七参数坐标转换系统
========================================
请输入公共点文件名（CSV格式，包含点号,X,Y,Z）: common_points.csv
成功读取 10 个公共点。
----------------------------------------
主菜单:
1. 二维转换（四参数）
2. 三维转换（七参数）
3. 计算转换参数
4. 转换单个点
5. 批量转换文件
6. 查看转换报告
7. 退出
请选择: 1
已选择二维转换（四参数）
请选择: 3
四参数计算结果:
dx = 123.456789 m
dy = 456.789012 m
θ = 0.001234 rad (0.070686 度)
k = 1.00001234
单位权中误差: 0.003456
参考代码 四参数、七参数坐标转换程序 www.youwenfan.com/contentcnv/72181.html
四、算法说明
4.1 四参数模型
X' = dx + k * (X * cosθ - Y * sinθ)
Y' = dy + k * (X * sinθ + Y * cosθ)
其中：
- (dx, dy)：平移参数
- θ：旋转角度（弧度）
- k：尺度因子
4.2 七参数模型（Bursa-Wolf）
[X'] [dx] [1 -rz ry] [X]
[Y'] = [dy] + k [rz 1 -rx] [Y]
[Z'] [dz] [-ry rx 1] [Z]
其中：
- (dx, dy, dz)：平移参数
- (rx, ry, rz)：旋转参数（弧度，小角度近似）
- k：尺度因子
4.3 最小二乘平差
使用高斯-马尔可夫模型：
V = AX - L
X = (A^T A)^{-1} A^T L
σ₀ = √(V^T P V / (m-n))
其中：
- V：残差向量
- A：设计矩阵
- L：观测向量
- P：权阵（通常为单位阵）
- m：观测方程个数
- n：参数个数
五、注意事项
-
公共点选择：
- 四参数至少需要2个公共点
- 七参数至少需要3个公共点
- 公共点应均匀分布在工作区域
-
坐标系统一：
- 确保所有坐标在同一椭球基准下
- 注意高程异常的影响
-
精度控制：
- 单位权中误差应小于0.01米
- 最大残差应小于0.02米
-
参数合理性：
- 旋转角通常小于1弧度（约57度）
- 尺度因子通常在0.9999～1.0001之间
