林木直径是林业生产、科研重要的测树因子之一，准确、高效地对林木直径进行测算，也是林业工作者的重要研究内容。从目前获取林木直径的途径看，可分为直接测量和间接测算两大类，前者简单、适用性强，但需要大量的财力支持，同时劳动强度大，效率低，且获得的信息也比较单一。除直接测量以外的所有其他方法叫间接测量，可以说，到目前为止还没有一种间接办法能取代直接测量，因此，寻找一种简单高效的树木测算方法一直是林业工作者的梦想。

计算机视觉(图像理解)是人工智能的一个分支，它是在图像处理技术上发展起来的，却与图像处理有本质的不同，直接目的是从图像中获取人们感兴趣的信息。这一高新技术在诸多领域(王雪峰等，2000)已经得到很好的应用，如果将这一技术成功的应用于林业，不仅能获得更为丰富的信息，而且还能大大减轻劳动强度，节省大量资金，这必将给森林经理调查带来质的飞跃。本文按着在林业中可能应用的步骤，给出从图像抽取3D点信息的具体算法，并用实例对其验证。

1 算法概述 1.1 计算摄像机内参数

已知空间点及对应图像点u的齐次坐标，则针孔模型(章毓晋，2000)可为：

(1)

其中：是空间点X的齐次表示，P是3×4的投影矩阵。令，将上式两边同乘h后消除Z_c，此时(1)变为下面的简单形式：

(2)

实际上，P是由如下几个矩阵合成的：

(3)

ρ是上三角阵，代表摄像机内参数，主要包括一个像素在横、纵向的物理长度(d_x, d_y)，像主点坐标(u₀, v₀)，扭曲系数等；R是3×3的旋转矩阵，在R的9个参数中只有3个独立参数，即空间旋转的3个欧拉角(α_x, α_y, α_z); S^T=[s_x  s_y  s_z]是平移向量。R、S的物理意义是将世界坐标系原点带回到摄像机坐标系原点所需要的外参数。把(3)代入(2)式并去掉齐次坐标表示有：

(4)

这是针孔模型的基本方程式。下面对(2)或(4)式中参数的求解过程，叫摄像机定标。

对(2)进行矩阵拉直运算，有：

(5)

⊗是Kronecker积。如果已知n个样本点(X_t, h_t)，按(5)构造系数矩阵A：

(6)

此时有：

(7)

将(6)式中的A分解成分解成，则(7)可改写为：

(8)

由(7)式知，等式两边同乘一不为0常数，方程仍然成立，所以我们可以设p=1，即，这样(7)的最小二乘解为：

(9)

是的估计值。

在三维重建过程中，需要两台或多台摄像机(或一台在两个或多个位置摄像)，如果各摄像机相对位置保持不变，则可由投影矩阵直接重建出空间信息，而不需要了解投影矩阵中各参数的具体含义。很明显，这对于林业野外操作是不适用的，就是说，需要分解投影矩阵，只有不限制摄像机外参数才更有实际意义。理论上，得到投影矩阵后，就可以根据正交矩阵性质从中分解出摄像机的内、外参数。需要注意的是，尽管在求解P时令p=1，但分解时必须考虑这一常数。这种从投影矩阵分解参数计算方法的最大优点是速度快，但由于实际数据误差的存在，但难保证分解出的旋转矩阵R正交，从而影响最终的重建结果。因此，仅把这种线性结果作为初始参数，然后由度量误差模型(唐守正等，1996)直接计算(4)式中参数。这样做由于考虑了所有约束条件且在求解时直接计算各独立参数，从而消除了矩阵分解所代来的误差，提高了参数估计精度。

1.2 重建前的摄像机外参数求解

把外参数与内参数求解分离且放在重建前来做，完全是出于对林业实际问题的考虑。我们知道，如果与是一对匹配点，则存在外极线方程(马颂德等，1998)：

(10)

ρ₁、ρ₂是摄像机一、二的内参数矩阵，R₂₁、S₂₁是第二个摄像机相对于第一个摄像机的旋转及平移，且‖S₂₁‖=L_b(基线长)。[S₂₁]_x是向量S₂₁的反对称矩阵。外极线方程的最原始作用是用于对匹配点的搜索，如果已知匹配点对，我们将其发展成为对参数的求解。

当已知n个匹配点对由(10)式就可解出外参数。容易看出，(10)式有5个独立外参数，它们是3个旋转参数(α₂₁)，2个平移参数(S₂₁中的2个)，就是说，至少要有5对匹配点才能求出外参数。实际应用，所选匹配点数要远大于5。

推广到k(k≥2)个匹配点对应同一空间的情况，此时形成C_k²个外极线方程组，可求出全部外参数。

1.3 三维重建

获得全部摄像机参数后，由匹配点重建其对应空间点仅仅是一个简单计算问题，本文不再给出计算公式。需要指出的是，这种计算得到的空间点，通常是使世界坐标系与某一个摄像机坐标系重合即某一摄像机坐标系中的点。这对于解决我们的问题已经足够了，因为我们关心的是树木间的相对位置。文中给出的三维点是把世界坐标系重合到第一个摄像机坐标系的结果。

2 树木直径抽取实例

为验证上述算法在林业应用中的可行性，我们对办公楼前的小树林进行摄影，得到以下2张图像，同时量测基线长度为4 299 mm。现在的任务是从这2张图像中重建树木胸径。需要说明的是，图 1a中编号是为本文的便于说明加进去的，摄影时并不存在。

2.1 参数计算

内参数获取是在室内完成的。首先由(2)计算P并分解出各参数作为初始值，然后由度量误差模型精确计算。在由度量误差模型计算过程中，需要注意3点：1)由于图像点对整个系统的影响要大于空间点，因此，假设图像点是有误差的；2)在计算中要不断求解F(X, u, c)=0中的u，此图像点坐标可由(4)直接计算；3)误差矩阵赋值为2×2单位阵，把计算结果中的内参数列于表 1。表 1中的外参数是从图 1a、b中由人机交互方式抽取22个明显的匹配特征点，利用度量误差模型解(15)式得到的，此时假设匹配点是有误差的，误差阵赋4×4单位阵，方程F(X, u, c)=0中X为空，u是匹配点。另一值得注意的是焦距因子，最为简单的处理办法是保持定标与重建时二者的统一，这种作法势必对实际摄像造成了一定限制，所以本文没有采用这种办法，而是把焦距单独考虑，仅在定标时给出d_x、d_y，在应用中取实际使用焦距，增加了灵活性。实际上，目前的大多数CCD像机在摄影同时，也记录了所使用的焦距及其他因子，可以直接从文件中读出。

2.2 三维重建

图像特征点匹配是重建中最为关键的问题之一。直接用野外获取的树木图像匹配是困难的，这一方面是由于问题本身存在难度，另一方面是由于处于自然背景中的树木图像复杂且含有大量的噪声，因此，对原始图像预处理是必要的。本试验中，分别用一阶梯度算子、二阶微分算子以及Susan算子进行树干边缘检测实验。结果表明，对于树干阶跃型的边缘，Canny算子不但能很好地检测出树干的边缘强度，而且能够检测出树干边缘的方向。对于不同背景下的树干，可以通过改变Canny算子的标准差σ，来获得最佳的树干的边缘。图 2是取σ=1时的Canny检测后图像。

利用经过去噪、Canny检测的树木图像，在左图像中用鼠标给出胸径位置，由程序自动找出2个边缘点，然后计算两个边缘点的外极线，最后在外极线上下各15个象素点的窄条带内，搜索灰度最相关的点，这个点即为左图像点的匹配点，整理成表 2。这种在右图像中外极线经过的窄条带内寻找匹配点的作法，不仅大大减少了搜索范围，更重要的是降低了误匹配的可能性。

至此我们具备了重构空间点的全部信息，由(19)容易计算胸径部位左、右点在摄像机坐标系中的三维坐标，进一步计算两点间欧式距离得到胸径见表 3。

由表 3可知，10号树相对误差最大为9.9%，这主要是由于匹配误差造成的，从图 1a也容易看出，10号树胸径部位边缘不很清晰。除此之外，其他相对误差均在5%以下。最小相对误差为0.3%，平均相对误差是2.7%，结果还是令人满意的。

3 结论及讨论

本文给出了如何从图像出发抽取树木直径的算法、步骤，同时也给出了一个实例。结果表明：1)从图像出发抽取树木胸径是可行的。本文是以胸径为例，实际上，如果我们能得到树木上任意直径部位的匹配图像，则得到该处直径与胸径没有任何区别；2)本文在参数求解时均采用了统计学中较新的研究成果——度量误差模型算法，并得到了很好的结果；3)本文采用了内、外参数分离、不固定焦距的做法，使实际操作更加容易，实例结果也表明了这种策略的可行性。

本文在确定匹配点时采用了人机交互方式，如何完全由机器来做或尽可能减少人的参与，寻找更好的算法是一个值得研究的问题；另外，如何处理遮挡也是需要解决的问题。