1、微波遥感编程实习报告一、 实习目的本次实习希望通过编程实现基于 R-D 方程的 SAR 影像几何校正，比较校正前后SAR 影像特征与区别，分析距离成像几何原理，掌握 R-D 校正步骤。希望通过这次实习了解 SAR 成像几何原理，熟悉距离式成像过程。掌握距离 -多普勒方程进行校正的原理和难点，认识校正所需参数意义与读取方法。理解并掌握斜视 SAR 成像中距离式成像的机理与方法，了解 SAR 图像特殊几何特征成因与校正必要性；掌握基于 R-D 构像方程的 SAR 图像像素坐标与大地坐标间的定位关系；在此基础上编程实现 SAR 图像的几何校正过程。二、 实习内容1. 数据介绍SAR 影像的头文件信息
2、中包括该幅影像上 5 个成像点的 XYZ 坐标，第一点的成像时间和 5 个点间的成像时间间隔，影像的第一行成像时间与成像间隔，影像距离向的成像时间与时间间隔，SAR 影像 4 个角点的经纬度以及投影椭球的长短半轴。这些数据用于确定任意时刻下成像点的位置矢量、速度矢量和加速度矢量，确定对应的 DEM 范围以及确定 DEM 上各点对应 SAR 影像的 ij 坐标和灰度值。这些参数已经存储在 parameter.txt 里，可直接读取使用，但注意单位转换问题。DEM 影像是 ers2dem.img，可以利用 DEM 提供的 H 结合经纬度坐标 BL，利用投影椭球信息，计算出成像点的大地坐标 XYZ，
3、但注意裁剪 SAR 对应到 DEM 的范围。SAR 影像灰度信息则是存储在 DAT_01.001 内，注意 SAR 影像的数据格式，这里存储的是复数形式，需要计算强度信息作为灰度值。满足公式：强度=实部+虚部。数据的预处理包括利用影像 5 个成像点的 XYZ 坐标和第一点成像时间、成像时间间隔，利用多项式拟合的方法，确定轨道系数 ABC，从而可求出任意时刻的卫星位置矢量，位置矢量求导可得到速度矢量，再求导可得到加速度矢量，这些都是在之后的 R-D 方程中有所利用。注意为解求矩阵方便，可将第一点成像时间当作 0。其余的数据预处理包括单位转换以及由于双程传输引起的时间转换。考虑到 DEM 的分辨率
4、是 90m，而 SAR 影像分辨率是 30m，所以还要进行DEM 加密，这里使用的方法是双线性内插方法实现 DEM 加密。2. R-D 方程已知 R-D 方程由 3 条方程组成，包括距离方程、多普勒频率方程和椭球方程，如下图。图 1 R-D 方程利用 R-D 方程进行几何校正时，使用的是间接校正法。即利用 DEM 中各点的经纬度坐标 BL 结合 DEM 影像上的高度信息 H，加上投影椭球参数，计算出各点的大地坐标 XYZ，再由该点的位置矢量结合当时卫星的位置矢量、速度矢量和加速度矢量代入距离-多普勒方程，确定成像时刻，从而找到 SAR 影像上对应的（i，j ） ，取出此处的灰度值，赋给 DEM
5、 对应的大地坐标处。最后再内插得到的灰度影像，即可得到纠正后的 SAR 影像。三、 实习步骤下面将根据上述实习内容一步步地编程实习 SAR 影像几何校正。几何过程可按概括为以下的流程图：开始卫星头文件参数读取多项式拟合求解拟合参数 A B CD E M 影像读取经纬度 B L 结合D E M 影像中的 H求解大地坐标X Y Z地面点坐标和初始迭代时间代入多普勒方程迭代解求 D E M 上各点对应到 S A R 影像上的行列号 （ i , j ）读取复数类型的S A R 影像求出强度值作为灰度将 S A R 影像上对应的灰度值赋给该地面点坐标得到几何校正后的 S A R 影像并保存结束图 2 几
6、何校正流程图1. 数据读取和预处理需要利用的 SAR 影像头文件数据已经保存在 parameter.txt 中，编程实现数据读取及保存。图 3 参数读取并存储结合 5 点的位置矢量和第一点的成像时间和成像时间间隔，利用多项式拟合求解拟合参数，由拟合参数可以求解任意时刻的卫星位置矢量、速度矢量和加速度矢量，这里对拟合参数的求解是利用矩阵最小二乘法求解的。在矩阵计算时，由于时间 t 较大，高次方运算可能会超过数据范围，所以把初始成像时间设为 0。注意一些数据需要进行转换，如投影椭球的长短半轴长单位是千米，需要转换到米为单位；经纬度是以度为单位的，在由经纬度坐标 BLH 计算到大地坐标 XYZ 时需
7、要把经纬度转为弧度为单位；考虑到双程传输，时间需要进行改正。需要利用到 DEM 数据和 SAR 影像，利用 GDAL 库的 RasterIO 函数将影像数据读入。注意 SAR 影像 DAT_01.001 是以复数形式存储，分为实部和虚部，数据格式是 CInt16，而我们需要的是它的灰度信息，灰度信息是影像的强度，符合公式：强度= 实部+ 虚部 ，利用公式求出其强度信息并读出来存储起来。编程实现 DEM 影像读取和 SAR 影像读取并存储。图4 DEM 的读取图 5 SAR 的读取2. 由 R-D 方程计算行列号利用 R-D 方程间接校正，即利用 DEM 的大地坐标，结合卫星位置矢量、速度矢量、
8、加速度矢量代入 R-D 方程，求解出该点对应到 SAR 影像上的行列号。对多普勒方程求微分，有 2()AspsVRta式中 表示卫星速度矢量和位置矢量， 表示 DEM 大地坐标，对 DEM、 的任一点来说，成像时间未知，多普勒方程不为 0，可针对方位向时间进行迭代。以中间行成像时间为初始迭代时间 ，求出当前时刻的0midt，结合 DEM 上该点的 ，代入上式可以求出 ，令、 、 At，再进行下一次迭代，随着迭代次数的增加， 会减小，成像10+Att时间会越趋于真实值。可以通过设置一个阈值，当 时，迭代停止，t此时求出的时间为真实方位向时间。为使循环不陷入死循环，还需设置一个迭代次数限制，所以最
9、后迭代次数到了或是满足小于阈值的条件时，迭代停止，求出真实方位向成像时间。再利用第一行成像时间和方位向成像时间间隔，求出行号。这里设置阈值 且 ，公式为-15=0N(_)/itaziti求出方位向时间后，可求出正确的卫星位置矢量 ，结合 求出斜距 ，距离向时间为 。再由距离向时间求出列数|spR /trngeRc0(_)jtat_range最后注意判定行列号在 SAR 影像包括的范围内时，给该点赋上灰度值，若不在区域内时，把灰度值赋为 0。行列号计算编程如下：图 6 R-D 方程间接几何校正3. 结果图像的存储与输出利用 GDAL 库实现上一步中保存下来的灰度矩阵存为图像并输出，可以将校正前后
10、的 SAR 影像对比。代码部分如下：图 7 结果的存储与输出图 8 运行结果四、 实习心得通过这次几何校正编程实习，我弄清楚了利用 R-D 方程间接几何校正的具体流程，需要结合 DEM 数据，迭代求出 SAR 影像上的行列号。并且了解到 GDAL 库如何读取影像已经存储影像，GDAL 库还可用于影像截取指定区域实现影像内插，内插方式可以选择最近邻内插和双线性内插等方式。在这次实习中，我也遇到了许多问题。首先利用GDAL 中的 RasterIO 可以读取出影像各个像元的灰度值，对 DEM 来说其灰度值就是它的高程，还可以读取出影像左上角点的地理坐标，可以利用这个求出每个像元的地理坐标，这里是经纬
11、度值。我原本希望将计算出的经纬度值都以矩阵的方式存储起来，但程序运行到这里时总是提示错误，错误理由检查后发现是内存没分配到，在询问老师后，老师说每个程序运行时内存都是一定的，我这样写内存不足没法实现分配存储，建议我在每次需要利用 BLH 求解大地坐标 XYZ 时再计算经纬度 BL，也不要把它存储起来。然后就是计算对应 SAR 影像的行列号处，每次运行完我编写的程序后得到的影像都是一片黑，设断点调试程序时发现每次的 SAR 影像上对应的行号 i 计算出来的都是负值，不在 SAR 影像包括的范围内。经检查后发现是自己的经纬度 BL 弄反了，并且由经纬度坐标 BLH 计算地理坐标 XYZ 的公式就写
12、错了，导致每次运行结果都是负值，修改了这两处错误后，再运行编写的程序即可得到 SAR 影像几何校正后的结果。对比纠正前后的 SAR 影像，发现不仅是影像的大小、形状有变化，主要差别在于几何校正后的 SAR 影像包含了地理坐标信息， ，几何校正后的影像反映的才是雷达成像时雷达波扫过地表的真实区域范围和回波分布情况。总的来说，这次编程实习，我收获很大。首先是真正动手编程去实习基于 R-D 方程的几何纠正，这让我对课本知识有了更深刻地理解，而不是仅有个概念；其次是学习了怎么利用 GDAL 读图像并实现内插、保存等操作，这十分有用。当然还是很感谢老师在编程过程中提供的帮助，老师十分耐心地予以指导，在我数组指针无法动态提供内存时，让我们不需要保存这些中间数据，只需在使用时再直接计算即可。所以很感谢老师这几天的耐心指导和帮助，让我们能真的有所收获，而不是虚度光阴。
