業界別最適化Example集: 3D点群データの平面フィッティングによる地形評価 (建設)
nAG数値計算ライブラリ > nAGライブラリ数理最適化 > 業界別最適化Example集 > 3D点群データの平面フィッティングによる地形評価 (建設)

Keyword: 建設, 最適化, 地形評価, 平面フィッティング, NLP

建設分野の最適化問題:3D点群データの平面フィッティングによる地形評価問題

問題の概要

建設予定地

建設現場では、土地の平坦度を評価することが重要です。ドローンや測量機器を用いて取得した3D点群データを解析し、地面の平面性を評価するために、点群データから最適な平面をフィッティングする必要があります。この問題の解決により、現場の土地が建設に適しているかどうかを迅速かつ正確に判断できます。

地形評価問題(具体例)

具体的な例として、建設予定地の一部で取得された15個の3D座標点があります。これらの点から地形の平面性を評価するため、各点から平面までの距離の二乗和を最小化する平面を見つける問題を考えます。

この問題を解くことで、以下の情報が得られます:

  1. 地形の傾斜角度と方向

  2. 地形の平坦度(最大偏差、平均偏差、標準偏差)

  3. 最適な平面の方程式

座標データ

点番号 x座標 (m) y座標 (m) z座標 (m)
1 1.2 2.5 3.8
2 2.3 3.1 4.1
3 3.8 1.9 4.5
4 3.2 2.7 4.0
5 4.1 2.3 4.3
6 3.7 3.5 4.4
7 2.1 3.8 4.2
8 3.5 2.6 4.6
9 2.9 3.0 4.0
10 3.0 2.8 4.5
11 3.1 3.4 4.7
12 1.8 2.9 4.1
13 2.5 3.2 4.3
14 3.4 2.1 4.2
15 4.0 3.3 4.6
データ

問題の定式化

パラメータ

パラメータ 説明
xi 各点のx座標 (m) 表参照
yi 各点のy座標 (m) 表参照
zi 各点のz座標 (m) 表参照
n 点の数 15

xiyizi は実数値で、座標データ表に示されています。

決定変数

変数 説明 範囲
a 平面の法線ベクトルのx成分 <a<
b 平面の法線ベクトルのy成分 <b<
c 平面の法線ベクトルのz成分 <c<
d 平面の定数項 <d<

これらの変数は実数値で、平面の方程式 ax+by+cz+d=0 を定義します。

目的関数

目的
最小化 i=1n(|axi+byi+czi+d|a2+b2+c2)2
この目的関数は、各点から平面までの距離の二乗和を表しています。分子の絶対値 ax_i + by_i + cz_i + d

制約条件

制約 説明
法線ベクトルの制約 a2+b2+c2=1 平面の法線ベクトルが単位ベクトルであることを保証

この制約条件は、法線ベクトルの長さを1に正規化し、解の一意性を確保します。

コード例

以下に、この最適化問題を解くPythonコード例を示します。

import numpy as np
from naginterfaces.library import opt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 座標データ
coordinates = np.array([
    [1.2, 2.5, 3.8], [2.3, 3.1, 4.1], [3.8, 1.9, 4.5],
    [3.2, 2.7, 4.0], [4.1, 2.3, 4.3], [3.7, 3.5, 4.4],
    [2.1, 3.8, 4.2], [3.5, 2.6, 4.6], [2.9, 3.0, 4.0],
    [3.0, 2.8, 4.5], [3.1, 3.4, 4.7], [1.8, 2.9, 4.1],
    [2.5, 3.2, 4.3], [3.4, 2.1, 4.2], [4.0, 3.3, 4.6]
])
n = len(coordinates)
epsilon = 1e-8  # 数値的安定性のための微小値

def objfun(mode, x, grad, _nstate):
    """ 目的関数:各点から平面までの距離の二乗和 """
    a, b, c, d = x
    norm = np.sqrt(a**2 + b**2 + c**2) + epsilon
    
    if mode in [0, 2]:
        plane_eq = np.dot(coordinates, [a, b, c]) + d
        distances = np.abs(plane_eq) / norm
        obj_val = np.sum(distances**2)
    else:
        obj_val = 0.0
    
    if mode == 2:
        sign = np.sign(plane_eq)
        common_term = 2 * sign * distances / norm**2
        grad_abc = np.sum(common_term[:, np.newaxis] * (coordinates * norm - [a, b, c] * plane_eq[:, np.newaxis] / norm), axis=0)
        grad_d = np.sum(2 * sign * distances / norm)
        grad[:] = np.concatenate([grad_abc, [grad_d]])
    
    return obj_val, grad

def confun(mode, needc, x, cjac, _nstate):
    """ 非線形制約:法線ベクトルの長さを1に制約 """
    a, b, c, _ = x
    
    if mode in [0, 2]:
        cons = np.array([a**2 + b**2 + c**2 - 1])
    else:
        cons = np.array([])
    
    if mode == 2:
        cjac[0, :] = [2*a, 2*b, 2*c, 0]
    
    return cons, cjac

# 初期値と境界条件の設定
x_init = [0, 0, 1, 0]
nvar, ncnln = 4, 1
lower_bounds = np.concatenate([-np.inf * np.ones(nvar), [0]])
upper_bounds = np.concatenate([np.inf * np.ones(nvar), [0]])

# ソルバーオプションの設定
comm = opt.nlp1_init('nlp1_solve')
opt.nlp1_option_string('Print Level = 0', comm)
opt.nlp1_option_string('Major Iteration Limit = 100', comm)

# ソルバーの呼び出し
result = opt.nlp1_solve(
    np.zeros((0, nvar)), lower_bounds, upper_bounds, objfun, x_init, comm,
    confun=confun, cjac=np.zeros((ncnln, nvar))
)

# 最適なパラメータと地形メトリクスの計算
a, b, c, d = result.x
normal = np.array([a, b, c])
normal /= np.linalg.norm(normal)
distances = np.abs(np.dot(coordinates, normal) + d)
slope_angle = np.arccos(abs(c)) * 180 / np.pi
slope_direction = np.array([a, b])
slope_direction /= np.linalg.norm(slope_direction)

# 目的関数の最適値を計算
obj_val, _ = objfun(0, result.x, None, None)

# 結果の出力
print("=== 建設現場地形評価 ===")
print(f"平面の方程式パラメータ: a={a:.6f}, b={b:.6f}, c={c:.6f}, d={d:.6f}")
print(f"目的関数の最適値: {obj_val:.6f}")
print(f"最大偏差: {np.max(distances):.6f} m")
print(f"平均偏差: {np.mean(distances):.6f} m")
print(f"標準偏差: {np.std(distances):.6f} m")
print(f"傾斜角度: {slope_angle:.6f}°")
print(f"傾斜の向き (X, Y): ({slope_direction[0]:.6f}, {slope_direction[1]:.6f})")
print("  注: 傾斜の向きは、XY平面上で最も急な下り勾配の方向を示します。")
print("      (1, 0)はX軸正方向、(0, 1)はY軸正方向を表します。")

# 3Dプロット
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')

scatter = ax.scatter(coordinates[:, 0], coordinates[:, 1], coordinates[:, 2], c=coordinates[:, 2], cmap='viridis', s=50)
plt.colorbar(scatter, label='Z value (m)')

x_range = np.linspace(coordinates[:, 0].min(), coordinates[:, 0].max(), 10)
y_range = np.linspace(coordinates[:, 1].min(), coordinates[:, 1].max(), 10)
X, Y = np.meshgrid(x_range, y_range)
Z = (-a * X - b * Y - d) / c
ax.plot_surface(X, Y, Z, alpha=0.5, cmap='coolwarm')

midpoint = coordinates[:, :2].mean(axis=0)
ax.quiver(midpoint[0], midpoint[1], Z.mean(), slope_direction[0], slope_direction[1], 0, 
          length=0.5, normalize=True, color='red', arrow_length_ratio=0.3)

ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.set_title('Terrain and Optimal Plane')

plt.tight_layout()
plt.savefig('terrain_evaluation.png', dpi=300, bbox_inches='tight')
print("\n3Dプロットを 'terrain_evaluation.png' として保存しました。")
print("赤い矢印は傾斜の向きを示しています。")

結果例

上記のコードを実行すると、以下のような結果が得られます。

=== 建設現場地形評価 ===
平面の方程式パラメータ: a=-0.227024, b=-0.161285, c=0.960441, d=-2.978646
目的関数の最適値: 0.440030
最大偏差: 0.298829 m
平均偏差: 0.139197 m
標準偏差: 0.099797 m
傾斜角度: 16.169698°
傾斜の向き (X, Y): (-0.815217, -0.579156)
  注: 傾斜の向きは、XY平面上で最も急な下り勾配の方向を示します。
      (1, 0)はX軸正方向、(0, 1)はY軸正方向を表します。

3Dプロットを 'terrain_evaluation.png' として保存しました。
赤い矢印は傾斜の向きを示しています。
3Dプロット

この結果から、以下のことが分かります:

  1. 地形の傾斜角度は約16.17°であり、これは中程度の傾斜を示しています。

  2. 最大偏差は約0.2988 m であり、これは最も平面から離れた点がおよそ30 cm の高低差があることを示しています。この値が許容範囲内かどうかは、プロジェクトの要件によって判断する必要があります。

  3. 平均偏差は約0.1392 m で、標準偏差は約0.0998 m です。これらの値は、地形の全体的な起伏の程度を示しています。

  4. 傾斜の向きは主に北西方向(X軸負方向とY軸負方向)を指しています。

まとめ

この例題を通じて、3D点群データから平面をフィッティングし、地形の平面性を評価する手法を示しました。この手法は建設現場において土地の評価に役立つだけでなく、他の様々な分野でも応用可能です。 例えば、製造業での部品の平坦性検査や、航空機の表面検査など、3Dデータを扱う多くの場面で利用できます。


関連情報
© 隴鯉ス・隴幢スャ郢昜ケ斟礼ケ晢スシ郢晢ス。郢晢スェ郢ァ�ォ郢晢スォ郢ァ�「郢晢スォ郢ァ�エ郢晢スェ郢ァ�コ郢晢ソス郢ァ�コ郢ァ�ー郢晢スォ郢晢スシ郢晉軸�ス�ェ陟台ク茨スシ螟ゑス、�セ 2024
Privacy Policy  /  Trademarks