# hsRPC **Repository Path**: iam002/hs-rpc ## Basic Information - **Project Name**: hsRPC - **Description**: 【SAR001】遥感影像RPC模型的拟合 - **Primary Language**: Unknown - **License**: MIT - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 0 - **Forks**: 0 - **Created**: 2025-02-28 - **Last Updated**: 2025-03-10 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README 遥感影像的**有理函数模型 (Rational Function Model, RFM)** 是一种广泛应用的几何校正模型，其通过比值多项式关联地面坐标与影像坐标，具有形式简洁、通用性强等特点。有理函数的系数简称为 RPC（Rational Polynomial Coefficients），故也称为 RPC 模型。 RFM的核心是用有理多项式替代传统传感器严密成像模型（如共线方程或距离多普勒模型），将地面点的大地坐标（纬度、经度、高程）与其对应的影像坐标（行、列）通过比值多项式关联。这种设计使模型独立于具体传感器参数，适用于光学、SAR等多种遥感数据。 ### RFM模型的优势 1. **通用性与独立性** RFM不依赖具体传感器参数（如姿态角、轨道数据），可适用于多源遥感数据，解决了传统严密模型参数保密或获取困难的问题。例如，SAR影像的距离多普勒模型常因参数保密性受限，而RFM可提供替代方案。 2. **高精度与计算效率** RFM通过高阶多项式可逼近严密模型的精度（误差通常在亚像素级别），同时计算复杂度低，适合大规模遥感数据的并行快速处理。 3. **标准化与兼容性** 其参数通常以**RPC文件（Rational Polynomial Coefficients）**形式存储，如`.rpb`或`.rpc`文件，包含正则化参数和多项式系数，被主流遥感软件（如ENVI、GDAL）广泛支持。 ### 数学形式与参数结构 RFM 的数学表达式为比值多项式，具体形式如下： $$ \begin{aligned} \text{Line} &= \frac{N_l(\text{Lat},\ \text{Lon},\ \text{Hgt})} {D_l(\text{Lat},\ \text{Lon},\ \text{Hgt})} \\ \text{Sample} &= \frac{N_s(\text{Lat},\ \text{Lon},\ \text{Hgt})} {D_s(\text{Lat},\ \text{Lon},\ \text{Hgt})} \\ \end{aligned} $$ 其中，$N_l$、$D_l$、$N_s$、$D_s$ 为三次多项式，每个多项式包含20项系数。三次多项式的形式约定如下： $$ \begin{aligned} N_l(L,P,H) & = \sum_{i=0}^{3} \sum_{j=0}^{3} \sum_{k=0}^{3} a_t L^i P^j H^k \\ & = a_0 + a_1 L + a_2 P + a_3 H + a_4 LP \\ & + a_5 LH + a_6 PH + a_7 L^2 + a_8 P^2 + a_9 H^2 \\ & + a_{10} PLH + a_{11} L^3 + a_{12} LP^2 + a_{13} LH^2 + a_{14}L^2P\\ & + a_{15} P^3 + a_{16} PH^2 + a_{17} L^2 H + a_{18} P^2 H + a_{19} H^3 \end{aligned} $$ $$ \begin{aligned} D_l(L,P,H) &= \sum_{i=0}^{3} \sum_{j=0}^{3} \sum_{k=0}^{3} b_t L^i P^j H^k \\ N_s(L,P,H) &= \sum_{i=0}^{3} \sum_{j=0}^{3} \sum_{k=0}^{3} c_t L^i P^j H^k \\ D_s(L,P,H) &= \sum_{i=0}^{3} \sum_{j=0}^{3} \sum_{k=0}^{3} d_t L^i P^j H^k \\ \end{aligned} $$ 特别地，为了使模型求解具有唯一性，通常将分母多项式的常数项设置为 1（即 $b_0 = 1$，$d_0 = 1$），则 RPC 模型共计有 20 + 19 + 20 + 19 = 78 个参数需要求解。 ### RPC模型求解步骤 #### 获取采样点 RPC 模型的求解一般采用**不依赖地形**的算法进行求解。具体的做法是：在覆盖遥感影像经纬高范围的三维空间（$x$ : 经度，$y: $ 纬度，$z$ : 高程）进行均匀采样，得到若干个虚拟控制点的物方坐标： $$ Q_i = [x_i,\ y_i,\ z_i], \ i \in [1,\ N] $$ 通过严密几何定位模型，获取每个虚拟控制点对应的像方坐标（$c$ : 列索引，$r$ ：行索引）： $$ q_i = [c_i,\ r_i], i \in [1,\ N] $$ ![](img/resamp.png) #### 计算偏移量和缩放系数 为提高参数求解稳定性，RFM对物方坐标和像方坐标进行**正则化处理**，具体的公式为： $$ x' = (x - \text{offset}) / \text{scale} $$ RPC 的偏移量 offset 和缩放系数 scale 可根据虚拟控制点来计算： $$ \begin{aligned} \text{longOffset} &= \left( 10^{-6} \cdot \left\lfloor \max(x_i) \cdot 10^6 \right\rfloor + 10^{-6} \cdot \left\lfloor \min(x_i) \cdot 10^6 \right\rfloor \right)/2 \\ \text{latOffset} &= \left( 10^{-6} \cdot \left\lfloor \max(y_i) \cdot 10^6 \right\rfloor + 10^{-6} \cdot \left\lfloor \min(y_i) \cdot 10^6 \right\rfloor \right)/2 \\ \text{heightOffset} &= \left( 1000\cdot \left\lfloor \frac{\max(z_i)}{1000} \right\rfloor + 1000\cdot \left\lfloor \frac{\min(z_i)}{1000} \right\rfloor \right)/2 \\ \text{lineOffset} &= \text{height} /2 \\ \text{sampOffset} &= \text{width} /2 \\ \end{aligned} $$ $$ \begin{aligned} \text{longScale} &= \left( 10^{-6} \cdot \left\lfloor \max(x_i) \cdot 10^6 \right\rfloor - 10^{-6} \cdot \left\lfloor \min(x_i) \cdot 10^6 \right\rfloor \right)/2 \\ \text{latScale} &= \left( 10^{-6} \cdot \left\lfloor \max(y_i) \cdot 10^6 \right\rfloor - 10^{-6} \cdot \left\lfloor \min(y_i) \cdot 10^6 \right\rfloor \right)/2 \\ \text{heightScale} &= \left( 1000\cdot \left\lfloor \frac{\max(z_i)}{1000} \right\rfloor -1000\cdot \left\lfloor \frac{\min(z_i)}{1000} \right\rfloor \right)/2 \\ \text{lineScale} &= \text{height} /2 \\ \text{sampScale} &= \text{width} /2 \\ \end{aligned} $$ 上式中，$\text{height}$ 和 $\text{width}$ 分别表示图像的高和宽，由于标准的 RPC 文件没有给图像大小预留参数，建议将图高和图宽存放到像方坐标的偏移量和缩放系数中。 后文，除特别说明外，默认像方坐标和物方坐标均进行了正则化处理。 另外，由于标准的RPC文件为文本格式，我个人建议将偏移量和缩放系数保留小数点前n位（如小数点后6位），舍弃多余的尾数，以免后续读取RPC文件引入额外的舍入误差；而对于求解得到的系数，应该尽可能的保存多的小数点位数。 #### 转换为线性模型 前文所述，RPC模型的分母多项式的常数项为1，将比值多项式的分母移到等式的右侧，稍加变形，可以将比值多项式的非线性模型转换为线性模型： $$ \begin{aligned} (a_0 + a_1 L + \cdots + a_{19} H^3) &= (1 + b_1 L + \cdots + b_{19} H^3) r \\ (a_0 + a_1 L + \cdots + a_{19} H^3) &- (b_1 L + \cdots + b_{19} H^3) r = r \end{aligned} $$ 表示为矩阵形式为： $$ \begin{bmatrix} 1 & L & \cdots & H^3 & -rL & \cdots & -rH^3 \end{bmatrix} \begin{bmatrix} a_0 \\ \vdots \\ a_{19} \\ b_1 \\ \vdots \\ b_{19} \end{bmatrix} = r $$ #### SVD + ICCV 求解 我们可以使用**最小二乘法求解** RPC 模型。为了更好的数值稳定性，强烈推荐使用**奇异值分解法**进行求解。 下面给出奇异值分解求解 $Ax=b$ 最小二乘解的公式推导。 对矩阵 $A$ 进行 SVD 分解： $$ \begin{aligned} A &= UDV^T \\ A^T &= VDU^T \\ A^TA &= VD^2V^T \end{aligned} $$ 其中，$A$ 和 $U$ 为 $m\times n$ 矩阵，$D$ 和 $V$ 为 $n \times n$ 矩阵。 $Ax = b$ 最小二乘解为： $$ x_{svd} = VD^{-1}U^T b $$ 相比使用伪逆求解，SVD 求解的数值稳定性更高。 接下来，可以使用 ICCV（iteration by correcting characteristic value）方法进一提高 SVD 的求解精度。这是一种迭代方法，将 SVD 的求解结果作为迭代初值： $$ x^{(k)} = (A^TA+I)^{-1}(A^Tb+x^{(k-1)}) $$ 其中，$I$ 为 $n\times n$ 单位矩阵。 使用 SVD 分解结果可进一步化简 ICCV 迭代公式： $$ \begin{aligned} x^{(k)} &= (A^TA+I)^{-1}(A^Tb+x^{(k-1)}) \\ &= (VD^2V^T +I)^{-1} (VDU^Tb+x^{(k-1)}) \\ &= (VD^2V^T +VV^T)^{-1} (VDU^Tb+x^{(k-1)}) \\ &= (V(D^2+I)V^T)^{-1} (VDU^Tb+x^{(k-1)}) \\ &= V(D^2+I)^{-1}V^T (VDU^Tb+x^{(k-1)}) \\ &= V(D^2+I)^{-1}DU^Tb + V(D^2+I)^{-1}V^T x^{(k-1)}\\ x^{(k)} &= M_0 + M_1 x^{(k-1)} \end{aligned} $$ $M_0$ 和 $M_1$ 是在迭代过程中是常数。迭代终止条件为 $|x^{(k)} - x^{(k-1)}| < \xi$。 **补充**：RFM 的各项系数的数量级可达到 1e-7 数量级，RPC和拟合点数据**必须使用 double 数据类型**才可以保证拟合精度（single 数据类型的拟合误差高达 0.5 个像素）。 ### Maltab 实现代码 [hsRPC](+hsRPC) 包的内容: - [hsRPC.fitting](+hsRPC/fitting.m) 使用 SVD + ICCV 拟合 RPC 模型 - [hsRPC.project](+hsRPC/project.m) 应用RPC模型进行投影, 将物方/像方坐标转换为像方/物方坐标 - [hsRPC.read](+hsRPC/read.m) 读取标准RPC文件 - [hsRPC.write](+hsRPC/write.m) 将RPC结构体保存为标准RPC文件 代码演示： ![](img/demo.gif) ### 实验结果记录 #### 实验1: 使用 single 数据拟合 RPC 模型 04：[demo_fitrpc_in_single.m](demo_fitrpc_in_single.m) ![](img/demo_fitrpc_in_single_04.png) 数据归一化范围为 [-1, 1] ，SVD 分解法拟合失败。相比之下，LC 算法更加稳健些。 #### 实验2 : 使用 double 数据拟合 RPC 模型 但截断为 single 04：[demo_userpc_in_single.m](demo_userpc_in_single.m) ![](img/demo_userpc_in_single_04.png) 精度存在损失，两算法结果无明显差异。 #### 实验3：简单神经网络 float32 拟合效果不行 双层全连接网络（参数量等同于RPC模型）： ```python class SimpleFCN(nn.Module): def __init__(self, input_channels, output_channels): super(SimpleFCN, self).__init__() self.m_fc = nn.Sequential( nn.Linear(input_channels, 16), nn.ReLU(), nn.Linear(16, output_channels) ) def forward(self, x): return self.m_fc(x) ``` 01：[train_rpc_net.py](train_rpc_net.py) 归一化后的误差 ```bash Epoch [096/100], Loss: 0.00000047 Epoch [097/100], Loss: 0.00000054 Epoch [098/100], Loss: 0.00000053 Epoch [099/100], Loss: 0.00000053 Epoch [100/100], Loss: 0.00000046 Test Loss: 0.00000033 ``` ![](img/test_lost.png) 实际误差： ![](img/net_error.png)