发明名称 一种基于星载可见光相机的空间碎片快速检测定位方法
摘要 本发明公开的一种基于星载可见光相机的空间碎片快速检测定位方法,涉及空间碎片天基可见光相机在轨快速识别、匹配及定位方法,属于航天器制导与控制领域。本发明通过连通域检测方法对星载可见光相机拍摄的图片中碎片轨迹进行编号;对目标几何特性进行识别,并保存碎片成像轨迹几何特征;根据图像采集的时间和碎片成像几何特征,制定帧间匹配准则,进行帧间匹配,使帧与帧之间碎片成像轨迹一一对应;根据帧间匹配结果,解算目标方位角α、β,完成碎片方位角定位。本发明可实现对空间碎片成像的识别、匹配、定位,且能够满足高匹配率、高定位精度。此外,根据帧间匹配结果进行ROI区域选择,提高连通域检测的计算效率,还可解决星上计算资源有限的问题。
申请公布号 CN105910583A 申请公布日期 2016.08.31
申请号 CN201610262261.X 申请日期 2016.04.25
申请人 北京理工大学 发明人 翟光;赵琪;张景瑞;梁斌
分类号 G01C11/00(2006.01)I;G01C11/04(2006.01)I 主分类号 G01C11/00(2006.01)I
代理机构 代理人
主权项 一种基于星载可见光相机的空间碎片快速检测定位方法,其特征在于:包括如下步骤,步骤一:通过连通域检测的方法,对星载可见光相机拍摄的图片中碎片轨迹进行编号;采用8连通域检测方法检测出星载可见光相机拍摄的图片中所有连通域;8连通域检测方法是指遍历灰度矩阵A中像素点A(i,j)周围的8个像素点A(i‑1,j‑1)、A(i‑1,j)、A(i‑1,j+1)、A(i,j‑1)、A(i,j+1)、A(i+1,j‑1)、A(i+1,j)、A(i+1,j+1),若周围的8个像素点中存在满足灰度阈值的像素点,则认为该像素点与A(i,j)连通;遍历图像中每个像素,进行8连通域检测,并对每个连通区域进行编号;步骤二:计算碎片成像几何特征;所述的碎片成像几何特征包括:成像面积m<sub>00</sub>、形心(u<sub>0</sub>,v<sub>0</sub>)、椭圆长半轴a、椭圆短半轴b、椭圆离心率e、椭圆长轴与横轴正向夹角θ;步骤三:制定帧间匹配准则,进行帧间匹配,使帧与帧之间碎片成像轨迹一一对应;步骤3.1:帧间匹配特征选取;由于碎片与相机存在相对运动,导致碎片成像在不同帧中的位置不同,所以在选择匹配特征时,选择在短时间内几乎不发生变化的特征作为匹配特征;步骤二中得到的六个碎片图像特征包括:成像面积m<sub>00</sub>、形心(u<sub>0</sub>,v<sub>0</sub>)、椭圆长半轴a、椭圆短半轴b、椭圆离心率e、椭圆长轴与横轴正向夹角θ,其中,短时间内几乎不变的特征包括:成像面积m<sub>00</sub>、椭圆长半轴a、椭圆短半轴b、椭圆离心率e、椭圆长轴与横轴正向夹角θ;上述五个特征中,选取成像面积m<sub>00</sub>、椭圆长半轴a、椭圆长轴与横轴正向夹角θ三个特征作为选取的帧间匹配特征;为了对碎片成像位置进行约束,结合形心(u<sub>0</sub>,v<sub>0</sub>)、椭圆半长轴a、椭圆长轴与横轴正向夹角θ三个特征计算相邻帧同一碎片的形心偏移距离阈值;步骤3.2:根据步骤3.1选取的特征设定用于步骤3.3帧间匹配的阈值;所述的阈值包括相邻帧之间碎片成像形心偏移距离阈值l<sub>thr</sub>、面积变化阈值s<sub>thr</sub>、椭圆长轴与横轴正向夹角变化阈值θ<sub>thr</sub>、半长轴变化阈值a<sub>thr</sub>;记相机曝光时间为t<sub>ep</sub>,拍照间隔时间为t<sub>in</sub>,相邻帧之间碎片成像形心偏移距离阈值为l<sub>thr</sub>,面积变化阈值为s<sub>thr</sub>,椭圆长轴与横轴正向夹角变化阈值为θ<sub>thr</sub>,半长轴变化阈值为a<sub>thr</sub>;碎片在曝光时间内移动的距离近似等于椭圆长轴2a,所以相邻帧碎片形心移动距离可近似表示为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mi>&Delta;</mi><mi>l</mi><mo>&ap;</mo><mfrac><mrow><mn>2</mn><mi>a</mi></mrow><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub></mfrac><mo>&CenterDot;</mo><mrow><mo>(</mo><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub><mo>+</mo><msub><mi>t</mi><mrow><mi>i</mi><mi>n</mi></mrow></msub><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000021.GIF" wi="1110" he="119" /></maths>因此,设定相邻帧之间碎片成像形心偏移距离阈值l<sub>thr</sub>:l<sub>thr</sub>=Δl+δ<sub>l</sub>   (2)式(2)中,δ<sub>l</sub>为形心偏移阈值冗余量;对于面积变化阈值s<sub>thr</sub>、椭圆长轴与横轴正向夹角变化阈值θ<sub>thr</sub>、椭圆半长轴变化阈值a<sub>thr</sub>,由于在短时间内碎片成像的几何特征变化很小,所以,近似取相邻两帧s<sub>thr</sub>、θ<sub>thr</sub>、a<sub>thr</sub>为:s<sub>str</sub>=δ<sub>s</sub>   (3)θ<sub>thr</sub>=δ<sub>θ</sub>   (4)a<sub>thr</sub>=δ<sub>a</sub>   (5)其中,δ<sub>s</sub>、δ<sub>θ</sub>、δ<sub>a</sub>分别为面积、夹角、半长轴冗余量;步骤3.3:进行帧间匹配,使帧与帧之间碎片成像轨迹一一对应;对于步骤3.2设定的阈值,按影响程度进行排序;首先,满足形心偏移阈值l<sub>thr</sub>是匹配的先决条件;对于满足该条件的,根据面积阈值s<sub>thr</sub>、椭圆半长轴阈值a<sub>thr</sub>、椭圆长轴与横轴夹角阈值进行匹配θ<sub>thr</sub>;因不同的特征对匹配的影响程度不同,所以引入匹配度的概念,对各个特征分配不同的匹配权重k<sub>i</sub>;在满足特征阈值的前提下,根据当前帧碎片成像特征与前一帧图像特征的误差,确定特征匹配度mch<sub>i</sub>;最后,计算得到所有特征的匹配度Mch:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>M</mi><mi>c</mi><mi>h</mi><mo>=</mo><munderover><mo>&Sigma;</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><msub><mi>k</mi><mi>i</mi></msub><mo>&CenterDot;</mo><msub><mi>mch</mi><mi>i</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000022.GIF" wi="1108" he="115" /></maths>若匹配度Mch大于设定的阈值,则认为相邻两帧中两个像为同一目标;步骤四:解算目标方位角α、β,完成碎片方位角定位;因碎片在短时间内的相对运动可视为匀速直线运动,所以根据相邻帧图像中同一碎片所成像的形心计算碎片在曝光结束时刻所在(u,v)坐标;上一帧碎片形心为(u<sub>c</sub>′,v<sub>c</sub>′),当前帧碎片形心为(u<sub>c</sub>,v<sub>c</sub>),相机曝光时间t<sub>ep</sub>、拍照间隔t<sub>in</sub>,得拍照结束时刻目标坐标(u<sub>t</sub>,v<sub>t</sub>)为:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><msub><mi>u</mi><mi>t</mi></msub><mo>=</mo><msub><mi>u</mi><mi>c</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>u</mi><mi>c</mi></msub><mo>-</mo><msubsup><mi>u</mi><mi>c</mi><mo>&prime;</mo></msubsup></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>n</mi></mrow></msub><mo>+</mo><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub></mrow></mfrac><mo>&CenterDot;</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>v</mi><mi>t</mi></msub><mo>=</mo><msub><mi>v</mi><mi>c</mi></msub><mo>+</mo><mfrac><mrow><msub><mi>v</mi><mi>c</mi></msub><mo>-</mo><msubsup><mi>v</mi><mi>c</mi><mo>&prime;</mo></msubsup></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>n</mi></mrow></msub><mo>+</mo><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub></mrow></mfrac><mo>&CenterDot;</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>t</mi><mrow><mi>e</mi><mi>p</mi></mrow></msub></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000031.GIF" wi="1158" he="271" /></maths>结合针孔成像坐标系,根据图像坐标系O<sub>p</sub>‑X<sub>p</sub>Y<sub>p</sub>相机坐标系O<sub>c</sub>‑X<sub>c</sub>Y<sub>c</sub>Z<sub>c</sub>、数字图像坐标系的转换关系,计算目标方位角α、β;相机坐标系O<sub>c</sub>‑X<sub>c</sub>Y<sub>c</sub>Z<sub>c</sub>中,O<sub>c</sub>为相机光心,Z<sub>c</sub>为光轴,X<sub>c</sub>、Y<sub>c</sub>与成像平面的X<sub>p</sub>、Y<sub>p</sub>轴平行;由针孔成像模型各坐标系间的比例关系:<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><msub><mi>x</mi><mi>p</mi></msub><mo>=</mo><mfrac><msub><mi>x</mi><mi>c</mi></msub><msub><mi>z</mi><mi>c</mi></msub></mfrac><mo>&CenterDot;</mo><mi>f</mi></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>y</mi><mi>p</mi></msub><mo>=</mo><mfrac><msub><mi>y</mi><mi>c</mi></msub><msub><mi>z</mi><mi>c</mi></msub></mfrac><mo>&CenterDot;</mo><mi>f</mi></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000032.GIF" wi="1054" he="262" /></maths>记每个像素在X<sub>p</sub>与Y<sub>p</sub>方向上物理尺寸为dX、dY,CCD相机面阵为m×m像素,取图像坐标系原点在数字图像坐标系(u,v)中的坐标(u<sub>0</sub>,v<sub>0</sub>)为:<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><msub><mi>u</mi><mn>0</mn></msub><mo>=</mo><mfrac><mi>m</mi><mn>2</mn></mfrac><mo>+</mo><mn>0.5</mn></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>v</mi><mn>0</mn></msub><mo>=</mo><mfrac><mi>m</mi><mn>2</mn></mfrac><mo>+</mo><mn>0.5</mn></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000033.GIF" wi="1070" he="239" /></maths>式(9)中,m一般为偶数,因此需加0.5后,图像坐标系原点(u<sub>0</sub>,v<sub>0</sub>)位于图像正中;根据图像坐标系O<sub>p</sub>‑X<sub>p</sub>Y<sub>p</sub>与数字图像坐标系o‑uv间的关系,得图像坐标系中的点(x<sub>p</sub>,y<sub>p</sub>)对应的(u,v)坐标为:<maths num="0006" id="cmaths0006"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>u</mi><mo>=</mo><msub><mi>u</mi><mn>0</mn></msub><mo>-</mo><mfrac><msub><mi>x</mi><mi>p</mi></msub><mrow><mi>d</mi><mi>X</mi></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mi>v</mi><mo>=</mo><msub><mi>v</mi><mn>0</mn></msub><mo>-</mo><mfrac><msub><mi>y</mi><mi>p</mi></msub><mrow><mi>d</mi><mi>Y</mi></mrow></mfrac></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000034.GIF" wi="1069" he="255" /></maths>偏航角α、俯仰角β为:<maths num="0007" id="cmaths0007"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>&alpha;</mi><mo>=</mo><mi>arctan</mi><mrow><mo>(</mo><mfrac><msub><mi>x</mi><mi>c</mi></msub><msub><mi>z</mi><mi>c</mi></msub></mfrac><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>&beta;</mi><mo>=</mo><mi>arctan</mi><mrow><mo>(</mo><mfrac><msub><mi>y</mi><mi>c</mi></msub><msqrt><mrow><msubsup><mi>x</mi><mi>c</mi><mn>2</mn></msubsup><mo>+</mo><msubsup><mi>z</mi><mi>c</mi><mn>2</mn></msubsup></mrow></msqrt></mfrac><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000035.GIF" wi="1149" he="278" /></maths>由针孔成像模型,相机坐标系与图像坐标系的转换关系为:<maths num="0008" id="cmaths0008"><math><![CDATA[<mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>x</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mi>p</mi></msub></mtd></mtr><mtr><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>=</mo><mfrac><mn>1</mn><msub><mi>z</mi><mi>c</mi></msub></mfrac><mfenced open = "[" close = "]"><mtable><mtr><mtd><mi>f</mi></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mi>f</mi></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd></mtr></mtable></mfenced><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>x</mi><mi>c</mi></msub></mtd></mtr><mtr><mtd><msub><mi>y</mi><mi>c</mi></msub></mtd></mtr><mtr><mtd><msub><mi>z</mi><mi>c</mi></msub></mtd></mtr><mtr><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000041.GIF" wi="1221" he="271" /></maths>结合式(10)、(11)、(12),得:<maths num="0009" id="cmaths0009"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>&alpha;</mi><mo>=</mo><mi>arctan</mi><mrow><mo>(</mo><mfrac><mrow><mrow><mo>(</mo><mrow><mi>u</mi><mo>-</mo><msub><mi>u</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>&CenterDot;</mo><mi>d</mi><mi>X</mi></mrow><mi>f</mi></mfrac><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>&beta;</mi><mo>=</mo><mi>arctan</mi><mrow><mo>(</mo><mfrac><mrow><mrow><mo>(</mo><mrow><mi>v</mi><mo>-</mo><msub><mi>v</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>&CenterDot;</mo><mi>d</mi><mi>Y</mi></mrow><msqrt><mrow><msup><mrow><mo>(</mo><mrow><mrow><mo>(</mo><mrow><mi>u</mi><mo>-</mo><msub><mi>u</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow><mo>&CenterDot;</mo><mi>d</mi><mi>X</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><mi>f</mi><mn>2</mn></msup></mrow></msqrt></mfrac><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000973350690000042.GIF" wi="1286" he="357" /></maths>解算出目标方位角α、β后,即完成碎片方位角定位。
地址 100081 北京市海淀区中关村南大街5号