RSD高分卫星数据处理能力提升——日正射处理数千景高分数据集

news2024/9/27 19:23:12

李国春

通常认为，能够单日处理几百景高分辨率对地观测卫星数据的系统就已经是非常优秀的卫星数据处理系统了。RSD此次优化将其处理能力提升超过了一个数量级，达到了单日正射处理数千景高分辨率卫星数据集的水平。

不仅如此，RSD达到如此高的处理能力是在成本万元以内的消费级个人电脑上实现的，而且还是单CPU的独立机器。感兴趣的读者可以试验一下把这些数据集从一块硬盘拷贝到另一块硬盘所需要的时间，就能比较出RSD的处理速度多么恐怖。如果你有价值几万几十万元的高档工作站，处理速度会更快。

多搞几台这种消费级的计算单元，组建一个计算机集群，将这个脚本程序改成多机并行的，花很少的钱就可以挑战任何先进的对地观测卫星数据处理系统了。你要是感兴趣就联系我。

上一帖《一个高效的通用光学卫星数据正射校正程序》里面各种数据集正射都是几秒的时间，有同学反馈说我是不是把分钟错当成秒了。不是的，就是这个速度。RSD处理速度快主要得益于算法先进，其次是数据处理优化准确。RSD绝大多数处理模块都没有优化，这部分正射优化了是因为今年病了住院做不了其它推进性工作，我就整理下思路优化了一下。核心函数我把它改成 Sw7WristbandsBuffer，7Wristbands不是什么7个波段而是医院的腕带，用来纪念我2022年大难不死。:)

跑题了，言归正传。算法就不说了，我们自己独有的，和任何流行算法都不一样。在工程实现上也进行了精准的优化。顺便说一句，有人认为使用了高效的编程工具理所当然的就快，其实不然。工具不是决定因素，就像建筑工程的质量根本不取决于你使用了多么先进的装备一样。

下面进行详细介绍。

一、数据来源

本次试验使用的是高分2号卫星的多光谱和全色数据共1454景L1A级原始数据集。将这些数据集名称收集写入到一个文本文件，再输入到本试验脚本程序。脚本逐个读取文件名，在对应位置找到该数据集，进行处理。文本文件是通过运行代码 通配符查找文件.c 自动生成的。详见《搜索本地文件》。

运行通配符查找文件.c程序代码，在本地计算机一共搜索到（MSS和PAN各）1454个数据集。MSS和PAN各自结果分别保存在一个文本文件中。

结果见图1。

图1 本机搜索到的1455个GF2 MSS L1A文件

图1中的GF2.txt文件只列出了搜索到的TIFF文件，正射校正还需要一个RPC参数文件，通常RPC参数文件与TIFF文件有一样的主文件名，可以在程序里自动搜索，这里无需列出。

与MSS类似的还有一个PAN文件。需要重新搜索一次，并单独保存为一个文本文件，结果参见图1。

二、处理过程

在RSD的脚本窗口打开程序演示批量正射校正.c 。代码如下：

main(){
	STRING name;
	name = OpenFileDialog(TRUE);//name 返回待处理数据文件名的文本文件。参数TRUE 表示打开文件，FALSE 表示保存文件
	if(name == ""){
		Print("已退出主程序");
		return FALSE;
	}
	int nID1 = Open(name,TRUE);	//参数TRUE表示文本。
	STRING outPath;
	BrowseFolder(outPath);//选择一个结果文件的输出路径
	Print("开始生成中国区 SRTM 90m DEM 数据...");
	DWORD tc = GetTickCount();
	//创建 DEM 数据，正射校正用///
	double		demGeoBox[4];	//排列依次为：最小纬度、最小经度、最大纬度和最大经度
	longlong	demHeight = 30000;
	longlong	demWidth = 42000;
	short		demBuffer[demHeight][demWidth];
	double		dbOrgLat = 46.0;	//设置左上角位置的经纬度
	double		dbOrgLon = 95.0;
	short		defaltElv = 0.0;	
	double		demPixsize = 5.0/6000.0;	//SRTM90m = 0.0008333。ASTGTM2、3 DEM 为 0.00027778度/像元.demPixsize = 1.0/3600.0;
	DemSrtm90Buffer(demBuffer,dbOrgLat,dbOrgLon,demHeight,demWidth,defaltElv,demGeoBox);
	double dmtm = (GetTickCount() - tc) / 1000.0;
	Print("生成中国区 SRTM 90m DEM 数据用时 %f 秒。",dmtm);
	STRING mainName,outName,tifName,rpcName,outName,projStr;
	double box[4];	//地理位置4角坐标，排列依次为：最小N、最小E、最大N和最大E
	double pixsize = 0.8;//GF2 MSS 3.2m， GF2 PAN 0.8m
	longlong	width,height;
	int			bands,datatype,rank,boo,n;;
	double datasetTime,outpuTime,totalTime;
	int dims[16];
	int count = 0;
	DWORD bgntc = GetTickCount();
	do{
		//读来一行文件名，并查找对应的 RPC 文件///
		tifName = "";
		n = Read(nID1,tifName);
		if(!PathFileExists(tifName))
		{
			continue;
		}
		rpcName = tifName - '.';
		rpcName = rpcName + ".rpb";
		if(!PathFileExists(rpcName))
		{
			rpcName = tifName - '.';
			rpcName = rpcName + "_rpc.txt";
			if(!PathFileExists(rpcName))
			{
				continue;
			}
		}
		Print("开始对%s进行正射校正...",tifName);//开始进行正射校正
		tc = GetTickCount();
		projStr = NULL;//投影参数串，指定 projStr=NULL; 时，自动探测带号，WGS84大地基准
		Alloc(WORD retBuffer[1][1][1] = {0};); //哑数组，接收处理结果。用 Alloc为每个文件单独申请和删除（每个文件大小不同）
		int boo = Sw7WristbandsBuffer(tifName,rpcName,retBuffer,box,demBuffer,demGeoBox,demPixsize,projStr,TRUE,CC,pixsize);
		datasetTime = (GetTickCount() - tc) / 1000.0;
		//正射校正结果导出为TIFF数据格式文件。
		mainName = FileNameMain(tifName);
		outName = outPath + mainName + ".tif"; //输出文件(.tif），输出路径 + 输入文件的主文件名+ 扩展名 .TIFF
		tc = GetTickCount();
		InqVariable(retBuffer,rank,dims,datatype);//读正射校正后PAN数据的数据范围、维数和数据类型信息
		height = dims[0];
		width = dims[1];
		bands = dims[2];		GdSaveTiffFile(outName,retBuffer,height,width,bands,datatype,BIP,projStr,box,pixsize);
		Free(retBuffer);
		outpuTime = (GetTickCount() - tc) / 1000.0;

		//处理计数、计时
		totalTime = (GetTickCount() - bgntc) / 1000.0;
		count++;
		Print("正射校正用时 %f 秒，导出结果用时 %f 秒。共处理了 %d 个数据集，累积用时 %f 分钟。",datasetTime,outpuTime,count,totalTime/60);
	}while(n);
	Close(nID1);
	Print("主程序结束");
	return TRUE;
}