后方交会检校(2)--误差方程及各项系数解求及精度

上文中明确了两个式子:
$$
\lbrace
\begin{matrix}\frac{f}{\overline{Z}}dXs+\frac{x}{\overline{Z}}dZs-(f+\frac{x^{2}}{f})d\varphi-\frac{xy}{f}d\omega+yd\kappa -\frac{f}{\overline{Z}}dX-\frac{x}{\overline{Z}}dZ+\frac{x}{f}df+dx_{0}-(x-x_{cal})=0 \\\ \frac{f}{\overline{Z}}dYs+\frac{y}{\overline{Z}}dZs -\frac{xy}{f}d\varphi-(f+\frac{x^{2}}{f})d\omega-xd\kappa -\frac{f}{\overline{Z}}dX-\frac{y}{\overline{Z}}dZ+\frac{y}{f}df+dy_{0}-(y-y_{cal})=0
\end{matrix}$$

上式即称为线性化的共线方程

$$\lbrace \begin{matrix}\frac{f}{\overline{Z}}dXs+\frac{x}{\overline{Z}}dZs-(f+\frac{x^{2}}{f})d\varphi-\frac{xy}{f}d\omega+yd\kappa -\frac{f}{\overline{Z}}dX-\frac{x}{\overline{Z}}dZ+\frac{x}{f}df+dx_{0}-(x-x_{cal})=v_{x}\\\\\frac{f}{\overline{Z}}dYs+\frac{y}{\overline{Z}}dZs -\frac{xy}{f}d\varphi-(f+\frac{x^{2}}{f})d\omega-xd\kappa -\frac{f}{\overline{Z}}dX-\frac{y}{\overline{Z}}dZ+\frac{y}{f}df+dy_{0}-(y-y_{cal})=v_{y}\end{matrix}$$

上式则称为误差方程

这里面$-(x-x_{cal})$即每一步计算后的值与实际测量值的差,在平差计算过程中,多余的观测数据就是为了进行平差计算获得更精确更均匀的结果。分别命名为$-l_{x}$和$-l_{y}$.一个控制点物方坐标和一个像点坐标组成的两个方程也就有两个误差方程,误差方程的基本公式为

$$V=AX-L$$

有的地方也写为

$$V=A\Delta -L$$

其中V代表误差矩阵,X是待求解的每个元素的改正数矩阵,L是x-x计 y-y计的矩阵。

通常,不会解算物方空间的坐标,同时加入对畸变系数的解算和矫正(将在第四篇提到),将X分为$X_{in}$、$X_{out}$、$X_{distortion}$,则上式变形为

$$V=AX_{out}+BX_{in}+CX_{distortion}-(xy)_{res}$$

其中,A B C都是各项未知数的系数,$xy_{res}$是像点坐标的误差

$$\lbrace
\begin{matrix}x_{res}=x-x_{cal}+x_{dis}\\\\y_{res}=y-y_{cal}+y_{dis}\end{matrix}$$

其中,右边三项分别为像点观测坐标,像点计算坐标,畸变误差。
在后方交会时,当只计算外方位元素时:

$$V_{2n\times 1}=A_{2n\times6}X_{6\times1}-L_{2n\times1}$$.

其中,式子每两行具体为

$$\begin{bmatrix}v_{x}\\\\v_{y}\end{bmatrix}=\begin{bmatrix}\frac{f}{\overline{Z}}&0&\frac{x}{\overline{Z}}&-(f+\frac{x^{2}}{f})&-\frac{xy}{f}&y\\\\0&\frac{f}{\overline{Z}}&\frac{y}{\overline{Z}}&-\frac{xy}{f}&-(f+\frac{y^{2}}{f})&-x\end{bmatrix}\begin{bmatrix}dXs\\\\dYs\\\\dZs\\\\d\varphi\\\\d\omega\\\\d\kappa\end{bmatrix}-\begin{bmatrix}l_{x}\\\\l_{y}\end{bmatrix}$$

若要求出V方程的最佳估计值,在摄影测量中,一般按照最小二乘原理进行解答,即$$V^{T}V=min$$这就转化为了对该式求极值的问题。

$$V^{T}V=(AX-L)^{T}(AX-L)=X^{T}A^{T}AX-2X^{T}A^{T}L+L^{T}L\Rightarrow \frac{\partial V^{T}V}{\partial X}=2A^{T}AX-2A^{T}L=0$$

上式中大X的含义是待求解的每个元素的残差矩阵包含了Xs Ys Zs Ω φ k,而不是实际测量三维点坐标X。此外,推导过程中对$V^{T}V$求导,则可得到对应6个残差的最佳估计值。这个过程称为误差方程式的法化

上式整理后可得$$A^{T}AX-A^{T}L=A^{T}(AX-L)=0$$

有的地方亦写作$$A^{T}P(A\Delta -L)=0$$Δ与X同含义,而P是测量平差中进行稳健估计而使用的权矩阵,即
$$\begin{bmatrix}P_{1}&&&&&\\&P_{2}&&&&\\&&.&&&\\&&&.&&\\&&&&.&\\&&&&&P_{n}\end{bmatrix}$$
非主要内容略微说一下测量平差的估值。主要目的是处理含有误差(其中一些误差可能较大,定义为比最大偶然误差还要大的误差,称为粗差)依赖一定的数学模型,包含函数模型和随即模型,按照某种估计准则,求出未知参数的最优估值,并评定其精度。当粗差归为函数模型,粗差表现为观测量误差绝对值较大且偏离群体,处理思想是在正式进行最小二乘平差之前探测和定位粗差,然后剔除含有粗差的观测值,得到一组比较净化的观测值,模型为均值漂移模型;将粗差归为随机模型,粗差即表现为先验随机模型和实际随机模型的差异过大,处理方法为根据逐次迭代平差的结果不断地改变观测值的权或者方差,最终使粗差观测值的权趋于零或者方差趋于无穷大,模型为方差膨胀模型

该方法的实践逻辑是由于不大准确知道观测数据中有效信息和有害信息所占的比例及其它们具体包含在哪些观测中,要冒着损失一些效率的风险,去获得较可靠的、具有实际意义的、较有效的估值。针对极大似然估计的稳健估计模型主要使用M估计。而M估计应用最广的是选权迭代法。而P可以使等权独立观测的P=I,也可以是权函数矩阵代替观测权阵P。这里不展开来讲。

法化后的方程称为法方程式,一般写作

$$NX=A^{T}L$$,其中$N=A^{T}A$也把N称为法方程的系数

上两式就是法方程,也是△Xs,△Ys,△Zs,△k,△φ,△ω的解。

答解改正数:$$X=(A^{T}A)^{-1}A^{T}L$$或$$\Delta =(A^{T}PA)^{-1}A^{T}PL$$

计算每一步迭代值:$$X_{s}^{k+1}=X_{s}^{k}+dXs^{k+1},Y_{s}^{k+1}=Y_{s}^{k}+dYs^{k+1},Z_{s}^{k+1}=Z_{s}^{k}+dZs^{k+1}$$
$$\varphi ^{k+1}=\varphi ^{k}+d\varphi ^{k+1},\omega ^{k+1}=\omega ^{k}+d\omega ^{k+1},\kappa ^{k+1}=\kappa ^{k}+d\kappa ^{k+1}$$

这里注意,平坦地区不适宜内外元素同时答解,因为dXs0与dx0相关,dYs0与dy0相关。

空间后方交会的平面二维坐标的选点,应尽可能的布满整个区域。或者覆盖中心和边缘的更全面的位置,如图。不均匀以及不足够数量的控制点参与迭代会极大的影响最终计算的精度。



具体的迭代流程如下:

①确定各未知数的初值$X_{s}^{0},Y_{s}^{0},Z_{s}^{0},\varphi ^{0},\omega ^{0},\kappa ^{0}$
②对每个控制点计算误差方程系数aij和lx、ly。按照最上面的公式列出误差方程。
③答解线性方程组,得到△Xs,△Ys,△Zs,△k,△ω,△φ


④将增量△和初值相加,得到新的一个外方位元素。

⑤各个增量△与一个给定的限差作比较,若小于则停止迭代运算;若不是则重复②到⑤

分别对上式各步骤进行解释:
①$X_{s}^{0},Y_{s}^{0}$可以利用X Y的物方点均值设置,$Z_{s}^{0}$一般为H即航摄高度,有时也用mf表示。m是航摄比例尺,f是相机焦距。这些值也可以直接给定某个初值。三个角元素则一般设为0,k有时可根据航迹图给定。

②该步骤可分解为计算旋转矩阵R(R的具体值参考下文)、计算得到的变换坐标$\overline{X},\overline{Y},\overline{Z}$(其实就是地面辅助坐标系的坐标)、计算x计、y计(即由相应地面坐标点外方位元素计算得到的像点坐标,$$x_{cal}=-f\frac{\overline{X}}{\overline{Z}},y_{cal}=-f\frac{\overline{Y}}{\overline{Z}}$$、计算误差方程各项系数、组成误差方程V=AX-L

③该步骤其实就是矩阵运算,求出$A^{T}PA$ ($A^{T}A$) 及其逆矩阵。求出$A^{T}L$ ($A^{T}PL$),两式相乘求取X[△Xs,△Ys,△Zs,△k,△ω,△φ]

④根据改正数的值加上上一步的值得到各项新值。

⑤检查是否收敛时,一般只对角元素改正数设定限差,为0.1’(按弧度值是0.0000291rad)当三个角改正数都小于限差时迭代结束。
平差理论中,单位权中的误差为

$$m_{0}=\sqrt{\frac{\sum v_{i}^{2}}{2n-t}}​$$其中m0为单位权中误差。vi是第i个方程的残差,n是控制点个数,t是未知数个数。只解求外方位精度则只有6个精度。

根据传播定律,$$m^{2}=N^{-1}m_{0}^{2}$$式中$N^{-1}$是法方程系数矩阵的逆矩阵。即

$$N^{-1}=\begin{bmatrix}Q_{11}&Q_{12}&.&.&.&Q_{1t}\\Q_{21}&Q_{22}&.&.&.&Q_{2t} \\\ .&.&.&.&.&.\\\ .&.&.&.&.&.\\\ .&.&.&.&.&.\\Q_{t1}&Q_{t2}&.&.&.&Q_{tt}\end{bmatrix}$$,则第i个未知数中的误差为:

$m_{i}=\sqrt{Q_{ii}}m_{0}$,其中$Q_{ii}$为$N^{-1}$的各项系数的对角元素,也称为第i个未知数的权倒数。

——一点说明

再次声明,无论半片还是全片,增加控制点都可以提高精度,但精度提高的不是很大;采用相同数量控制点,全片比半片精度成倍的提高,分布良好更有利于精度改善。

热评文章