基于分形布朗运动的干涉SAR影像仿真方法与流程

本发明涉及雷达图像处理技术领域，特别涉及一种基于分形布朗运动的干涉sar影像仿真方法。

背景技术：

合成孔径雷达干涉(insar)技术已经成为了最为有效的全球测图手段之一。现阶段，srtm提供了全球56°s～60°n范围内的数字高程模型(dem)数据，其标称精度达到了16m，格网大小30米，在全球尺度的科研以及工程应用中发挥了不可估量的作用。在2016年，德宇航宣布完成了全球dem的测绘，90％地区的精度为3.5米，如果不考虑南极地区，全球90％地区的精度则达到了1.3米，数据已经可以用于精细的军事化定位和军事目标识别。insar地形测绘的巨大潜力，使得insar技术在工程应用上迎来了一个新的时代。然而，我国干涉sar卫星尚不具备业务化运行能力，数据匮乏。现阶段我国可用的sar卫星只有高分三号，但高分三号并没有针对干涉性能进行专门的设计，因此干涉性能较差。先期实验结果表明，在接近300景的高分三号影像中只能够形成不到10对干涉对。因此，使用模拟仿真的方法获取干涉对，并基于此进行insar地形测绘相关技术研究，是进行国产sar卫星指标论证、数据处理以及测量检校研究的重要途径。

现阶段干涉影像仿真过程中多采用真实的dem数据，例如srtm数据，完成干涉对数据的生成。这种仿真方法虽然能够在一定程度上满足干涉测量的基本要求，然而srtm的dem数据最高空间分辨率只有30米，无法满足高空间分辨率的仿真需求，而更高分辨率的dem数据无法免费获取，仿真成本极高。

分形布朗运动(fractionalbrownianmotion,fbm)是一种二维随机运动，这种随机运动具有自相似性以及长距离相关性。同样的，在地学研究过程中，所有的地学信息在一定的距离内都有较好的自相似性。特别是对于地形来说，虽然较小尺度内的地形差异较大，然而在宏观尺度上来看，地形的相似性极高，例如山地地区一般可绵延数十乃至数百公里，而我国东部平原地区则可以绵延数千公里。地形的这种分布特征与fbm的特征极为一致。

fbm在sar影像中主要用于后向散射系数的仿真，或者进行sar影像的海波分析，而极少用于干涉影像仿真过程。fbm是纯数学模型，能够以定量的方式表征insar处理过程中的高程在方位向和距离向的统计特征，现有技术在干涉影像仿真过程中尚未用到fbm。本发明则从仿真层面出发，以fbm模型为基础，开展地形仿真，几何参数仿真，干涉相位仿真，后向散射系数仿真，干涉影像仿真，提供可用于干涉处理的干涉对，完成影像的全链路仿真分析。本发明拟为国产sar卫星的立项发射提供仿真依据。

技术实现要素：

本发明提供了一种基于分形布朗运动的干涉sar影像仿真方法，所述方法包括以下步骤：

步骤s1，基于分形布朗运动(fbm)仿真sar影像的外接dem；

步骤s2，使用在步骤s1仿真的所述sar影像的外接dem仿真所述sar影像的相位；

步骤s3，使用在步骤s1仿真的所述sar影像的外接dem仿真所述sar影像的振幅；

步骤s4，基于步骤s2仿真得到的sar影像的相位和步骤s3仿真得到的sar影像的振幅，仿真所述sar影像。

优选的，步骤s1具体包括以下子步骤：

步骤s1.1，从现有sar影像中选取一个干涉影像对，指定其中一景影像为主影像，另一景影像为从影像，并提取所述主影像和所述从影像的成像几何参数；

步骤s1.2，构建斜距-多普勒(r-d)模型，并基于所述主影像的成像几何参数和构建的r-d模型，确定所述主影像对应的外接dem参数；

步骤s1.3，使用分形布朗运动(fbm)仿真所述主影像的外接dem。

优选的，步骤s1.2具体包括以下子步骤：

步骤s1.2.1，构建斜距方程。如下所示：

其中，r0是近地点斜距，△r是斜距分辨率，(xs,ys,zs)是成像时刻的卫星坐标，(xp,yp,zp)是地面点的三维坐标，x为地面点在所述主影像中的距离向像素坐标；

步骤s1.2.2，构建多普勒方程，如下所示：

其中，(vxsvysvzs)为成像时刻的卫星速度，(xs,ys,zs)为成像时刻的卫星坐标，(xp,yp,zp)是地面点的三维坐标，fd为多普勒效应值，λ为雷达波长，r0是近地点斜距，△r是斜距分辨率，x为地面点在所述主影像中的距离向像素坐标；

步骤s1.2.3，构建椭球方程，如下所示：

其中，(xp,yp,zp)是地面点的三维坐标，h为地面点的高程，re与rp分别为椭球的赤道半径与极半径；

步骤s1.2.4，利用步骤1.2.1-1.2.3构建的斜距方程、多普勒方程和椭球方程，解算所述主影像的角点坐标；

步骤s1.2.5，根据所述主影像的角点坐标，确定所述主影像的外接dem参数。

优选的，步骤s1.3具体包括以下子步骤：

步骤s1.3.1，计算能够提供满足fbm特征的二维矩阵的循环分块矩阵；

步骤s1.3.2，步骤s1.3.2，计算与所述循环分块矩阵有相同协方差的随机场；

步骤s1.3.3，对在s1.3.2中获取的随机场进行改正，得到最终的sar影像的外接dem。

优选的，步骤s1.3.2具体包括以下子步骤：

步骤s1.3.2.1，计算所述循环分块矩阵的特征值；

步骤s1.3.2.2，生成二维复随机场；

步骤s1.3.2.3，将步骤s1.3.2.1获取的特征值构成特征值矩阵，并与步骤s1.3.2.2生成的所述二维复随机场相乘，计算得到与所述循环分块矩阵有相同协方差的随机场；

步骤s1.3.2.4，对s1.3.2.3获取的随机场进行傅里叶变换，并抽取傅里叶变换后的随机场的正频部分的数据；

步骤s1.3.2.5，使用s1.3.2.4中得到的傅里叶变换的正频数据解算独立随机场，得到与所述正频数据的实部和虚部分别对应的两个独立随机场；

步骤s1.3.2.6，将所述两个独立随机场的数据进行拉伸变换，得到拉伸变换后的随机场。

优选的，步骤s1.3.3具体包括以下子步骤：

步骤s1.3.3.1，为所述随机场添加增量平稳特性；

步骤s1.3.3.2，将步骤s1.3.3.1获得的具有增量平稳特性的随机场的数据进行拉伸变换；

步骤s1.3.3.3，将步骤s1.3.2.6获得的随机场与步骤s1.3.3.2获得的随机场相加；

步骤s1.3.3.4，将步骤s1.3.3.3获得的随机场拉伸到指定高程范围。

优选的，步骤s2具体包括以下子步骤：

步骤s2.1，使用所述sar影像的外接dem仿真所述主影像和所述从影像的干涉相位；

步骤s2.2，使用随机噪声模型仿真从影像的相位；

步骤s2.3，使用步骤s2.1得到的所述从影像的干涉相位和步骤s2.2得到的所述从影像的相位，计算所述主影像的相位；其中：

所述主影像的相位的表达式为

其中，为步骤s2.2获得的从影像的相位，为主影像像素空间中每个像素点的干涉相位，为噪声相位。

优选的，步骤s2.1具体包括以下子步骤：

步骤s2.1.1，将步骤s1仿真得到的dem转换到sar影像空间；

步骤s2.1.2，使用通过步骤s2.1.1转换后的dem，计算sar影像空间中每个像素点的垂直基线和平行基线；

步骤s2.1.3，使用步骤s2.1.2得到的垂直基线和平行基线计算干涉相位。

优选的，步骤s3具体包括以下子步骤：

步骤s3.1，仿真所述主影像的振幅；

步骤s3.2，仿真所述从影像的振幅；

其中，步骤s3.1具体包括以下子步骤：

步骤s3.1.1，计算所述主影像像素空间每个像素点的本地法线向量。

本地法线向量使用如下公式计算：

其中，为像素点(1,y)对应的视线向向量，为像素点(0,y)对应的视线向向量，为像素点(1,y-1)对应的视线向向量；

步骤s3.1.2，基于所述主影像像素空间中像素点的本地法线向量，判断所述像素点是否是平地像素；具体包括：

本地法线向量与本地坐标向量的夹角可表示为χ，其表达式为：

其中，为本地法线向量的单位向量，为像素点坐标矢量的单位向量；优选的，当χ<1°时判断所述像素点为平地像素。

步骤s3.1.3，计算每个像素点的名义入射角β，具体如下：

式中，rh为卫星到地心的距离，r1为主影像斜距，表达为r1＝r0+x△r，其中r0为主影像近地点斜距，x为像素点的x方向坐标，△r为影像的斜距分辨率，re为地球曲率半径，h为地面点高程。

步骤s3.1.4，基于对所述像素点是否为平地像素的判断和像素点的名义入射角，计算每个像素点的本地入射角θ，具体如下：

其中，β为步骤s3.1.3中得到的名义入射角，为卫星速度的单位矢量，为本地法线向量的单位矢量；

步骤s3.1.5，基于对所述像素点是否为平地像素的判断和像素点的本地入射角，计算每个像素点的振幅值am，具体如下：

其中，θ为步骤s3.1.4中得到的每个像素点的本地入射角，μ为像素尺寸归一化参数，表达式如下：

其中，χ为本地法线向量与本地坐标向量的夹角，θ为步骤s3.1.4中得到的每个像素点的本地入射角。

优选的，步骤s3优选的还包括以下子步骤：

步骤s3.1.6，计算叠掩区振幅值；

步骤s3.1.7，计算阴影区振幅值。

优选的，步骤s4包括以下子步骤：

步骤s4.1，仿真所述主影像；

对主影像每个像素来说，其像素值表达为一个复数，即

f(x,y,m)＝am+bmi

其中，f(x,y,m)为主影像每个像素的像素值，(x,y)为主影像像素的坐标，m为主影像像素的复数值，am为步骤s3.1仿真的主影像振幅，为步骤2.3仿真的主影像的相位；

步骤s4.2，仿真所述从影像；

对从影像每个像素来说，其像素值表达为一个复数，即

f(x,y,s)＝as+bsi

其中，f(x,y,s)为从影像每个像素的像素值，(x,y)为为从影像像素的坐标，s为从影像像素的复数值，as为步骤s3.2仿真的从影像振幅，为步骤2.2仿真的从影像相位。

本发明仿真过程中采用的参数为真实的tandem-x卫星干涉参数，参数具有极好的工程应用性，能够确保在干涉处理过程中，所有的误差来源不会与干涉几何状态有太大关联，从而使得研究人员能够基于本发明得到更为精确的误差仿真结论。

附图说明

图1为根据本发明实施例的基于分形布朗运动的sar干涉影像仿真方法的流程图。

图2(a)和图2(b)为根据本发明实施例的基于分形布朗运动仿真得到的干涉sar影像的外接数字高程模型fbm结果，其中，图2(a)为实部随机场改正之后得到的结果，图2(b)为虚部随机场改正之后得到的结果。

图3为根据本发明实施例的insar在地形测绘中的观测几何示意图。

图4.(a)为根据本发明实施例的基于分形布朗运动仿真得到的干涉sar影像的外接dem的主影像的四个角点连接得到的矩形；图4(b)为投影后的sar主影像坐标空间高程图。

图5(a)为根据本发明实施例的仿真的总相位；图5(b)为根据本发明实施例的仿真的地形相位。

图6(a)为根据本发明实施例的仿真的从影像相位，图6(b)为根据本发明实施例的仿真的主影像相位。

图7为根据本发明实施例的平地像素判定方法的示意图。

图8(a)为所述主影像振幅仿真结果，图8(b)为从影像振幅仿真结果。

图9(a)为主影像的相位和振幅叠加的结果，图9(b)为从影像的相位和振幅叠加的结果。

具体实施方式

下面对本发明进行更为详细的描述，描述过程中参照了上述附图，并使用示例性实施例进行了详解。

针对背景技术中提出的问题，本发明提出了一种使用分形布朗运动(fractionalbrownianmotion,fbm)的干涉sar影像仿真方法。这种方法采用分形布朗运动模型确定dem，并采用干涉几何参数获取影像干涉对，对于影像模拟仿真来说具有极强的实用性。

本发明的主要目的是使用分形布朗运动(fractionalbrownianmotion,fbm)仿真dem数据，并利用干涉几何参数，完成干涉影像的仿真。仿真结果包含两景可干涉的sar影像。

如图1所示，第一，本发明首先仿真dem，仿真过程中，使用主影像几何参数确定dem覆盖范围，以及dem的格网大小，从而确定主影像的外接dem的基本属性，并使用fbm仿真外接dem。随后再次使用主影像的几何参数，将外接dem进行精确的编码，投影到slc影像空间中，使得dem数据与slc数据一一对应。第二，使用dem仿真影像相位。仿真过程中，认为从影像中只包含相位噪声，从而使用高斯白噪声确定从影像相位。并联合从影像相位、噪声、平地相位以及地形相位仿真主影像相位。第三，仿真slc影像振幅。仿真过程中，需要用到fbm生成的dem数据以及主、从影像的几何参数。第四，将仿真的相位和振幅合成主从slc影像。第五，将主从slc影像进行常规干涉处理，获取dem数据，与仿真的dem进行交叉验证，并给出对应的精度报告。

本发明提出的基于分形布朗运动的干涉sar影像仿真方法具体包括以下步骤。

步骤s1，基于分形布朗运动(fbm)仿真干涉sar影像的外接数字高程模型(dem)。

首先，从现有sar影像中选取一对干涉数据，随意指定其中一景影像作为参考影像，即主影像；另一景影像作为从影像。其次，确定所述主影像的外接dem参数，包括所述主影像的外接dem的左上角经纬度、经纬度分辨率、行列数、高程范围等。最后，使用fbm仿真外接dem。

具体来说，步骤s1包括以下步骤：

步骤s1.1，从现有sar影像中选取一个干涉影像对，指定其中一景影像为主影像，另一景影像为从影像，并提取所述主影像和所述从影像的成像几何参数。需要说明的是，主、从影像可以由本领域技术人员根据实际需要随意指定。

根据本发明的优选实施例，选取四川甘孜藏族自治州理塘县的tandem-x的干涉影像对，干涉对的成像时间为2013年9月12日，指定其中的tdx-1传感器对应的影像为主影像，tsx-1传感器对应的影像为从影像。提取两景影像的几何参数，主、从影像的成像几何参数分别如表1和表2所示。

表1.主影像的成像几何参数

表2.从影像的成像几何参数

步骤s1.2，构建斜距-多普勒(r-d)模型，并基于所述主影像的成像几何参数和构建的r-d模型，确定所述主影像对应的外接dem参数。

具体的，采用步骤s1.1中给出的主影像的成像几何参数(见表1)，结合sar成像过程中常用的r-d模型确定主影像四个角点的坐标。r-d模型是进行sar影像定位的重要方程，它通过使用s1.1中提供的斜距、多普勒参数，结合地面点的椭球参数，建立某一点p的像素坐标(x,y)与地面坐标(xp,yp,zp)之间的对应关系。r-d模型共包含三个子模型，即斜距方程、多普勒方程以及椭球方程。步骤s1.2具体包括以下步骤：

步骤s1.2.1，构建斜距方程。

斜距方程表达的是，成像时刻的卫星坐标与地面坐标之间的距离等于头文件中获取的斜距。斜距方程如下所示：

其中，r0是近地点斜距，△r是斜距分辨率，(xs,ys,zs)是成像时刻的卫星坐标，(xp,yp,zp)是地面点的三维坐标，x为地面点在所述主影像中的距离向像素坐标。

各参数的解算具体如下：

其中，近地点斜距r0和斜距分辨率△r是已知量，分别是表1中的主影像的近地点斜距和距离分辨率。

卫星的三维坐标(xs,ys,zs)具体计算过程如下：

卫星的三维坐标是采用卫星的n个状态矢量插值而成。首先使用n个状态矢量构建卫星位置与卫星成像时刻的多项式关系，优选的，采用4阶多项式，以成像时刻的卫星坐标x轴的状态矢量为例，其与卫星成像时刻之间的关系可表达为：

其中，(t1t2…tn)为n个状态矢量对应的时刻，(a1a2a3a4a5)为多项式系数，为卫星成像时刻的n个状态矢量的x轴分量，在s1.1中，表1给出了19个主影像的状态矢量，因此n＝19。那么在方程(2)中，状态矢量的x轴分量、19个状态矢量对应的时间均为已知量，仅有(a1a2a3a4a5)为未知量。将表1的参数带入公式(2)中，令

那么依据最小二乘原理，即可得到

同样的，将成像时刻的卫星坐标y轴的状态矢量带入公式(2)–(5)，可以得到y方向的拟合参数

以及z方向拟合参数

y为地面点在所述主影像中的方位向像素坐标，y与成像时刻t之间存在下述关系

t＝t0+△t·(y-1)(8)

其中，t0为方位向成像起始时间，即步骤s1.1中，表1给出的方位向起始时间，△t为方位向成像时间间隔，即步骤s1.1中，表1给出的方位向单行时间。

由此，构建的斜距方程为：

式(7)即公式(1)、(2)、(6)的实施例，式中各变量的含义与公式(1)、(2)、(6)相同。至此，斜距方程中未知量只包含两类，即地面坐标(xp,yp,zp)以及像素坐标(x,y)。

步骤s1.2.2，构建多普勒方程。

多普勒方程描述的是卫星与地面的相对运动带来的多普勒效应。成像时刻的卫星速度(vxsvysvzs)与卫星侧视方向矢量((xs-xp)(ys-yp)(zs-zp))之间的叉乘代表了对应的多普勒效应。

构建的多普勒方程的具体表达式为：

其中，(vxsvysvzs)为成像时刻的卫星速度，(xs,ys,zs)为成像时刻的卫星坐标，(xp,yp,zp)是地面点的三维坐标，fd为多普勒效应值，λ为雷达波长，r0是近地点斜距，△r是斜距分辨率，x为地面点在所述主影像中的距离向像素坐标。

在星载sar处理过程中，会通过星上偏航导引以及地面的成像，保证多普勒效应为0，称为零多普勒效应。因此在本发明中，方程右侧的值恒为0。在上述方程中，未知量仅包含成像时刻的卫星速度。

成像时刻的卫星速度解算过程如下。

步骤s1.2.2.1，计算卫星的速度矢量。

卫星速度矢量是采用卫星的n个速度矢量插值而成。首先使用n个速度矢量构建卫星速度与卫星成像时间的多项式关系，优选的，采用4阶多项式。其构建方式与步骤s1.2.1中计算卫星的三维坐标的过程完全一致，此处不再赘述。以主影像为例，采用与公式(2)–(7)中成像时刻的卫星位置相同的解算方法，可得到成像时刻的卫星速度与成像时间之间的关系：

式中，(vxsvysvzs)为成像时刻的卫星速度，y是地面点在所述主影像中的方位向像素坐标，t是成像时刻。

步骤s1.2.2.2，构建多普勒方程。

式中，(vxsvysvzs)为成像时刻的卫星速度，(xs,ys,zs)为成像时刻的卫星坐标，(xp,yp,zp)是地面点的三维坐标。

至此，多普勒方程中未知量只包含两类，即地面坐标(xp,yp,zp)，以及地面点在所述主影像中的方位向和距离向像素坐标(x,y)。

步骤s1.2.3，构建椭球方程。

构建的椭球方程为：

其中，(xp,yp,zp)是地面点的三维坐标，h为地面点的高程，re与rp分别为椭球的赤道半径与极半径，在wgs84坐标系统中，赤道半径为6378137.00m，极半径为6356752.31m。因此严格来讲，如果想要获取地面点的精确坐标，需确保地面点高程为已知量。在本实施例中，精确的高程信息是使用步骤s1.3中的仿真过程获取的。使用仿真获取精确的高程信息之后，椭球方程中的未知参数只有地面坐标。

步骤s1.2.4，利用步骤1.2.1-1.2.3构建的斜距方程、多普勒方程和椭球方程，解算主影像的角点坐标。

影像四个角点的像素坐标对应四个不同的地面坐标。坐标解算过程中需要用到步骤s1.2.1、s1.2.2以及s1.2.3中的三个方程。在方程中，未知量共包含三类，即像素坐标、地面点坐标、地面点高程。在此步骤中，需要确定像素坐标、地面点高程，以解算地面点坐标。

仿真所用数据在四川省理塘县，此处地势起伏较大，高程在3500-6000米范围。精确的高程需要在步骤s1.3中通过仿真获得，此处使用概要的高程值，优选的，取3500米，左上角坐标为(0,0)。将坐标与高程代入方程(9)(12)(13)中，即可解得对应的左上角的坐标为(-947867.6,5436136.9,3194792.1)，将此坐标转换为对应的经纬度，为(99.89089°e,30.23596°n)。

同样的，设定高程为3500米，陆续解得其余三个角点坐标对应的经纬度坐标为(100.05818°e,30.21120°n)(100.027203°e,30.05569°n)(99.86016°e,30.08055°n)。

步骤s1.2.5，根据所述主影像的角点坐标，确定所述主影像的外接dem参数。

由步骤s1.2.4可知，主影像覆盖的范围为99.86016°e–100.05818°e，30.05569°n–30.23596°n。外接dem范围应比影像范围稍大，优选的，设定如下外接dem参数：

表3.外接dem参数

步骤s1.3，使用分形布朗运动(fbm)仿真所述主影像的外接dem。

一个标准的分形布朗运动定义在某概率空间上，指数为h(0<h<1)的一个随机过程b，满足如下条件：

1、以概率1满足b(0,0)＝0，且b(x,y)为(x,y)的连续函数；

2、对任意变量(x,y)，(h,k)∈r²，其二维增量b(x+h,y+k)-b(x,y)服从期望为0，方差为(h²+k²)^h的正态分布，其概率满足

其中，(x,y)为概率密度分布函数的自变量，(h,k)为自变量的增量，p为概率密度函数，b为布朗运动过程，h为hurst参数，exp为自然指数，r为积分参数，代表从负无穷到z的积分。

反之亦然，凡是符合以上条件的，均可认为是fbm。实现fbm的方法包括hosking算法，cholesky算法，daviesandharte算法等。优选的，本实施例采用daviesandharte算法。

步骤s1.3具体包括以下子步骤：

步骤s1.3.1，计算能够提供满足以上提到的2个fbm条件的二维矩阵的循环分块矩阵。

计算循环分块矩阵的目的是寻找一个方阵，使得方阵可提供满足fbm条件的二维矩阵。一般来说，循环分块矩阵c的表达形式为：

其中，γ(·)是方阵的协方差函数，m为要仿真的dem的行数s和列数l之间的较大值，在本实施例中，表3给出的dem的行数和列数均为8000。每个循环分块矩阵c均可以进行如下分解：

c＝qλq^*(16)

其中，q为单位矩阵，q^*为单位矩阵q的转置的复共轭，λ为循环分块矩阵c的特征值构成的对角阵，单位矩阵q中的任意元素定义为：

其中，qj,k为单位矩阵q中坐标为(j,k)的元素，对角阵λ中的任意元素定义为：

其中，λk为循环分块矩阵c的第k个特征值，rj为循环分块矩阵c第一行的第j+1个元素，exp为自然指数，m为要仿真的dem的行数s和列数l之间的较大值，在本实施例中，表3给出的dem行数和列数均为8000。由于循环分块矩阵c为正定对称矩阵，因此特征值为正实数。如果某矩阵的特征值为特征向量与循环分块矩阵c相同，那么此矩阵依然为正定实矩阵。这样的矩阵可定义为

s＝qλ^1/2q^*(19)

此矩阵特征值将满足fbm的第二条特性，即二维增量符合正态分布，其中q为单位矩阵，q*为单位矩阵q的转置的复共轭，λ为为循环分块矩阵c的特征值构成的对角阵。

因此本步骤的核心是获取符合条件的循环分块矩阵。矩阵大小为预期大小的一半。由s1.2.5可知，此循环分块矩阵大小为4000×4000。为了构建循环分块矩阵，需首先进行行列坐标归一化，随后计算距离权重矩阵的参数，并以此构建单块矩阵、分块循环矩阵，最终解算循环分块矩阵的特征值。步骤s1.3.1具体包括以下子步骤：

步骤s1.3.1.1，将待仿真的单块矩阵行列坐标归一化。

本步骤的目的是避免解算过程中的数值过大，造成数据溢出。坐标归一化公式为：

其中，(tx,ty)是归一化之后的待仿真单块矩阵的行列坐标，r为拉伸因子，优选的，可令r＝2。s为待仿真的单块矩阵影像列数，l为影像行数，x为待仿真的单块矩阵的列坐标，y为待仿真的单块矩阵的行坐标。在本实施例中，s＝4000，l＝4000。

步骤s1.3.1.2，计算距离权重矩阵的参数。

距离权重矩阵参数，是进行步骤s1.3.2.3中单块距离权重矩阵的输入参数，用以表现fbm的距离自相似性，详细的距离权重矩阵参数描述如下。

首先，确定fbm的hurst参数。hurst参数的范围为0<h<1，h越小，地形差异越小，h越大，地形差异越大，当h为0.5时，为普通的布朗运动。优选的，令h＝0.7。

其次，确定fbm的其他参数，

当h≤0.75时，

当h>0.75时，

其中，h为hurst参数，β、α、c1、c2为中间参数，无实际的数学意义，只是为了方便进行后续的公式描述。在本实施例中，h＝0.7，对应的β＝0,α＝0.7,c1＝0.3,c2＝0.7。这4个参数将陆续在步骤s1.3.1.3中得到具体使用。

步骤s1.3.1.3，使用步骤s1.3.1.2中的参数计算单块矩阵。单块矩阵是循环分块矩阵的基本构成，是步骤s1.3.1.4的输入参数。

计算单块矩阵的目的是提供具有距离自相似性的循环分块矩阵。具体过程描述如下。

引入各向同性函数r，其表述为

其中，其中(x,y)为待仿真的分块矩阵中每个元素的坐标，(x0,y0)为待仿真的单块矩阵中同性函数的起算点坐标，即(0,0)。待仿真单块矩阵中每个元素的距离权重函数表述为

其中c1、c2、α以及β为步骤s1.3.1.2中确定的距离权重参数。

步骤s1.3.1.4，计算循环分块矩阵。

将距离权重函数进行左右镜像和上下镜像，扩展之后的结果即为循环分块矩阵。扩展结果见下式右侧。

其中，ρsl为距离权重矩阵的第s列第l行的元素，s为距离权重矩阵的总列数，l为距离权重矩阵的总行数。在本实施例中，s＝4000，l＝4000。最终的循环分块矩阵为8000×8000的矩阵。

步骤s1.3.2，计算与所述循环分块矩阵有相同协方差的随机场。

这一步骤是为了求解公式(19)中的单位矩阵q，使得单位矩阵q满足公式(17)的定义。

步骤s1.3.2具体包括以下子步骤：

步骤s1.3.2.1，计算循环分块矩阵的特征值。

循环分块矩阵的特征值设定为λ，根据特征值的定义，λ满足如下公式

其中e为单位矩阵，c为循环分块矩阵，解算上述方程中的λ，即可得到矩阵的特征值，对于本实施例来说，矩阵具备8000个特征值，表示为(λ1λ2……λ8000)，特征值为步骤s1.3.2中的公式(17)中的矩阵q的创立提供前提条件。

步骤s1.3.2.2，生成二维复随机场。

二维复随机场fcr中的每个元素的表达式为：

fcr(j,k)＝a+bi(27)

其中，fcr(j,k)为二维复随机场中第j列第k行的元素，a为元素的实部，b为元素的虚部，在本实施例中，随机数矩阵大小为(8000×8000)。因此实部和虚部均可组成一个8000×8000的矩阵。

步骤s1.3.2.3，将步骤s1.3.2.1获取的特征值构成特征值矩阵，并与步骤s1.3.2.2生成的所述二维复随机场相乘，计算得到与所述循环分块矩阵有相同协方差的随机场。

将步骤s1.3.2.1获取的特征值构成特征值矩阵，并与所述二维复随机场相乘，得到具有循环分块矩阵特征的随机场，即

其中，fr为与循环分块矩阵有相同协方差的随机场，fcr为步骤s1.3.2.1中生成的二维复随机场，(λ1λ2……λ8000)为循环分块矩阵的特征值。

步骤s1.3.2.4，对s1.3.2.3获取的随机场fr进行傅里叶变换，并抽取傅里叶变换后的随机场fr的正频部分的数据。

在复数域进行fft变换，其表达式为

f(j,k)＝fft(fr)(j,k),j＝1,2,...,8000,k＝1,2,...,8000(29)

其中，f(j,k)为傅里叶变换之后，矩阵第j列第k行的元素，fft为快速傅里叶变换函数，fr为与循环分块矩阵有相同协方差的随机场，来自于步骤s1.3.2.2。将以上公式的矩阵形式拆分为元素形式，那么每个元素(j,k)的变换结果可表达为：

其中，f(j,k)为傅里叶变换之后，矩阵第j列第k行的元素，fr为与循环分块矩阵有相同协方差的随机场，来自于步骤s1.3.2.2。m为fr的矩阵大小，m与n是积分自变量，exp[]为指数函数，对于本实施例来说，m＝8000，对于上面的等式来说，傅里叶变换前后的数值均为复数，变换之后的矩阵，左下角为正频部分，右上角为负频部分，两者的频谱特征相同。我们只取左下角正频部分的数据进行后续的仿真分析。

步骤s1.3.2.5，使用s1.3.2.4中得到的傅里叶变换的正频数据解算独立随机场，得到与所述正频数据的实部和虚部分别对应的两个独立随机场；

在复数的fft变换之后，实频空间与复频空间是正交的，两者形成的随机场相互独立，因此s1.3.2.4中的正频部分数据的实部与虚部，分别构成两个独立随机场。

步骤s1.3.2.6，将步骤s1.3.2.5中得到的所述两个独立随机场的数据进行拉伸变换，得到拉伸变换后的随机场。

步骤s1.3.2.5获取独立随机场数据是频率域数据，值域范围极大，在计算过程中易造成数据溢出。因此需要将数据拉伸到[0,1]之间。步骤s1.3.2.6具体步骤包括：

s1.3.2.6.1，减掉所述独立随机场原点处的值，确保所述独立随机场的原点值为0，以符合步骤s1.3中描述的fbm的第一个特性。

s1.3.2.6.2，对所述两个独立随机场的数据分别进行数据拉伸。

拉伸公式为：

其中，f′(j,k)为拉伸之后第j列第k行的值，f(e,k)为拉伸之前第j列第k行的值，fmin为s1.3.2.6.1中经修正的独立随机场矩阵的最小值，fmax为矩阵的最大值。

步骤s1.3.3，对在s1.3.2中获取的随机场进行改正，得到最终的sar影像的外接dem。

步骤s1.3.3具体包括以下子步骤：

步骤s1.3.3.1，为所述随机场添加增量平稳特性。

随机场改正的目的是，使得fbm在噪声分形的基础上，具备增量平稳性。随机场改正的方法，是在随机场中的每个像素值中，添加具备增量平稳特性的数值，此数值可表达为：

其中，g为增量平稳矩阵，g(j,k)为增量平稳矩阵g的第j列第k行的元素，n1与n2分别为两个随机数，txj为步骤s1.3.1.1中待仿真的单块矩阵行列坐标归一化后的第j列的x坐标，tyk为步骤s1.3.1.1中待仿真的单块矩阵行列坐标归一化后的第k列的y坐标，c2为步骤s1.3.1.2中定义的距离权重矩阵参数。

步骤s1.3.3.2，将步骤s1.3.3.1获得的具有增量平稳特性的随机场的数据进行拉伸变换。具体的，将所述数据拉伸到[0,1]。

拉伸过程与步骤s1.3.2.6相同。此处不再赘述。

步骤s1.3.3.3，将步骤s1.3.2.5获得的随机场与步骤s1.3.3.2获得的随机场相加。

步骤s1.3.3.4，将步骤s1.3.3.3获得的随机场拉伸到指定高程范围。

由步骤s1.2.5可知，dem高程范围为3500–6000m。因此，将fbm结果拉伸到3500–6000m。即可得到主影像的外接dem。

优选的，选用虚部生成的fbm结果进行后续仿真。

仿真结果如图2(a)和图2(b)所示。其中图2(a)为步骤s1.3.2.5中得到的实部随机场，经步骤s1.3.3改正之后得到的外接dem结果，图2(b)为步骤s1.3.2.5中得到的虚部随机场，经步骤s1.3.3改正之后得到的外接dem结果。

步骤s2，使用在步骤s1仿真的所述sar影像的外接dem仿真所述sar影像的相位。

按照干涉的理论，主影像与从影像的相位差构成干涉相位。因此，为了完成主从影像的相位仿真，需要首先完成主从影像的干涉相位仿真，仿真过程中需要用到步骤s1.1中的主、从影像的成像几何参数，以及步骤s1.3中的仿真的外接dem。随后完成从影像相位的仿真。最终完成主影像相位的仿真。

各步骤详细如下。

步骤s2.1，使用所述sar影像的外接dem，仿真所述主影像和所述从影像的干涉相位。

计算过程中用到的干涉几何如图3所示。后续的讨论将基于此干涉几何进行。其中，s1为主影像对应的天线相位中心位置，s2为从影像对应的天线相位中心位置。b为基线长度，α为基线倾角，θ为侧视角，ρ为基线矢量与los向的夹角，b⊥为垂直基线，r和r+△r分别为主从影像天线相位中心到地面点p的距离，h为卫星高度，rh为卫星到地心的距离，rh为地面点到地心的距离，re为地球曲率半径，h为地面点高程。本图是雷达干涉的基本原理图。

步骤s2.1具体包括以下子步骤：

步骤s2.1.1，将步骤s1仿真得到的dem转换到sar影像空间。

采用步骤s1.2.1、s1.2.2以及s1.2.3中的r-d模型，将步骤s2得到的外接dem中每个坐标点投影到sar影像像素空间中，并保留像素坐标(x,y)满足1≤x≤s，且1≤y≤l的点，其中s为仿真得到的dem的总列数，l为总行数。投影之后的高程图如图4所示，其中图4(a)为步骤s1.3仿真得到的外接dem，黑色框线为主影像对应的影像范围，图4(b)为转换到sar影像空间的dem。

步骤s2.1.2，使用通过步骤s2.1.1转换后的dem，计算sar影像空间中每个像素点的垂直基线和平行基线。

干涉相位主要包括两类，即平地效应相位，以及地形相位。其中平地效应相位由平行基线表达，地形相位由垂直基线表达。因此本步骤需解算这两个基线参数。本步骤的输入值为步骤s2.1.1转换的dem，输出值用以进行步骤s2.1.3的相位解算。

步骤s2.1.2具体包括以下子步骤：

步骤s2.1.2.1，计算每个像素点的视线向向量。

每个像素点的视线向的向量表达如下：

其中，为像素坐标(x,y)处的地物对应的视线向向量，为主影像的相位中心坐标矢量，可通过步骤s1.2.1中卫星的三维坐标的解算方法解算获得。为每个像素点的坐标矢量，坐标矢量使用r-d模型解算获得，解算过程参见步骤s1.2.1-s1.2.4，其中的高程信息由步骤s2.1.1转换的dem提供。

步骤s2.1.2.2，计算每个像素点垂直视线向的向量

其中，为卫星飞行方向矢量，采用卫星位置矢量与步骤s2.1.2.1中解算得到的像素点视线向向量解算获得，即

卫星位置矢量的解算过程参见步骤s1.2.1。

步骤s2.1.2.3，计算sar影像空间中每个像素点的基线矢量。

步骤s2.1.2.3具体包括以下子步骤：

步骤s2.1.2.3.1，计算每个点的像素矢量在从影像中的像素坐标。

计算过程中，使用的是步骤s1.2.1中的斜距方程，以及步骤s1.2.2中的多普勒方程。即使用两个方程解算两个未知数。解算过程中使用最小二乘或者牛顿迭代法均可，此处不再赘述。

步骤s2.1.2.3.2，计算从影像中像素坐标对应的位置矢量

卫星位置矢量的解算过程参见步骤s1.2.1。

步骤s2.1.2.3.3，计算每个像素点的基线矢量

其中，为主影像中卫星的位置矢量，为从影像中卫星的位置矢量。两卫星矢量均可通过步骤s1.2.1中卫星的三维坐标的解算方法解算获得。

步骤s2.1.2.4，根据每个像素点的基线矢量，计算每个像素点的垂直基线和平行基线。

每个像素点的垂直基线b⊥以及平行基线b||可通过下式进行计算：

其中，为步骤s2.1.2.1计算得到的像素点(x,y)的视线向量，为步骤s2.1.2.2计算得到的像素点垂直视线向的向量。

步骤s2.1.3，使用步骤s2.1.2得到的垂直基线和平行基线计算干涉相位。

步骤s2.1.3具体包括以下子步骤：

步骤s2.1.3.1，计算从影像每个像素点的斜距r2。

其中，ρ为基线矢量与雷达视线向的夹角，表达为cos^-1为反余弦函数，为基线向量的单位向量，即为步骤s2.1.2.3中计算的基线向量，为基线向量的模长，为像素点视线向量的单位向量，即为步骤步骤s2.1.2.1中计算的像素点的视线向量，(x,y)为像素点的像素坐标，为视线向量的模长。r1为主影像中的斜距，表达为r1＝r0+x△r，r0为主影像近地点斜距，x为像素点的x方向坐标，△r为影像的斜距分辨率，近地点斜距和斜距分辨率在表1中均有描述。

步骤s2.1.3.2，计算每个像素点的总相位

其中，λ为雷达波长，本实施例采用的是x波段数据，其波长为3cm，其中r1为主影像中的斜距，表达为r1＝r0+x△r，其中r0为主影像近地点斜距，x为像素点的x方向坐标，△r为影像的斜距分辨率，近地点斜距和斜距分辨率在表1中均有描述。r2为步骤s2.1.3.1中解算得到的从影像每个像素点的斜距。本实施例中，计算得到的总相位如图5(a)所示。

步骤s2.1.3.3，计算每个像素点的地形相位

地形相位的表达式如下：

其中，为步骤s2.1.3.2中获得的每个像素点的总相位，λ为雷达波长，本实施例采用的是x波段数据，其波长为3cm，b||为步骤s2.1.2.4解算得到的每个像素点的平行基线。计算得到的地形相位如图5(b)所示。

步骤s2.2，使用随机噪声模型仿真从影像的相位。

从影像的相位完全由随机噪声组成。依据干涉相位噪声与相干性的关系可知

其中，为相位噪声的方差，γ为相干性，地形测绘过程中，一般令γ>0.9。nl为从影像的多视数，多视数可降低影像噪声，在仿真过程中，为了保持相位的真实性，一般令多视数为1。这种情况下，可解得相位噪声为

依据误差传播定律，单景影像的相位噪声表达为

其中，为公式42得到的相位噪声。

因此，生成均值为0，标准差为0.242的随机噪声，其大小为8000×8000，作为从影像相位。即

其中为从影像相位，为服从正态分布的相位，n(0,0.242)代表均值为0，方差为0.242的标准正态分布，结果如图6(a)所示。

步骤s2.3，使用步骤s2.1得到的所述从影像的干涉相位和步骤s2.2得到的所述从影像的相位，计算所述主影像的相位。

主影像相位的表达式为

其中，为步骤s2.2获得的从影像相位，为主影像像素空间中每个像素点的干涉相位，由步骤s2.1.3.2获得，为噪声相位，由其生成过程与步骤s2.2相同，不再赘述。仿真得到的主影像相位如图6(b)所示。

步骤s3，使用在步骤s1仿真的所述sar影像的外接dem仿真所述sar影像的振幅。

sar影像为侧视成像，影像中存在叠掩和阴影。因此影像由三类主要特征构成，第一类即平地和非平地像素构成的一般像素，第二类为叠掩区的像素，第三类为阴影区的像素。计算过程中，会使用到所述主影像以及从影像像素空间中每个像素点的本地法线向量、像平面法线、名义入射角、本地入射角等。具体步骤如下。

步骤s3.1，仿真所述主影像的振幅。

步骤s3.1.1，计算所述主影像像素空间每个像素点的本地法线向量。

本地法线向量使用如下公式计算：

其中，为像素点(1,y)对应的视线向向量，为像素点(0,y)对应的视线向向量，为像素点(1,y-1)对应的视线向向量，代表向量的叉乘。公式中的三个视线向向量均可由步骤s2.1.2.1获得。本步骤为步骤s3.1.2提供本地法线向量的单位向量，进而判断某一像素是否隶属于平地、非平地、叠掩、阴影区域。

步骤s3.1.2，基于所述主影像像素空间中像素点的本地法线向量，判断所述像素点是否是平地像素。

本地法线向量与本地坐标向量的夹角可表示为χ，其表达式为

其中，cos^-1为反余弦函数，为本地法线向量的单位向量，其表达式为

为步骤s3.1.1计算得到的每个像素点的本地法线向量，为本地法线向量的模长，为像素点坐标矢量的单位向量，其表达式与上式相同，需采用每个限速点的坐标矢量解算得来，坐标矢量使用r-d模型解算获得，解算过程参见步骤s1.2.1-s1.2.4，其中的高程信息由步骤s2.1.1转换的dem提供。当χ小于一定角度时，可认为是平地像素(如图7所示)。优选的，令χ<1°时为平地像素。

步骤s3.1.3，计算每个像素点的名义入射角β。

式中，各参数的含义与步骤s2.1相同，cos^-1为反余弦函数，rh为卫星到地心的距离，r1为主影像斜距，表达为r1＝r0+x△r，其中r0为主影像近地点斜距，x为像素点的x方向坐标，△r为影像的斜距分辨率，近地点斜距和斜距分辨率在表1中均有描述。re为地球曲率半径，h为地面点高程。

步骤s3.1.4，基于对所述像素点是否为平地像素的判断和像素点的名义入射角，计算每个像素点的本地入射角θ。

其中，β为步骤s3.1.3中得到的名义入射角，cos^-1为反余弦函数，为卫星速度的单位矢量，由步骤s1.2.2.1获得，为本地法线向量的单位矢量，由步骤s3.1.3获得。

步骤s3.1.5，基于对所述像素点是否为平地像素的判断和像素点的本地入射角，计算每个像素点的振幅值am。

公式表达式如下：

其中，θ为步骤s3.1.4中得到的每个像素点的本地入射角，μ为像素尺寸归一化参数，表达式如下：

其中，χ为本地法线向量与本地坐标向量的夹角，θ为步骤s3.1.4中得到的每个像素点的本地入射角

至此，本发明已经完成了主影像像素空间中大多数像素的振幅计算。然而，sar影像为侧视成像，在成像过程中会生成特有的叠掩、阴影等现象。因此需要在下面的步骤中计算叠掩区和阴影区的振幅值。

步骤s3.1.6，计算叠掩区振幅值。

步骤s3.1.6.1，确定叠掩点。

计算每个点的本地坡向角ψ，其表达式为

其中，cos^-1为反余弦函数，为本地法线向量的单位向量，通过公式(48)计算得来，为像平面法线的单位向量，将ψ大于90°的定义为叠掩点。像平面法线的单位向量表达为

其中，|nimg|为像平面法线的模长，nimg为像平面法线，其表达式为

为主影像成像时刻的卫星速度矢量，其解算方法参见步骤s1.2.2.1。为像素点(1,y-1)对应的视线向向量。

步骤s3.1.6.2，计算叠掩点振幅值。

寻找叠掩点左侧的最近有效像素及其振幅值(xl,yl,aml)，右侧的最近有效像素及其振幅值(xr,yr,amr)，使用线性插值获取此点(x,y)的振幅am，即

步骤s3.1.7，计算阴影区灰度值。

步骤s3.1.7.1，确定阴影点。

本地入射角大于90°的像素点定义为阴影点。本地入射角由步骤s3.1.4计算获得。

步骤s3.1.7.2，计算阴影点振幅值。

寻找阴影点左侧的最近有效像素及其振幅值(xl,yl,aml)，右侧的最近有效像素及其振幅值(xr,yr,amr)，使用线性插值获取此点(x,y)的振幅am，即

步骤s3.2，仿真所述从影像的振幅。

从影像振幅仿真过程中，用到了表2中提供的参数。具体步骤与s3.1类似，此处不再赘述。图8给出了使用主影像以及从影像的仿真结果。

步骤s4，基于步骤s2仿真得到的sar影像的相位和步骤s3仿真得到的sar影像的振幅，仿真所述sar影像。

步骤s4.1，仿真所述主影像。

对主影像每个像素来说，其像素值表达为一个复数，即

f(x,y,m)＝am+bmi(58)

其中，f(x,y,m)为主影像每个像素的像素值，(x,y)为主影像像素的坐标，m为主影像像素的复数值，am为步骤s3.1仿真的主影像振幅，为步骤2.3仿真的主影像的相位。

步骤s4.2，仿真所述从影像。

对从影像每个像素来说，其像素值表达为一个复数，即

f(x,y,s)＝as+bsi(59)

其中，f(x,y,s)为从影像每个像素的像素值，(x,y)为为从影像像素的坐标，s为从影像像素的复数值，as为步骤s3.2仿真的从影像振幅，为步骤2.2仿真的从影像相位。

图9给出了相位和振幅叠加的结果，其中图9(a)为主影像结果，图9(b)为从影像结果。

现有技术多采用srtm数据或已知的dem数据进行仿真，仿真过程中无法从纯数学角度探讨地形测绘过程中的误差来源，本发明提出了基于fbm的仿真方法，一方面，fbm能够表达数据的自相似性，这与地形本身一定范围内的相似性一致；另一方面，fbm能够表达数据的随机性，这与地形本身极小范围内的随机起伏一致。因此，fbm的数学特征能够很好的描述地形信息，这使得从理论层面进行地形特性的研究成为可能。此外，传统仿真过程中，参数为仿真系统生成，而非真实参数，本发明仿真过程中采用的参数为真实的tandem-x卫星干涉参数，参数具有极好的工程应用性，能够确保在干涉处理过程中，所有的误差来源不会与干涉几何状态有太大关联，从而使得研究人员能够基于本发明得到更为精确的误差仿真结论。本发明为基于干涉sar技术的卫星测绘提供了良好的参考。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。