数字高程模型（DEM）内插

QuietForest_L

已于 2024-11-28 15:46:53 修改

阅读量716

点赞数 9

文章标签：算法

于 2024-11-22 20:41:20 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_74942631/article/details/143974992

版权

一、实习目的

掌握移动曲面法数字高程模型内插原理及其内插子程序的设计方法，了解其它逐点高程内插方法的基本原理。

二、实习内容

根据提供的10个数据点的坐标（Xn,Yn,Zn）和待求点的平面坐标（Xp,Yp），要求利用移动二次曲面拟合法，由格网点P（Xp,Yp）周围的10个已知点内插出待求格网点P的高程，编制相应的程序并进行调试，最后解算出格网点P的高程并提交源程序代码。

三、资料准备

已知数据点坐标

点号	X	Y	Z
1	102	110	15
2	109	113	18
3	105	115	19
4	103	103	17
5	108	105	21
6	105	108	15
7	115	104	20
8	118	108	15
9	116	113	17
10	113	118	22

编程计算点（110,110）上的高程。

四. 操作步骤

、读入已知点的坐标，建立以待定点为原点的局部坐标系；
、建立误差方程式：
组成法方程，解算六个系数：
整理计算结果，以实习报告的形式递交成果。

五. 代码过程

#define _CRT_SECURE_NO_WARNINGS 1
#include <stdio.h>
#include <math.h>


typedef struct
{
	double index,X,Y,Z;
}data;

data Date[10];//共有10个点
#pragma region InverseMatrix

int InverseMatrix(double B[6][6], double inverse[6][6]) {
	// 创建增广矩阵，将 B 扩展成 [B | I]
	double aug[6][2 * 6];
	for (int i = 0; i < 6; i++) {
		for (int j = 0; j < 6; j++) {
			aug[i][j] = B[i][j];         // 左边是矩阵 B
			aug[i][j + 6] = (i == j) ? 1.0 : 0.0;  // 右边是单位矩阵 I
		}
	}

	// 高斯-约当消元法
	for (int i = 0; i < 6; i++) {
		// 找到当前列中绝对值最大的元素并交换行
		double maxElement = fabs(aug[i][i]);
		int maxRow = i;
		for (int k = i + 1; k < 6; k++) {
			if (fabs(aug[k][i]) > maxElement) {
				maxElement = fabs(aug[k][i]);
				maxRow = k;
			}
		}

		// 如果无法求逆（矩阵不可逆），返回错误代码
		if (maxElement == 0.0) {
			printf("矩阵不可逆\n");
			return 0;
		}

		// 交换当前行与最大元素所在行
		if (maxRow != i) {
			for (int k = 0; k < 2 * 6; k++) {
				double temp = aug[i][k];
				aug[i][k] = aug[maxRow][k];
				aug[maxRow][k] = temp;
			}
		}

		// 归一化当前行，使对角线元素变为 1
		double diag = aug[i][i];
		for (int k = 0; k < 2 * 6; k++) {
			aug[i][k] /= diag;
		}

		// 消去其他行的当前列
		for (int j = 0; j < 6; j++) {
			if (j != i) {
				double factor = aug[j][i];
				for (int k = 0; k < 2 * 6; k++) {
					aug[j][k] -= factor * aug[i][k];
				}
			}
		}
	}

	// 提取右侧的单位矩阵部分作为逆矩阵
	for (int i = 0; i < 6; i++) {
		for (int j = 0; j < 6; j++) {
			inverse[i][j] = aug[i][j + 6];
		}
	}
	return 1;  // 成功返回 1
}

#pragma endregion

#pragma region Import data
int read_data(const char* filename) {
	FILE* file = fopen(filename, "r");
	if (!file) {
		printf("无法打开文件\n");
		return 1;
	}

	char buffer[100];
	fgets(buffer, sizeof(buffer), file); // 跳过标题行
	for (int i = 0; fscanf(file, "%lf,%lf,%lf,%lf", &Date[i].index, &Date[i].X, &Date[i].Y, &Date[i].Z) == 4; i++) {
		Date[i].X -= 110;
		Date[i].Y -= 110;
	}

	fclose(file);
	return 0;
}

#pragma endregion


double calculate_distance(double x, double y) {
	return sqrt(x * x + y * y);
}


int main()
{
	read_data("./内插点数据.csv");
	printf("****************导入点坐标*******************\n");
	for (int i = 0; i < 10; i++) {	 //输出读取的数据验证结果
		printf("点号: %lf, X: %lf, Y: %lf, Z: %lf\n",
			Date[i].index, Date[i].X, Date[i].Y, Date[i].Z);
	}

	double P[10]; // 存储每个点的权值
	for (int i = 0; i < 10; i++) {
		double distance = calculate_distance(Date[i].X, Date[i].Y);
		if (distance == 0) {
			printf("待定点与已知点 %d 重合，权值设置为最大\n", i);
			P[i] = 1e10; // 防止除以零，设为一个很大的值
		}
		else {
			P[i] = 1.0 / distance; // 权值为距离的倒数
		}
	}

	//建立误差方程&解算法方程
	double L[10];
	double result[6] = {0};
	double A[10][6];//误差方程式系数
	double AT[6][10];
	double A1[6][10];
	double ATA[6][6];
	double ATA1[6][6];//逆矩阵

	for (int i = 0; i < 10; i++)
	{
		L[i] = Date[i].Z * P[i];;
		A[i][0] = Date[i].X * Date[i].X * P[i];
		A[i][1] = Date[i].X * Date[i].Y * P[i];
		A[i][2] = Date[i].Y * Date[i].Y * P[i];
		A[i][3] = Date[i].X * P[i];
		A[i][4] = Date[i].Y * P[i]; 
		A[i][5] = 1 * P[i];
	}

	for (int i = 0; i < 10; i++)//转置
	{
		for (int j = 0; j < 6; j++)
		{
			AT[j][i] = A[i][j];
		}
	}
	
	for (int i = 0; i < 6; i++)
	{
		for (int m = 0; m < 6; m++)
		{
			ATA[i][m] = 0;
			for (int j = 0; j < 10; j++)
			{
				ATA[i][m] += AT[i][j] * A[j][m];
			}
		}
	}

	InverseMatrix(ATA, ATA1);

	for (int i = 0; i < 6; i++)
	{
		for (int m = 0; m < 10; m++)
		{
			A1[i][m] = 0;
			for (int j = 0; j < 6; j++)
			{
				A1[i][m] += ATA1[i][j] * AT[j][m];
			}
		}
	}

	for (int i = 0; i < 6; i++)
	{
		for (int j = 0; j < 10; j++)
		{
			result[i] += A1[i][j] * L[j];
		}
	}

// 计算拟合值并输出误差
	double fitZ[10]; // 用于存储拟合值
	double error[10]; // 用于存储误差值
	double mse = 0.0; // 初始化 MSE

	printf("\n****************拟合坐标*******************\n");
	for (int i = 0; i < 10; i++) {
		// 计算拟合值
		fitZ[i] = result[0] * Date[i].X * Date[i].X
			+ result[1] * Date[i].X * Date[i].Y
			+ result[2] * Date[i].Y * Date[i].Y
			+ result[3] * Date[i].X
			+ result[4] * Date[i].Y
			+ result[5];

		// 计算误差
		error[i] = Date[i].Z - fitZ[i];

		// 累加平方误差
		mse += error[i] * error[i];

		// 输出每个点的误差
		printf("点号: %lf, 实际高程: %lf, 拟合高程: %lf, 误差: %lf\n",
			Date[i].index, Date[i].Z, fitZ[i], error[i]);
	}

	// 计算均方误差
	mse /= 10.0; // 将总平方误差除以点的数量

	// 输出均方误差
	printf("均方误差 (MSE): %lf\n", mse);



	printf("\n六个参数为：%f %f %f %f %f %f \n", result[0], result[1], result[2], result[3], result[4], result[5]);

	printf("\nP的高程值 = %f ", result[5]);

	FILE* fp2;
	fp2 = fopen("result.txt", "w");
	fprintf(fp2, "****************数字高程模型内插*******************\n");
	fprintf(fp2, "\n");
	fprintf(fp2, " 六个参数最终结果:\n ");
	fprintf(fp2, "% f % f % f % f % f % f \n", result[0], result[1], result[2], result[3], result[4], result[5]);
	fprintf(fp2, "\n");
	fprintf(fp2, "P点的高程值：\n");
	fprintf(fp2, "%f\n", result[5]);

	return;
}

注意：上述代码的部分来源于网络，但大致流程如上，本代码里计算P时采用距离倒数权重的方式，并没有设置为E。如何是设置距离倒数为权重，需要注意如果是求已知点的高程，可能为NAN，但是本程序里面预防了这种情况，当距离为0时，直接将该点的高程权重拉为最大，即不考虑其他点，直接采用改点的高程值。

读取数据时采用的csv表格的形式，首先将提供的数据导入到csv中，然后再进行读取。