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

微波遥感编程实习报告微波遥感编程实习报告 一 实习目的 本次实习希望通过编程实现基于 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 影像几何校正 几何过程可按概 括为以下的流程图 开始 卫星头文件参 数读取 多项式拟合 求解拟合参 数ABC DEM影像读取 经纬度BL结合 DEM影像中的H 求解大地坐标 XYZ 地面点坐标和 初始迭代时间 代入多普勒方 程 迭代解求DEM上各 点对应到SAR影像上 的行列号 i j 读取复数类型的 SAR影像求出强度 值作为灰度 将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 影像上的行列 号 对多普勒方程求微分 有 2 A ssp ssps V RR t a RRV 式中表示卫星速度矢量和位置矢量 表示 DEM 大地坐标 对 DEM 的任一点来说 成像时间未知 多普勒方程不为 0 可针对方位向时间进行迭 代 以中间行成像时间为初始迭代时间 求出当前时刻的 0 mid tt 结合 DEM 上该点的 代入上式可以求出 令 At 再进行下一次迭代 随着迭代次数的增加 会减小 成像 10 AtttAt 时间会越趋于真实值 可以通过设置一个阈值 当时 迭代停止 At 此时求出的时间为真实方位向时间 为使循环不陷入死循环 还需设置一个 迭代次数限制 所以最后迭代次数到了或是满足小于阈值的条件时 迭代停 止 求出真实方位向成像时间 再利用第一行成像时间和方位向成像时间间 隔 求出行号 这里设置阈值且 公式为 15 10 10N 0 N it tazitazi 求出方位向时间后 可求出正确的卫星位置矢量 结合求出斜距 距离向时间为 再由距离向时间求出列数 sp RRR trangeR c 0 jtrangetranget range 最后注意判定行列号在 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 方 程