CN113514035A - 一种全球数字高程模型约束的影像区域网平差方法 - Google Patents
一种全球数字高程模型约束的影像区域网平差方法 Download PDFInfo
- Publication number
- CN113514035A CN113514035A CN202110785058.1A CN202110785058A CN113514035A CN 113514035 A CN113514035 A CN 113514035A CN 202110785058 A CN202110785058 A CN 202110785058A CN 113514035 A CN113514035 A CN 113514035A
- Authority
- CN
- China
- Prior art keywords
- image
- elevation
- ground
- rfm
- adjustment
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 45
- 238000004364 calculation method Methods 0.000 claims abstract description 40
- 238000012937 correction Methods 0.000 claims abstract description 29
- PXFBZOLANLWPMH-UHFFFAOYSA-N 16-Epiaffinine Natural products C1C(C2=CC=CC=C2N2)=C2C(=O)CC2C(=CC)CN(C)C1C2CO PXFBZOLANLWPMH-UHFFFAOYSA-N 0.000 claims description 47
- 230000009466 transformation Effects 0.000 claims description 46
- 239000011159 matrix material Substances 0.000 claims description 22
- 238000013507 mapping Methods 0.000 claims description 21
- 238000005259 measurement Methods 0.000 claims description 13
- 238000012892 rational function Methods 0.000 claims description 6
- 238000011478 gradient descent method Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 2
- 230000008901 benefit Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 241000132092 Aster Species 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000000691 measurement method Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 210000000988 bone and bone Anatomy 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000002542 deteriorative effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C11/00—Photogrammetry or videogrammetry, e.g. stereogrammetry; Photographic surveying
- G01C11/04—Interpretation of pictures
Landscapes
- Engineering & Computer Science (AREA)
- Multimedia (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种全球数字高程模型约束的影像区域网平差方法，首先构建立体影像平差区域网，布设密集连接点。随后获取测区内全球数字高程模型数据，将测区划分为平地、丘陵和山地三种地形区域，设置利用全球数字高程模型数据针对平差计算获得的不同地形区域内连接点对应地面点高程值改正数的约束规则；最后开展区域网平差计算，在每次迭代计算时均利用该约束规则对连接点对应地面点的高程值改正数进行约束，直至平差结束，并根据平差结果为每张影像生成新的成像几何模型参数文件。本发明通过在不同地形区域采用不同规则使用全球数字高程模型对区域网平差进行高程约束，实现无地面控制条件下卫星遥感影像高程精度的大幅度提升。
Description
技术领域
本发明涉及摄影测量技术领域，尤其涉及一种全球数字高程模型约束的影像区域网平差方法。
背景技术
数字高程模型是用一组有序数值阵列形式表示地面高程的一种实体地面模型。目前，通过互联网可以公开免费获取的全球范围数字高程模型产品主要包括：①SRTM-DEM(全称Space Shuttle Radar Topography Mission)：由美国太空总署等制作并免费提供，其数据覆盖范围为北纬60°到南纬56°之间的所有陆地区域，占地球陆地表面的80％，全球最高分辨率为30米，官方标称高程中误差10米。②ASTER GDEM(全称Advanced SpaceborneThermal Emission and Reflection Radiometer Global Digital Elevation Model)：由美国NASA和日本联合制作并免费提供，其数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，占地球陆地表面的99％，全球最高空间分辨率为30米，官方标称高程中误差12米。③ALOS World 3D-30m(全称ALOS Global Digital Surface Model"ALOS World 3D-30m(AW3D30)")：由日本航空航天勘探局制作并免费提供，数据覆盖全球所有陆地区域，全球最高分辨率为5米，目前可免费获取的最高空间分辨率为30米，官方标称高程中误差3米。根据大量学者研究表明，经过多年来不断编辑和更新，前述各种免费全球数字高程模型产品的高程精度一般要高于其官方标称精度，且高程精度与地形类型存在较强的相关性，在高程起伏较小的平坦地区其高程精度极高，如ALOS World 3D-30m和SRTM-DEM在平地地形区域的高程中误差甚至可以优于2～5米，完全可以满足我国1:5万比例尺立体测图的高程精度要求。
在无地面控制条件下，提升国产卫星影像几何精度(尤其是高程精度)，使得国产卫星影像能够满足1:50000甚至更大比例尺测图精度要求，一直是当前国产卫星追求的理想目标。大量机构和学者的研究表明，在无地面控制条件下主要的国产遥感卫星(如资源三号系列卫星、天绘一号系列卫星)影像的平面中误差一般为8～25米，可以满足我国1:50000比例尺测图的平面精度要求(注：我国1:5万测图平面精度要求在平地和丘陵为25米、山地和高山地为37.5米)；影像高程中误差一般为6～15米，达不到1:50000比例尺测图的高程精度要求(注：1:5万测图高程精度要求在平地、丘陵、山地、高山地分别为3米、5米、8米、14米)。因此高程精度是制约国产卫星影像满足1:5万比例尺测图精度的关键因素，如何充分利用丰富且免费的全球数字高程模型数据资源，在不再借助其他地面控制资料条件下，通过区域网平差方法实现卫星影像高程精度的较大幅度提升并满足1:50000甚至更大比例尺的测图要求，这对基于国产卫星的全球1:50000比例尺测图以及国内控制获取困难地区的高精度测图都拥有非常重要的实际意义和研究价值。
在目前已有的利用全球数字高程模型产品辅助卫星影像区域网平差的研究和技术方面，主要集中在利用全球数字高程模型约束的弱交会条件下的平面影像区域网平差方面，仅为了保障弱交会影像区域网平差的正确性，并不能提升卫星立体影像高程精度。还有部分研究学者先通过卫星立体影像对制作整景DEM，再通过与全球数字高程模型产品进行配准，修正DEM的定向参数，最后再反向改正立体影像的成像几何模型参数。该思路虽然也能达到提升卫星立体影像的高程精度的目的，但是没有顾及全球数字高程模型产品在不同地形区域的高程精度差异，更没有充分利用全球数字高程模型产品在平坦地形区域高程精度极高的特性，因此对立体影像提升效果还存在较大提升空间，且工序比较繁杂，算法也比较复杂。
综上所述，目前急需寻找一种能有效区分全球数字高程模型产品在不同地形区域高程精度差异，并充分利用其在平地地形区域高程精度极高的特性，简单有效的提升卫星立体影像高程精度的技术方法。
发明内容
本发明的目的在于提供一种全球数字高程模型约束的影像区域网平差方法，从而解决现有技术中存在的前述问题。
为了实现上述目的，本发明采用的技术方案如下：
一种全球数字高程模型约束的影像区域网平差方法，包括如下步骤，
S1、将测区内连续覆盖的遥感卫星立体影像构建一个平差区域网，并采用像方仿射变换补偿的基于有理函数模型的区域网平差模型；有理函数模型即RFM；
S2、在平差区域网内立体影像上布设密集的连接点；
S3、选取并从互联网上下载测区内一种或多种免费全球数字高程模型产品，并根据地表高程起伏情况将测区概略地划分为平地、丘陵和山地三种地形区域，结合全球数字高程模型产品在不同地形区域的高程精度经验值，分别为平差计算时求解的位于三种地形区域的连接点对应地面点高程值改正数设置不同的高程约束规则；
S4、确定连接点对应地面点的平面坐标初值、高程值初值以及影像RFM仿射变换模型参数初值；
S5、针对每个连接点，将影像RFM的仿射变换参数和连接点对应地面点的三维坐标作为未知数，逐点构建连接点的误差方程；
S6、对误差方程进行法化，获取法方程；采用最小二乘平差原理对法方程进行整体平差计算，计算每张影像RFM的仿射变换参数以及连接点对应地面点的三维坐标的改正数；
S7、每完成一次平差计算，判断是否满足平差收敛条件，若是，则平差计算结束；若否，则利用平差计算获取的改正数更新影像RFM的仿射变换参数以及连接点对应地面点的三维坐标，同时利用高程约束规则对更新后的不同地形区域连接点对应地面点的高程值进行修改，重复步骤S5到S6，再次进行平差迭代计算，直到满足收敛条件为止；
S8、利用最后一次平差迭代计算得到的影像RFM的仿射变换参数，对影像RFM进行像方仿射变换补偿，为每张影像生成新的RFM。
优选的，步骤S1中的RFM是利用有理多项式建立影像的像方与其对应的物方之间的数学映射关系，RFM方程定义为，
其中，(P,L,H)为归一化后的地面坐标；(x,y)为归一化后的影像坐标；NumS(P,L,H)、DenS(P,L,H)、NumL(P,L,H)和DenL(P,L,H)分别表示一般多项式，多项式中每一个变量P,L,H的幂均不超过3次，所有变量的幂之和也不超过3次；
影像RFM的误差补偿采用像方仿射变换模型，则RFM描述的像点坐标(x,y)和地面点坐标(P,L,H)之间的关系被修正为，
其中，(△x,△y)为影像坐标(x,y)的系统误差像方补偿值，其值为，
△y＝a0+a1·y+a2·x
△x＝b0+b1·y+b2·x
其中，(a0,a1,a2,b0,b1,b2)表示仿射变换模型的参数。
优选的，步骤S2中先采用影像自动匹配算法匹配连接点，当影像自动匹配算法匹配的连接点分布不均匀、数量不足或维度不够时，采用人工判读方法补测连接点。
优选的，步骤S3中将测区划分为平地、丘陵和山地三种地形区域时，利用全球数字高程模型计算各立体影像覆盖范围的地面坡度角和高程差，当立体影像覆盖范围内地面坡度角小于2度或高程差小于80米，则该立体影像覆盖区域为平地地形区域；当立体影像覆盖范围内地面坡度角小于6度或高程差小于300米，则该立体影像覆盖区域为丘陵地形区域；否则则为山地地形区域。
优选的，步骤S3中，根据所采用的全球数字高程模型产品在不同地形区域的高程精度经验值，分别设置每次平差迭代时计算平地、丘陵和山地地形区域连接点对应地面点高程改正值的高程约束规则，具体如下，
hplan_corr∈[(hinit-σplan)，(hinit+σplan)]
hhill_corr∈[(hinit-σhill)，(hinit+σhill)]
hmount_corr∈[(hinit-σmount)，(hinit+σmount)]
其中，hplan_corr、hhill_corr和hmount_corr分别为平差计算获取的平地、丘陵和山地地形区域的连接点对应地面点的高程改正值，hinit为根据连接点对应地面点平面经纬度坐标从全球数字高程模型产品上直接读取的高程值，σplan、σhill和σmount分别为全球数字高程模型产品在平地、丘陵和山地地形区域的高程中误差经验值。
优选的，步骤S4中，连接点对应地面点平面坐标的初值是通过影像RFM对连接点进行直接前方交会获取；影像RFM仿射变换模型的初值由人为预设；连接点对应地面点高程初值是根据连接点对应地面点平面坐标初值从全球高程模型产品上读取获得。
优选的，步骤S5中，逐点构建的连接点的误差方程如下，
V＝AX+BY-L,P
式中，V＝[vx vy]T表示连接点行和列坐标观测值的残差向量；X＝[△a0 △a1 △a2△b0 △b1 △b2]T为影像RFM的像方坐标系统误差仿射变换补偿模型参数的改正数向量，Y＝[△lat △lon △h]T为连接点对应地面点三维坐标的改正数向量；
为未知数X的系数矩阵；
为未知数Y的系数矩阵；
为初值计算得到的常数项；P为权矩阵。
优选的，步骤S6具体包括如下内容，
S61、根据最小二乘平差原理，对误差方程进行法化，得到的法方程如下，
记为，
S62、对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后利用前方交会的方式更新连接点对应地面点的三维坐标，将上式变换为，
Nt＝G
其中，
N中
为2阶对角矩阵，在计算时通过对每个小矩阵求逆的方式得到；
S63、求解出
之后便能够分别得到N、G，针对公式Nt＝G，利用共轭梯度下降法进行迭代求解，当两次求解得到的t的差值小于设定的阈值或者求解次数超过预设次数时结束迭代，输出得到最终的t，即仿射变换参数的改正数。
优选的，步骤S7中，利用高程约束规则对每次平差迭代计算后获取的改正后的连接点对应地面点高程值进行修改，具体为，
针对平地地形区域的连接点，采用下式进行高程约束，
针对丘陵地形区域的连接点，采用下式进行高程约束，
针对山地地形区域的连接点，采用下式进行高程约束，
优选的，步骤S8具体包括如下内容，
S81、采用影像原始RFM计算影像四个角点对应的地面范围，利用全球数字高程模型计算该影像覆盖区域的最大和最小高程值；将影像覆盖区域在高程方向以一定的间隔分层，在平面上以一定的格网大小建立地面规则格网，构建一个虚拟的地面点立体空间格网，作为控制点格网；利用影像原始RFM以及平差获取的仿射变换参数，计算所有控制点在影像上的影像坐标；加密控制点格网和层，建立独立的检查点；
S82、采用最小二乘平差原理求解RFM参数；首先将RFM方程变形为
Fy＝NumL(P,L,H)-y*DenL(P,L,H)＝0
Fx＝NumS(P,L,H)-x*DenS(P,L,H)＝0
则获取线性化误差方程为，
V＝Bx-l,W
其中，
x＝[a1i a2ja3i a4j]T，W为权矩阵；
根据最小二乘平差原理，求解获取影像新RFM参数为，
x＝(BTB)-1BTl
S83、利用影像新RFM计算检查点对应的影像坐标，比较与利用影像原始RFM以及仿射变换参数计算的检查点对应影像坐标之间的差值，评定影像新RFM精度，在达到精度要求时输出影像新的RFM参数，为每张影像生成新RFM。
本发明的有益效果是：1、能够将国产卫星立体影像的高程精度有效提升至满足我国1:50000甚至1:25000比例尺测图精度要求，使得基于国产卫星影像的全球范围无地面控制条件下高精度测图成为了可能。2、使用了可免费获取全球数字高程模型数据，相对传统的使用地面控制点的区域网平差技术方案，减少了外业控制测量或高精度控制点收集环节，大大降低了作业成本，提升了作业效率，尤其是保障了很多无法开展外业控制测量和收集高精度控制资料的区域(如境外区域、环境恶劣或交通闭塞的测图困难地区等)高精度测图。3、相对其他利用免费全球数字高程模型辅助区域网平差的技术方法，本发明既充分利用了全球数字高程模型在地形平坦区域高程精度极高的优势，又有效避免了全球数字高程模型在局部区域高程错误或精度低下的影响，因此对立体影像高程精度的提升效果更加有效和优越。同时作业步骤更加简单，因此作业效率更高，可操作性和可靠性也更强。
附图说明
图1是本发明实施例中区域网平差方法的原理流程图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。
如图1所示，本实施例中，提供了一种全球数字高程模型约束的影像区域网平差方法，包括如下步骤，
S1、将测区内连续覆盖的遥感卫星立体影像构建一个平差区域网，并采用像方仿射变换补偿的基于有理函数模型(RFM)的区域网平差模型；
S2、在平差区域网内立体影像上布设密集的连接点；
S3、选取并从互联网上下载测区内一种或多种免费全球数字高程模型产品，并根据地表高程起伏情况将测区概略地划分为平地、丘陵和山地三种地形区域，结合全球数字高程模型产品在不同地形区域的高程精度经验值，分别为平差计算时求解的位于三种地形区域的连接点对应地面点高程值改正数设置不同的高程约束规则；
S4、确定连接点对应地面点的平面坐标初值、高程值初值以及影像RFM仿射变换模型参数初值；
S5、针对每个连接点，将影像RFM的仿射变换参数和连接点对应地面点的三维坐标作为未知数，逐点构建连接点的误差方程；
S6、对误差方程进行法化，获取法方程；采用最小二乘平差原理对法方程进行整体平差计算，计算每张影像RFM的仿射变换参数以及连接点对应地面点的三维坐标的改正数；
S7、每完成一次平差计算，判断是否满足平差收敛条件，若是，则平差计算结束；若否，则利用平差计算获取的改正数更新影像RFM的仿射变换参数以及连接点对应地面点的三维坐标，同时利用高程约束规则对更新后的不同地形区域连接点对应地面点的高程值进行修改，重复步骤S5到S6，再次进行平差迭代计算，直到满足收敛条件为止；
S8、利用最后一次平差迭代计算得到的影像RFM的仿射变换参数，对影像RFM进行像方仿射变换补偿，为每张影像生成新的RFM。
本实施例中，参见附图1，区域网平差方法主要包括八部分内容，分别为建立区域网和平差模型、布设连接点、利用全球数字高程模型设置高程约束规则、确定平差中各未知变量的初值、构建误差方程、平差计算、平差迭代和高程约束、生成影像新RFM。下面分别对这八部分内容进行详细的解释说明。
一、建立区域网和平差模型
该部分对应步骤S1。步骤S1中RFM为有理函数模型，其是利用有理多项式建立影像的像方(影像像素坐标)与其对应的物方(地面大地坐标)之间的数学映射关系，RFM方程定义为，
其中，为了避免计算过程中参数数值量级差别过大而引起入舍误差，需要将影像像点坐标和对应的地面点坐标均归一化到-1～1之间，以增强参数求解的稳定性；(P,L,H)为归一化后的地面坐标；(x,y)为归一化后的影像坐标；NumS(P,L,H)、DenS(P,L,H)、NumL(P,L,H)和DenL(P,L,H)分别表示一般多项式，多项式中每一个变量P,L,H的幂均不超过3次，所有变量的幂之和也不超过3次；其形式如下，
NL(P,L,H)＝a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2+a18L2H+a19P2H+a20H3
DL(P,L,H)＝b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2+b18L2H+b19P2H+b20H3
Ns(P,L,H)＝c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2+c18L2H+c19P2H+c20H3
Ds(P,L,H)＝d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2+d18L2H+d19P2H+d20H3
其中，aj,bj,cj,dj(j＝0,1,…,19)为有理多项式系数。
影像RFM的误差补偿采用像方仿射变换模型，则RFM描述的像点坐标(x,y)和地面点坐标(P,L,H)之间的关系被修正为，
其中，(△x,△y)为影像坐标(x,y)的系统误差像方补偿值，其值为，
△y＝a0+a1·y+a2·x
△x＝b0+b1·y+b2·x
其中，(a0,a1,a2,b0,b1,b2)表示仿射变换模型的参数。
二、布设连接点
该部分内容对应步骤S2。步骤S2中可采用影像自动匹配算法或人工判读/量测方法布设连接点；为了避免全球数字高程模型产品局部区域高程错误或精度低下等偶然因素对后续区域网平差整体高程精度的影响，在区域网中连接点数量应尽量充足且点位分布应尽量均匀。
先采用影像自动匹配算法匹配连接点，常见的影像自动匹配算法有基于像方灰度的影像匹配算法(如最小二乘法)、基于物方的影像匹配算法、基于像方特征的跨接法影像匹配、金字塔多级影像匹配、SIFT等。要求区域网中每个立体影像范围内的连接点应分布均匀，根据实践经验，本实施例中每一标准景立体影像上连接点数量大于100个；相邻立体影像之间的连接点数量大于10个；区域内90％的连接点维度(即该连接点连接影像的数量)应等于重叠区影像的总数量。
当影像自动匹配算法匹配的连接点分布不均匀、数量不足或维度不够时，采用人工判读方法/量测方法补测连接点，人工判读/量测的连接点应位于影像清晰、特征明显、反差较大、易于转刺和量测的固定目标上。
三、利用全球数字高程模型设置高程约束规则
该部分内容对应步骤S3。步骤S3中，选取并从互联网上下载测区内一种或多种免费全球数字高程模型产品，免费全球数字高程模型产品主要包括：
①SRTM-DEM：由美国美国太空总署(NASA)等制作并免费提供，全称Space ShuttleRadar Topography Mission数字表面模型。其数据覆盖范围为北纬60°到南纬56°之间的所有陆地区域，占地球陆地表面的80％。全球最高分辨率为30米。
②ASTER GDEM：由美国NASA和日本联合制作并免费提供，全称AdvancedSpaceborne Thermal Emission and Reflection Radiometer Global DigitalElevation Model，即先进星载热发射和反射辐射仪全球数字高程模型。其数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，占地球陆地表面的99％。全球最高空间分辨率为30米。
③ALOS World 3D-30m：由日本航空航天勘探局制作并免费提供，全称ALOSGlobal Digital Surface Model"ALOS World 3D-30m(AW3D30)"。数据覆盖全球所有陆地区域。全球最高分辨率为5米，目前可免费获取的最高空间分辨率为30米。
将测区概略划分为平地、丘陵和山地三种地形区域时，利用全球数字高程模型计算各立体影像覆盖范围的地面坡度角和高程差，当立体影像覆盖范围内地面坡度角小于2度或高程差小于80米，则该立体影像覆盖区域为平地地形区域；当立体影像覆盖范围内地面坡度角小于6度或高程差小于300米，则该立体影像覆盖区域为丘陵地形区域；否则则为山地地形区域。
根据所采用的全球数字高程模型产品在不同地形区域的高程精度经验值，分别设置每次平差迭代时计算平地、丘陵和山地地形区域连接点对应地面点高程改正值的高程约束规则，具体如下，
hplan_corr∈[(hinit-σplan)，(hinit+σplan)]
hhill_corr∈[(hinit-σhill)，(hinit+σhill)]
hmount_corr∈[(hinit-σmount)，(hinit+σmount)]
其中，hplan_corr、hhill_corr和hmount_corr分别为平差计算获取的平地、丘陵和山地地形区域的连接点对应地面点的高程改正值，hinit为根据连接点对应地面点平面经纬度坐标从全球数字高程模型产品上直接读取的高程值，σplan、σhill和σmount分别为全球数字高程模型产品在平地、丘陵和山地地形区域的高程中误差经验值。
四、确定平差中各未知变量的初值
该部分对应步骤S4。连接点对应地面点平面坐标的初值是通过影像RFM对连接点进行直接前方交会获取；影像RFM仿射变换模型的初值由人为预设；连接点对应地面点高程初值是根据连接点对应地面点平面坐标初值从全球高程模型产品上读取获得。
五、构建误差方程
该部分对应步骤S5。步骤S5中，逐点构建的连接点的误差方程如下，
V＝AX+BY-L,P
式中，V＝[vx vy]T表示连接点行和列坐标观测值的残差向量；X＝[△a0 △a1 △a2△b0 △b1 △b2]T为影像RFM的像方坐标系统误差仿射变换补偿模型参数的改正数向量，Y＝[△lat △lon △h]T为连接点对应地面点三维坐标的改正数向量；
为未知数X的系数矩阵；
为未知数Y的系数矩阵；
为初值计算得到的常数项；P为权矩阵，通常各类观测值的权值由他们的先验信息确定，一般可取观测值10倍先验标准差，并在每次平差迭代计算后重新计算各类观测值的权。
六、平差计算
该部分对应步骤S6。步骤S6具体包括如下内容，
S61、根据最小二乘平差原理，对误差方程进行法化，得到的法方程如下，
记为，
S62、对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后利用前方交会的方式更新连接点对应地面点的三维坐标，将上式变换为，
Nt＝G
其中，
N中
为2阶对角矩阵，在计算时通过对每个小矩阵求逆的方式得到；这一步可以多线程并行处理，能够大大节省计算时间，提升计算效率。
S63、求解出
之后便能够分别得到N、G，针对公式Nt＝G，利用数学中的共轭梯度下降法进行迭代求解，当两次求解得到的t的差值小于设定的阈值(本实施例中设定的阈值为0.1个像元pixel，具体可以根据实际情况进行确定)或者求解次数超过预设次数(本实施例中预设次数为20，具体可以根据实际情况记性确定)时结束迭代，输出得到最终的t，即仿射变换参数的改正数。
七、平差迭代和高程约束
该部分对应步骤S7。步骤S7中每完成依次平差计算，需要判断其是否满足平差迭代收敛条件，收敛条件即为所有影像RFM的仿射变换参数中平移参数小于预设阈值。预设阈值可以根据实际情况进行选择，以便更好的满足实际需求。
如果满足迭代收敛条件，则平差结束；反之，则利用平差计算获取的改正数更新影像RFM的仿射变换参数以及连接点所对应地面点的三维坐标，同时利用高程约束规则对更新后的不同地形区域连接点对应地面点的高程值进行修改，具体为：
针对平地地形区域的连接点，采用下式进行高程约束，
针对丘陵地形区域的连接点，采用下式进行高程约束，
针对山地地形区域的连接点，采用下式进行高程约束，
约束修改后，重复步骤S5-步骤S6，再次进行平差迭代计算，直到满足迭代收敛条件为止；最后输出每张影像RFM的仿射变换参数。
八、生成影像新RFM
该部分对应步骤S8。步骤S8具体包括如下内容，
S81、采用影像原始RFM计算影像四个角点对应的地面范围；利用全球数字高程模型数据计算该影像覆盖区域的最大最小高程；然后将影像覆盖区域在高程方向以一定的间隔分层，为了防止设计矩阵状态恶化，一般高程方向分层的层数超过2，本实施例中高程分为10层，层数可以根据实际情况进行设置；在平面上以一定的格网大小建立地面规则格网，构建一个虚拟的地面点立体空间格网，作为控制点格网。本实施例中平面分为15×15格网，就是将该影像对应影像范围分成15×15的格子，共有16×16个格网点，格网和格网点都可以根据实际情况进行设置；最后利用影像原始RFM以及平差获取的仿射变换参数，计算所有控制点在影像上的影像坐标；加密控制点格网和层，建立独立的检查点。
S82、采用最小二乘平差原理求解RFM参数；首先将RFM方程变形为
Fy＝NumL(P,L,H)-y*DenL(P,L,H)＝0
Fx＝NumS(P,L,H)-x*DenS(P,L,H)＝0
则获取线性化误差方程为，
V＝Bx-l,W
其中，
x＝[a1i a2ja3i a4j]T，W为权矩阵；
根据最小二乘平差原理，求解获取影像新RFM参数为，
x＝(BTB)-1BTl
S83、利用影像新RFM计算检查点对应的影像坐标，比较与利用影像原始RFM以及仿射变换参数计算的检查点对应影像坐标之间的差值，评定影像新RFM精度，在达到精度要求时输出影像新的RFM参数，为每张影像生成新RFM。
通过采用本发明公开的上述技术方案，得到了如下有益的效果：
本发明提供的一种全球数字高程模型约束的影像区域网平差方法，能够将国产卫星立体影像的高程精度有效提升至满足我国1:50000甚至1:25000比例尺测图精度要求，使得基于国产卫星影像的全球范围无地面控制条件下高精度测图成为了可能。使用了可免费获取全球数字高程模型数据，相对传统的使用地面控制点的区域网平差技术方案，减少了外业控制测量或高精度控制点收集环节，大大降低了作业成本，提升了作业效率，尤其是保障了很多无法开展外业控制测量和收集高精度控制资料的区域(如境外区域、环境恶劣或交通闭塞的测图困难地区等)高精度测图。相对其他利用免费全球数字高程模型辅助区域网平差的技术方法，本发明既充分利用了全球数字高程模型在地形平坦区域高程精度极高的优势，又有效避免了全球数字高程模型在局部区域高程错误或精度低下的影响，因此对立体影像高程精度的提升效果更加有效和优越。同时作业步骤更加简单，因此作业效率更高，可操作性和可靠性也更强。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。
Claims (10)
1.一种全球数字高程模型约束的影像区域网平差方法，其特征在于：包括如下步骤，
S1、将测区内连续覆盖的遥感卫星立体影像构建一个平差区域网，并采用像方仿射变换补偿的基于有理函数模型的区域网平差模型；有理函数模型即RFM；
S2、在平差区域网内立体影像上布设密集的连接点；
S3、选取并从互联网上下载测区内一种或多种免费全球数字高程模型产品，并根据地表高程起伏情况将测区概略地划分为平地、丘陵和山地三种地形区域，结合全球数字高程模型产品在不同地形区域的高程精度经验值，分别为平差计算时求解的位于三种地形区域的连接点对应地面点高程值改正数设置不同的高程约束规则；
S4、确定连接点对应地面点的平面坐标初值、高程值初值以及影像RFM仿射变换模型参数初值；
S5、针对每个连接点，将影像RFM的仿射变换参数和连接点对应地面点的三维坐标作为未知数，逐点构建连接点的误差方程；
S6、对误差方程进行法化，获取法方程；采用最小二乘平差原理对法方程进行整体平差计算，计算每张影像RFM的仿射变换参数以及连接点对应地面点的三维坐标的改正数；
S7、每完成一次平差计算，判断是否满足平差收敛条件，若是，则平差计算结束；若否，则利用平差计算获取的改正数更新影像RFM的仿射变换参数以及连接点对应地面点的三维坐标，同时利用高程约束规则对更新后的不同地形区域连接点对应地面点的高程值进行修改，重复步骤S5到S6，再次进行平差迭代计算，直到满足收敛条件为止；
S8、利用最后一次平差迭代计算得到的影像RFM的仿射变换参数，对影像RFM进行像方仿射变换补偿，为每张影像生成新的RFM。
2.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S1中的RFM是利用有理多项式建立影像的像方与其对应的物方之间的数学映射关系，RFM方程定义为，
其中，(P,L,H)为归一化后的地面坐标；(x,y)为归一化后的影像坐标；NumS(P,L,H)、DenS(P,L,H)、NumL(P,L,H)和DenL(P,L,H)分别表示一般多项式，多项式中每一个变量P,L,H的幂均不超过3次，所有变量的幂之和也不超过3次；
影像RFM的误差补偿采用像方仿射变换模型，则RFM描述的像点坐标(x,y)和地面点坐标(P,L,H)之间的关系被修正为，
其中，(△x,△y)为影像坐标(x,y)的系统误差像方补偿值，其值为，
△y＝a0+a1·y+a2·x
△x＝b0+b1·y+b2·x
其中，(a0,a1,a2,b0,b1,b2)表示仿射变换模型的参数。
3.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S2中先采用影像自动匹配算法匹配连接点，当影像自动匹配算法匹配的连接点分布不均匀、数量不足或维度不够时，采用人工判读方法补测连接点。
4.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S3中将测区划分为平地、丘陵和山地三种地形区域时，利用全球数字高程模型计算各立体影像覆盖范围的地面坡度角和高程差，当立体影像覆盖范围内地面坡度角小于2度或高程差小于80米，则该立体影像覆盖区域为平地地形区域；当立体影像覆盖范围内地面坡度角小于6度或高程差小于300米，则该立体影像覆盖区域为丘陵地形区域；否则则为山地地形区域。
5.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S3中，根据所采用的全球数字高程模型产品在不同地形区域的高程精度经验值，分别设置每次平差迭代时计算平地、丘陵和山地地形区域连接点对应地面点高程改正值的高程约束规则，具体如下，
hplan_corr∈[(hinit-σplan)，(hinit+σplan)]
hhill_corr∈[(hinit-σhill)，(hinit+σhill)]
hmount_corr∈[(hinit-σmount)，(hinit+σmount)]
其中，hplan_corr、hhill_corr和hmount_corr分别为平差计算获取的平地、丘陵和山地地形区域的连接点对应地面点的高程改正值，hinit为根据连接点对应地面点平面经纬度坐标从全球数字高程模型产品上直接读取的高程值，σplan、σhill和σmount分别为全球数字高程模型产品在平地、丘陵和山地地形区域的高程中误差经验值。
6.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S4中，连接点对应地面点平面坐标的初值是通过影像RFM对连接点进行直接前方交会获取；影像RFM仿射变换模型的初值由人为预设；连接点对应地面点高程初值是根据连接点对应地面点平面坐标初值从全球高程模型产品上读取获得。
7.根据权利要求1所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S5中，逐点构建的连接点的误差方程如下，
V＝AX+BY-L,P
式中，V＝[vx vy]T表示连接点行和列坐标观测值的残差向量；X＝[△a0 △a1 △a2 △b0△b1 △b2]T为影像RFM的像方坐标系统误差仿射变换补偿模型参数的改正数向量，Y＝[△lat △lon △h]T为连接点对应地面点三维坐标的改正数向量；
为未知数X的系数矩阵；
为未知数Y的系数矩阵；
为初值计算得到的常数项；P为权矩阵。
8.根据权利要求7所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S6具体包括如下内容，
S61、根据最小二乘平差原理，对误差方程进行法化，得到的法方程如下，
记为，
S62、对误差方程进行变换消去X，只求解其中的仿射变换参数的改正数，然后利用前方交会的方式更新连接点对应地面点的三维坐标，将上式变换为，
Nt＝G
其中，
N中
为2阶对角矩阵，在计算时通过对每个小矩阵求逆的方式得到；
S63、求解出
之后便能够分别得到N、G，针对公式Nt＝G，利用共轭梯度下降法进行迭代求解，当两次求解得到的t的差值小于设定的阈值或者求解次数超过预设次数时结束迭代，输出得到最终的t，即仿射变换参数的改正数。
9.根据权利要求5所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S7中，利用高程约束规则对每次平差迭代计算后获取的改正后的连接点对应地面点高程值进行修改，具体为，
针对平地地形区域的连接点，采用下式进行高程约束，
针对丘陵地形区域的连接点，采用下式进行高程约束，
针对山地地形区域的连接点，采用下式进行高程约束，
10.根据权利要求2所述的全球数字高程模型约束的影像区域网平差方法，其特征在于：步骤S8具体包括如下内容，
S81、采用影像原始RFM计算影像四个角点对应的地面范围，利用全球数字高程模型计算该影像覆盖区域的最大和最小高程值；将影像覆盖区域在高程方向以一定的间隔分层，在平面上以一定的格网大小建立地面规则格网，构建一个虚拟的地面点立体空间格网，作为控制点格网；利用影像原始RFM以及平差获取的仿射变换参数，计算所有控制点在影像上的影像坐标；加密控制点格网和层，建立独立的检查点；
S82、采用最小二乘平差原理求解RFM参数；首先将RFM方程变形为
Fy＝NumL(P,L,H)-y*DenL(P,L,H)＝0
Fx＝NumS(P,L,H)-x*DenS(P,L,H)＝0
则获取线性化误差方程为，
V＝Bx-l,W
其中，
x＝[a1i a2j a3ia4j]T，W为权矩阵；
根据最小二乘平差原理，求解获取影像新RFM参数为，
x＝(BTB)-1BTl
S83、利用影像新RFM计算检查点对应的影像坐标，比较与利用影像原始RFM以及仿射变换参数计算的检查点对应影像坐标之间的差值，评定影像新RFM精度，在达到精度要求时输出影像新的RFM参数，为每张影像生成新RFM。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110785058.1A CN113514035B (zh) | 2021-07-12 | 2021-07-12 | 一种全球数字高程模型约束的影像区域网平差方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110785058.1A CN113514035B (zh) | 2021-07-12 | 2021-07-12 | 一种全球数字高程模型约束的影像区域网平差方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113514035A true CN113514035A (zh) | 2021-10-19 |
| CN113514035B CN113514035B (zh) | 2022-03-01 |
Family
ID=78067473
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110785058.1A Active CN113514035B (zh) | 2021-07-12 | 2021-07-12 | 一种全球数字高程模型约束的影像区域网平差方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113514035B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114964169A (zh) * | 2022-05-13 | 2022-08-30 | 中国科学院空天信息创新研究院 | 像方物方协同改正的遥感影像平差方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101604018A (zh) * | 2009-07-24 | 2009-12-16 | 中国测绘科学研究院 | 高分辨率遥感影像数据处理方法及其系统 |
| US20110249860A1 (en) * | 2010-04-12 | 2011-10-13 | Liang-Chien Chen | Integrating and positioning method for high resolution multi-satellite images |
| CN103823981A (zh) * | 2014-02-28 | 2014-05-28 | 武汉大学 | 一种数字高程模型辅助的卫星影像区域网平差方法 |
| CN110388898A (zh) * | 2019-06-27 | 2019-10-29 | 中国科学院遥感与数字地球研究所 | 构建虚拟控制点约束的多源多重覆盖遥感影像平差方法 |
-
2021
- 2021-07-12 CN CN202110785058.1A patent/CN113514035B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101604018A (zh) * | 2009-07-24 | 2009-12-16 | 中国测绘科学研究院 | 高分辨率遥感影像数据处理方法及其系统 |
| US20110249860A1 (en) * | 2010-04-12 | 2011-10-13 | Liang-Chien Chen | Integrating and positioning method for high resolution multi-satellite images |
| CN103823981A (zh) * | 2014-02-28 | 2014-05-28 | 武汉大学 | 一种数字高程模型辅助的卫星影像区域网平差方法 |
| CN110388898A (zh) * | 2019-06-27 | 2019-10-29 | 中国科学院遥感与数字地球研究所 | 构建虚拟控制点约束的多源多重覆盖遥感影像平差方法 |
Non-Patent Citations (1)
| Title |
|---|
| 刘佳 等: "基于有理多项式模型区域网平差的GF-1影像几何校正", 《农业工程学报》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114964169A (zh) * | 2022-05-13 | 2022-08-30 | 中国科学院空天信息创新研究院 | 像方物方协同改正的遥感影像平差方法 |
| CN114964169B (zh) * | 2022-05-13 | 2023-05-30 | 中国科学院空天信息创新研究院 | 像方物方协同改正的遥感影像平差方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113514035B (zh) | 2022-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110388898B (zh) | 构建虚拟控制点约束的多源多重覆盖遥感影像平差方法 | |
| KR101965965B1 (ko) | 위성영상과 제공 rpc로부터 제작된 수치표고모델의 자동 오차보정 방법 | |
| Yang et al. | Large-scale block adjustment without use of ground control points based on the compensation of geometric calibration for ZY-3 images | |
| CN112017224B (zh) | Sar数据区域网平差处理方法和系统 | |
| CN102506824B (zh) | 一种城市低空无人机系统生成数字正射影像图的方法 | |
| Li et al. | Rigorous photogrammetric processing of HiRISE stereo imagery for Mars topographic mapping | |
| CN105528500B (zh) | 一种分米级星载tdi ccd立体测绘相机成像仿真方法和系统 | |
| CN113358091B (zh) | 一种利用三线阵立体卫星影像生产数字高程模型dem的方法 | |
| CN111724465B (zh) | 基于平面约束优选虚拟控制点的卫星影像平差方法及装置 | |
| CN112529946B (zh) | 高分立体模型优化方法、系统、电子设备及可读存储介质 | |
| CN109919835B (zh) | 基于多源卫星遥感影像联合平差的境外电力选线方法 | |
| CN106780321A (zh) | 一种cbers‑02卫星hr传感器影像整体严密定向与纠正拼接方法 | |
| US8855439B2 (en) | Method for determining a localization error in a georeferenced image and related device | |
| CN109709551B (zh) | 一种星载合成孔径雷达影像的区域网平面平差方法 | |
| CN108961199A (zh) | 多源遥感数据时空融合方法及装置 | |
| Zhang et al. | Block adjustment for satellite imagery based on the strip constraint | |
| CN114838740B (zh) | 一种考虑不同经纬度区域的卫星图像几何定标方法 | |
| CN110660099B (zh) | 基于神经网络的遥感影像处理的有理函数模型拟合方法 | |
| CN113514035B (zh) | 一种全球数字高程模型约束的影像区域网平差方法 | |
| Zhang et al. | Simplified SEBAL method for estimating vast areal evapotranspiration with MODIS data | |
| CN113720351A (zh) | 星载激光测高数据与遥感立体影像的联合区域网平差方法 | |
| CN113324527A (zh) | 一种同轨激光测高点与三线阵立体影像联合测绘处理方法 | |
| Zhou et al. | DEM generation using Ziyuan-3 mapping satellite imagery without ground control points | |
| KR102015817B1 (ko) | 입체 위성영상의 제공 rpc 자동 보정 방법 | |
| CN109579796A (zh) | 一种投影后影像的区域网平差方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |