CN112883325A - 一种单像空间后方交会数值计算方法 - Google Patents
一种单像空间后方交会数值计算方法 Download PDFInfo
- Publication number
- CN112883325A CN112883325A CN202110053783.XA CN202110053783A CN112883325A CN 112883325 A CN112883325 A CN 112883325A CN 202110053783 A CN202110053783 A CN 202110053783A CN 112883325 A CN112883325 A CN 112883325A
- Authority
- CN
- China
- Prior art keywords
- equation
- matrix
- formula
- image
- calculating
- 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
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Multimedia (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Image Analysis (AREA)
Abstract
本发明提供了一种单像空间后方交会数值计算方法，从像片中选取三个不在同一条直线上的控制点，根据三条摄影光线在像方空间和物方空间相应光线间顶角相等的原理，构建三元二次方程组；将六个外方位元素的待定值分成两个步骤来求解，从而可以确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。该方法无需提供任何初始值，即可准确估计出像片的外方位元素，且针对竖直摄影和大倾角摄影等情况均适用，具有计算过程相对简便、计算速度快、计算结果准确、稳定可靠等优点。
Description
技术领域
本发明涉及摄影测量技术领域，具体涉及一种单像空间后方交会数值计算方法。
背景技术
空间后方交会，是指利用航摄像片上三个以上不在一条直线上的控制点按共线方程计算该像片外方位元素的方法，是单幅影像解析过程中的一个步骤，是摄影测量领域的基础理论之一，其目标是通过拍摄空间中的全局控制点确定相机的外方位元素。空间后方交会的原理为影像的共线条件方程，利用内参数已知的相机拍摄空间中若干个坐标已知的全局控制点，然后结合像点坐标和控制点坐标建立共线条件约束方程，即可求解相机在当前位姿状态下的外部参数，实现对相机的定位。空间后方交会对布点空间要求不高，接收成本较低，可以快速确定受限空间内的姿态信息，具有结构简单、方便灵活等优点，在摄影测量领域得到了广泛的应用。
单幅影像的空间后方交会，如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统(GPS)、惯性导航系统(INS)以及星相摄影机来获取影像的外方位元素；也可利用影像覆盖范围内一定数量的控制点的空间坐标与影像坐标，根据共线条件方程反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。
目前，解决单像空间后方交会问题主要有共线条件方程解法、角锥体法和直接线性变换法。角锥体法是通过确定摄影光束的外方位元素，将六个外方位元素的待定值分成两个步骤来解求，从而确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。但传统的单像空间后方交会计算依赖于初始值，并且像片的外方位元素计算过程相对复杂，存在计算速度慢、计算结果不够准确、稳定性差等问题。
发明内容
针对现有技术存在的上述不足，本发明的目的在于：提供一种单像空间后方交会数值计算方法，基于摄影光束角锥体在像方空间和物方空间相应光线间顶角相等的原理，将六个外方位元素的待定值分成两个步骤来求解，从而可以确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。该方法无需提供任何初始值，即可准确估计出像片的外方位元素，且针对竖直摄影和大倾角摄影等情况均适用，具有计算过程相对简便、计算速度快、计算结果准确、稳定可靠等优点。
一种单像空间后方交会数值计算方法，包括以下步骤：
从像片中选取三个不在同一条直线上的控制点，分别为A点、B点和C点；三个控制点和摄站点S构成角锥体和三条摄影光线，三个控制点在三条摄影光线上的像点分别为a、b和c；根据三条摄影光线在像方空间和物方空间相应光线间顶角相等的原理，构建三元二次方程组；
针对三元二次方程组，采用奇异值分解得到三个奇异值，若空间后方交会存在解，则三个奇异值中有一个为零，其余两个符号相反，得到自由变量t构成的两组解；根据自由变量t构成的两组解计算角锥体的实际棱长，基于角锥体的实际棱长和边长公式将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程；根据线性方程组计算得到线性相关的无穷多组解，并由约束方程约束得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，线元素解和角元素解一一对应；
基于共线条件方程，根据三个像点构建旋转矩阵；根据旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；根据旋转矩阵R计算旋转角度，得到外方位角元素；其中，RT为旋转矩阵R的转置矩阵。
进一步地，所述摄站点S的坐标为(XS,YS,ZS)，三个控制点的坐标分别为A(XA,YA,ZA)、B(XB,YB,ZB)和C(XC,YC,ZC)，三个像点的坐标分别为a(xa,ya)、b(xb,yb)和c(xc,yc)，三条摄影光线间所构成角锥体内顶角的余弦为：
设定各边长AB＝l3，BC＝l1，CA＝l2，SA＝u，SB＝v，SC＝w，从而构建三元二次方程组(公式1)：
l1 2＝v2+w2-2vwcos(α)
l2 2＝w2+u2-2wucos(β)
l3 2＝v2+u2-2vucos(γ)
其中，l1、l2和l3根据三个控制点的坐标计算，cos(α)、cos(β)和cos(γ)根据三个像点坐标和相机焦距f计算。
进一步地，针对所述三元二次方程组进行奇异值分解之前，将三元二次方程组进行等价变换，构建齐次的三元二次方程，并以二次型表示：
XTPX＝0；XTQX＝0 (公式2)
其中，X、P和Q为定义符号，XT表示X的转置，t1和t2表示自由变量。
进一步地，针对所述三元二次方程组，采用奇异值分解得到三个奇异值，若空间后方交会存在解，则三个奇异值中有一个为零，其余两个符号相反，得到自由变量t构成的两组解，具体包括：
根据线性方程组的高斯消元法，用公式2中的第一个方程减去第二个方程的λ倍，得到方程：
XT(P-λQ)X＝0 (公式3)
转换为矩阵(P-λQ)的行列式值等于0；对矩阵(P-λQ)进行奇异值分解SVD，得到关于矩阵(P-λQ)的奇异值(s1,s2,s3)的方程：
s1y1 2+s2y2 2+s3y3 2＝0 (公式4)
其中，y1、y2和y3为定义符号；若空间后方交会存在解，则三个奇异值中必须有一个为零，其余两个符号相反，得到自由变量t构成的两组解；假定s3＝0，则得到：
其中，t为定义符号，是自由变量，代表任意实数；Y表示y1、y2和y3的集合。
进一步地，所述根据自由变量t构成的两组解计算角锥体的实际棱长，具体包括：
基于自由变量t构成的两组解和(公式2)，得到比例尺为λ的4组解(v',u',w')；(v',u',w')的值与实际的(v,u,w)存在一个比例系数λ，根据(v',u',w')和(公式1)得到：
l1 2+l2 2+l3 2＝λ·2(v'2+u'2+w'2-v'w'cos(α)-w'u'cos(β)-u'v'cos(γ))
从而计算出比例系数λ的值，再根据比例系数λ计算实际棱长：
其中，v｀、u｀和w｀均为定义符号，表示任意比例尺下的棱长距离。
进一步地，所述基于角锥体的实际棱长和边长公式将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程，具体包括：
根据角锥体的实际棱长和边长公式，反解外方位线元素，边长公式为：
将物方坐标原点移至控制点A(XA,YA,ZA)处，构建两个线性方程组(公式7)和一个约束方程(公式8)：
XS'2+YS'2+ZS'2＝u2 (公式8)
m＝((XB-XA)2+(YB-YA)2+(ZB-ZA)2-v2+u2)/2
n＝((XC-XA)2+(YC-YA)2+(ZC-ZA)2-w2+u2)/2
其中，XS'、YS'、ZS'、m和n均表示定义符号，XS'、YS'和ZS'表示物方坐标原点平移后的外方位线元素。
进一步地，所述根据线性方程组计算得到线性相关的无穷多组解，并由约束方程约束得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，线元素解和角元素解一一对应；具体包括：
采用奇异值分解SVD计算线性方程组(公式7)，得到线性相关的无穷多组解x；由于(公式7)为两个三元一次方程，无法唯一确定3个未知数，但可获取到3个未知数的线性关系：
根据(公式9)和约束方程(公式8)得到关于t的一元二次方程，计算出t后基于(公式9)并考虑原点偏移后得到：
其中，con1、con2、con3、l1、l2和l3均表示线性系数。
进一步地，所述基于共线条件方程，根据三个像点构建旋转矩阵，具体包括：
像点a、像点b和像点c根据共线条件方程有如下关系：
三个像点联列方程构建旋转矩阵:
Q＝RTN (公式10)
其中，R为旋转矩阵，RT为旋转矩阵R的转置矩阵，Q和N为定义符号，R表示为：
进一步地，所述根据旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；具体包括：
计算(公式10)时需要考虑旋转矩阵的正交性RRT＝I，I为单位阵，将R的计算转化为计算与R相关的可变量最小值，RT为矩阵R的转置矩阵；
考虑矩阵迹的性质tr(AB)＝tr(BA)；对矩阵S进行奇异值分解SVD，S＝UΣVT，tr(RUΣVT)＝tr(ΣVTRU)＝tr(ΣM)；其中，U、Σ和V为定义符号，VT为矩阵V转置矩阵，M矩阵为正交矩阵，M矩阵表示为：
其中，d表示矩阵阶数，σi为奇异值，mii代表M矩阵的主对角线元素，I为单矩阵；为了求得最大迹，R＝VUT，需要mii＝1(i＝1,2,…,d)的其他元素为0，即M为单位阵：
根据(公式12)计算R矩阵，如果其行列式为-1，则不是需要的旋转矩阵；同样，为了求得(公式11)最大的迹，需要mii＝1(i＝1,2,…,d-1)、mdd＝-1且其他元素为0。
进一步地，所述根据旋转矩阵R计算旋转角度，得到外方位角元素，具体包括：
根据(公式11)和(公式12)计算R矩阵，得到：
其中，UT为矩阵U转置矩阵；根据(公式12)和(公式13)得到：
根据式(公式14)中的R，计算出旋转角度(phi,omega,kappa)：
omega＝-arcsin(R[5])
从而计算得到外方位角元素。
相比于现有技术，本发明具有以下优点：
本发明提供了一种单像空间后方交会数值计算方法，基于摄影光束角锥体在像方空间和物方空间相应光线间顶角相等的原理，将六个外方位元素的待定值分成两个步骤来求解，先计算三个外方位线元素，然后计算三个外方位角元素，线元素解和角元素解一一对应，从而可以确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。该方法无需提供任何初始值，即可准确估计出像片的外方位元素，且针对竖直摄影和大倾角摄影等情况均适用，具有计算过程相对简便、计算速度快、计算结果准确、稳定可靠等优点。
附图说明
图1为本发明实施例中一种单像空间后方交会数值计算方法的计算流程图；
图2为本发明实施例中三个控制点和摄站点构成的角锥体示意图；
图3为本发明实施例中采用的具体像片点号列表；
图4为本发明实施例中采用的具体像片示意图。
具体实施方式
下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例仅用于更加清楚地说明本发明的技术方案，因此只是作为示例，而不能以此来限制本发明的保护范围。
实施例：
参照图1和图2，一种单像空间后方交会数值计算方法，包括以下步骤：
从像片中选取三个不在同一条直线上的控制点，分别为A点、B点和C点；三个控制点和摄站点S构成角锥体和三条摄影光线，三个控制点在三条摄影光线上的像点分别为a、b和c；根据三条摄影光线在像方空间和物方空间相应光线间顶角相等的原理，构建三元二次方程组；
针对三元二次方程组，采用奇异值分解得到三个奇异值，若空间后方交会存在解，则三个奇异值中有一个为零，其余两个符号相反，得到自由变量t构成的两组解；根据自由变量t构成的两组解计算角锥体的实际棱长，基于角锥体的实际棱长和边长公式将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程；根据线性方程组计算得到线性相关的无穷多组解，并由约束方程约束得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，线元素解和角元素解一一对应；
基于共线条件方程，根据三个像点构建旋转矩阵；根据旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；根据旋转矩阵R计算旋转角度，得到外方位角元素；其中，RT为旋转矩阵R的转置矩阵。
其中，所述摄站点S的坐标为(XS,YS,ZS)，三个控制点的坐标分别为A(XA,YA,ZA)、B(XB,YB,ZB)和C(XC,YC,ZC)，三个像点的坐标分别为a(xa,ya)、b(xb,yb)和c(xc,yc)，三条摄影光线间所构成角锥体内顶角的余弦为：
设定各边长AB＝l3，BC＝l1，CA＝l2，SA＝u，SB＝v，SC＝w，从而构建三元二次方程组(公式1)：
l1 2＝v2+w2-2vwcos(α)
l2 2＝w2+u2-2wucos(β)
l3 2＝v2+u2-2vucos(γ)
其中，l1、l2和l3根据三个控制点的坐标计算，cos(α)、cos(β)和cos(γ)根据三个像点坐标和相机焦距f计算。
上述方法是基于摄影光束角锥体在像方空间和物方空间相应光线间顶角相等的原理，将六个外方位元素的待定值分成两个步骤来求解，先计算三个外方位线元素，然后计算三个外方位角元素，线元素解和角元素解一一对应，从而可以确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。该方法无需提供任何初始值，即可准确估计出像片的外方位元素，且针对竖直摄影和大倾角摄影等情况均适用，具有计算过程相对简便、计算速度快、计算结果准确、稳定可靠等优点。
参照图3和图4，下面具体以一个像片作为示例来展示上述方法的计算过程：
首选从像片中选取三个控制点，分别为控制点4714、控制点4718和控制点4788，控制点4714的坐标为(3.47687，-0.32870，8.35018)，对应的像点坐标为(2208.95，37.53)；控制点4718的坐标为(2.32350，0.36637，8.51494)，对应的像点坐标为(2926.66，1053.96)；控制点4788的坐标为(2.15176，-0.23152，8.43674)，对应的像点坐标为(2381.19，1258.15)。相机型号为Cannon，相机内方位元素：像幅大小(wx h)：3680x 2456像素，焦距7343像素，像主点(x0,y0):(0,0)。
具体计算步骤如下：
①利用三条摄影光线(束)在像方空间和物方空间相应光线间顶角相等的原理，构建三元二次方程组(公式1)：
l1 2＝v2+w2-2vwcos(α)
l2 2＝w2+u2-2wucos(β)
l3 2＝v2+u2-2vucos(γ)
采用上面的控制点4714、控制点4718和控制点4788计算得到：l1＝1.356664，l2＝0.626958，l3＝1.3314926，cos(α)＝0.988109，cos(β)＝0.997469，cos(γ)＝0.987527。
②将三元二次方程组进行等价变换，构建齐次的三元二次方程，并以二次型表示；具体为，消去公式1中的第一式，即用第一式分别乘以
从而得到齐次的三元二次方程组，以二次型表示为：
XTPX＝0；XTQX＝0 (公式2)
其中，X、P和Q为定义符号，XT表示X的转置，t1和t2表示自由变量。
采用上面的控制点4714、控制点4718和控制点4788计算后得到：
P[0]＝0.2135660704480973；P[1]＝-0.2110266086638881；
P[2]＝0.0000000000000000；P[3]＝-0.2110266086638881；
P[4]＝-0.7864339295519026；P[5]＝0.9974692555440020；
P[6]＝0.0000000000000000；P[7]＝0.9974692555440020；
P[8]＝-1.0000000000000000；
Q[0]＝-0.0367640204593520；Q[1]＝-0.9517823766622202；
Q[2]＝0.9875270219198669；Q[3]＝-0.9517823766622202；
Q[4]＝0.9632359795406479；Q[5]＝0.0000000000000000；
Q[6]＝0.9875270219198669；Q[7]＝0.0000000000000000；
Q[8]＝-1.0000000000000000。
③将公式2变换为二次型的公式3，奇异值分解矩阵(P-λQ)，再带入到公式3中，得到关于矩阵(P-λQ)奇异值(s1,s2,s3)的公式4，若空间后方交会存在解，三个奇异值中必须有一个为零，其余两个符号相反，公式4的解描述为公式5的形式，得到自由变量t构成的两组解。具体如下：
公式2的三元二次齐次方程组有两个方程、三个未知数，由于是齐次的，所以如果Y是方程组的解，那么kY必定也是方程组的解。根据线性方程组的高斯消元法，用公式2中的第一个方程减去第二个方程的λ倍，得到方程：
XT(P-λQ)X＝0 (公式3)
转换为使得矩阵(P-λQ)的行列式值等于0，即：
|P-λQ|＝0
采用奇异值分解SVD：P-λQ＝UTdiag(s1,s2,s3)U，假定中间未知数Y＝UX，则X＝U- 1Y，将公式3转换为：
s1y1 2+s2y2 2+s3y3 2＝0 (公式4)
其中，U、UT、y1、y2和y3均为定义符号；假定Q可逆，则|P-λQ|＝0等价于|PQ-1-λE|＝0，E为单位阵，根据特征值与特征向量的定义，|P-λQ|＝0转换为计算PQ-1的特征值。由于|P-λQ|＝0，则diag(s1,s2,s3)中，某一个特征值必须为0，其余两个符号应该相反，假定s3＝0，则得到：
其中，t为定义符号，是自由变量，代表任意实数；Y表示y1、y2和y3的集合。
采用上面的控制点4714、控制点4718和控制点4788计算，三个奇异值分别为：
s1＝-1.9058299423570109；
s2＝0.33296367871096866；
s3＝-2.1441812215520123e-006；
s3值过小，考虑到计算舍与误差，解有效。y1＝0.41798082386239727或-0.41798082386239727。
④将自由变量t构成的两组解带入到公式2中，得到比例尺为λ的共计4组解，再利用公式1计算角锥体的实际棱长，具体地：
将X＝U-1Y带入XTQX＝0中，整理出关于t的一元二次方程：
(U-1Y)TQ(U-1Y)＝0
YT(UQUT)Y＝0
其中，UQUT为对阵矩阵，XT为X的转置矩阵，YT为Y的转置矩阵。计算出两个t根后再利用X＝U-1Y即可计算出最多4组(v',u',w')值，由于(y1,y2,y3)是为简化计算而设置的，此时(v',u',w')值与实际的(v,u,w)在一个比例系数λ，根据(v',u',w')和(公式1)得到：
l1 2+l2 2+l3 2＝λ·2(v'2+u'2+w'2-v'w'cos(α)-w'u'cos(β)-u'v'cos(γ))
除λ其他值均已知，从而计算出比例系数λ的值，再根据比例系数λ计算实际棱长：
其中，v｀、u｀和w｀均为定义符号，表示任意比例尺下的棱长距离。
采用上面的控制点4714、控制点4718和控制点4788计算，实际数据仅能得到2组(v'w'u')解和比例尺λ，分别为：
v'＝21.076562272540261；u'＝19.789720821225259；w'＝20.443214812124037；
λ＝0.39875789和v'＝38.744740993377377；u'＝40.752064051444385；
w'＝39.945473955792686；λ＝0.2103959087。
对应的实际棱长为：
u＝8.4044455193169636；v＝7.8913073362620194；w＝8.1518932217916245；
和u＝8.1517349898338285；v＝8.5740675487517475；w＝8.4043642926039972。
⑤基于角锥体的实际棱长和边长公式，将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程，具体地：
根据角锥体的实际棱长和边长公式，反解外方位线元素，边长公式为：
为了简化计算，将物方坐标原点移至控制点A(XA,YA,ZA)处，即消去公式6中的第一式，构建两个线性方程组(公式7)和一个约束方程(公式8)：
XS'2+YS'2+ZS'2＝u2 (公式8)
m＝((XB-XA)2+(YB-YA)2+(ZB-ZA)2-v2+u2)/2
n＝((XC-XA)2+(YC-YA)2+(ZC-ZA)2-w2+u2)/2
其中，XS'、YS'、ZS'、m和n均表示定义符号，XS'、YS'和ZS'表示物方坐标原点平移后的外方位线元素。
⑥根据线性方程组计算得到线性相关的无穷多组解，即公式9，再将其带入到约束方程(公式8)后得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，考虑到所有可能的计算舍余误差，共计最大可能24组解，剔除掉重复解后剩余的线元素解用于计算旋转角，每一个线元素解对应一个唯一的角元素解。具体地：
利用奇异值分解SVD计算线性方程组(公式7)，Ax＝b，即min||Ax＝b||2，对A矩阵进行奇异值分解：
考虑到对于一个正交矩阵M来说，||Mv||＝||v||和MMT＝E，则有如下推导：
min||Ax-b||2＝min||UΣVTx-b||2＝min||ΣVTx-UTb||2
令y＝VTxb'＝UTb，则该公示的线性优化问题转换为：min||Σy-b'||2，将y和b'带入到Ax＝b中，Σ为对角矩阵，求得y后再利用y＝VTx即可计算x。另外，由于由于(公式7)为两个三元一次方程，无法唯一确定3个未知数，但可获取到3个未知数的线性关系：
根据(公式9)和约束方程(公式8)得到关于t的一元二次方程，计算出t后基于(公式9)并考虑原点偏移后得到：
其中，con1、con2、con3、l1、l2和l3均表示线性系数，x、x'、b、b'、y、Σ和V为定义符号，VT为矩阵V转置矩阵，MT为矩阵M转置矩阵。
采用上面的控制点4714、控制点4718和控制点4788计算，共计得到4组解，分别为：
第一组：xs1＝1.9526876186281710；ys1＝2.6535657388948155；
zs1＝16.0588928869197680；
第二组：xs2＝1.1867144211787379；ys2＝4.7090544870209907；
zs2＝2.0253236080444870；
第三组：xs3＝4.5924755116375975；ys3＝-3.9462040445349937；
zs3＝15.570104563126696；
第四组：xs4＝3.7602266099948993；ys4＝-1.7128644725196600；
zs4＝0.3222830239902752。
⑦基于共线条件方程，根据三个像点构建旋转矩阵，具体包括：
其中共线条件方程为：
其中，(x，y)为像点，(X,Y,Z)为物点，f为相机焦距，通过旋转矩阵计算6个外方位元素。旋转矩阵为：
像点a、像点b和像点c根据共线条件方程有如下关系：
三个像点联列方程构建旋转矩阵：
Q＝RTN (公式10)
其中，R为旋转矩阵，RT为旋转矩阵R的转置矩阵，Q和N为定义符号，R表示为：
采用上面的控制点4714、控制点4718和控制点4788计算，得到：
Q[0]＝0.28806821859406445；Q[1]＝0.36699367471235789；
Q[2]＝0.30445030346938218；Q[3]＝0.0048947944147075565；
Q[4]＝0.13216329119659970；Q[5]＝0.16086272966024284；
Q[6]＝0.95759738012573947；Q[7]＝0.92078689563937810；
Q[8]＝0.93884673665282381；
基于4组结果，公式10中的N分别为：
N[0]＝0.18137798295197818；N[1]＝0.047006234724193727
N[2]＝0.024434748325130502；N[3]＝-0.35484619503113352；
N[4]＝-0.28983537286491395；N[5]＝-0.35391762552811284；
N[6]＝-0.91716203866721335；N[7]＝-0.95592147717964271；
N[8]＝-0.93495735593385765；
N[0]＝0.27251703498002183；N[1]＝0.14407167202279558；
N[2]＝0.11839736248257809；N[3]＝-0.59941780221055296；
N[4]＝-0.55031043238591426；N[5]＝-0.60606626047546175；
N[6]＝0.75261727593695860；N[7]＝0.82243649075644987；
N[8]＝0.78655308305953431；
N[0]＝-0.13683056653349099；N[1]＝-0.26461737503042015；
N[2]＝-0.29039655455122010；N[3]＝0.44376863012247131；
N[4]＝0.50298061622650381；N[5]＝0.44199310756790833；
N[6]＝-0.88563355795800180；N[7]＝-0.82279289285482815；
N[8]＝-0.84871192637269077；
N[0]＝-0.034735871025922335；N[1]＝-0.16755155453002249；
N[2]＝-0.19137075406588944；N[3]＝0.16979755834638624；
N[4]＝0.24250452838364825；N[5]＝0.17625741427285804；
N[6]＝0.98486659423684197；N[7]＝0.95557209580857738；
N[8]＝0.96556230166784096。
⑧考虑到旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；具体地：
计算时需考虑旋转矩阵的正交性RRT＝E。解算线性方程组Ax＝b，等价于最小化||Ax＝b2||，即min||Ax＝b2||。对于公式10中的一个线性方程组：
||Rxi-yi||2＝(Rxi-yi)T(Rxi-yi)＝xi Txi-yi TRxi-xi TRyi+yi Tyi
对于标量yi TRxi和xi TRyi来说，标量的转置任然等于标量本身，因而
||Rxi-yi||2＝xi Txi-2yi TRxi+yi Tyi
该公式中xi Txi和yi Tyi为常量，问题转换为求与R相关的可变量最小值。
考虑矩阵迹的性质tr(AB)＝tr(BA)，tr(YTRX)＝tr(RYTX)＝tr(RS)；对矩阵S进行奇异值分解SVD，S＝UΣVT，tr(RUΣVT)＝tr(ΣVTRU)＝tr(ΣM)；M矩阵为正交矩阵：
其中，d表示矩阵阶数，σi为奇异值，mii代表M矩阵的主对角线元素，N、xi、yi为定义符号，RT为矩阵R的转置矩阵，xi T为xi的装置矩阵，yi T为yi的装置矩阵。
为了求得最大迹，R＝VUT，需要mii＝1(i＝1,2,…,d)的其他元素为0，即M为单位阵：
根据(公式12)计算R矩阵，如果其行列式为-1，则不是需要的旋转矩阵；同样，为了求得(公式11)最大的迹，需要mii＝1(i＝1,2,…,d-1)、mdd＝-1且其他元素为0。
采用上面的控制点4714、控制点4718和控制点4788计算，S矩阵分别为：
S[0]＝-0.32083331916948793；S[1]＝-0.29728026276011521；
S[2]＝-0.89927216496673890；S[3]＝-0.099706595171477966；
S[4]＝-0.93358833135117070；S[5]＝0.34419677866628434；
S[6]＝-0.94187290872495122；S[7]＝0.20009316064826682；
S[8]＝0.26988544027380629；
S[0]＝0.16742300207862593；S[1]＝-0.55915072321452453；
S[2]＝0.75810043275477623；S[3]＝0.039420624116983355；
S[4]＝-0.17315833784934129；S[5]＝0.23890689622085931；
S[6]＝0.50477788376301791；S[7]＝-1.6497228825149117；
S[8]＝2.2164458700934393；
S[0]＝-0.22494075956525333；S[1]＝0.44699127913578440；
S[2]＝-0.81547327218633914；S[3]＝-0.082356443133693530；
S[4]＝0.13974794763512155；S[5]＝-0.24960412805734980；
S[6]＝-0.64732266089573010；S[7]＝1.3030534244324219；
S[8]＝-2.4025077108512929；
S[0]＝-0.12975954533901612；S[1]＝0.19157253142446784；
S[2]＝0.92836341600935879；S[3]＝-0.053098611715459994；
S[4]＝0.061234569524275244；S[5]＝0.28643526016571241；
S[6]＝-0.36721006279161061；S[7]＝0.55137138709405542；
S[8]＝2.7294989500295959；
计算出R反带入公式10中，提出第1、第3组解，得两组旋转矩阵：
第一组：
R[0]＝0.049384753562629391；R[1]＝-0.96028169927131035；
R[2]＝0.27462739149648602；R[3]＝0.57088223126319548；
R[4]＝-0.19848056898269367；R[5]＝-0.79667994939264208；
R[6]＝0.81954537650060322；R[7]＝0.19612374099254803；
R[8]＝0.53840584513503709；
第二组：
R[0]＝-0.078584720664683638；R[1]＝-0.99687886486717947；
R[2]＝-0.0075478777929717522；R[3]＝0.98977022431484907；
R[4]＝-0.077115504405732338；R[5]＝-0.12003375375278946；
R[6]＝0.11907705378362887；R[7]＝-0.016903483705247080；
R[8]＝0.99274111806695975；
通过剩余得控制点剔除第一组旋转矩阵，得到最后R t结果为：(第4组线元素)
t[0]＝3.7602266099948993；
t[1]＝-1.7128644725196600；
t[2]＝0.3222830239902752；
⑨根据旋转矩阵R计算旋转角度，得到外方位角元素，具体地：
根据(公式11)和(公式12)计算R矩阵，得到：
其中，UT为矩阵U转置矩阵；根据(公式12)和(公式13)得到：
根据式(公式14)中的R，计算出旋转角度(phi,omega,kappa)：
omega＝-arcsin(R[5])
从而计算得到外方位角元素。
采用上面的控制点4714、控制点4718和控制点4788计算，计算公式为：
phi＝atan2(-R[2],R[8])；
omega＝-asin(R[5])；
kappa＝atan2(R[3],R[4])；
结果为：
phi＝0.007602921065232400；
omega＝0.12032388190105718；
kappa＝1.6485517782118468。
通过上面的示例，可以看出，本申请提供的方法是基于摄影光束角锥体在像方空间和物方空间相应光线间顶角相等的原理，将六个外方位元素的待定值分成两个步骤来求解，从而可以确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。每个像点、摄影中心、物点构成一条空间光线，当存在2条以上光线时，利用像点坐标和外方位元素，通过前方交会计算出物点坐标，当确定出被拍摄物体多个关键点在影像上的像点后，即可获取每个关键点的物方坐标，通过非接触的摄影方式实现量测的目的。该方法无需提供任何初始值，即可准确估计出像片的外方位元素，且针对竖直摄影和大倾角摄影等情况均适用，具有计算过程相对简便、计算速度快、计算结果准确、稳定可靠等优点。
上述数值计算方法的应用场景如下：根据影像覆盖范围内至少3个已知三维坐标(X Y Z)的点和同一影像上对应的像点，以共线条件方程为解析基础(或任何描述共线关系的模型)，确定影像的6个外方位元素，包括3个线元素和3个角元素，进而确定拍摄影像的摄影中心在拍摄瞬间的位置和姿态信息。
具体实施时，还可以将上述方法应用在其他地方，比如摄影测量中利用至少3个控制点计算影像的外方位元素。运动恢复结构(Structure-from-Motion，SfM)中估计相片的初始位姿信息或投影矩阵。同步定位与地图构建(Simultaneous Localization andMapping，SLAM)中建造增量式地图时估计机器人的位置，实现自主导航与定位。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的保护范围当中。
Claims (10)
1.一种单像空间后方交会数值计算方法，其特征在于，包括以下步骤：
从像片中选取三个不在同一条直线上的控制点，分别为A点、B点和C点；三个控制点和摄站点S构成角锥体和三条摄影光线，三个控制点在三条摄影光线上的像点分别为a、b和c；根据三条摄影光线在像方空间和物方空间相应光线间顶角相等的原理，构建三元二次方程组；
针对三元二次方程组，采用奇异值分解得到三个奇异值，若空间后方交会存在解，则三个奇异值中有一个为零，其余两个符号相反，得到自由变量t构成的两组解；根据自由变量t构成的两组解计算角锥体的实际棱长，基于角锥体的实际棱长和边长公式将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程；根据线性方程组计算得到线性相关的无穷多组解，并由约束方程约束得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，线元素解和角元素解一一对应；
基于共线条件方程，根据三个像点构建旋转矩阵；根据旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；根据旋转矩阵R计算旋转角度，得到外方位角元素；其中，RT为旋转矩阵R的转置矩阵。
2.根据权利要求1所述的单像空间后方交会数值计算方法，其特征在于，所述摄站点S的坐标为(XS,YS,ZS)，三个控制点的坐标分别为A(XA,YA,ZA)、B(XB,YB,ZB)和C(XC,YC,ZC)，三个像点的坐标分别为a(xa,ya)、b(xb,yb)和c(xc,yc)，三条摄影光线间所构成角锥体内顶角的余弦为：
设定各边长AB＝l3，BC＝l1，CA＝l2，SA＝u，SB＝v，SC＝w，从而构建三元二次方程组(公式1)：
l1 2＝v2+w2-2vwcos(α)
l2 2＝w2+u2-2wucos(β)
l3 2＝v2+u2-2vucos(γ)
其中，l1、l2和l3根据三个控制点的坐标计算，cos(α)、cos(β)和cos(γ)根据三个像点坐标和相机焦距f计算。
3.根据权利要求2所述的单像空间后方交会数值计算方法，其特征在于，针对所述三元二次方程组进行奇异值分解之前，将三元二次方程组进行等价变换，构建齐次的三元二次方程，并以二次型表示：
XTPX＝0；XTQX＝0 (公式2)
其中，X、P和Q为定义符号，XT表示X的转置，t1和t2表示自由变量。
4.根据权利要求3所述的单像空间后方交会数值计算方法，其特征在于，针对所述三元二次方程组，采用奇异值分解得到三个奇异值，若空间后方交会存在解，则三个奇异值中有一个为零，其余两个符号相反，得到自由变量t构成的两组解，具体包括：
根据线性方程组的高斯消元法，用公式2中的第一个方程减去第二个方程的λ倍，得到方程：
XT(P-λQ)X＝0 (公式3)
转换为矩阵(P-λQ)的行列式值等于0；对矩阵(P-λQ)进行奇异值分解SVD，得到关于矩阵(P-λQ)的奇异值(s1,s2,s3)的方程：
s1y1 2+s2y2 2+s3y3 2＝0 (公式4)
其中，y1、y2和y3为定义符号；若空间后方交会存在解，则三个奇异值中必须有一个为零，其余两个符号相反，得到自由变量t构成的两组解；假定s3＝0，则得到：
其中，t为定义符号，是自由变量，代表任意实数；Y表示y1、y2和y3的集合。
5.根据权利要求4所述的单像空间后方交会数值计算方法，其特征在于，所述根据自由变量t构成的两组解计算角锥体的实际棱长，具体包括：
基于自由变量t构成的两组解和(公式2)，得到比例尺为λ的4组解(v',u',w')；(v',u',w')的值与实际的(v,u,w)存在一个比例系数λ，根据(v',u',w')和(公式1)得到：
l1 2+l2 2+l3 2＝λ·2(v'2+u'2+w'2-v'w'cos(α)-w'u'cos(β)-u'v'cos(γ))
从而计算出比例系数λ的值，再根据比例系数λ计算实际棱长：
其中，v’、u’和w’均为定义符号，表示任意比例尺下的棱长距离。
6.根据权利要求5所述的单像空间后方交会数值计算方法，其特征在于，所述基于角锥体的实际棱长和边长公式将物方坐标原点移至控制点A(XA,YA,ZA)处，构建线性方程组和约束方程，具体包括：
根据角锥体的实际棱长和边长公式，反解外方位线元素，边长公式为：
将物方坐标原点移至控制点A(XA,YA,ZA)处，构建两个线性方程组(公式7)和一个约束方程(公式8)：
XS'2+YS'2+ZS'2＝u2 (公式8)
m＝((XB-XA)2+(YB-YA)2+(ZB-ZA)2-v2+u2)/2
n＝((XC-XA)2+(YC-YA)2+(ZC-ZA)2-w2+u2)/2
其中，XS'、YS'、ZS'、m和n均表示定义符号，XS'、YS'和ZS'表示物方坐标原点平移后的外方位线元素。
7.根据权利要求6所述的单像空间后方交会数值计算方法，其特征在于，所述根据线性方程组计算得到线性相关的无穷多组解，并由约束方程约束得到唯一解；物方坐标原点偏移到控制点A(XA,YA,ZA)后计算得到外方位线元素，线元素解和角元素解一一对应；具体包括：
采用奇异值分解SVD计算线性方程组(公式7)，得到线性相关的无穷多组解x；由于(公式7)为两个三元一次方程，无法唯一确定3个未知数，但可获取到3个未知数的线性关系：
根据(公式9)和约束方程(公式8)得到关于t的一元二次方程，计算出t后基于(公式9)并考虑原点偏移后得到：
其中，con1、con2、con3、l1、l2和l3均表示线性系数。
8.根据权利要求7所述的单像空间后方交会数值计算方法，其特征在于，所述基于共线条件方程，根据三个像点构建旋转矩阵，具体包括：
像点a、像点b和像点c根据共线条件方程有如下关系：
三个像点联列方程构建旋转矩阵：
Q＝RTN (公式10)
其中，R为旋转矩阵，RT为旋转矩阵R的转置矩阵，Q和N为定义符号，R表示为：
9.根据权利要求8所述的单像空间后方交会数值计算方法，其特征在于，所述根据旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为与R相关的可变量计算；根据矩阵迹的性质以及最大迹的条件，可知R即为求解的旋转矩阵；具体包括：
计算(公式10)时需要考虑旋转矩阵的正交性RRT＝E，E为单位阵，将R的计算转化为计算与R相关的可变量最小值，RT为矩阵R的转置矩阵；
考虑矩阵迹的性质tr(AB)＝tr(BA)；对矩阵S进行奇异值分解SVD，S＝UΣVT，tr(RUΣVT)＝tr(ΣVTRU)＝tr(ΣM)；其中，U、Σ和V为定义符号，VT为矩阵V转置矩阵，M矩阵为正交矩阵，M矩阵表示为：
其中，d表示矩阵阶数，σi为奇异值，mii代表M矩阵的主对角线元素；为了求得最大迹，R＝VUT，需要mii＝1(i＝1,2,…,d)的其他元素为0，即M为单位阵：
根据(公式12)计算R矩阵，如果其行列式为-1，则不是需要的旋转矩阵；同样，为了求得(公式11)最大的迹，需要mii＝1(i＝1,2,…,d-1)、mdd＝-1且其他元素为0。
10.根据权利要求9所述的单像空间后方交会数值计算方法，其特征在于，所述根据旋转矩阵R计算旋转角度，得到外方位角元素，具体包括：
根据(公式11)和(公式12)计算R矩阵，得到：
其中，UT为矩阵U转置矩阵；根据(公式12)和(公式13)得到：
根据式(公式14)中的R，计算出旋转角度(phi,omega,kappa)：
omega＝-arcsin(R[5])
从而计算得到外方位角元素。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110053783.XA CN112883325B (zh) | 2021-01-15 | 2021-01-15 | 一种单像空间后方交会数值计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110053783.XA CN112883325B (zh) | 2021-01-15 | 2021-01-15 | 一种单像空间后方交会数值计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112883325A true CN112883325A (zh) | 2021-06-01 |
| CN112883325B CN112883325B (zh) | 2023-09-22 |
Family
ID=76048077
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110053783.XA Active CN112883325B (zh) | 2021-01-15 | 2021-01-15 | 一种单像空间后方交会数值计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112883325B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102889882A (zh) * | 2012-09-03 | 2013-01-23 | 北京信息科技大学 | 一种基于光束平差的三维重建方法 |
| KR20130121290A (ko) * | 2012-04-27 | 2013-11-06 | 서울시립대학교 산학협력단 | 회전식 라인 카메라로 획득한 실내 전방위 영상의 지오레퍼런싱 방법 |
| CN105910586A (zh) * | 2016-04-06 | 2016-08-31 | 华北理工大学 | 基于具有时间属性的像片获取实际地理信息的方法 |
| CN108594255A (zh) * | 2018-04-20 | 2018-09-28 | 武汉大学 | 一种激光测距辅助光学影像联合平差方法及系统 |
| CN109883399A (zh) * | 2019-03-15 | 2019-06-14 | 武汉理工大学 | 一种基于焦距修正的重叠多片摄影测量交替趋近算法 |
| CN111457896A (zh) * | 2020-04-20 | 2020-07-28 | 中国人民解放军空军航空大学 | 一种单像空间后方交会非迭代方法 |
-
2021
- 2021-01-15 CN CN202110053783.XA patent/CN112883325B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20130121290A (ko) * | 2012-04-27 | 2013-11-06 | 서울시립대학교 산학협력단 | 회전식 라인 카메라로 획득한 실내 전방위 영상의 지오레퍼런싱 방법 |
| CN102889882A (zh) * | 2012-09-03 | 2013-01-23 | 北京信息科技大学 | 一种基于光束平差的三维重建方法 |
| CN105910586A (zh) * | 2016-04-06 | 2016-08-31 | 华北理工大学 | 基于具有时间属性的像片获取实际地理信息的方法 |
| CN108594255A (zh) * | 2018-04-20 | 2018-09-28 | 武汉大学 | 一种激光测距辅助光学影像联合平差方法及系统 |
| CN109883399A (zh) * | 2019-03-15 | 2019-06-14 | 武汉理工大学 | 一种基于焦距修正的重叠多片摄影测量交替趋近算法 |
| CN111457896A (zh) * | 2020-04-20 | 2020-07-28 | 中国人民解放军空军航空大学 | 一种单像空间后方交会非迭代方法 |
Non-Patent Citations (2)
| Title |
|---|
| 张一平;龚志辉;王勃: "基于大倾角的角锥体空间后方交会研究", 测绘工程, vol. 21, no. 2, pages 17 - 20 * |
| 陈继华;孙荣旭;张斌: "单像空间后方交会算法的研究", 激光与光电子学进展, no. 002, pages 272 - 277 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112883325B (zh) | 2023-09-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110296691B (zh) | 融合imu标定的双目立体视觉测量方法与系统 | |
| CN106780321B (zh) | 一种cbers-02卫星hr传感器影像整体严密定向与纠正拼接方法 | |
| US8165728B2 (en) | Method and system for providing a GPS-based position | |
| US9466143B1 (en) | Geoaccurate three-dimensional reconstruction via image-based geometry | |
| CN110264528B (zh) | 一种鱼眼镜头双目像机快速自标定方法 | |
| CN106885585B (zh) | 一种基于光束法平差的星载摄影测量系统一体化检校方法 | |
| CN102519433B (zh) | 一种利用rpc反演星载线阵传感器几何定标参数方法 | |
| CN108845335A (zh) | 一种基于图像和导航信息的无人机地面目标定位方法 | |
| US20090141043A1 (en) | Image mosaicing apparatus for mitigating curling effect | |
| JP2013187862A (ja) | 画像データ処理装置、画像データ処理方法および画像データ処理用のプログラム | |
| CN108663043B (zh) | 基于单个相机辅助的分布式pos主子节点相对位姿测量方法 | |
| CN105809706B (zh) | 一种分布式多像机系统的全局标定方法 | |
| CN108594255B (zh) | 一种激光测距辅助光学影像联合平差方法及系统 | |
| CN110969665A (zh) | 一种外参标定方法、装置、系统及机器人 | |
| CN114972545B (zh) | 一种高光谱卫星的在轨数据快速预处理方法 | |
| CN115423863A (zh) | 相机位姿估计方法、装置及计算机可读存储介质 | |
| Liu et al. | A new approach to fast mosaic UAV images | |
| CN110660099B (zh) | 基于神经网络的遥感影像处理的有理函数模型拟合方法 | |
| CN112857328B (zh) | 一种无标定摄影测量方法 | |
| CN109099888A (zh) | 一种位姿测量方法、设备及存储介质 | |
| El Habchi et al. | CGA: A new approach to estimate the geolocation of a ground target from drone aerial imagery | |
| CN112883325A (zh) | 一种单像空间后方交会数值计算方法 | |
| CN116124094A (zh) | 基于无人机侦察图像与组合导航信息的多目标协同定位方法 | |
| Castanheiro et al. | Modeling Hyperhemispherical Points and Calibrating a Dual-Fish-Eye System for Close-Range Applications | |
| CN104019800A (zh) | 大侧摆线阵ccd遥感图像对地定位的方法 |
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 |