微波遥感集中编程实习报告二.doc

微波遥感编程实习报告
一、 实习目的
本次实习希望通过编程实现基于R-D方程的SAR影像几何校正，比较校正前后SAR影像特征与区别，分析距离成像几何原理，掌握R-D校正步骤。希望通过这次实习了解SAR成像几何原理，熟悉距离式成像过程。掌握距离-多普勒方程进行校正的原理和难点，认识校正所需参数意义与读取方法。理解并掌握斜视SAR成像中距离式成像的机理与方法，了解SAR图像特殊几何特征成因与校正必要性；掌握基于R-D构像方程的SAR图像像素坐标与大地坐标间的定位关系；在此基础上编程实现SAR图像的几何校正过程。
二、 实习内容
1. 数据介绍
SAR影像的头文件信息中包括该幅影像上5个成像点的XYZ坐标，第一点的成像时间和5个点间的成像时间间隔，影像的第一行成像时间与成像间隔，影像距离向的成像时间与时间间隔，SAR影像4个角点的经纬度以及投影椭球的长短半轴。这些数据用于确定任意时刻下成像点的位置矢量、速度矢量和加速度矢量，确定对应的DEM范围以及确定DEM上各点对应SAR影像的ij坐标和灰度值。这些参数已经存储在parameter.txt里，可直接读取使用，但注意单位转换问题。
DEM影像是ers2dem.img，可以利用DEM提供的H结合经纬度坐标BL，利用投影椭球信息，计算出成像点的大地坐标XYZ，但注意裁剪SAR对应到DEM的范围。
SAR影像灰度信息则是存储在DAT_01.001内，注意SAR影像的数据格式，这里存储的是复数形式，需要计算强度信息作为灰度值。满足公式：强度²=实部²+虚部²。
数据的预处理包括利用影像5个成像点的XYZ坐标和第一点成像时间、成像时间间隔，利用多项式拟合的方法，确定轨道系数ABC，从而可求出任意时刻的卫星位置矢量，位置矢量求导可得到速度矢量，再求导可得到加速度矢量，这些都是在之后的R-D方程中有所利用。注意为解求矩阵方便，可将第一点成像时间当作0。其余的数据预处理包括单位转换以及由于双程传输引起的时间转换。
考虑到DEM的分辨率是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对应的大地坐标处。最后再内插得到的灰度影像，即可得到纠正后的SAR影像。
三、 实习步骤
下面将根据上述实习内容一步步地编程实习SAR影像几何校正。几何过程可按概括为以下的流程图：

图2 几何校正流程图
1. 数据读取和预处理
需要利用的SAR影像头文件数据已经保存在parameter.txt中，编程实现数据读取及保存。

图3 参数读取并存储
结合5点的位置矢量和第一点的成像时间和成像时间间隔，利用多项式拟合求解拟合参数，由拟合参数可以求解任意时刻的卫星位置矢量、速度矢量和加速度矢量，这里对拟合参数的求解是利用矩阵最小二乘法求解的。在矩阵计算时，由于时间t较大，高次方运算可能会超过数据范围，所以把初始成像时间设为0。
注意一些数据需要进行转换，如投影椭球的长短半轴长单位是千米，需要转换到米为单位；经纬度是以度为单位的，在由经纬度坐标BLH计算到大地坐标XYZ时需要把经纬度转为弧度为单位；考虑到双程传输，时间需要进行改正。
需要利用到DEM数据和SAR影像，利用GDAL库的RasterIO函数将影像数据读入。注意SAR影像DAT_01.001是以复数形式存储，分为实部和虚部，数据格式是CInt16，而我们需要的是它的灰度信息，灰度信息是影像的强度，符合公式：强度²=实部²+虚部²，利用公式求出其强度信息并读出来存储起来。编程实现DEM影像读取和SAR影像读取并存储。

图4 DEM的读取
图5 SAR的读取
2. 由R-D方程计算行列号
利用R-D方程间接校正，即利用DEM的大地坐标，结合卫星位置矢量、速度矢量、加速度矢量代入R-D方程，求解出该点对应到SAR影像上的行列号

。对多普勒方程求微分，有

式中表示卫星速度矢量和位置矢量，表示DEM大地坐标，对DEM的任一点来说，成像时间未知，多普勒方程不为0，可针对方位向时间进行迭代。以中间行成像时间为初始迭代时间，求出当前时刻的，结合DEM上该点的，代入上式可以求出，令，再进行下一次迭代，随着迭代次数的增加，会减小，成像时间会越趋于真实值。可以通过设置一个阈值，当时，迭代停止，此时求出的时间为真实方位向时间。为使循环不陷入死循环，还需设置一个迭代次数限制，所以最后迭代次数到了或是满足小于阈值的条件时，迭代停止，求出真实方位向成像时间。再利用第一行成像时间和方位向成像时间间隔，求出行号。这里设置阈值且，公式为

求出方位向时间后，可求出正确的卫星位置矢量，结合求出斜距，距离向时间为。再由距离向时间求出列数
最后注意判定行列号在SAR影像包括的范围内时，给该点赋上灰度值，若不在区域内时，把灰度值赋为0。行列号计算编程如下：

图6 R-D方程间接几何校正
3. 结果图像的存储与输出
利用GDAL库实现上一步中保存下来的灰度矩阵存为图像并输出，可以将校正前后的SAR影像对比。代码部分如下：

图7 结果的存储与输出

图8 运行结果
四、 实习心得
通过这次几何校正编程实习，我弄清楚了利用R-D方程间接几何校正的具体流程，需要结合DEM数据，迭代求出SAR影像上的行列号。并且了解到GDAL库如何读取影像已经存储影像，GDAL库还可用于影像截取指定区域实现影像内插，内插方式可以选择最近邻内插和双线性内插等方式。在这次实习中，我也遇到了许多问题。首先利用GDAL中的RasterIO可以读取出影像各个像元的灰度值，对DEM来说其灰度值就是它的高程，还可以读取出影像左上角点的地理坐标，可以利用这个求出每个像元的地理坐标，这里是经纬度值。我原本希望将计算出的经纬度值都以矩阵的方式存储起来，但程序运行到这里时总是提示错误，错误理由检查后发现是内存没分配到，在询问老师后，老师说每个程序运行时内存都是一定的，我这样写内存不足没法实现分配存储，建议我在每次需要利用BLH求解大地坐标XYZ时再计算经纬度BL，也不要把它存储起来。然后就是计算对应SAR影像的行列号处，每次运行完我编写的程序后得到的影像都是一片黑，设断点调试程序时发现每次的SAR影像上对应的行号i计算出来的都是负值，不在SAR影像包括的范围内。经检查后发现是自己的经纬度BL弄反了，并且由经纬度坐标BLH计算地理坐标XYZ的公式就写错了，导致每次运行结果都是负值，修改了这两处错误后，再运行编写的程序即可得到SAR影像几何校正后的结果。
对比纠正前后的SAR影像，发现不仅是影像的大小、形状有变化，主要差别在于几何校正后的SAR影像包含了地理坐标信息，，几何校正后的影像反映的才是雷达成像时雷达波扫过地表的真实区域范围和回波分布情况。

总的来说，这次编程实习，我收获很大。首先是真正动手编程去实习基于R-D方程的几何纠正，这让我对课本知识有了更深刻地理解，而不是仅有个概念；其次是学习了怎么利用GDAL读图像并实现内插、保存等操作，这十分有用。当然还是很感谢老师在编程过程中提供的帮助，老师十分耐心地予以指导，在我数组指针无法动态提供内存时，让我们不需要保存这些中间数据，只需在使用时再直接计算即可。所以很感谢老师这几天的耐心指导和帮助，让我们能真的有所收获，而不是虚度光阴。

（注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。可复制、编制，期待你的好评与关注）