储量估算完整管线

从钻孔数据库 .dmd 和矿体约束 .dmf 出发，依次完成格式转换、特高品位处理、样长组合、块段模型创建和距离幂估值，最终输出估值后的 .dmb 块段模型。

前置条件

uv add pandas numpy

• 环境变量 DIMINE_HOME 指向包含 DmPyBindInterface.pyd 的目录
• Python >= 3.12
• 输入文件：钻孔数据库 .dmd + 矿体约束 .dmf

流程概览

钻孔数据库(.dmd)          矿体模型(.dmf)
      │                        │
      ▼                        │
① DMD→DMG 格式转换             │
      │                        │
      ▼                        │
② 特高品位处理（迭代收敛）      │
      │                        │
      ▼                        │
③ 样长组合                     │
      │                        │
      ▼                        ▼
④ 创建空块段模型  ←── 矿体包围盒
      │
      ▼
⑤ 搜索椭球体参数 ←── MVEE 算法
      │
      ▼
⑥ 距离幂估值（IDW，迭代扩张）

---

步骤一：DMD → DMG 格式转换 + 品位字段检测

将钻孔数据库从 .dmd 转换为 .dmg 分段文件格式，并自动检测有效的品位字段。

from dimine_python_sdk.lib.io import DmdFile, DmgFile, convert_file

dmd_path = "01_地质勘探数据库.dmd"
dmg_raw_path = "01_raw.dmg"

# 1.1 格式转换
convert_file(dmd_path, dmg_raw_path)
print(f"已转换: {dmd_path} → {dmg_raw_path}")

# 1.2 加载 DMG 分段数据，自动检测品位字段
dmg = DmgFile(dmg_raw_path)
segments_df = dmg.to_segments_dataframe()
dmg.close()

print(f"DMG 分段数据: {len(segments_df)} 行")
print(f"列名: {list(segments_df.columns)}")

# 1.3 从数值列中排除非品位字段，自动筛选有效品位
numeric_cols = segments_df.select_dtypes(include="number").columns.tolist()
excluded = {
    "起始", "结束", "样长", "X", "Y", "Z", "覆盖长度", "覆盖率",
    "测斜深度", "方位角", "倾角", "样本编号", "钻孔编号",
    "start_x", "start_y", "start_z", "end_x", "end_y", "end_z",
    "SAMPLE-ID", "SAMFROM", "SAMTO",
}
candidate_fields = [c for c in numeric_cols if c not in excluded]

# 1.4 过滤无有效数据的字段（全为 ≤0、NaN、Inf）
valid_grade_fields = []
for field in candidate_fields:
    col = segments_df[field]
    is_invalid = (
        col.isna() | col.isin([float("inf"), float("-inf")]) | (col <= 0)
    ).all()
    if not is_invalid:
        valid_grade_fields.append(field)

print(f"检测到的有效品位字段: {valid_grade_fields}")

# 1.5 读取钻孔元信息（仅用于显示）
db = DmdFile(dmd_path)
print(f"钻孔数: {len(db.collar)}")
print(f"样品数: {len(db.sample)}")

---

步骤二：特高品位处理（迭代收敛）

特高品位处理需要迭代执行：一轮处理会降低矿体平均品位，原本不超限的样品在新平均值下可能成为新的特高品位。需要反复迭代直到再无样品超过阈值。

import shutil
from dimine_python_sdk.lib.io import DmgFile
from dimine_python_sdk.lib.prospecting import DrillFunctionWrapper
from dimine_python_sdk.lib.prospecting.models import HighGradeProcessParam

# 2.1 拷贝原始文件（特高品位处理原地修改）
dmg_hg_path = "02_high_grade.dmg"
shutil.copy2(dmg_raw_path, dmg_hg_path)

# 2.2 参数配置
HG_AVERAGE_MULTIPLE = 8   # 平均值倍数（铁矿石=8，金/有色金属=6，煤矿/非金属=7）
HG_PROCESS_MODE = 0        # 0=国内规范（均值倍数），1=国外规范（累积频率）
HG_REPLACE_METHOD = 3      # 0=剔除 1=给定值 2=相邻样品均值 3=矿体平均品位 4=单一工程 5=截止品位
MAX_HG_ITERATIONS = 5      # 安全上限

REPLACE_NAMES = {
    0: "剔除法", 1: "给定值法", 2: "相邻样品平均值法",
    3: "矿体平均品位法", 4: "单一工程法", 5: "截止品位法",
}

for grade_field in valid_grade_fields:
    print(f"\n处理元素: {grade_field}")
    print(f"  模式: {'国内规范' if HG_PROCESS_MODE == 0 else '国外规范'}")
    print(f"  倍数: {HG_AVERAGE_MULTIPLE}x, 替换方法: {REPLACE_NAMES[HG_REPLACE_METHOD]}")

    for iteration in range(1, MAX_HG_ITERATIONS + 1):
        # 读取当前 DMG 品位数据
        dmg = DmgFile(dmg_hg_path)
        segments_df = dmg.to_segments_dataframe()
        dmg.close()

        if grade_field not in segments_df.columns:
            print(f"  [SKIP] 字段不存在")
            break

        # 过滤有效品位（排除 ≤0、NaN、Inf）
        grades = segments_df[grade_field].dropna()
        grades = grades[~grades.isin([float("inf"), float("-inf")])]
        grades = grades[grades > 0]

        if len(grades) == 0:
            print(f"  [SKIP] 无有效品位数据")
            break

        # 计算阈值
        avg_grade = float(grades.mean())
        threshold = avg_grade * HG_AVERAGE_MULTIPLE
        exceeding_count = int((grades > threshold).sum())

        print(f"  迭代 #{iteration}: 样品数={len(grades)}, "
              f"均值={avg_grade:.4f}, 阈值={threshold:.4f}, "
              f"超限={exceeding_count}, 最大值={float(grades.max()):.4f}")

        if exceeding_count == 0:
            print(f"  [收敛] 已无超限样品")
            break

        # 调用 C++ 特高品位处理
        hg_params = HighGradeProcessParam(
            input_file=dmg_hg_path,
            grade_field=grade_field,
            process_mode=HG_PROCESS_MODE,
            average_multiple=HG_AVERAGE_MULTIPLE,
            frequency=0.95,
            replace_method=HG_REPLACE_METHOD,
            assign_value=2.6,
            adjoin_average=3.2,
            average_method=0,
            contain_mode=0,
        )
        result = DrillFunctionWrapper.extra_high_grade_process(hg_params)
        print(f"  C++ 处理: {result.get('message', result)}")
    else:
        print(f"  [警告] 达到最大迭代次数 {MAX_HG_ITERATIONS}，可能未完全收敛")

replace_method 枚举：

值	方式	适用场景
0	剔除法	样品充足时
1	给定值法	有经验阈值时
2	相邻样品平均值法	通用推荐
3	矿体平均品位法	矿体品位稳定时
4	单一工程法	钻孔独立性强时
5	截止品位法	有工业指标时

---

步骤三：样长组合

将不等长样品合并为指定长度的等长组合样，覆盖率低于阈值的残段将被丢弃。

from dimine_python_sdk.lib.prospecting import DrillFunctionWrapper
from dimine_python_sdk.lib.prospecting.models import SampleLengthCombineParam

dmg_combined_path = "03_combined.dmg"

combine_params = SampleLengthCombineParam(
    input_file=dmg_hg_path,
    combine_length=2.0,       # 组合样长度（m）
    combine_percent=0.75,      # 最小覆盖率
    output_file=dmg_combined_path,
)

print(f"组合长度: {combine_params.combine_length}m")
print(f"最小覆盖率: {combine_params.combine_percent}")

result = DrillFunctionWrapper.sample_length_combine(combine_params)
print(f"组合结果: {result}")

---

步骤四：从矿体包围盒创建空块段模型

从矿体 .dmf 模型中提取所有几何顶点，计算包围盒，加上边距后创建空块段模型。

import math
import numpy as np
from dimine_python_sdk.lib.io import DmfFile, DmbFile
from dimine_python_sdk.lib.prospecting import create_block_model
from dimine_python_sdk.lib.prospecting.models import (
    CreateBlockModelParams,
    FieldDefinition,
)
from dimine_python_sdk.models.types import Shell, Point as DmfPoint, Line as DmfLine

dmf_constraint_path = "约束_矿体.dmf"
dmb_output_path = "04_block_model.dmb"

# 4.1 从 DMF 计算矿体包围盒
dmf = DmfFile(dmf_constraint_path)
all_points = []

for layer in dmf.layers:
    for entity in layer.entities:
        if isinstance(entity, Shell):
            pts = entity.geometry.points
            if pts is not None and len(pts) > 0:
                all_points.append(pts)
        elif isinstance(entity, DmfLine):
            geom = entity.geometry
            if geom is not None and len(geom) > 0:
                all_points.append(geom)
        elif isinstance(entity, DmfPoint):
            geom = entity.geometry
            if geom is not None and len(geom) > 0:
                all_points.append(geom.reshape(1, 3))

pts = np.vstack(all_points)
dmf.close()

min_x, max_x = float(pts[:, 0].min()), float(pts[:, 0].max())
min_y, max_y = float(pts[:, 1].min()), float(pts[:, 1].max())
min_z, max_z = float(pts[:, 2].min()), float(pts[:, 2].max())

print(f"矿体包围盒:")
print(f"  X: {min_x:.2f} ~ {max_x:.2f}")
print(f"  Y: {min_y:.2f} ~ {max_y:.2f}")
print(f"  Z: {min_z:.2f} ~ {max_z:.2f}")

# 4.2 计算块段模型参数（加边距并取整）
BLOCK_SIZE_X, BLOCK_SIZE_Y, BLOCK_SIZE_Z = 5.0, 5.0, 5.0
MARGIN_XY, MARGIN_Z = 20.0, 10.0

origin_x = min_x - MARGIN_XY
origin_y = min_y - MARGIN_XY
origin_z = min_z - MARGIN_Z

length_x_raw = max_x - min_x + 2 * MARGIN_XY
length_y_raw = max_y - min_y + 2 * MARGIN_XY
length_z_raw = max_z - min_z + 2 * MARGIN_Z

nx = math.ceil(length_x_raw / BLOCK_SIZE_X)
ny = math.ceil(length_y_raw / BLOCK_SIZE_Y)
nz = math.ceil(length_z_raw / BLOCK_SIZE_Z)

length_x = nx * BLOCK_SIZE_X
length_y = ny * BLOCK_SIZE_Y
length_z = nz * BLOCK_SIZE_Z

print(f"\n块段模型参数:")
print(f"  原点: ({origin_x:.2f}, {origin_y:.2f}, {origin_z:.2f})")
print(f"  单元块: {BLOCK_SIZE_X}m × {BLOCK_SIZE_Y}m × {BLOCK_SIZE_Z}m")
print(f"  范围: {length_x:.1f}m × {length_y:.1f}m × {length_z:.1f}m")
print(f"  块数: {nx} × {ny} × {nz} = {nx * ny * nz:,}")

# 4.3 创建空块段模型
block_params = CreateBlockModelParams(
    file_name=dmb_output_path,
    origin_x=origin_x,
    origin_y=origin_y,
    origin_z=origin_z,
    size_x=BLOCK_SIZE_X,
    size_y=BLOCK_SIZE_Y,
    size_z=BLOCK_SIZE_Z,
    length_x=length_x,
    length_y=length_y,
    length_z=length_z,
    field_definition=[
        FieldDefinition(field_name="岩石类型", field_type="字节型"),
        FieldDefinition(field_name="体重", field_type="浮点型"),
    ],
)

create_result = create_block_model(block_params)
print(f"创建结果: {create_result}")

# 4.4 验证块段模型属性
with DmbFile(dmb_output_path) as block_model:
    print(f"\n块段模型属性:")
    print(f"  原点: {block_model.origin}")
    print(f"  尺寸: X={block_model.xyz_length[0]:.1f}m, "
          f"Y={block_model.xyz_length[1]:.1f}m, "
          f"Z={block_model.xyz_length[2]:.1f}m")
    print(f"  最大精度层级: {block_model.max_level}")
    print(f"  最小块段: {block_model.min_size}")
    print(f"  块段总数: {block_model.total_block_count:,}")
    print(f"  字段定义: {block_model.field_definitions}")

---

步骤五：距离幂估值（IDW，迭代扩张）

这是最关键的步骤，分四个阶段：

1. 5a — 从矿体模型用 MVEE 算法计算搜索椭球体角度和矿体几何尺寸
2. 5b — 从钻孔孔口数据计算勘探线间距
3. 5c — 基于勘探线间距 + 矿体几何比例推导搜索椭球体三轴半径
4. 5d — 迭代执行 IDW 估值，未估值块段 > 0 时主轴半径翻倍重试

5a. MVEE 搜索椭球体参数

from dimine_python_sdk.lib.prospecting import compute_search_ellipsoid_from_dmf

print("[5a] 计算搜索椭球体参数（MVEE 算法）...")

ellipsoid = compute_search_ellipsoid_from_dmf(
    dmf_constraint_path, voxel_size=5.0
)

print(f"  MVEE 椭球体（矿体几何）:")
print(f"    主轴半径（半走向长）: {ellipsoid.ellipsoid_main_radius:.1f}m")
print(f"    次轴半径（半延伸长）: {ellipsoid.second_radius:.1f}m")
print(f"    短轴半径（半厚度）:   {ellipsoid.short_radius:.1f}m")
print(f"    方位角: {ellipsoid.ellipsoid_azimuth:.1f}°")
print(f"    倾伏角: {ellipsoid.ellipsoid_plunge:.1f}°")
print(f"    倾角:   {ellipsoid.ellipsoid_dip:.1f}°")
print(f"    次轴/主轴: {ellipsoid.ellipsoid_second_main_rate:.4f}")
print(f"    短轴/主轴: {ellipsoid.ellipsoid_short_main_rate:.4f}")
print(f"    椭球中心: {ellipsoid.center}")
print(f"    主轴方向: {ellipsoid.pc1_direction}")

5b. 勘探线间距计算

import math
import numpy as np

def find_xy_columns(collar_df):
    """从孔口表识别 X/Y 坐标列名"""
    x_candidates = ["横坐标", "EAST", "X", "Easting", "easting", "x"]
    y_candidates = ["纵坐标", "NORTH", "Y", "Northing", "northing", "y"]
    x_col = next((c for c in x_candidates if c in collar_df.columns), None)
    y_col = next((c for c in y_candidates if c in collar_df.columns), None)
    if x_col is None or y_col is None:
        numeric_cols = collar_df.select_dtypes(include="number").columns.tolist()
        if len(numeric_cols) >= 2:
            x_col, y_col = numeric_cols[0], numeric_cols[1]
        else:
            raise ValueError("无法识别X/Y坐标列")
    return x_col, y_col


def compute_exploration_line_spacing(collar_df, strike_direction,
                                      line_field=None, same_line_threshold=5.0):
    """计算勘探线间距

    优先按 line_field 分组计算；否则将所有孔口投影到走向方向取相邻间距中位数。
    """
    if len(collar_df) < 2:
        return 50.0

    # 方法1：按勘探线字段分组
    if line_field is not None and line_field in collar_df.columns:
        lines = collar_df[line_field].dropna().unique()
        if len(lines) >= 2:
            x_col, y_col = find_xy_columns(collar_df)
            centers = []
            for line in sorted(lines):
                line_df = collar_df[collar_df[line_field] == line]
                centers.append((line_df[x_col].mean(), line_df[y_col].mean()))

            dx_s, dy_s = strike_direction[0], strike_direction[1]
            norm_xy = math.sqrt(dx_s**2 + dy_s**2)
            if norm_xy < 1e-6:
                return 50.0
            dx_s, dy_s = dx_s / norm_xy, dy_s / norm_xy

            projections = sorted(cx * dx_s + cy * dy_s for cx, cy in centers)
            gaps = [projections[i] - projections[i-1]
                    for i in range(1, len(projections))
                    if projections[i] - projections[i-1] > same_line_threshold]
            if gaps:
                return float(np.median(gaps))

    # 方法2：投影法（无勘探线字段时）
    x_col, y_col = find_xy_columns(collar_df)
    dx_s, dy_s = strike_direction[0], strike_direction[1]
    norm_xy = math.sqrt(dx_s**2 + dy_s**2)
    if norm_xy < 1e-6:
        return 50.0
    dx_s, dy_s = dx_s / norm_xy, dy_s / norm_xy

    projections = sorted(collar_df[x_col] * dx_s + collar_df[y_col] * dy_s)
    gaps = [projections[i] - projections[i-1]
            for i in range(1, len(projections))
            if projections[i] - projections[i-1] > same_line_threshold]
    if not gaps:
        return 50.0
    return float(np.median(gaps))


print("[5b] 计算勘探线间距...")

# 自动检测勘探线字段
line_field = None
for candidate in ["勘探线", "勘探线号"]:
    if candidate in db.collar.columns:
        line_field = candidate
        break

spacing = compute_exploration_line_spacing(
    db.collar, ellipsoid.pc1_direction, line_field=line_field
)
print(f"  勘探线间距: {spacing:.1f}m")
print(f"  钻孔数: {len(db.collar)}")

5c. 搜索椭球体半径计算

基于勘探线间距和矿体几何比例推导三轴搜索半径：

def compute_search_radii(ellipsoid_params, line_spacing, composite_length,
                          spacing_multiplier=1.0, short_sample_multiplier=3):
    """根据勘探线间距和矿体几何计算搜索椭球体三轴半径

    算法：
      (1) 长半轴 = 勘探线间距 × spacing_multiplier（建议 1.0~1.2）
      (2) 次半轴 = 长半轴 × (延伸长度 / 矿体走向长度)
      (3) 短半轴 = max(长半轴 × (厚度 / 矿体走向长度),
                        组合样长 × short_sample_multiplier)
    """
    strike_len = 2.0 * ellipsoid_params.ellipsoid_main_radius
    extension_len = 2.0 * ellipsoid_params.second_radius
    thickness = 2.0 * ellipsoid_params.short_radius

    main_radius = line_spacing * spacing_multiplier

    second_rate = (extension_len / strike_len) if strike_len > 0 else 0.8
    second_radius = main_radius * second_rate

    short_by_ratio = main_radius * (thickness / strike_len) if strike_len > 0 else main_radius * 0.2
    short_min = composite_length * short_sample_multiplier
    short_radius = max(short_by_ratio, short_min)
    short_rate = short_radius / main_radius if main_radius > 0 else 0.2

    return {
        "main_radius": round(main_radius, 2),
        "second_main_rate": round(second_rate, 4),
        "short_main_rate": round(short_rate, 4),
    }


print("[5c] 计算搜索椭球体半径...")

SPACING_MULTIPLIER = 1.2      # 勘探线间距倍数（1.0~1.2）
SHORT_SAMPLE_MULTIPLIER = 3   # 短轴组合样长倍数（2~4）

radii = compute_search_radii(
    ellipsoid, spacing, composite_length=2.0,
    spacing_multiplier=SPACING_MULTIPLIER,
    short_sample_multiplier=SHORT_SAMPLE_MULTIPLIER,
)

print(f"  长半轴 = 勘探线间距({spacing:.1f}m) × {SPACING_MULTIPLIER} = {radii['main_radius']}m")
print(f"  次半轴/主轴 = {radii['second_main_rate']:.4f}")
print(f"  短半轴/主轴 = {radii['short_main_rate']:.4f}")

5d. 迭代 IDW 估值

from dimine_python_sdk.lib.prospecting import BlockModelEvaluator
from dimine_python_sdk.lib.prospecting.models import (
    BlockModelDistancePowerParams,
    EntityConstraint,
)

print("[5d] 迭代距离幂估值...")

# 约束条件：矿体实体内部
constraints = [
    EntityConstraint(
        file=dmf_constraint_path,
        range=0,            # 0=内部，1=外部
        bool_operate="and",
    ),
]

# IDW 参数
IDW_POWER = 2
MAX_IDW_ITERATIONS = 10
angle_kwargs = {
    "angle_main": ellipsoid.ellipsoid_azimuth,
    "angle_second": ellipsoid.ellipsoid_plunge,
    "angle_short": ellipsoid.ellipsoid_dip,
}

current_main_radius = radii["main_radius"]
evaluator = BlockModelEvaluator()

for iteration in range(1, MAX_IDW_ITERATIONS + 1):
    sec_r = current_main_radius * radii["second_main_rate"]
    short_r = current_main_radius * radii["short_main_rate"]

    print(f"\n  --- IDW 迭代 #{iteration}/{MAX_IDW_ITERATIONS} ---")
    print(f"  搜索椭球体: 主轴={current_main_radius:.1f}m, "
          f"次轴={sec_r:.1f}m, 短轴={short_r:.1f}m")

    idw_params = BlockModelDistancePowerParams(
        sample_file=dmg_combined_path,
        variable=valid_grade_fields,
        extra_attribute="",
        power=IDW_POWER,
        min_value=0.0001,
        max_value=9999.0,
        single_block_min=4,
        single_block_max=8,
        sub_block_main=1,
        sub_block_second=1,
        sub_block_short=1,
        main_radius=current_main_radius,
        second_main_rate=radii["second_main_rate"],
        short_main_rate=radii["short_main_rate"],
        **angle_kwargs,
        octant_max=0,
        project_count_min=0,
        single_project_sample_max=0,
    )

    result = evaluator.distance_power_evaluation(
        block_model_file=dmb_output_path,
        constraint_params=constraints,
        evaluation_params=idw_params,
        overwrite_result=True,
    )

    print(f"  估值结果: 成功={result.success}")
    if result.elements:
        print(f"  {'元素':<8} {'单元块总数':>12} {'已估值':>12} "
              f"{'本次估值':>12} {'未估值':>12}")
        print(f"  {'-'*8} {'-'*12} {'-'*12} {'-'*12} {'-'*12}")
        for e in result.elements:
            print(f"  {e.element:<8} {e.total_blocks:>12,} "
                  f"{e.evaluated_blocks:>12,} "
                  f"{e.current_evaluated_blocks:>12,} "
                  f"{e.unevaluated_blocks:>12,}")

    # 检查是否还有未估值块段
    has_unevaluated = any(e.unevaluated_blocks > 0 for e in result.elements)

    if not has_unevaluated:
        print(f"\n  [收敛] 第 {iteration} 轮后所有块段已估值完成！")
        break

    # 仍有未估值块段，主轴半径翻倍
    current_main_radius *= 2
    print(f"  仍有未估值块段 → 主轴半径扩大至 {current_main_radius:.1f}m")
else:
    remaining = []
    if result.elements:
        remaining = [f"{e.element}={e.unevaluated_blocks}"
                     for e in result.elements if e.unevaluated_blocks > 0]
    print(f"\n  [警告] 达到最大迭代次数，剩余未估值: {', '.join(remaining) if remaining else '无'}")

# 验证估值后字段
with DmbFile(dmb_output_path) as block_after:
    print(f"\n估值后字段定义: {block_after.field_definitions}")

IDW 估值参数说明

参数	默认值	说明
power	2	距离幂次，常用 1–3
single_block_min	4	单块最小样品数
single_block_max	8	单块最大样品数
octant_max	0	八分圆最大样品数，0=不限制
project_count_min	0	最少工程数（保证空间代表性）
single_project_sample_max	0	单工程最大样品数（避免单孔主导）
min_value / max_value	0.0001 / 9999	估值范围截断
sub_block_main/second/short	1	子块离散数，>1 可提高精度

---

完整管线脚本

以下是把上述所有步骤整合在一起的端到端脚本：

"""
储量估算完整管线
输入: 01_地质勘探数据库.dmd + 约束_矿体.dmf
输出: 估值后的块段模型 .dmb + 储量统计结果

Usage:
    python reserve_estimation_pipeline.py
"""
import math
import os
import shutil
import sys
from pathlib import Path

import numpy as np

# ==============================================================================
# 路径与参数配置
# ==============================================================================
# 输入文件
DMD_FILE = Path("C:/Users/76081/Desktop/yanshi/储量估值/01_地质勘探数据库.dmd")
DMF_CONSTRAINT = Path("C:/Users/76081/Desktop/yanshi/储量估值/约束_矿体.dmf")

# 输出目录
OUTPUT_DIR = Path("./reserves_output")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# 中间/输出文件
DMG_RAW = OUTPUT_DIR / "01_raw.dmg"
DMG_HG = OUTPUT_DIR / "02_high_grade.dmg"
DMG_COMBINED = OUTPUT_DIR / "03_combined.dmg"
DMB_OUTPUT = OUTPUT_DIR / "04_block_model.dmb"

# -- 特高品位处理 --
HG_AVERAGE_MULTIPLE = 8
HG_PROCESS_MODE = 0
HG_REPLACE_METHOD = 3

# -- 样长组合 --
COMPOSITE_LENGTH = 2.0
COMBINE_PERCENT = 0.75

# -- 块段模型 --
BLOCK_SIZE_X, BLOCK_SIZE_Y, BLOCK_SIZE_Z = 5.0, 5.0, 5.0
MARGIN_XY, MARGIN_Z = 20.0, 10.0

# -- IDW 估值 --
IDW_POWER = 2
SPACING_MULTIPLIER = 1.2
SHORT_SAMPLE_MULTIPLIER = 3
MAX_IDW_ITERATIONS = 10

# -- DIMINE_HOME --
if not os.environ.get("DIMINE_HOME"):
    os.environ["DIMINE_HOME"] = r"E:\数采软件\x64\x64"

# ==============================================================================
# 导入 SDK
# ==============================================================================
from dimine_python_sdk.lib.io import DmdFile, DmgFile, DmbFile, convert_file, DmfFile
from dimine_python_sdk.lib.prospecting import (
    DrillFunctionWrapper,
    BlockModelEvaluator,
    create_block_model,
    compute_search_ellipsoid_from_dmf,
)
from dimine_python_sdk.lib.prospecting.models import (
    HighGradeProcessParam,
    SampleLengthCombineParam,
    CreateBlockModelParams,
    FieldDefinition,
    BlockModelDistancePowerParams,
    EntityConstraint,
)
from dimine_python_sdk.models.types import Shell, Point as DmfPoint, Line as DmfLine


# ==============================================================================
# 工具函数
# ==============================================================================
def find_xy_columns(collar_df):
    x_candidates = ["横坐标", "EAST", "X", "Easting", "easting", "x"]
    y_candidates = ["纵坐标", "NORTH", "Y", "Northing", "northing", "y"]
    x_col = next((c for c in x_candidates if c in collar_df.columns), None)
    y_col = next((c for c in y_candidates if c in collar_df.columns), None)
    if x_col is None or y_col is None:
        numeric_cols = collar_df.select_dtypes(include="number").columns.tolist()
        if len(numeric_cols) >= 2:
            x_col, y_col = numeric_cols[0], numeric_cols[1]
        else:
            raise ValueError("无法识别X/Y坐标列")
    return x_col, y_col


def compute_exploration_line_spacing(collar_df, strike_direction,
                                      line_field=None, same_line_threshold=5.0):
    if len(collar_df) < 2:
        return 50.0
    if line_field is not None and line_field in collar_df.columns:
        lines = collar_df[line_field].dropna().unique()
        if len(lines) >= 2:
            x_col, y_col = find_xy_columns(collar_df)
            centers = [(collar_df[collar_df[line_field] == l][x_col].mean(),
                         collar_df[collar_df[line_field] == l][y_col].mean())
                        for l in sorted(lines)]
            dx_s, dy_s = strike_direction[0], strike_direction[1]
            norm_xy = math.sqrt(dx_s**2 + dy_s**2)
            if norm_xy < 1e-6:
                return 50.0
            dx_s, dy_s = dx_s / norm_xy, dy_s / norm_xy
            proj = sorted(cx * dx_s + cy * dy_s for cx, cy in centers)
            gaps = [proj[i] - proj[i-1] for i in range(1, len(proj))
                    if proj[i] - proj[i-1] > same_line_threshold]
            if gaps:
                return float(np.median(gaps))
    x_col, y_col = find_xy_columns(collar_df)
    dx_s, dy_s = strike_direction[0], strike_direction[1]
    norm_xy = math.sqrt(dx_s**2 + dy_s**2)
    if norm_xy < 1e-6:
        return 50.0
    dx_s, dy_s = dx_s / norm_xy, dy_s / norm_xy
    proj = sorted(collar_df[x_col] * dx_s + collar_df[y_col] * dy_s)
    gaps = [proj[i] - proj[i-1] for i in range(1, len(proj))
            if proj[i] - proj[i-1] > same_line_threshold]
    return float(np.median(gaps)) if gaps else 50.0


def compute_search_radii(ellipsoid_params, line_spacing, composite_length,
                          spacing_multiplier=1.0, short_sample_multiplier=3):
    strike_len = 2.0 * ellipsoid_params.ellipsoid_main_radius
    extension_len = 2.0 * ellipsoid_params.second_radius
    thickness = 2.0 * ellipsoid_params.short_radius
    main_radius = line_spacing * spacing_multiplier
    second_rate = (extension_len / strike_len) if strike_len > 0 else 0.8
    short_by_ratio = main_radius * (thickness / strike_len) if strike_len > 0 else main_radius * 0.2
    short_min = composite_length * short_sample_multiplier
    short_radius = max(short_by_ratio, short_min)
    short_rate = short_radius / main_radius if main_radius > 0 else 0.2
    return {
        "main_radius": round(main_radius, 2),
        "second_main_rate": round(second_rate, 4),
        "short_main_rate": round(short_rate, 4),
    }


def compute_dmf_bounding_box(dmf_path):
    dmf = DmfFile(dmf_path)
    all_points = []
    for layer in dmf.layers:
        for entity in layer.entities:
            if isinstance(entity, Shell):
                pts = entity.geometry.points
                if pts is not None and len(pts) > 0:
                    all_points.append(pts)
            elif isinstance(entity, DmfLine):
                geom = entity.geometry
                if geom is not None and len(geom) > 0:
                    all_points.append(geom)
            elif isinstance(entity, DmfPoint):
                geom = entity.geometry
                if geom is not None and len(geom) > 0:
                    all_points.append(geom.reshape(1, 3))
    pts = np.vstack(all_points)
    dmf.close()
    return {
        "min_x": float(pts[:, 0].min()), "max_x": float(pts[:, 0].max()),
        "min_y": float(pts[:, 1].min()), "max_y": float(pts[:, 1].max()),
        "min_z": float(pts[:, 2].min()), "max_z": float(pts[:, 2].max()),
    }


# ==============================================================================
# 主流程
# ==============================================================================
def main():
    print("=" * 60)
    print("  储量估算完整管线")
    print("=" * 60)

    # --- [1] DMD → DMG 格式转换 ---
    print("\n[1/7] DMD → DMG 格式转换...")
    convert_file(str(DMD_FILE), str(DMG_RAW))

    dmg = DmgFile(str(DMG_RAW))
    segments_df = dmg.to_segments_dataframe()
    dmg.close()

    # 自动检测品位字段
    numeric_cols = segments_df.select_dtypes(include="number").columns.tolist()
    excluded = {"起始", "结束", "样长", "X", "Y", "Z", "覆盖长度", "覆盖率",
                "测斜深度", "方位角", "倾角", "样本编号", "钻孔编号",
                "start_x", "start_y", "start_z", "end_x", "end_y", "end_z",
                "SAMPLE-ID", "SAMFROM", "SAMTO"}
    candidate_fields = [c for c in numeric_cols if c not in excluded]
    valid_grade_fields = []
    for field in candidate_fields:
        col = segments_df[field]
        if not (col.isna() | col.isin([float("inf"), float("-inf")]) | (col <= 0)).all():
            valid_grade_fields.append(field)
    print(f"  有效品位字段: {valid_grade_fields}")

    db = DmdFile(str(DMD_FILE))
    print(f"  钻孔数: {len(db.collar)}, 样品数: {len(db.sample)}")

    # --- [2] 特高品位处理（迭代收敛） ---
    print(f"\n[2/7] 特高品位处理（{HG_AVERAGE_MULTIPLE}x 均值倍数）...")
    shutil.copy2(str(DMG_RAW), str(DMG_HG))

    for grade_field in valid_grade_fields:
        print(f"  处理元素: {grade_field}")
        for iteration in range(1, 6):
            dmg = DmgFile(str(DMG_HG))
            seg = dmg.to_segments_dataframe()
            dmg.close()

            grades = seg[grade_field].dropna()
            grades = grades[~grades.isin([float("inf"), float("-inf")])]
            grades = grades[grades > 0]

            avg = float(grades.mean())
            threshold = avg * HG_AVERAGE_MULTIPLE
            exceeding = int((grades > threshold).sum())

            print(f"    迭代#{iteration}: 均值={avg:.4f}, 阈值={threshold:.4f}, 超限={exceeding}")
            if exceeding == 0:
                print(f"    [收敛]")
                break

            hg_params = HighGradeProcessParam(
                input_file=str(DMG_HG), grade_field=grade_field,
                process_mode=HG_PROCESS_MODE, average_multiple=HG_AVERAGE_MULTIPLE,
                frequency=0.95, replace_method=HG_REPLACE_METHOD,
            )
            DrillFunctionWrapper.extra_high_grade_process(hg_params)

    # --- [3] 样长组合 ---
    print(f"\n[3/7] 样长组合（{COMPOSITE_LENGTH}m）...")
    combine_params = SampleLengthCombineParam(
        input_file=str(DMG_HG), combine_length=COMPOSITE_LENGTH,
        combine_percent=COMBINE_PERCENT, output_file=str(DMG_COMBINED),
    )
    DrillFunctionWrapper.sample_length_combine(combine_params)

    # --- [4] 创建空块段模型 ---
    print("\n[4/7] 创建空块段模型...")
    bbox = compute_dmf_bounding_box(str(DMF_CONSTRAINT))
    origin_x = bbox["min_x"] - MARGIN_XY
    origin_y = bbox["min_y"] - MARGIN_XY
    origin_z = bbox["min_z"] - MARGIN_Z
    nx = math.ceil((bbox["max_x"] - bbox["min_x"] + 2 * MARGIN_XY) / BLOCK_SIZE_X)
    ny = math.ceil((bbox["max_y"] - bbox["min_y"] + 2 * MARGIN_XY) / BLOCK_SIZE_Y)
    nz = math.ceil((bbox["max_z"] - bbox["min_z"] + 2 * MARGIN_Z) / BLOCK_SIZE_Z)

    block_params = CreateBlockModelParams(
        file_name=str(DMB_OUTPUT),
        origin_x=origin_x, origin_y=origin_y, origin_z=origin_z,
        size_x=BLOCK_SIZE_X, size_y=BLOCK_SIZE_Y, size_z=BLOCK_SIZE_Z,
        length_x=nx * BLOCK_SIZE_X, length_y=ny * BLOCK_SIZE_Y, length_z=nz * BLOCK_SIZE_Z,
        field_definition=[
            FieldDefinition(field_name="岩石类型", field_type="字节型"),
            FieldDefinition(field_name="体重", field_type="浮点型"),
        ],
    )
    create_block_model(block_params)
    print(f"  块数: {nx} × {ny} × {nz} = {nx * ny * nz:,}")

    # --- [5] MVEE 椭球体 + 勘探线间距 ---
    print("\n[5/7] 计算搜索椭球体...")
    ellipsoid = compute_search_ellipsoid_from_dmf(str(DMF_CONSTRAINT), voxel_size=5.0)
    print(f"  方位角={ellipsoid.ellipsoid_azimuth:.1f}°, "
          f"倾伏角={ellipsoid.ellipsoid_plunge:.1f}°, "
          f"倾角={ellipsoid.ellipsoid_dip:.1f}°")

    line_field = next((c for c in ["勘探线", "勘探线号"] if c in db.collar.columns), None)
    spacing = compute_exploration_line_spacing(
        db.collar, ellipsoid.pc1_direction, line_field=line_field
    )
    print(f"  勘探线间距: {spacing:.1f}m")

    radii = compute_search_radii(ellipsoid, spacing, COMPOSITE_LENGTH,
                                  spacing_multiplier=SPACING_MULTIPLIER,
                                  short_sample_multiplier=SHORT_SAMPLE_MULTIPLIER)
    print(f"  搜索半径: 主轴={radii['main_radius']}m")

    # --- [6] 迭代 IDW 估值 ---
    print(f"\n[6/6] 距离幂估值（IDW, 幂次={IDW_POWER}）...")
    constraints = [EntityConstraint(file=str(DMF_CONSTRAINT), range=0, bool_operate="and")]
    angle_kwargs = {
        "angle_main": ellipsoid.ellipsoid_azimuth,
        "angle_second": ellipsoid.ellipsoid_plunge,
        "angle_short": ellipsoid.ellipsoid_dip,
    }
    current_main_radius = radii["main_radius"]
    evaluator = BlockModelEvaluator()

    for iteration in range(1, MAX_IDW_ITERATIONS + 1):
        idw_params = BlockModelDistancePowerParams(
            sample_file=str(DMG_COMBINED), variable=valid_grade_fields,
            extra_attribute="", power=IDW_POWER,
            min_value=0.0001, max_value=9999.0,
            single_block_min=4, single_block_max=8,
            sub_block_main=1, sub_block_second=1, sub_block_short=1,
            main_radius=current_main_radius,
            second_main_rate=radii["second_main_rate"],
            short_main_rate=radii["short_main_rate"],
            octant_max=0, project_count_min=0, single_project_sample_max=0,
            **angle_kwargs,
        )
        result = evaluator.distance_power_evaluation(
            block_model_file=str(DMB_OUTPUT),
            constraint_params=constraints,
            evaluation_params=idw_params,
            overwrite_result=True,
        )
        print(f"  迭代#{iteration}: 主轴={current_main_radius:.1f}m, 成功={result.success}")
        if result.elements:
            for e in result.elements:
                print(f"    {e.element}: 总块={e.total_blocks}, 已估值={e.evaluated_blocks}, "
                      f"本次={e.current_evaluated_blocks}, 未估值={e.unevaluated_blocks}")

        if not any(e.unevaluated_blocks > 0 for e in result.elements):
            print(f"  [收敛] 所有块段已估值完成！")
            break
        current_main_radius *= 2
        print(f"  主轴半径扩大至 {current_main_radius:.1f}m")
    else:
        print(f"  [警告] 达到最大迭代次数")

    # 清理
    db.close()

    print(f"\n{'=' * 60}")
    print(f"  流程完成！输出目录: {OUTPUT_DIR.resolve()}")
    print(f"{'=' * 60}")
    for label, path in [
        ("原始 DMG", DMG_RAW), ("特高品位处理后 DMG", DMG_HG),
        ("样长组合后 DMG", DMG_COMBINED), ("块段模型（估值后）", DMB_OUTPUT),
    ]:
        size_kb = path.stat().st_size / 1024 if path.exists() else 0
        print(f"  {label}: {path.name} ({size_kb:.1f} KB)")


if __name__ == "__main__":
    main()

---

输出说明

输出文件清单

序号	文件	说明
①	01_raw.dmg	原始 DMG 分段文件（DMD→DMG 转换结果）
②	02_high_grade.dmg	特高品位处理后的 DMG
③	03_combined.dmg	样长组合后的 DMG
④	04_block_model.dmb	估值后的块段模型（含品位字段）

---

参数速查

特高品位处理

参数	建议值	说明
average_multiple	8（铁）/ 6（有色、金）/ 7（煤、非金属）	均值倍数阈值
process_mode	0	0=国内规范，1=国外规范（累积频率）
replace_method	3	3=矿体平均品位法（推荐）

块段模型

参数	建议值	说明
size_x/y/z	5m / 5m / 5m	单元块尺寸，根据矿体规模和开采方式调整
margin_xy	20m	XY 方向边距
margin_z	10m	Z 方向边距

搜索椭球体

参数	建议值	说明
spacing_multiplier	1.0–1.2	长半轴 = 勘探线间距 × 倍数
short_sample_multiplier	2–4	短半轴 ≥ 组合样长 × 倍数

---

异常处理

from dimine_python_sdk.lib.prospecting import (
    BlockDataError,
    BlockModelCreateError,
    BlockModelEvaluationError,
)

try:
    evaluator.distance_power_evaluation(...)
except BlockModelEvaluationError as e:
    print(f"估值失败: {e}")
    if e.response:
        print(f"响应详情: {e.response}")
except BlockModelCreateError as e:
    print(f"块段模型创建失败: {e}")
except BlockDataError as e:
    print(f"块段模型数据异常: {e}")

---

相关参考

• API 文档：dimine_python_sdk.lib.prospecting
• 指南：储量估算
• 参数模型：dimine_python_sdk.lib.prospecting.models
• 钻孔数据库导入示例