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

./微波遥感编程实习报告实习目的本次实习希望通过编程实现基于R-D方程的SAR影像几何校正,比较校正前后SAR影像特征与区别,分析距离成像几何原理,掌握R-D校正步骤。希望通过这次实习了解SAR成像几何原理,熟悉距离式成像过程。掌握距离-多普勒方程进行校正的原理和难点,认识校正所需参数意义与读取方法。理解并掌握斜视SAR成像中距离式成像的机理与方法,了解SAR图像特殊几何特征成因与校正必要性；掌握基于R-D构像方程的SAR图像像素坐标与坐标间的定位关系；在此基础上编程实现SAR图像的几何校正过程。实习容数据介绍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加密。R-D方程已知R-D方程由3条方程组成,包括距离方程、多普勒频率方程和椭球方程,如下图。图1R-D方程利用R-D方程进行几何校正时,使用的是间接校正法。即利用DEM中各点的经纬度坐标BL结合DEM影像上的高度信息H,加上投影椭球参数,计算出各点的坐标XYZ,再由该点的位置矢量结合当时卫星的位置矢量、速度矢量和加速度矢量代入距离-多普勒方程,确定成像时刻,从而找到SAR影像上对应的〔i,j,取出此处的灰度值,赋给DEM对应的坐标处。最后再插得到的灰度影像,即可得到纠正后的SAR影像。实习步骤下面将根据上述实习容一步步地编程实习SAR影像几何校正。几何过程可按概括为以下的流程图：图2几何校正流程图数据读取和预处理需要利用的SAR影像头文件数据已经保存在parameter.txt中,编程实现数据读取及保存。图3参数读取并存储结合5点的位置矢量和第一点的成像时间和成像时间间隔,利用多项式拟合求解拟合参数,由拟合参数可以求解任意时刻的卫星位置矢量、速度矢量和加速度矢量,这里对拟合参数的求解是利用矩阵最小二乘法求解的。在矩阵计算时,由于时间t较大,高次方运算可能会超过数据围,所以把初始成像时间设为0。注意一些数据需要进行转换,如投影椭球的长短半轴长单位是千米,需要转换到米为单位；经纬度是以度为单位的,在由经纬度坐标BLH计算到坐标XYZ时需要把经纬度转为弧度为单位；考虑到双程传输,时间需要进行改正。需要利用到DEM数据和SAR影像,利用GDAL库的RasterIO函数将影像数据读入。注意SAR影像DAT_01.001是以复数形式存储,分为实部和虚部,数据格式是CInt16,而我们需要的是它的灰度信息,灰度信息是影像的强度,符合公式：强度²=实部²+虚部²,利用公式求出其强度信息并读出来存储起来。编程实现DEM影像读取和SAR影像读取并存储。图4DEM的读取图5SAR的读取由R-D方程计算行列号利用R-D方程间接校正,即利用DEM的坐标,结合卫星位置矢量、速度矢量、加速度矢量代入R-D方程,求解出该点对应到SAR影像上的行列号。对多普勒方程求微分,有式中表示卫星速度矢量和位置矢量,表示DEM坐标,对DEM的任一点来说,成像时间未知,多普勒方程不为0,可针对方位向时间进行迭代。以中间行成像时间为初始迭代时间,求出当前时刻的,结合DEM上该点的,代入上式可以求出,令,再进行下一次迭代,随着迭代次数的增加,会减小,成像时间会越趋于真实值。可以通过设置一个阈值,当时,迭代停止,此时求出的时间为真实方位向时间。为使循环不陷入死循环,还需设置一个迭代次数限制,所以最后迭代次数到了或是满足小于阈值的条件时,迭代停止,求出真实方位向成像时间。再利用第一行成像时间和方位向成像时间间隔,求出行号。这里设置阈值且,公式为求出方位向时间后,可求出正确的卫星位置矢量,结合求出斜距,距离向时间为。再由距离向时间求出列数最后注意判定行列号在SAR影像包括的围时,给该点赋上灰度值,若不在区域时,把灰度值赋为0。行列号计算编程如下：图6R-D方程间接几何校正结果图像的存储与输出利用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方程的几何纠正,这让我对课本知识有了更深刻地理解,而