一个高效的通用光学卫星数据正射校正程序
李国春
随着高分辨率对地观测卫星发射的日益增多，对数据处理软件的要求也越来越高。通常每个系列卫星都有自己的数据特点并需要专门的处理软件，但卫星数量的增加为每种卫星单独设计软件的压力越来越大。本文介绍的一种处理方案旨在能够正射校正处理大多数的光学遥感卫星的观测数据。
将卫星运营部门发布的卫星数据处理成用户数据最重要的处理成本是正射校正部分。这里提供的处理方案可以快速无人工干预处理由有理式（RPC）参数描述的一级数据。
本文演示处理了几种卫星数据，更多类别的卫星（尤其是近年大量发射的新卫星）数据用户可以用本文程序试验处理，遇有问题可以及时反馈以便修正处理。
一、99脚本语言处理代码
1. 准备正射校正用DEM数据
Print("开始生成中国区 SRTM 90m DEM 数据...");
DWORD tc = GetTickCount();
//创建 DEM 数据，正射校正用///
double demGeoBox[4]; //排列依次为：最小纬度、最小经度、最大纬度和最大经度
longlong demHeight = 36000;
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度/像元.pixsize = 1.0/3600.0;
DemSrtm90Buffer(demBuffer,dbOrgLat,dbOrgLon,demHeight,demWidth,defaltElv,demGeoBox);
double dmtm = (GetTickCount() - tc) / 1000.0;
Print("生成中国区 SRTM 90m DEM 数据用时 %f 秒。",dmtm);
这段代码从SRTM90读来高程数据，存放在DemSrtm90Buffer指向的缓存。数据左上角纬度46.0°N，经度95°E。东西向42000个像元，南北向36000个像元。每个像元5/6000°。这样根据经纬度位置能够定位到对应的高程。
正射校正用DEM数据可以使用SRTM90、ASTGDEM2、3等通用DEM，也可以使用用户自行定义的DEM数据，只需指定开始经纬度位置、像元大小和数据范围，并将数据预存到一个数据缓存中即可。
2. 选择待校正的数据集
//选择待进行正射校正的卫星数据文件，自动查找对应的RPC数据文件 ///
longlong width,height;
int bands,datatype;
STRING tifName = OpenFileDialog(TRUE,"*.tif;*.tiff");//打开一个tif文件
STRING rpcName = tifName - '.';
rpcName = rpcName + ".rpb";
if(!PathFileExists(rpcName))
{
rpcName = tifName - '.';
rpcName = rpcName + "_rpc.txt";
if(!PathFileExists(rpcName))
{
Print("没有找到多光谱数据 %s 对应的 RPC 文件。",tifName);
return FALSE;
}
}
这段代码用OpenFileDialog函数通过对话框选择一个待处理的文件tifName，并通过文件名查找对应的有理式系数（RPC）文件rpcName。先查找.rpb文件，如果该文件不存在再查找.txt文件。
3. 选择结果输出目录并生成一个输出文件名
//选择一个结果文件的输出路径，并指定导出文件的类型（.img或者.tif）/
STRING outPath;
BrowseFolder(outPath);
STRING mainName = FileNameMain(tifName);
STRING outName = outPath + mainName + ".tif";//正射校正后数据的导出文件名
输出文件名指定.tif扩展名自动保存为TIFF文件，指定为.img时自动保存为带.hdr的IMG格式。
4. 开始计时
Print("开始对%s进行正射校正...",tifName);//开始进行正射校正/
DWORD tc = GetTickCount();
5. 获取数据行数、列数、波段数、数据类型和四角经纬度信息
//检测数据集基本信息///
longlong width,height;
int bands,datatype;
double quad[8]; //数据四角经纬度，每角一对，书写顺序
int b = TifInqDataInfo(tifName,height,width,bands,datatype); //读取 TIFF 文件的数据范围、数据类型和波段数
int bo2 = SwRpcGeoQuad(rpcName,quad,height,width,demBuffer,demHeight,demWidth,SINT16BL,demGeoBox,demPixSize);//四角坐标
6. 指定投影参数、四角经纬度转换成Northing、Easting坐标
STRING projstr = "";
double box[4];
ProjGeoCoordTransForEx(projstr,2,2,quad,box);//四角经纬度转换成投影坐标 N、E
投影参数 是RSD投影串，格式为：
投影串定义
/ ID = <投影ID>, //投影ID，例如UTM、ALBERS_CONIC_EQAREA等
/ a = <椭球体半长轴>, //椭球体半长轴，如WGS84 为6378137
/ f = <椭球体扁率>, //椭球体扁率(实为扁率的倒数1/f),如WGS84 为298.2572,
/ SP1 = <第1标准纬线>, //圆锥投影的第1标准纬线，麦卡托的割线，UTM、高斯投影的带号。
/ SP2 = <第2标准纬线>, //圆锥投影的第2标准纬线。
/ org_lat = <原点纬度>, //投影的纬度原点。
/ org_lon = <中央经线>, //投影的中央经线。
/ FN = <伪北>, //False Northing，伪北坐标。
/ FE = <伪东>, //False Easting，伪东坐标。
/ K = <系数>, //变换系数。
/ DATUM = <大地基准>, //缺省为WGS84。
这里STRING projstr = "";定义了一个空串，为缺省的UTM投影，投影参数由数据范围自动探测得到。
7. 估算像元尺度和正射后图像大小
//估算像元大小，南北方向和东西方向都估算，取平均值，然后取整。（或者根据卫星已知信息直接指定像元尺度）
double nsize,esize;
double pixsize = SwEvalPixelSize(quad,height,width,nsize,esize); //估算像元尺度的函数
//正射校正后的图像大小/
height = (box[2]-box[0]) / pixsize;
width = (box[3]-box[1]) / pixsize;
8. 正射校正
//正射校正，///
int precis = TRUE; //精度选项，TRUE时全程使用 8 字节浮点精度。
int resamp = CC;//重采样像元重定位方式 FW=前向重采样，NN=最近邻点法，BL=二次线性插值法，CC=三次褶积法。
WORD retBuffer[height][width][bands]; //哑数组，接收处理结果
int boo = Sw7WristbandsBuffer3(tifName,rpcName,retBuffer,box,projstr,pixsize,demBuffer,demGeoBox,demPixSize,precis,resamp);
正射校正结果保存在retBuffer缓存中。
9. 统计处理时间
double tm = (GetTickCount() - tc) / 1000.0;
Print("正射校正结束，并保存在数组retBuffer中。总计用时 %f 秒。",tm);
10. 输出结果、结束程序
//正射校正结果导出为TIFF数据格式文件。
Print("开始导出正射校正结果...");
tc = GetTickCount();
int rank,dims,datatype;
int dims[16];
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);
double outpuTime = (GetTickCount() - tc) / 1000.0;
Print("数据保存为 %s 用时 %f 秒。",outName,outpuTime);
Print("主程序结束");
return 1;
二、运行代码
1. 准备代码
将上述代码放入一个main{}内，保存为一个扩展名为.c的文本文件。启动RSD，在 脚本编辑 窗口打开这个文件，如图1。
图1 运行脚本代码
2. 调试运行
在脚本的第20行的行号后面点击，出现一个小方块，这个小方块是程序断点。点击工具栏里面的小三角，程序开始运行，DEM数据生成以后，程序执行到第20行发生中断，该行代码变为粉色。程序此时停在这里等待。
3. 观察变量（可视化计算）
打开 监视 窗口，双击点亮20行之前的任意一个变量，然后右击该变量，在弹出菜单添加监视，则该变量及其对应的值出现在监视窗口。用户在这里可以监视变量值的正确与否。
对于二维以上数组的观察可以使用预览图的方式。例如上述代码第19行读来全国DEM数据在demBuffer中，观察一下读来的数据是否正确，可以双击点亮该变量，右击，在弹出菜单选 SeeSee。出现图2的DEM图像，说明数据正确。
图2 观察变量合成的图像
4. 继续运行程序并输出结果
点击工具条上的小三角，程序继续运行。先输入一个待正射校正的高分数据集。再指定一个输出路径。运行结束后在指定的输出目录里面就会找到这个正射校正后的 .tif数据文件。
在RSD中打开这个TIFF文件，图3为GF2正射校正结果。
图3 GF2多光谱数据正射校正结果
从图3中红框可见，GF2多光谱数据正射校正用时6.187秒。
三、测试多种数据正射校正
1. GF1 多光谱数据
图4 GF1多光谱数据正射校正结果
GF1多光谱数据正射校正用时4.891秒。
2. GF多模多光谱数据
图4 GF多模多光谱数据正射校正结果
GF多模多光谱数据正射校正用时17.453秒。
3. 吉林1号多光谱数据
图5 吉林1号多光谱数据正射校正结果
吉林1号多光谱数据正射校正用时1.453秒。
4. 资源3多光谱数据
图6 资源3号多光谱数据正射校正结果
资源3号多光谱数据正射校正用时15.719秒。
其它还试验了多种卫星数据，这里就不一一列举了。希望同学们用该脚本测试更多的卫星数据，如遇到问题，将结果及时反馈给我，以便及时解决。