CN112835043B - 一种任意方向的二维形变的监测方法 - Google Patents

Abstract

本发明公开了一种任意方向的二维形变的监测方法，包括以下步骤：步骤1：对数据进行数据预处理；步骤2：对预处理的数据进行时空基准统一；步骤3：对时空基准同一的数据进行建立数学模型；步骤4：对建立数学模型的数据进行模型解算及精度验证；该方法不仅有利于对地质灾害结果进行直观的解译，还可以减少必要观测量，更有利于对解算结果进行检核。

Description

一种任意方向的二维形变的监测方法

技术领域

本发明涉及合成孔径雷达干涉测量技术，特别涉及一种任意方向的二维形变的监测方法。

背景技术

随着InSAR技术的使用越来越广泛，其固有的缺陷也随之显现出来。由于InSAR技术为侧视成像技术，所以其解算结果为一维的LOS向形变结果，然而仅凭一维LOS向形变可能会对形变结果造成误判，错判的现象。目前，恢复三维形变的方法有：

(1)多轨DInSAR方法：通过联合多轨道多平台的InSAR数据构建线性方程组进行三维形变的解算，但是目前SAR卫星均为极轨卫星，运行方向为南北向，使得LOS形变对南北向形变分量不敏感，从而导致南北向形变解算精度较低；

(2)多轨DInSAR与MAI联合：此方法从子孔径SAR影像生成的前、后视影像中提取方位向形变信息，将其和LOS向形变组合，因此只需升降轨数据即可完成三维形变解算，但是MAI提取的方位向形变精度较低；

(3)多轨DInSAR和GPS联合：该方法通过融合DInSAR数据和GPS数据，得到高时空分辨率的三维形变结果，但是该方法需要将GPS测量得到的稀疏数据插值到InSAR结果中，为了保证必要的精度，必须增加GPS站的数量和分布，这无疑增加了成本和工作量；

(4)多轨DInSAR与OT联合：Offset Tracking(OT)可以求出斜距向和方位向的形变信息，因此只需升降轨数据即可求出三维地表形变场，同时方位向形变对南北向形变最为敏感，可以弥补LOS向对南北向不敏感的不足，但是该方法计算量大且精度低；

(5)附加先验信息的单轨DInSAR或者单轨OT技术：单轨D-InSAR技术通过三维形变间关系的先验信息和单轨数据，结合D-InSAR技术或者OT技术，即可解算出地表三维形变场，但是该方法存在一些限制影响高精度的三维形变，如：没有冗余测量，LOS形变测量值的误差无法削减；先验信息构建的模型存在误差等。

但是，由于目前SAR卫星均为接近南北向飞行的极轨卫星，卫星飞行方向与LOS向互相垂直，故InSAR技术对南北方向形变信号难以捕获，因此南北方向形变解算的精度较低；对于某些地质灾害(如：滑坡等)，我们并不关心其三维方向的形变结果，如：当滑坡灾害坡面方向为东西向时，传统的三维形变将导致解算的冗余，因此分解任意方向的二维形变可以减少必要观测量，有利于对结果进行检核；当滑坡灾害坡面方向为任意方向时，那么我们更关心的是沿坡向的形变，传统的三维形变结果反而不利于我们对结果进行直观解译。因此，我们提出一种任意方向的二维形变的监测方法，该方法不仅有利于对地质灾害结果进行直观的解译，还可以减少必要观测量，更有利于对解算结果进行检核。

发明内容

本发明目的在于，克服目前InSAR三维形变解算方法的缺陷，提出一种任意方向的二维形变的监测方法，该方法易于实施，应用广泛，有利于对地质灾害结果进行直观的解译，可以减少必要观测量，更有利于对解算结果进行检核。

发明的技术解决方案如下：

一种任意方向的二维形变的监测方法，其特征在于，包括以下步骤：

步骤1：对数据进行数据预处理：即对获取的多种数据或者升降轨数据分别进行影像配准，干涉对构网，利用外部DEM数据去除模拟地形相位，得到差分干涉相位后，对差分干涉相位进行滤波，解缠，得到各数据LOS向形变序列结果；

步骤2：对预处理的数据进行时空基准统一；

步骤3：对时空基准统一的数据进行建立数学模型：即根据LOS向数据的形变结果建立观测方程，通过严密的几何关系建立平差的函数模型和随机模型；

步骤4：对建立数学模型的数据进行模型解算及精度验证；即采用迭代最小二乘的方法进行求解，通过残差进行重新定权，直至结果符合精度要求为止。

优选地，步骤1的干涉对构网具体为：

在对所有影像配准后，采用二连接的干涉对构网方式，即：1和2组成干涉对，1和3组成干涉对，2和3组成干涉对，用于任何一景影像都处于一个三角形闭合环中，有利于克服了时间失相干，同时又可以检验其中的误差。

优选地，步骤2包括以下子步骤：

步骤21：利用外部DEM数据将LOS向形变结果地理编码至地理坐标系下；

步骤22：将不同数据重叠区域进行裁剪，进而重采样至统一地理格网中；

步骤23：对不同数据选取相同时间基准作为起始的时间参考，将形变结果进行平移，达到时间基准的统一。

优选地，步骤3包括以下子步骤：

步骤31：根据InSAR的理论知识，建立严密的数学关系式，从而对LOS向形变进行任意方向二维形变分解；

其中，任意方向可以是两部分：第一部分：通过外部DEM数据得到坡向，根据地形计算该区域的形变，有利于得到沿坡面形变结果；

第二部分：人为指定某一固定方向，计算该区域相对于指定方向的形变，有利于对形变进行直观解译；考虑到升降轨数据成像的几何差异，建立升轨数学模型和降轨数学模型。

优选地，升轨数学模型的求解步骤：

步骤311：假定卫星以升轨方向飞行，由于LOS向形变结果为竖直向形变结果和水平向形变结果的投影之和，可得：

d_los1＝d_u×cosθ₁+(-d_h×sinθ₁)

式中，d_u为竖直向形变在LOS向投影，d_h为水平向形变在LOS向投影，θ为雷达的局部入射角；同时假定形变靠近卫星为正，远离卫星为负；

步骤312：将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献，即x方向形变为0，那么：

d_h＝d_y×sinγ_asc

式中，γ_asc＝270°-(α₁-α')，其中，α₁为卫星飞行方向与北方向的夹角，α'为x方向与北方向的夹角；

步骤313：将步骤412中的公式式代入步骤411式中，可得：

d_los1＝d_u×cosθ₁+d_y×cos(α₁-α')sinθ₁

降轨数学模型：

步骤321：假定卫星以降轨方向飞行，可得：

d_los2＝d_u×cosθ₂+(-d_h×sinθ₂)

步骤322：将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献即x方向形变为0，那么：

d_h＝d_y×sinγ_des

式中，γ_des＝(α₂-α')-90°。

步骤323：将步骤322式代入步骤321式中，可得：

d_los2＝d_u×cosθ₂+d_y×cos(α₂-α')sinθ₂

利用多个轨道数学模型或升降轨数学模型，得到竖直向和沿任意方向水平向的形变结果：

优选地，步骤3中

对于不同数据时间基准不统一的情况，提出保留不同数据时间间隔不同，其中，针对方程秩亏的现象，假设相邻位移序列的形变速率相等；

增加虚拟观测的方式对其进行约束，进而对形变序列直接进行分解；

假设三个不同几何的数据分别有N₁、N₂、N₃期，N＝N₁+N₂+N₃，数据获取的时间分别为

T＝T₁∪T₂∪T₃；

通过时序InSAR的方法得到相应的LOS向形变序列

d＝d₁∪d₂∪d₃；

通过多个方向的LOS形变分解出二维形变序列；

根据竖直向和沿任意方向水平向的形变结果，得到：

Gx＝d

其中，设计矩阵:

G＝diag(A⁽¹⁾A⁽²⁾A⁽³⁾)，

其中：

A⁽¹⁾＝(cos(α⁽¹⁾-α')sinθ⁽¹⁾cosθ⁽¹⁾)；A⁽²⁾＝(cos(α⁽²⁾-α')sinθ⁽²⁾cosθ⁽²⁾)，A⁽³⁾＝(cos(α⁽³⁾-α')sinθ⁽³⁾cosθ⁽³⁾)；

对于G矩阵某行中，某列是否为0与该轨道是否有观测值有关；

x＝(p₁ p₂ ... p_N)^T，其中

y,u分别是任意水平方向形变值和竖直方向形变值；

假设时间均不重叠且d₁ ⁽¹⁾，d₁ ⁽²⁾，d₁ ⁽³⁾≠0，观测量个数为N，待求参数个数为2N；

显然，秩亏，因为第一期LOS向观测值为0，那么分解得到的水平方向和竖直方向也应为0，为此，

假设y₁＝0，u₁＝0，在方程解算中，以增加虚拟观测的方式对其进行约束；

Cx＝0

其中，

为了减轻方程病态并使其满秩，假设相邻位移序列的形变速率相等，即

i＝1,2..N-2；那么，求解：

Rx＝0

其中，

把三个式子合并观测方程

利用最小二乘估计，解算方程。

优选地，步骤4中通过迭代最小二乘求解具体过程为：

步骤61：若

则满足解算的精度要求，迭代结束；

其中，

分别为两次的残差平方和；

步骤62：利用观测方程得到的参数估计值

将其代入

中得到整个模型的解算残差，利用残差的均方根误差

来评定模型参数解算的内符合精度。

有益效果

1.本发明公开了一种任意方向二维形变解算方法，该方法通过利用多个轨道数据或升降轨数据，将InSAR数据得到的LOS向形变结果分解为竖直向形变和某一任意方向，可以是人为指定；或给定数据计算，如地形坡向，的水平形变。

2.本发明的方法原理直观，易于编程实现和应用扩展，是一种稳健可靠的InSAR数据二维形变分解的方法，不仅有利于对地质灾害结果进行直观的解译，还可以减少必要观测量，有利于对解算结果进行检核。

附图说明

图1为二连接构网示意图。

图2为升轨卫星LOS向形变分解示意图。

图3为升轨卫星水平形变分解示意图。

图4为降轨卫星LOS向形变分解示意图。

图5为降轨卫星水平形变分解示意图。

图6为研究区域SRTM数据。

图7为研究区域坡向图。

图8为竖直向形变速率图。

图9为任意水平形变速率图。

图10为总流程图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图10所示，一种任意方向的二维形变的监测方法，包括以下步骤：

步骤1：数据预处理：对获取的多种数据或者升降轨数据分别进行影像配准，干涉对构网，利用外部DEM数据去除模拟地形相位，得到差分干涉相位后，对其进行滤波，解缠，得到各数据LOS向形变序列结果。

步骤2：时空基准统一：(1)首先利用外部DEM数据将LOS向形变结果地理编码至地理坐标系下，然后将不同数据重叠区域进行裁剪，进而重采样至统一地理格网中；(2)对不同数据选取相同时间基准作为起始的时间参考，然后将形变结果进行平移，从而达到时间基准的统一。

步骤3：建立数学模型：根据LOS向数据的形变结果建立观测方程，然后通过严密的几何关系建立平差的函数模型和随机模型。

步骤4：模型解算及精度验证：采用迭代最小二乘的方法进行求解，通过残差进行重新定权，直至结果符合精度要求为止。

步骤1包括如下几个子步骤：

步骤a：干涉对构网

首先，在对所有影像配准后，采用二连接的干涉对构网方式(如附图1所示，即：1和2组成干涉对，1和3组成干涉对，2和3组成干涉对)，这样任何一景影像都处于一个三角形闭合环中，既克服了时间失相干，又可以检验其中的误差；

步骤b：InSAR干涉相位计算得到LOS向形变序列结果

首先，用外部DEM数据通过地理编码得到SAR坐标系下的DEM，从而得到模拟地形相位，然后用生成的干涉图减去模拟干涉相位去除地形相位。为了消除相位残差，提高干涉图的信噪比，对差分干涉相位进行滤波后相位解缠，对解缠后的差分干涉图的基线进行精化，从而有效的去除了轨道相位。然后，通过不同滤波方式去除大气相位和噪声，最后将形变相位转化为LOS向形变序列结果。

步骤3包括如下几个子步骤：

该步骤建立数学模型，其包括函数模型和随机模型。

步骤a：建立二维形变序列分解函数模型

考虑到升降轨数据成像的几何差异，我们分两种情况来建立数学关系。

针对升轨数据：

如图2所示，假定卫星以升轨方向飞行，由于LOS向形变结果为竖值向形变结果和水平向形变结果的投影之和，故可得：

d_los1＝d_u×cosθ₁+(-d_h×sinθ₁) (1)

式中，d_u为竖值向形变在LOS向投影，d_h为水平向形变在LOS向投影，θ为雷达的局部入射角；同时假定形变靠近卫星为正，远离卫星为负。

如附图3所示，将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献(即x方向形变为0)，那么：

d_h＝d_y×sinγ_asc (2)

式中，γ_asc＝270°-(α₁-α')，其中，α₁为卫星飞行方向与北方向的夹角，α'为x方向与北方向的夹角。将(2)式代入(1)式中，可得：

d_los1＝d_u×cosθ₁+d_y×cos(α₁-α')sinθ₁ (3)

针对降轨数据：

如图4所示，假定卫星以降轨方向飞行，可得：

d_los2＝d_u×cosθ₂+(-d_h×sinθ₂) (4)

同理，如图5所示，将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献(即x方向形变为0)，那么：

d_h＝d_y×sinγ_des (5)

式中，γ_des＝(α₂-α')-90°。

将(5)式代入(4)式中，可得：

d_los2＝d_u×cosθ₂+d_y×cos(α₂-α')sinθ₂ (6)

因此，我们可以利用多个轨道或升降轨数据，就可以得到竖值向和沿任意方向水平向的形变结果：

步骤b：建立二维形变序列分解随机模型

假设不同角度的LOS向观测值对应的观测时间相同或者对应的时间间隔可以忽略，那么就可以利用步骤a中式子得到二维形变，但是实际上，忽略不同数据采集的时间间隔显然会带来误差，从而不利于我们对形变的解译。所以，我们将保留不同数据时间间隔不同，进而对二维形变序列进行分解。

假设三个不同几何的数据分别有N₁、N₂、N₃期，N＝N₁+N₂+N₃。数据获取的时间分别为

T＝T₁∪T₂∪T₃。通过时序InSAR的方法(如SBAS-InSAR)得到相应的LOS向形变序列

d＝d₁∪d₂∪d₃。通过多个方向的LOS形变分解出二维形变序列。

根据步骤a中公式，我们可以得到：

Gx＝d (8)

其中，设计矩阵:

G＝diag(A⁽¹⁾ A⁽²⁾ A⁽³⁾)，

其中：

A⁽¹⁾＝(cos(α⁽¹⁾-α')sinθ⁽¹⁾cosθ⁽¹⁾)；A⁽²⁾＝(cos(α⁽²⁾-α')sinθ⁽²⁾cosθ⁽²⁾)，A⁽³⁾＝(cos(α⁽³⁾-α')sinθ⁽³⁾cosθ⁽³⁾)；对于G矩阵某行中，某列是否为0与该轨道是否有观测值有关；x＝(p₁ p₂ ... p_N)^T，其中

y,u分别是任意水平方向形变值和竖直方向形变值；我们假设时间均不重叠且

那么观测量个数为N，待求参数个数为2N。

显然，式(8)秩亏。因为第一期LOS向观测值为0，那么分解得到的水平方向和竖直方向也应为0，为此，我们假设y₁＝0，u₁＝0，在方程解算中，我们以增加虚拟观测的方式对其进行约束；

Cx＝0

其中，

同时为了减轻方程病态并使其满秩，我们假设相邻位移序列的形变速率相等，即

i＝1,2..N-2。那么，我们就可以通过式(9)进行求解：

Rx＝0 (9)

其中，

把三个式子合并观测方程

从而利用最小二乘估计，解算方程。

步骤4包括如下几个子步骤：

步骤a：模型解算

通过迭代最小二乘求解：若

则满足解算的精度要求，迭代结束。

其中，

分别为两次的残差平方和。

步骤b：精度验证

利用式(11)得到的参数估计值

将其代入

中得到整个模型的解算残差，利用残差的均方根误差

来评定模型参数解算的内符合精度。

本发明在具体实施时，步骤如下：

步骤1：数据预处理

选取第一景影像为主影像，利用强度互相关配准的方法，选取1024个同名点，将所有的从影像与主影像进行配准，配准精度优于0.001像元；为了降低斑点噪声的影响，对数据进行多视处理，其中，数据多视为5：1，进而得到了15m分辨率的结果；为了克服时间去相干的影响，我们采用二连接的干涉对构网方式(如图1所示，即：1和2组成干涉对，1和3组成干涉对，2和3组成干涉对)，这样任何一景影像都处于一个三角形闭合环中，既克服了时间失相干，又可以检验其中的误差；我们用30mx30m的SRTM数据(图6)通过地理编码得到SAR坐标系下的DEM，从而得到模拟地形相位，然后用生成的干涉图减去模拟干涉相位去除地形相位。为了消除相位残差，提高干涉图的信噪比，我们对差分干涉相位进行了Goldstein滤波，滤波因子为0.5，滤波窗口大小为32×32。低相干区域可能会受到相位噪声的影响，所以我们在解缠时先生成平均的相干系数图，然后让相干性低于0.25的区域不参与解缠，从而保证了解缠的可靠性。在进行相位解缠时，我们在一块稳定区域相同位置处选择解缠的参考点，在利用最小费用流进行解缠后，我们对解缠后的差分干涉图的基线进行精化，从而有效的去除了轨道相位。通过上述过程我们得到了精确的解缠后的差分干涉相位，进而将其转化为LOS向形变序列结果。

步骤2：时空基准统一

为了统一空间基准，我们对得到的LOS向形变序列结果进行反地理编码，并且重采样到相同的地理坐标范围，从而得到15m分辨率的数据LOS向的形变序列结果；为了统一时间基准，我们选择升轨道数据的第一期影像的获取日期为统一的时间参考，然后将不同数据时间进行平移；

步骤3：建立数学模型

根据外部DEM数据得到研究区域的坡向(图7，规定北方向为0°)，然后根据式(1)到(11)沿坡向和竖直方向进行LOS向形变序列分解，进而得到任意水平方向和竖直方向的形变速率图，结果见图8和图9；

步骤4：模型解算及精度验证

利用残差的RMSE对结果进行精度验证，结果表明，本发明精度满足实际生产生活要求。

Claims (6)

1.一种任意方向的二维形变的监测方法，其特征在于，包括以下步骤：

步骤1：对数据进行数据预处理：即对获取的多种数据或者升降轨数据分别进行影像配准，干涉对构网，利用外部DEM数据去除模拟地形相位，得到差分干涉相位后，对差分干涉相位进行滤波，解缠，得到各数据LOS向形变序列结果；

步骤2：对预处理的数据进行时空基准统一；

步骤3：对时空基准统一的数据进行建立数学模型：即根据LOS向数据的形变结果建立观测方程，

通过严密的几何关系建立平差的函数模型和随机模型；

步骤4：对建立数学模型的数据进行模型解算及精度验证；即采用迭代最小二乘的方法进行求解，通过残差进行重新定权，直至结果符合精度要求为止；

所述步骤3中

对于不同数据时间基准不统一的情况，提出保留不同数据时间间隔不同，其中，针对方程秩亏的现象，假设相邻位移序列的形变速率相等；

增加虚拟观测的方式对其进行约束，进而对形变序列直接进行分解；

假设三个不同几何的数据分别有N₁、N₂、N₃期，N＝N₁+N₂+N₃，数据获取的时间分别为

T＝T₁∪T₂∪T₃；

通过时序InSAR的方法得到相应的LOS向形变序列

d＝d₁∪d₂∪d₃；

通过多个方向的LOS形变分解出二维形变序列；

根据竖直向和沿任意方向水平向的形变结果，得到：

Gx＝d

其中，设计矩阵:

G＝diag(A⁽¹⁾ A⁽²⁾ A⁽³⁾)，

其中：

A⁽¹⁾＝(cos(α⁽¹⁾-α')sinθ⁽¹⁾ cosθ⁽¹⁾)；A⁽²⁾＝(cos(α⁽²⁾-α')sinθ⁽²⁾ cosθ⁽²⁾)，A⁽³⁾＝(cos(α⁽³⁾-α')sinθ⁽³⁾cosθ⁽³⁾)；

x＝(p₁ p₂ ... p_N)^T，其中

y,u分别是任意水平方向形变值和竖直方向形变值；

假设时间均不重叠且d₁ ⁽¹⁾，d₁ ⁽²⁾，d₁ ⁽³⁾≠0，观测量个数为N，待求参数个数为2N；

假设y₁＝0，u₁＝0，在方程解算中，以增加虚拟观测的方式对其进行约束；

Cx＝0

其中，

为了减轻方程病态并使其满秩，假设相邻位移序列的形变速率相等，即

i＝1,2..N-2；

求解：

Rx＝0

其中，

把三个式子合并观测方程

利用最小二乘估计，解算方程。

2.根据权利要求1所述的一种任意方向的二维形变的监测方法，其特征在于，

所述步骤1的干涉对构网具体为：

在对所有影像配准后，采用二连接的干涉对构网方式，即：1和2组成干涉对，1和3组成干涉对，2和3组成干涉对，用于任何一景影像都处于一个三角形闭合环中，有利于克服了时间失相干，同时又可以检验其中的误差。

3.根据权利要求1所述的一种任意方向的二维形变的监测方法，其特征在于，所述步骤2包括以下子步骤：

步骤21：利用外部DEM数据将LOS向形变结果地理编码至地理坐标系下；

步骤22：将不同数据重叠区域进行裁剪，进而重采样至统一地理格网中；

步骤23：对不同数据选取相同时间基准作为起始的时间参考，将形变结果进行平移，达到时间基准的统一。

4.根据权利要求1所述的一种任意方向的二维形变的监测方法，其特征在于，所述步骤3包括以下子步骤：

步骤31：根据InSAR的理论知识，建立严密的数学关系式，从而对LOS向形变进行任意方向二维形变分解；

5.根据权利要求4所述的一种任意方向的二维形变的监测方法，其特征在于，

升轨数学模型的求解步骤：

步骤311：假定卫星以升轨方向飞行，由于LOS向形变结果为竖直向形变结果和水平向形变结果的投影之和，可得：

d_los1＝d_u×cosθ₁+(-d_h×sinθ₁)

式中，d_u为竖直向形变在LOS向投影，d_h为水平向形变在LOS向投影，θ₁为升轨方向的雷达局部入射角；同时假定形变靠近卫星为正，远离卫星为负；

步骤312：将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献，即x方向形变为0，那么：

d_h＝d_y×sinγ_asc

式中，γ_asc＝270°-(α₁-α')，其中，α₁为卫星飞行方向与北方向的夹角，α'为x方向与北方向的夹角；

步骤313：将步骤312中的公式代入步骤311式中，可得：

d_los1＝d_u×cosθ₁+d_y×cos(α₁-α')sinθ₁

降轨数学模型：

步骤321：假定卫星以降轨方向飞行，可得：

d_los2＝d_u×cosθ₂+(-d_h×sinθ₂)

式中，θ₂为降轨方向的雷达局部入射角；

步骤322：将水平形变投影至NE平面，定义新的坐标系为xoy，假定水平方向的形变仅由y方向贡献即x方向形变为0，那么：

d_h＝d_y×sinγ_des

式中，γ_des＝(α₂-α')-90°

步骤323：将步骤322式代入步骤321式中，可得：

d_los2＝d_u×cosθ₂+d_y×cos(α₂-α')sinθ₂

利用多个轨道数学模型或升降轨数学模型，得到竖直向和沿任意方向水平向的形变结果：

6.根据权利要求1所述的一种任意方向的二维形变的监测方法，其特征在于，所述步骤4中通过迭代最小二乘求解具体过程为：

步骤41：若

则满足解算的精度要求，迭代结束；

其中，

分别为两次的残差平方和；

步骤42：利用观测方程得到的参数估计值

将其代入

中得到整个模型的解算残差，利用残差的均方根误差

来评定模型参数解算的内符合精度。

Images